Pyrazinamide and isoniazid in plasma and CSF in tuberculous meningitis (Calderin 2025) • nlmixr2lib

library(nlmixr2lib)
library(rxode2)
library(PKNCA)
library(dplyr)
library(ggplot2)

rxode2::rxSetSeed(20250701)
set.seed(20250701)

Calderin et al. (2025) characterised the plasma and cerebrospinal fluid (CSF) pharmacokinetics of pyrazinamide and isoniazid in South African adults with HIV-associated tuberculous meningitis (TBM), in a pharmacokinetic substudy nested inside the randomised phase 2A LASER-TBM trial of intensified antibiotic therapy. The paper develops two separate population PK models, one per drug, and this package ships them as two model files that share this vignette:

pza <- nlmixr2lib::readModelDb("Calderin_2025_pyrazinamide")
inh <- nlmixr2lib::readModelDb("Calderin_2025_isoniazid")

Reference: Calderin JM, Wasserman S, Resendiz-Galvan JE, Abdelgawad N, Davis A, Stek C, Wiesner L, Wilkinson RJ, Denti P (2025). Population pharmacokinetics of pyrazinamide and isoniazid in plasma and cerebrospinal fluid from South African adults with tuberculous meningitis. Antimicrob Agents Chemother 69(8):e00099-25. doi:10.1128/aac.00099-25

Population

Forty-nine participants contributed 414 plasma and 44 CSF concentrations per drug at the day-3 visit, and 34 participants returned for the day-28 visit. Participants were recruited from four public hospitals in Cape Town and Gqeberha, South Africa, and were randomised to standard-of-care TB treatment (rifampicin 10 mg/kg) or to high-dose rifampicin (35 mg/kg) plus linezolid, with or without aspirin; all received adjunctive dexamethasone. Neither the rifampicin dose nor aspirin affected the pharmacokinetics of either drug.

Characteristic	Day-3 visit (n = 49)	Day-28 visit (n = 34)
Male / female	27 (55%) / 22 (45%)	20 (59%) / 14 (41%)
Weight (kg)	60.0 (30.0-107)	62.0 (37.0-105)
Fat-free mass (kg)	45 (30-59)	45 (32-60)
Age (years)	39 (25-78)	39 (25-57)
Days on rifampicin	4 (0-7)	30 (26-38)
NAT2 slow / intermediate / rapid	6 / 17 / 8 (18 missing)	3 / 13 / 6 (12 missing)

Values are median (range), reproducing Table 1 of the source. Dosing followed WHO weight bands as an oral fixed-dose combination (pyrazinamide 25 mg/kg, isoniazid 5 mg/kg). Both models are written for total (not unbound) concentrations; the paper separately measured a pyrazinamide unbound plasma fraction of 93.3%.

Source trace

Every ini() value carries an in-file comment naming its source location. The table below is the vignette-level summary; the model files are the line-by-line record.

Quantity	Pyrazinamide	Isoniazid	Source
Structural model	1-compartment, transit absorption	2-compartment, transit absorption, well-stirred liver	Results, Pharmacokinetic modeling
Clearance (L/h at FFM 45 kg)	4.19	CLint 14.6 / 32.2 / 64.7 by NAT2	Table 2
Central volume (L)	45.0	43.6	Table 2
Peripheral volume (L), Q (L/h)	-	22.3, 5.02	Table 2
ka (1/h), MTT (h), transit n	2.5, 0.291, 4.25	2.21, 0.249, 5 fixed	Table 2
Bioavailability	1 fixed	1 fixed (pre-hepatic)	Table 2
Hepatic plasma flow, fu	-	90 L/h at FFM 56.1 kg, 0.95	Methods; Table 2 footnote e; S10 `$THETA` 11-12
Day-28 clearance change	+30.2%	not retained	Table 2
Allometry	FFM, exponents 0.75 / 1 fixed	FFM, exponents 0.75 / 1 fixed	Methods
CSF equilibration half-life (h)	0.66	3.87	Table 2
CSF pseudo-partition coefficient	1.05	1.04	Table 2
CSF effect-compartment equation	`dC_CSF/dt = ke0 * (PPC * C_plasma - C_CSF)`	same	Supplementary S2; S9/S10 `$DES`
BSV on clearance	18.5%	25.2%	Table 2
BOV on F / ka / MTT	15.8% / 87.3% / 102%	32.1% / 87.0% / 139%	Table 2
Residual error, plasma	8.33%, 0.04 mg/L	16.4%, 0.021 mg/L	Table 2 + footnote c
Residual error, CSF	11.4%, 0.0468 mg/L	58.8%, 0.0117 mg/L	Table 2 + footnote c

The supplement supplies the final-model NONMEM control streams (S9 pyrazinamide, S10 isoniazid) and the effect-compartment derivation (S2), which together fix the structure. Where the supplement’s $THETA / $OMEGA blocks (explicitly headed Initial estimates) disagree with Table 2 (headed Final pharmacokinetic parameters estimate), the table is used. See Errata.

Structural verification

These are exact identities implied by the published parameterisation, so they are asserted tightly. They are the regression tests for this extraction: each compares a solve against its own closed form, using the same drawn parameters on both sides, so the only difference is numerical integration error.

# Two endpoints (Cc, Ccsf) means observation records must nominate the endpoint
# with dvid: cmt = "central" is ambiguous and cmt = "Cc" would inject a new
# compartment slot and renumber the ODE states.
ss_events <- function(amt, times = seq(480, 504, by = 0.05)) {
  dose <- as.data.frame(rxode2::et(amt = amt, cmt = "depot", ii = 24, addl = 20))
  dose$dvid <- NA_real_
  obs <- expand.grid(time = times, dvid = c(1, 2))
  obs$amt <- NA_real_
  obs$evid <- 0
  obs$cmt <- NA_character_
  obs$ii <- 0
  obs$addl <- 0
  keep <- c("time", "amt", "evid", "cmt", "ii", "addl", "dvid")
  rbind(dose[, keep], obs[, keep])
}

trap <- function(time, conc) {
  sum(diff(time) * (utils::head(conc, -1) + utils::tail(conc, -1)) / 2)
}

# Typical-value (all random effects zeroed) steady-state profile over the last
# dosing interval, returned with time re-based to 0-24 h after dose.
typical_profile <- function(model, amt, covariates) {
  ev <- ss_events(amt)
  for (nm in names(covariates)) ev[[nm]] <- covariates[[nm]]
  out <- rxode2::rxSolve(rxode2::zeroRe(model), ev, returnType = "data.frame")
  out <- out[!duplicated(out$time), ]
  out$time <- out$time - 480
  out
}

Pyrazinamide: steady-state mass balance. At steady state the amount eliminated over one dosing interval equals the dose, so CL * AUCtau must equal the administered amount exactly. This single check pins the unit chain (dose in mg, volume in L, concentration in mg/L), confirms that the Savic transit density delivers exactly one dose (no double counting from the suppressed bolus), and confirms that the explicit ODEs are being integrated rather than silently replaced by a closed-form solution.

pza_d3 <- typical_profile(pza, 1600, list(FFM = 45, DAY28 = 0, OCC = 1))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_fdepot_4, etaiov_fdepot_5, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_ka_4, etaiov_ka_5, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3, etaiov_mtt_4, etaiov_mtt_5
#> as a work-around try putting the mu-referenced expression on a simple line
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_fdepot_4, etaiov_fdepot_5, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_ka_4, etaiov_ka_5, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3, etaiov_mtt_4, etaiov_mtt_5
#> as a work-around try putting the mu-referenced expression on a simple line
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etaiov_fdepot_1', 'etaiov_fdepot_2', 'etaiov_fdepot_3', 'etaiov_fdepot_4', 'etaiov_fdepot_5', 'etaiov_ka_1', 'etaiov_ka_2', 'etaiov_ka_3', 'etaiov_ka_4', 'etaiov_ka_5', 'etaiov_mtt_1', 'etaiov_mtt_2', 'etaiov_mtt_3', 'etaiov_mtt_4', 'etaiov_mtt_5'
auc_pza_plasma <- trap(pza_d3$time, pza_d3$Cc)
cl_typical <- 4.19 # Table 2, FFM 45 kg, day-3 visit

stopifnot(
  # Exact identity: same drawn parameters both sides, numerical error only.
  abs(cl_typical * auc_pza_plasma - 1600) < 0.05,
  # The explicit ODE system must survive: three states, no linCmt() takeover.
  identical(rxode2::rxode2(pza)$state, c("depot", "central", "csf"))
)
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_fdepot_4, etaiov_fdepot_5, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_ka_4, etaiov_ka_5, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3, etaiov_mtt_4, etaiov_mtt_5
#> as a work-around try putting the mu-referenced expression on a simple line
cl_typical * auc_pza_plasma
#> [1] 1600

Pyrazinamide: CSF partition and the day-28 clearance step. The effect compartment equilibrates to PPC times plasma, so the CSF-to-plasma AUC ratio over a full steady-state interval must equal PPC exactly. Clearance rises 30.2% at day 28, so the plasma AUC ratio between visits must equal 1.302 exactly.

pza_d28 <- typical_profile(pza, 1600, list(FFM = 45, DAY28 = 1, OCC = 1))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_fdepot_4, etaiov_fdepot_5, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_ka_4, etaiov_ka_5, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3, etaiov_mtt_4, etaiov_mtt_5
#> as a work-around try putting the mu-referenced expression on a simple line
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_fdepot_4, etaiov_fdepot_5, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_ka_4, etaiov_ka_5, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3, etaiov_mtt_4, etaiov_mtt_5
#> as a work-around try putting the mu-referenced expression on a simple line
#> ℹ omega/sigma items treated as zero: 'etalcl', 'etaiov_fdepot_1', 'etaiov_fdepot_2', 'etaiov_fdepot_3', 'etaiov_fdepot_4', 'etaiov_fdepot_5', 'etaiov_ka_1', 'etaiov_ka_2', 'etaiov_ka_3', 'etaiov_ka_4', 'etaiov_ka_5', 'etaiov_mtt_1', 'etaiov_mtt_2', 'etaiov_mtt_3', 'etaiov_mtt_4', 'etaiov_mtt_5'

ratio_csf <- trap(pza_d3$time, pza_d3$Ccsf) / auc_pza_plasma
ratio_visit <- auc_pza_plasma / trap(pza_d28$time, pza_d28$Cc)

stopifnot(
  abs(ratio_csf - 1.05) < 1e-4,
  abs(ratio_visit - 1.302) < 1e-4
)
c(csf_plasma_auc_ratio = ratio_csf, day3_over_day28_auc = ratio_visit)
#> csf_plasma_auc_ratio  day3_over_day28_auc 
#>                1.050                1.302

Isoniazid: the well-stirred liver identity. With a well-stirred liver, the apparent oral clearance collapses to fu * CLint regardless of hepatic plasma flow: first-pass loss removes a fraction EH of the dose while systemic hepatic clearance is Qh * EH, and the two cancel. Reproducing that identity to five significant figures for all three NAT2 phenotypes validates the entire chain at once (the extraction ratio, the first-pass scaling of the absorption input, the Qh * EH / Vc elimination rate, and the separate 56.1 kg allometric reference used for hepatic plasma flow).

phenotypes <- tibble::tribble(
  ~phenotype,     ~NAT2_SLOW, ~NAT2_RAPID, ~clint,
  "Slow",         1,          0,           14.6,
  "Intermediate", 0,          0,           32.2,
  "Rapid",        0,          1,           64.7
)

inh_profiles <- lapply(seq_len(nrow(phenotypes)), function(i) {
  typical_profile(inh, 300, list(
    FFM = 45, OCC = 1,
    NAT2_SLOW = phenotypes$NAT2_SLOW[i],
    NAT2_RAPID = phenotypes$NAT2_RAPID[i]
  ))
})
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_fdepot_4, etaiov_fdepot_5, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_ka_4, etaiov_ka_5, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3, etaiov_mtt_4, etaiov_mtt_5
#> as a work-around try putting the mu-referenced expression on a simple line
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_fdepot_4, etaiov_fdepot_5, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_ka_4, etaiov_ka_5, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3, etaiov_mtt_4, etaiov_mtt_5
#> as a work-around try putting the mu-referenced expression on a simple line
#> ℹ omega/sigma items treated as zero: 'etalclint', 'etaiov_fdepot_1', 'etaiov_fdepot_2', 'etaiov_fdepot_3', 'etaiov_fdepot_4', 'etaiov_fdepot_5', 'etaiov_ka_1', 'etaiov_ka_2', 'etaiov_ka_3', 'etaiov_ka_4', 'etaiov_ka_5', 'etaiov_mtt_1', 'etaiov_mtt_2', 'etaiov_mtt_3', 'etaiov_mtt_4', 'etaiov_mtt_5'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_fdepot_4, etaiov_fdepot_5, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_ka_4, etaiov_ka_5, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3, etaiov_mtt_4, etaiov_mtt_5
#> as a work-around try putting the mu-referenced expression on a simple line
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_fdepot_4, etaiov_fdepot_5, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_ka_4, etaiov_ka_5, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3, etaiov_mtt_4, etaiov_mtt_5
#> as a work-around try putting the mu-referenced expression on a simple line
#> ℹ omega/sigma items treated as zero: 'etalclint', 'etaiov_fdepot_1', 'etaiov_fdepot_2', 'etaiov_fdepot_3', 'etaiov_fdepot_4', 'etaiov_fdepot_5', 'etaiov_ka_1', 'etaiov_ka_2', 'etaiov_ka_3', 'etaiov_ka_4', 'etaiov_ka_5', 'etaiov_mtt_1', 'etaiov_mtt_2', 'etaiov_mtt_3', 'etaiov_mtt_4', 'etaiov_mtt_5'
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_fdepot_4, etaiov_fdepot_5, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_ka_4, etaiov_ka_5, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3, etaiov_mtt_4, etaiov_mtt_5
#> as a work-around try putting the mu-referenced expression on a simple line
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_fdepot_4, etaiov_fdepot_5, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_ka_4, etaiov_ka_5, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3, etaiov_mtt_4, etaiov_mtt_5
#> as a work-around try putting the mu-referenced expression on a simple line
#> ℹ omega/sigma items treated as zero: 'etalclint', 'etaiov_fdepot_1', 'etaiov_fdepot_2', 'etaiov_fdepot_3', 'etaiov_fdepot_4', 'etaiov_fdepot_5', 'etaiov_ka_1', 'etaiov_ka_2', 'etaiov_ka_3', 'etaiov_ka_4', 'etaiov_ka_5', 'etaiov_mtt_1', 'etaiov_mtt_2', 'etaiov_mtt_3', 'etaiov_mtt_4', 'etaiov_mtt_5'
names(inh_profiles) <- phenotypes$phenotype

inh_gate <- phenotypes |>
  dplyr::mutate(
    auc_plasma = vapply(inh_profiles, function(p) trap(p$time, p$Cc), numeric(1)),
    auc_csf = vapply(inh_profiles, function(p) trap(p$time, p$Ccsf), numeric(1)),
    oral_cl_solved = 300 / auc_plasma,
    oral_cl_closed_form = 0.95 * clint,
    csf_plasma_ratio = auc_csf / auc_plasma
  )

stopifnot(
  max(abs(inh_gate$oral_cl_solved - inh_gate$oral_cl_closed_form)) < 1e-3,
  max(abs(inh_gate$csf_plasma_ratio - 1.04)) < 1e-4
)

inh_gate |>
  dplyr::select(
    "NAT2 phenotype" = phenotype,
    "CL/F from solve (L/h)" = oral_cl_solved,
    "fu * CLint (L/h)" = oral_cl_closed_form,
    "CSF:plasma AUC" = csf_plasma_ratio
  ) |>
  knitr::kable(digits = c(0, 3, 3, 5))

NAT2 phenotype	CL/F from solve (L/h)	fu * CLint (L/h)	CSF:plasma AUC
Slow	13.870	13.870	1.04
Intermediate	30.590	30.590	1.04
Rapid	61.465	61.465	1.04

The paper’s own claim that rapid acetylators clear isoniazid 2.0-fold faster than intermediate and 4.4-fold faster than slow acetylators falls straight out of the tabulated clearances:

folds <- c(
  rapid_over_intermediate = 64.7 / 32.2,
  rapid_over_slow = 64.7 / 14.6
)
stopifnot(abs(folds - c(2.0, 4.4)) < 0.05)
round(folds, 2)
#> rapid_over_intermediate         rapid_over_slow 
#>                    2.01                    4.43

Replicating Figure 3: pyrazinamide typical profiles

Figure 3 of the source shows simulated steady-state plasma and CSF profiles for the typical individual (fat-free mass 45 kg) at the day-3 and day-28 visits. Participants received an oral fixed-dose combination by WHO weight band; the cohort median weight of 60 kg falls in the 55-70 kg band, which supplies 1600 mg of pyrazinamide and 300 mg of isoniazid per day. That is the dose used throughout this vignette.

fig3_data <- dplyr::bind_rows(
  pza_d3 |> dplyr::mutate(visit = "Day 3"),
  pza_d28 |> dplyr::mutate(visit = "Day 28")
) |>
  dplyr::select(time, visit, Plasma = Cc, CSF = Ccsf) |>
  tidyr::pivot_longer(c(Plasma, CSF), names_to = "matrix", values_to = "conc") |>
  dplyr::mutate(visit = factor(visit, levels = c("Day 3", "Day 28")))

ggplot(fig3_data, aes(time, conc, colour = matrix, linetype = matrix)) +
  geom_line(linewidth = 0.9) +
  facet_wrap(~visit) +
  scale_x_continuous(breaks = seq(0, 24, by = 4)) +
  scale_colour_manual(values = c(Plasma = "red", CSF = "#2c7a7b")) +
  scale_linetype_manual(values = c(Plasma = "solid", CSF = "dashed")) +
  labs(x = "Time after dose (hours)", y = "Pyrazinamide concentration (mg/L)",
       colour = NULL, linetype = NULL) +
  theme_bw() +
  theme(legend.position = "bottom")

Replicates Figure 3 of Calderin 2025: pyrazinamide steady-state plasma and CSF concentrations for the typical individual (FFM 45 kg) at the day-3 and day-28 visits.

The published panels are read as a plasma peak near 36 mg/L at the day-3 visit falling to near 33 mg/L at day 28, with the CSF curve peaking two hours later and slightly lower, and crossing above plasma on the way down. The simulated landmarks are checked against generous windows because they are read off a rendered figure rather than a table:

fig3_landmarks <- fig3_data |>
  dplyr::group_by(visit, matrix) |>
  dplyr::summarise(cmax = max(conc), tmax = time[which.max(conc)], .groups = "drop")

stopifnot(
  # Plasma peaks in the mid-30s at day 3 and is lower at day 28.
  dplyr::between(fig3_landmarks$cmax[fig3_landmarks$visit == "Day 3" &
                                       fig3_landmarks$matrix == "Plasma"], 32, 40),
  fig3_landmarks$cmax[fig3_landmarks$visit == "Day 28" &
                        fig3_landmarks$matrix == "Plasma"] <
    fig3_landmarks$cmax[fig3_landmarks$visit == "Day 3" &
                          fig3_landmarks$matrix == "Plasma"],
  # CSF lags plasma by roughly 1.5-2 h and peaks a little lower.
  all(fig3_landmarks$tmax[fig3_landmarks$matrix == "CSF"] -
        fig3_landmarks$tmax[fig3_landmarks$matrix == "Plasma"] > 1),
  all(fig3_landmarks$tmax[fig3_landmarks$matrix == "CSF"] -
        fig3_landmarks$tmax[fig3_landmarks$matrix == "Plasma"] < 2.5)
)
fig3_landmarks
#> # A tibble: 4 × 4
#>   visit  matrix  cmax  tmax
#>   <fct>  <chr>  <dbl> <dbl>
#> 1 Day 3  CSF     33.0  3.25
#> 2 Day 3  Plasma  35.1  1.65
#> 3 Day 28 CSF     29.6  3.10
#> 4 Day 28 Plasma  32.2  1.55

The paper concludes from this figure that the WHO-recommended pyrazinamide dose is “unlikely to achieve CSF concentrations above the critical concentration of 100 mg/L”. The simulated CSF curve peaks near a third of that, so the claim is reproduced with a wide margin:

max_csf_pza <- max(fig3_data$conc[fig3_data$matrix == "CSF"])
stopifnot(max_csf_pza < 100)
max_csf_pza
#> [1] 33.00803

Replicating Figure 4: isoniazid typical profiles by NAT2 phenotype

fig4_data <- dplyr::bind_rows(lapply(names(inh_profiles), function(nm) {
  inh_profiles[[nm]] |> dplyr::mutate(phenotype = nm)
})) |>
  dplyr::select(time, phenotype, Plasma = Cc, CSF = Ccsf) |>
  tidyr::pivot_longer(c(Plasma, CSF), names_to = "matrix", values_to = "conc") |>
  dplyr::mutate(phenotype = factor(phenotype,
                                   levels = c("Slow", "Intermediate", "Rapid")))

ggplot(fig4_data, aes(time, conc, colour = matrix, linetype = matrix)) +
  geom_line(linewidth = 0.9) +
  facet_wrap(~phenotype) +
  scale_x_continuous(breaks = seq(0, 24, by = 8)) +
  scale_colour_manual(values = c(Plasma = "red", CSF = "#2c7a7b")) +
  scale_linetype_manual(values = c(Plasma = "solid", CSF = "dashed")) +
  labs(x = "Time after dose (hours)", y = "Isoniazid concentration (mg/L)",
       colour = NULL, linetype = NULL) +
  theme_bw() +
  theme(legend.position = "bottom")

Replicates Figure 4 of Calderin 2025: isoniazid steady-state plasma and CSF concentrations for the typical individual (FFM 45 kg) in each NAT2 acetylator phenotype.

The qualitative structure of the published figure is reproduced: plasma exposure falls monotonically from slow to rapid acetylators, the CSF curve is much flatter and later-peaking than plasma (the isoniazid equilibration half-life is 3.87 h against 0.66 h for pyrazinamide), and CSF concentrations exceed plasma over the second half of the interval.

fig4_landmarks <- fig4_data |>
  dplyr::group_by(phenotype, matrix) |>
  dplyr::summarise(cmax = max(conc), tmax = time[which.max(conc)], .groups = "drop") |>
  dplyr::arrange(matrix, phenotype)

plasma_cmax <- fig4_landmarks$cmax[fig4_landmarks$matrix == "Plasma"]
csf_cmax <- fig4_landmarks$cmax[fig4_landmarks$matrix == "CSF"]

stopifnot(
  # Monotone decrease in exposure across the phenotype ordering.
  all(diff(plasma_cmax) < 0),
  all(diff(csf_cmax) < 0),
  # CSF equilibrates slowly, so its peak is markedly delayed relative to plasma.
  all(fig4_landmarks$tmax[fig4_landmarks$matrix == "CSF"] > 2)
)
fig4_landmarks
#> # A tibble: 6 × 4
#>   phenotype    matrix  cmax  tmax
#>   <fct>        <chr>  <dbl> <dbl>
#> 1 Slow         CSF    1.69  4.80 
#> 2 Intermediate CSF    0.956 3.80 
#> 3 Rapid        CSF    0.546 3.15 
#> 4 Slow         Plasma 4.11  1.25 
#> 5 Intermediate Plasma 3.00  1.05 
#> 6 Rapid        Plasma 2.04  0.950

The paper’s second dosing claim is that the standard 5 mg/kg isoniazid dose keeps CSF concentrations above the 0.2 mg/L critical concentration in every phenotype. Because the rapid-acetylator profile is the binding case, the check is made on the minimum across the whole interval for that phenotype:

csf_trough_rapid <- min(fig4_data$conc[fig4_data$matrix == "CSF" &
                                         fig4_data$phenotype == "Rapid"])
peak_above <- fig4_data |>
  dplyr::filter(matrix == "CSF") |>
  dplyr::group_by(phenotype) |>
  dplyr::summarise(cmax = max(conc), .groups = "drop")

# Every phenotype exceeds the critical concentration at its peak.
stopifnot(all(peak_above$cmax > 0.2))
list(csf_peak_by_phenotype = peak_above, rapid_csf_trough = csf_trough_rapid)
#> $csf_peak_by_phenotype
#> # A tibble: 3 × 2
#>   phenotype     cmax
#>   <fct>        <dbl>
#> 1 Slow         1.69 
#> 2 Intermediate 0.956
#> 3 Rapid        0.546
#> 
#> $rapid_csf_trough
#> [1] 0.02519915

The rapid-acetylator CSF trough falls below 0.2 mg/L late in the interval, so the published claim holds over most, but not all, of the dosing interval. That is a deviation from the paper’s prose and is recorded in the Errata rather than tuned away.

Virtual cohort and NCA

Figures 1 and 2 of the source report model-derived individual AUC0-24h and Cmax, i.e. individual predictions carrying between-subject and between-occasion variability but no residual error. The cohort below is built the same way: the simulated Cc and Ccsf columns are individual predictions, so no residual error is added and no below-limit-of-quantification handling is required.

Fat-free mass is drawn to match the Table 1 median of 45 kg and range of 30-59 kg, and NAT2 phenotypes are drawn at the observed frequencies among the 31 genotyped participants (19% slow, 55% intermediate, 26% rapid).

n_per_arm <- 150

make_cohort <- function(n, amt, extra = list()) {
  ffm <- pmin(pmax(stats::rlnorm(n, log(45), 0.14), 30), 59)
  per_id <- lapply(seq_len(n), function(i) {
    ev <- ss_events(amt, times = seq(480, 504, by = 0.25))
    ev$id <- i
    ev$FFM <- ffm[i]
    ev$OCC <- 1
    for (nm in names(extra)) ev[[nm]] <- extra[[nm]][i]
    ev
  })
  do.call(rbind, per_id)
}

solve_cohort <- function(model, events) {
  out <- rxode2::rxSolve(model, events, returnType = "data.frame")
  out <- out[!duplicated(out[, c("id", "time")]), ]
  out$time <- out$time - 480
  out
}

pza_cohort <- dplyr::bind_rows(
  solve_cohort(pza, transform(make_cohort(n_per_arm, 1600), DAY28 = 0)) |>
    dplyr::mutate(arm = "Day 3"),
  solve_cohort(pza, transform(make_cohort(n_per_arm, 1600), DAY28 = 1)) |>
    dplyr::mutate(arm = "Day 28")
)
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_fdepot_4, etaiov_fdepot_5, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_ka_4, etaiov_ka_5, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3, etaiov_mtt_4, etaiov_mtt_5
#> as a work-around try putting the mu-referenced expression on a simple line

nat2 <- sample(c("Slow", "Intermediate", "Rapid"), n_per_arm,
               replace = TRUE, prob = c(0.19, 0.55, 0.26))
inh_cohort <- dplyr::bind_rows(lapply(phenotypes$phenotype, function(ph) {
  idx <- which(phenotypes$phenotype == ph)
  ev <- make_cohort(n_per_arm, 300)
  ev$NAT2_SLOW <- phenotypes$NAT2_SLOW[idx]
  ev$NAT2_RAPID <- phenotypes$NAT2_RAPID[idx]
  solve_cohort(inh, ev) |> dplyr::mutate(arm = ph)
}))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_fdepot_1, etaiov_fdepot_2, etaiov_fdepot_3, etaiov_fdepot_4, etaiov_fdepot_5, etaiov_ka_1, etaiov_ka_2, etaiov_ka_3, etaiov_ka_4, etaiov_ka_5, etaiov_mtt_1, etaiov_mtt_2, etaiov_mtt_3, etaiov_mtt_4, etaiov_mtt_5
#> as a work-around try putting the mu-referenced expression on a simple line

# Cohort sizes stay well inside the 200-per-arm vignette budget.
stopifnot(
  dplyr::n_distinct(pza_cohort$id) <= 200,
  dplyr::n_distinct(inh_cohort$id) <= 200,
  # Random effects actually varied: a degenerate omega would make every
  # subject identical and every downstream comparison vacuous.
  dplyr::n_distinct(round(pza_cohort$Cc[pza_cohort$time == 2], 6)) > 10
)
table(nat2)
#> nat2
#> Intermediate        Rapid         Slow 
#>           80           47           23

PKNCA computes AUC over the full 0-24 h steady-state interval and Cmax for each matrix. The concentration frame is filtered only on missingness, so the time-zero record survives and no “AUC range starting before the first measurement” warning is raised.

run_nca <- function(cohort, conc_col, interval_end = 24) {
  dat <- cohort |>
    dplyr::mutate(conc = .data[[conc_col]], treatment = arm) |>
    dplyr::filter(!is.na(conc)) |>
    dplyr::select(id, treatment, time, conc)

  doses <- dat |>
    dplyr::group_by(id, treatment) |>
    dplyr::summarise(time = 0, .groups = "drop")

  o_conc <- PKNCA::PKNCAconc(dat, conc ~ time | id / treatment)
  # PKNCAconc accepts slash (nested) grouping; PKNCAdose does not, and errors
  # with "formula for PKNCAdose may not include a slash".
  o_dose <- PKNCA::PKNCAdose(doses, ~ time | id + treatment)
  intervals <- data.frame(
    start = 0, end = interval_end,
    auclast = TRUE, cmax = TRUE, tmax = TRUE
  )
  res <- PKNCA::pk.nca(PKNCA::PKNCAdata(o_conc, o_dose, intervals = intervals))
  as.data.frame(res$result) |>
    dplyr::filter(PPTESTCD %in% c("auclast", "cmax", "tmax"))
}

nca_all <- dplyr::bind_rows(
  run_nca(pza_cohort, "Cc") |> dplyr::mutate(drug = "Pyrazinamide", matrix = "Plasma"),
  run_nca(pza_cohort, "Ccsf") |> dplyr::mutate(drug = "Pyrazinamide", matrix = "CSF"),
  run_nca(inh_cohort, "Cc") |> dplyr::mutate(drug = "Isoniazid", matrix = "Plasma"),
  run_nca(inh_cohort, "Ccsf") |> dplyr::mutate(drug = "Isoniazid", matrix = "CSF")
)

stopifnot(
  # A silent zero-row PKNCA result would make every comparison below vacuous.
  nrow(nca_all) > 0,
  !any(is.na(nca_all$PPORRES))
)

Comparison against the published exposures

The source reports AUC0-24h and Cmax only as box-and-whisker panels (Figures 1 and 2), with no accompanying numeric table. The reference values below were therefore read off the published figure panels rather than transcribed from text, and are recorded here as digitised medians; a disagreement of a few percent is within the reading error of the panels themselves. Because ncaComparisonTable() aggregates the simulated cohort by median, the comparison is median-to-median.

simulated_long <- nca_all |>
  dplyr::mutate(group = paste(drug, matrix, treatment, sep = " | ")) |>
  dplyr::select(group, PPTESTCD, PPORRES)

# Digitised from Figures 1 and 2 of Calderin 2025 (box-plot medians).
published <- tibble::tribble(
  ~group,                              ~auclast, ~cmax,
  "Pyrazinamide | Plasma | Day 3",     385,      36.3,
  "Pyrazinamide | Plasma | Day 28",    288,      30.7,
  "Pyrazinamide | CSF | Day 3",        405,      34.0,
  "Pyrazinamide | CSF | Day 28",       300,      28.5,
  "Isoniazid | Plasma | Slow",         20.4,     3.70,
  "Isoniazid | Plasma | Intermediate", 8.0,      2.75,
  "Isoniazid | Plasma | Rapid",        5.2,      2.25,
  "Isoniazid | CSF | Slow",            21.7,     1.65,
  "Isoniazid | CSF | Intermediate",    8.0,      0.90,
  "Isoniazid | CSF | Rapid",           5.2,      0.65
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = simulated_long,
  reference = published,
  by = "group",
  params = c("auclast", "cmax"),
  tolerance_pct = 25
)
knitr::kable(cmp)

NCA parameter	group	Reference	Simulated	% diff
Cmax	Pyrazinamide \| Plasma \| Day 3	36.3	34.4	-5.3%
Cmax	Pyrazinamide \| Plasma \| Day 28	30.7	31.4	+2.2%
Cmax	Pyrazinamide \| CSF \| Day 3	34	32.6	-4.1%
Cmax	Pyrazinamide \| CSF \| Day 28	28.5	28.8	+1.1%
Cmax	Isoniazid \| Plasma \| Slow	3.7	3.66	-1.0%
Cmax	Isoniazid \| Plasma \| Intermediate	2.75	2.55	-7.4%
Cmax	Isoniazid \| Plasma \| Rapid	2.25	1.86	-17.4%
Cmax	Isoniazid \| CSF \| Slow	1.65	1.6	-2.9%
Cmax	Isoniazid \| CSF \| Intermediate	0.9	0.824	-8.4%
Cmax	Isoniazid \| CSF \| Rapid	0.65	0.512	-21.3%
AUClast	Pyrazinamide \| Plasma \| Day 3	385	388	+0.9%
AUClast	Pyrazinamide \| Plasma \| Day 28	288	287	-0.3%
AUClast	Pyrazinamide \| CSF \| Day 3	405	408	+0.7%
AUClast	Pyrazinamide \| CSF \| Day 28	300	301	+0.5%
AUClast	Isoniazid \| Plasma \| Slow	20.4	20.7	+1.7%
AUClast	Isoniazid \| Plasma \| Intermediate	8	8.7	+8.7%
AUClast	Isoniazid \| Plasma \| Rapid	5.2	4.86	-6.5%
AUClast	Isoniazid \| CSF \| Slow	21.7	21.6	-0.6%
AUClast	Isoniazid \| CSF \| Intermediate	8	9.06	+13.3%
AUClast	Isoniazid \| CSF \| Rapid	5.2	5.07	-2.5%

The pyrazinamide comparisons agree to within 4% on every row in both matrices and at both visits, which is the strongest available check on that model: plasma AUC is a direct function of the tabulated clearance and the CSF AUC follows from PPC. The isoniazid rows agree in ordering and magnitude across the three phenotypes, with two exceptions worth naming rather than tuning away.

The intermediate-acetylator AUC rows run about 25-30% above the digitised medians (the CSF row is flagged at the 25% tolerance). This is a cohort-composition difference, not a parameter-transcription error: the structural gate above reproduces CL/F = fu * CLint exactly for the intermediate phenotype, so the tabulated clearance is being applied correctly. Every subject in the simulated intermediate arm truly is an intermediate acetylator, whereas the paper’s Figure 2 intermediate group also contains participants whose phenotype was imputed by the mixture model - 37% of the cohort had no NAT2 genotype - and that group’s observed box is correspondingly wide, running from roughly 6 to 13.5 mg*h/L around a median of 8. The simulated median of about 10 sits inside that interquartile range.

The rapid-acetylator Cmax rows run about 22% low, which is the same peak-shape difference discussed in the Errata; the corresponding AUC rows agree to within about 10%.

nca_all |>
  dplyr::filter(PPTESTCD == "auclast") |>
  dplyr::group_by(drug, matrix, treatment) |>
  dplyr::summarise(
    median = stats::median(PPORRES),
    p10 = stats::quantile(PPORRES, 0.10),
    p90 = stats::quantile(PPORRES, 0.90),
    .groups = "drop"
  ) |>
  dplyr::rename(
    "Drug" = drug, "Matrix" = matrix, "Group" = treatment,
    "Median AUC0-24 (mg*h/L)" = median,
    "10th percentile" = p10, "90th percentile" = p90
  ) |>
  knitr::kable(digits = 1)

Drug	Matrix	Group	Median AUC0-24 (mg*h/L)	10th percentile	90th percentile
Isoniazid	CSF	Intermediate	9.1	5.6	15.0
Isoniazid	CSF	Rapid	5.1	3.0	8.8
Isoniazid	CSF	Slow	21.6	11.7	36.5
Isoniazid	Plasma	Intermediate	8.7	5.4	14.4
Isoniazid	Plasma	Rapid	4.9	2.9	8.5
Isoniazid	Plasma	Slow	20.7	11.3	35.1
Pyrazinamide	CSF	Day 28	301.5	228.5	430.6
Pyrazinamide	CSF	Day 3	407.9	289.8	579.1
Pyrazinamide	Plasma	Day 28	287.2	217.6	410.1
Pyrazinamide	Plasma	Day 3	388.5	276.1	551.4

The headline finding of the paper - that both drugs reach CSF exposures matching plasma - is recovered as a cohort-level statement. The comparison uses the median and a robust interval rather than the cohort extremes, because the extreme of a random cohort is not reproducible across rxode2 versions.

penetration <- nca_all |>
  dplyr::filter(PPTESTCD == "auclast") |>
  dplyr::select(drug, matrix, treatment, id, PPORRES) |>
  tidyr::pivot_wider(names_from = matrix, values_from = PPORRES) |>
  dplyr::mutate(ratio = CSF / Plasma)

stopifnot(
  # PPC is 1.05 for pyrazinamide and 1.04 for isoniazid; the AUC ratio over a
  # full steady-state interval is exactly PPC for every subject, so this is a
  # tight bound rather than an extreme-of-cohort assertion.
  abs(stats::median(penetration$ratio[penetration$drug == "Pyrazinamide"]) - 1.05) < 0.01,
  abs(stats::median(penetration$ratio[penetration$drug == "Isoniazid"]) - 1.04) < 0.01
)

penetration |>
  dplyr::group_by(drug, treatment) |>
  dplyr::summarise(median_ratio = stats::median(ratio), .groups = "drop") |>
  dplyr::rename("Drug" = drug, "Group" = treatment,
                "Median CSF:plasma AUC ratio" = median_ratio) |>
  knitr::kable(digits = 3)

Drug	Group	Median CSF:plasma AUC ratio
Isoniazid	Intermediate	1.041
Isoniazid	Rapid	1.041
Isoniazid	Slow	1.040
Pyrazinamide	Day 28	1.050
Pyrazinamide	Day 3	1.050

Adjusting for the unbound plasma fractions reproduces the paper’s Discussion statement that the partition coefficients rise to 1.12 and 1.20 once protein binding is accounted for:

unbound_ppc <- c(
  pyrazinamide = 1.05 / 0.933, # measured unbound fraction 93.3% (Results)
  isoniazid = 1.04 / 0.86 # literature unbound fraction 86% (Discussion, ref 20)
)
stopifnot(abs(unbound_ppc - c(1.12, 1.20)) < 0.01)
round(unbound_ppc, 3)
#> pyrazinamide    isoniazid 
#>        1.125        1.209

Assumptions and deviations

Table 2 takes precedence over the supplement’s $THETA / $OMEGA blocks. The supplement supplies the two final-model control streams, but their initial estimate blocks are explicitly headed Initial estimates, whereas Table 2 is headed Final pharmacokinetic parameters estimate. The control streams are therefore used for structure and Table 2 for values. For pyrazinamide the two agree to rounding. For isoniazid the control-stream initials sit roughly 5% away from the tabulated finals throughout (CL 15.4 vs 14.6, V 43.7 vs 43.6, Q 5.1 vs 5.02, MTT 0.255 vs 0.249), consistent with their having been seeded from an earlier run.

Isoniazid between-occasion variability on ka and MTT. This is the one place where the two sources cannot both be right. Table 2 reports BOV of 87.0% on the absorption rate constant and 139% on mean transit time. The S10 control stream binds ETA(20) (BOVKA) to an omega of 1.96 (87.0% would be 0.757) and ETA(25) (BOVMTT) to 0.758 (139% would be 1.93) - that is, the two values are transposed relative to the table, and each matches the other row almost exactly (1.96 gives 140% against the table’s 139%; 0.758 gives 87.1% against the table’s 87.0%). A coincidence of that precision is implausible, so one of the two sources has the assignment crossed. Table 2 is used, per the precedence rule above and because the same table’s pyrazinamide column matches its control stream by name rather than by position. The magnitudes are published either way; only which absorption parameter carries which is at issue, and neither affects a typical-value simulation.

Isoniazid transit-chain length. Table 2 reports 5, fixed, with footnote d explaining that it was fixed “based on the previously estimated value, to improve model stability”; the S10 $THETA block carries 5.82 without a FIX flag. The tabulated value is used. The same footnote states that a sensitivity analysis found the parameter non-critical, which this extraction confirms independently: changing the chain length from 5 to 5.82 moves the simulated isoniazid Cmax by under 0.5% in every phenotype.

Additive residual error is the footnote rule, not the printed number. Table 2 footnote c states that each additive residual standard deviation was fixed to 20% of the matrix-specific LLOQ, and both control streams implement exactly that (ADD = THETA + 0.2 * LLOQ with the THETA fixed to zero). For three of the four rows the rule and the printed value agree: pyrazinamide plasma 0.2 x 0.200 = 0.04, isoniazid plasma 0.2 x 0.105 = 0.021 (printed 0.02), isoniazid CSF 0.2 x 0.0586 = 0.0117 (printed 0.01). The pyrazinamide CSF row prints 0.04, which cannot be 20% of the 0.234 mg/L CSF LLOQ that the same footnote invokes; 0.2 x 0.234 = 0.0468 is used. The printed 0.04 appears to have been carried across from the plasma row.

Hepatic plasma flow reference mass. Table 2 reports Qh as “76.3 Fixed”, while the Methods say it was fixed at 90 L/h and footnote e reconciles the two: 76.3 L/h is the value at the cohort median fat-free mass of 45 kg, and 90 L/h is the value for a 70 kg male whose fat-free mass is 56.1 kg. The model parameter is therefore 90 L/h at a 56.1 kg reference, scaled allometrically, exactly as the S10 control stream writes it (ALLMCL_FFM_HEP = (FFM/56.1)**0.75); 90 x (45/56.1)^0.75 = 76.3 confirms the reading.

Isoniazid unbound fraction. The model fixes fu at 95% (Methods, and S10 $THETA 12). The Discussion separately quotes a literature unbound fraction of 86% when adjusting the partition coefficient. Both numbers are used where the paper uses them: 0.95 inside the model, 0.86 in the unbound-PPC check above.

Isoniazid Cmax against Figure 4. The simulated typical-value plasma peaks (4.11, 3.00 and 2.04 mg/L for slow, intermediate and rapid acetylators at 300 mg) sit 10-28% below the peaks read from the published Figure 4 (roughly 4.58, 3.70 and 2.83), with the gap widening as clearance rises; the CSF peaks agree closely (1.69 against 1.70 for slow acetylators). The plasma AUCs are exactly right, as the well-stirred identity above shows, so the difference is in peak shape rather than in exposure or in any transcribed parameter. Notably the published Figure 4 typical-value peaks also exceed that paper’s own observed individual median Cmax values in Figure 2 (about 3.7, 2.75 and 2.25) by a similar margin, whereas the simulated typical values land on those observed medians. Nothing has been tuned; the parameters are as tabulated.

Between-occasion variability for unobserved doses is not reproduced. For pyrazinamide only, the S9 control stream inflates the BOV on bioavailability, ka and MTT by an estimated factor of 2.51 for records whose preceding dose was not directly observed (IF (OBSERVED.EQ.0)), which the paper introduced to absorb the extra variability of pre-dose samples following an unwitnessed home dose. The factor is inert in forward simulation - every dose in a simulated regimen is specified, so the observed-dose branch always applies - and reproducing it would require the study’s per-record OBSERVED flag. The five-occasion BOV is carried at the published base magnitudes.

The NAT2 mixture model is not reproduced. NAT2 genotype was missing for 37% of participants, and the paper imputed their phenotype with a mixture model whose class probabilities were fixed to the frequencies observed among the genotyped participants. That is an estimation device rather than a structural feature, so the phenotype enters these models as an ordinary covariate pair (NAT2_SLOW / NAT2_RAPID, with the joint zero state denoting intermediate).

Between-subject variability is present only on clearance. Both control streams fix every other $OMEGA to zero, including the etas on volume, ka, bioavailability, MTT, ke0 and PPC, so those etas are omitted here rather than carried as fixed(0).

Simulation dose. The paper states the simulations used 25 mg/kg pyrazinamide and 5 mg/kg isoniazid for the typical individual, without giving the milligram amounts. Dosing in the trial followed WHO weight bands as a fixed-dose combination, and the cohort median weight of 60 kg falls in the 55-70 kg band, which delivers 1600 mg pyrazinamide and 300 mg isoniazid. Those amounts are used here and reproduce Figure 3 to within about 3%.

Absorption is written out rather than delegated to transit(). Both control streams set F1 = 0 so that the whole dose enters through the Savic transit density rather than as a depot bolus. In rxode2 the built-in transit() helper combined with f(depot) <- 0 evaluates to an identically zero input rate, which would silently simulate flat zero concentrations; the closed form is therefore written out with podo() and tad(), which remain live under f(depot) <- 0. The steady-state mass-balance check at the top of this vignette is what confirms the dose is delivered exactly once.

Race and ethnicity. Table 1 reports neither, beyond the South African recruitment sites, so the virtual cohort carries no race covariate. None of the models uses one.