Skip to contents

Salmeterol and alpha-hydroxysalmeterol in plasma and urine (Thoueille 2026)

Source: vignettes/articles/Thoueille_2026_salmeterol.Rmd
Thoueille_2026_salmeterol.Rmd
ui <- rxode2::rxode(readModelDb("Thoueille_2026_salmeterol"))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etalurprod
#> as a work-around try putting the mu-referenced expression on a simple line
mod <- readModelDb("Thoueille_2026_salmeterol")

Model and source

  • Citation: Thoueille P, Danion A, Hostrup M, Petrou M, Deventer K, Buclin T, Girardin FR, Mazzoni I, Rabin O, Guidi M. Pharmacometric-Based Evaluation of Salmeterol and Its Metabolite alpha-Hydroxysalmeterol in Plasma and Urine: Practical Implications for Doping Control. CPT Pharmacometrics Syst Pharmacol. 2026. doi:10.1002/psp4.70187. Final NONMEM control stream and data-formatting description in Data S1.
  • Description: Joint plasma and urine population PK model for inhaled salmeterol and its major metabolite alpha-hydroxysalmeterol in healthy participants, chronic asthmatics and athletes / endurance-trained individuals (Thoueille 2026, WADA doping-control analysis). Intravenous-like bolus absorption into a two-compartment parent disposition; complete irreversible systemic conversion to alpha-hydroxysalmeterol in a metabolite plasma compartment whose volume is constrained equal to the parent central volume; separate cumulative urine compartments for parent and metabolite fed by first-order urinary excretion rate constants; and a cumulative urine-volume compartment driven by a constant estimated urine-production rate that approximates physiologic micturition. Athletes / endurance-trained individuals carry a 63% higher salmeterol plasma clearance and a 191% higher salmeterol urinary excretion rate constant.
  • Article: https://doi.org/10.1002/psp4.70187
  • Supplement (Data S1, final NONMEM control stream and data-formatting description): https://doi.org/10.1002/psp4.70187

Salmeterol is a long-acting beta-2 agonist on the World Anti-Doping Agency (WADA) Prohibited List, permitted by inhalation up to 200 ug in any 24 h period. WADA-accredited laboratories apply a minimum reporting level (MRL) of 10 ng/mL in urine to separate permitted therapeutic use from prohibited use. Thoueille 2026 pooled six studies to build a joint plasma-and-urine population PK model of salmeterol and its major metabolite alpha-hydroxysalmeterol, and used it to ask whether that MRL can actually discriminate the two.

Structurally the model is (Figure 2 and Data S1 $MODEL / $DES):

  • Intravenous-like bolus absorption into a two-compartment parent disposition (central, peripheral1). Maximum plasma concentrations occur within 20 min of inhalation, so no absorption compartment was retained.
  • Systemic conversion of parent into alpha-hydroxysalmeterol at rate kmet (k13), feeding a one-compartment metabolite disposition (central_ohsal) whose volume is constrained equal to the parent central volume for identifiability.
  • First-order urinary excretion of parent (kurine, k14) and metabolite (kurine_ohsal, k35) into cumulative urine compartments (urine, urine_ohsal).
  • A cumulative urine-volume compartment (urine_vol) driven by a constant estimated urine-production rate urprod (UR_PROD), which approximates physiologic micturition for the studies that did not record urine volumes. A urinary concentration is amount / volume accumulated since the last bladder voiding.
cat("ODE states:", paste(ui$state, collapse = ", "), "\n")
#> ODE states: central, peripheral1, central_ohsal, urine, urine_ohsal, urine_vol
cat("Endpoints:", paste(ui$predDf$var, collapse = ", "), "\n")
#> Endpoints: Cc, Cc_ohsal, urineSal, urineOhsal
cat("Covariates required:", paste(ui$allCovs, collapse = ", "), "\n")
#> Covariates required: ATHLETE, USG_CORRECTED, STUDY_SALM

Population

85 individuals from 6 pooled studies contributed the analysis dataset. 92 individuals provided 1175 concentrations (275 and 398 for salmeterol in plasma and urine, 185 and 317 for alpha-hydroxysalmeterol in plasma and urine); after M1 handling of 209 below-quantification-limit values, 966 concentrations from 85 individuals were analysed (Results 3). The median was 6 observations per individual (range 1 to 58).

The pooled studies (Table 1) are Jacobson & Hostrup 2022 (7 healthy + 14 athletes, DPI), Jessen 2021 (11 endurance-trained, DPI, the only source of plasma alpha-hydroxysalmeterol data), Jacobson 2017 (10 healthy, MDI and DPI), Hostrup 2012 (10 healthy + 10 asthmatics, DPI), Deventer 2011 (6 healthy, MDI) and Petrou et al. (24 healthy, DPI, unpublished). Bozzolino 2019 (10 asthmatics under chronic therapy) was excluded for unreliable dose and collection timing. Inhaled doses spanned 50 to 400 ug.

Demographic covariates were not available across the pooled studies, so only two categorical covariates could be screened: individual type and inhalation device. All participants were assumed to have normal renal function.

str(ui$population)
#> List of 10
#>  $ species       : chr "human"
#>  $ n_subjects    : num 85
#>  $ n_studies     : num 6
#>  $ age_range     : chr "Adult (per-study demographics not reported in the pooled analysis)"
#>  $ weight_range  : chr "Not reported"
#>  $ sex_female_pct: num NA
#>  $ disease_state : chr "Pooled healthy participants, chronic asthmatics, and healthy endurance-trained individuals / athletes. All assu"| __truncated__
#>  $ dose_range    : chr "Inhaled salmeterol 50-400 ug: 50 ug and 100 ug MDI, 100 ug / 200 ug / 400 ug DPI, single dose, plus a 200 ug 7-"| __truncated__
#>  $ regions       : chr "Australia, Belgium, Cyprus, Denmark (study locations of the pooled sources)"
#>  $ notes         : chr "Six studies pooled (Table 1): Jacobson & Hostrup 2022 (7 healthy + 14 athletes), Jessen 2021 (11 endurance-trai"| __truncated__

Source trace

Every ini() entry carries an in-file comment naming its origin in inst/modeldb/specificDrugs/Thoueille_2026_salmeterol.R. The table collects them for review. $THETA / $OMEGA / $SIGMA refer to the final-estimate blocks of the Data S1 NONMEM control stream, which carry more significant digits than the rounded Table 2 column.

Equation / parameter Value Source location
lvc (V1/F) 446 L Data S1 $THETA 1; Table 2 V1/F 446 (RSE 17%)
lq (Q/F) 1490 L/h Data S1 $THETA 2; Table 2 Q/F 1490 (RSE 22%)
lvp (V2/F) 871 L Data S1 $THETA 3; Table 2 V2/F 871 (RSE 12%)
lcl (CL_S/F) 193 L/h Data S1 $THETA 4; Table 2 CL_S/F 193 (RSE 10%)
lcl_ohsal (CL_alpha/F) 233 L/h Data S1 $THETA 5; Table 2 CL_alpha/F 233 (RSE 22%)
lkmet (k13) 0.3 1/h Data S1 $THETA 6; Table 2 k13 0.30 (RSE 16%)
lkurine (k14) 0.000943 1/h Data S1 $THETA 7; Table 2 k14 0.00094 (RSE 22%)
lkurine_ohsal (k35) 0.0147 1/h Data S1 $THETA 8; Table 2 k35 0.015 (RSE 22%)
lurprod (UR_PROD) 0.079 L/h Data S1 $THETA 9; Table 2 UR_PROD 0.079 (RSE 18%)
e_athlete_cl 1.63 Data S1 $THETA 10; Table 2 theta Athletes on CL_S (RSE 21%)
e_athlete_kurine 2.91 Data S1 $THETA 11; Table 2 theta Athletes on k14 (RSE 15%)
etalvc 0.0262 Data S1 $OMEGA 1; Table 2 omega V1 16% CV
etalq 0.554 Data S1 $OMEGA 2; Table 2 omega Q 86% CV
etalvp 0.171 Data S1 $OMEGA 3; Table 2 omega V2 43% CV
etalcl 0.103 Data S1 $OMEGA 4; Table 2 omega CL_S 33% CV
etalkmet 0.157 Data S1 $OMEGA 6; Table 2 omega k13 41% CV
etalkurine 0.081 Data S1 $OMEGA 7; Table 2 omega k14 29% CV
etalurprod 0.429 Data S1 $OMEGA 9; Table 2 omega UR_PROD USG-uncorrected 73% CV
expSd sqrt(0.0461) Data S1 $SIGMA BLOCK(2) 1,1; Table 2 sigma plasma salmeterol 22%
expSd_ohsal sqrt(0.145) Data S1 $SIGMA BLOCK(2) 2,2; Table 2 sigma plasma alpha-hydroxysalmeterol 38%
expSdUrineSalLow sqrt(0.0847) Data S1 $SIGMA 3; Table 2 sigma urine salmeterol (low) 29%
expSdUrineSalMid sqrt(0.166) Data S1 $SIGMA BLOCK(2) 1,1; Table 2 sigma urine salmeterol (mid) 41%
expSdUrineSalHigh sqrt(0.322) Data S1 $SIGMA 7; Table 2 sigma urine salmeterol (high) 57%
expSdUrineOhsalMid sqrt(0.144) Data S1 $SIGMA 4; Table 2 sigma urine alpha-hydroxysalmeterol (mid) 38%
expSdUrineOhsalHigh sqrt(0.264) Data S1 $SIGMA BLOCK(2) 2,2; Table 2 sigma urine alpha-hydroxysalmeterol (high) 51%
d/dt(central), d/dt(peripheral1), d/dt(central_ohsal), d/dt(urine), d/dt(urine_ohsal), d/dt(urine_vol) n/a Data S1 $DES DADT(1) to DADT(6)
k10 = cl/vc - kmet - kurine, k30 = cl_ohsal/vc_ohsal - kurine_ohsal n/a Data S1 $PK
vc_ohsal = vc (V3 = V1) n/a Data S1 $PK V3 = V1; Methods 2.2.1 and Discussion (identifiability)
USG-gated IIV on urprod n/a Data S1 $PK IF(USG.EQ.0) / IF(USG.EQ.1) branches, $OMEGA 10 is 0 FIX

Table 2 reports IIV as CV% = sqrt(exp(omega^2) - 1) * 100 and proportional residual errors as CV% = sqrt(sigma^2) * 100 (Table 2 note). The chunk below re-derives the published Table 2 column from the packaged variances, which checks the transcription of every variance term at once.

ini_df <- ui$iniDf
om <- ini_df[!is.na(ini_df$neta1) & ini_df$neta1 == ini_df$neta2, c("name", "est")]
om$`Table 2 CV%` <- round(sqrt(exp(om$est) - 1) * 100)
published_cv <- c(etalvc = 16, etalq = 86, etalvp = 43, etalcl = 33,
                  etalkmet = 41, etalkurine = 29, etalurprod = 73)
om$Published <- published_cv[om$name]
stopifnot(identical(om$`Table 2 CV%`, unname(om$Published)))

sd_names <- c(expSd = 22, expSd_ohsal = 38, expSdUrineSalLow = 29,
              expSdUrineSalMid = 41, expSdUrineSalHigh = 57,
              expSdUrineOhsalMid = 38, expSdUrineOhsalHigh = 51)
sig <- ini_df[ini_df$name %in% names(sd_names), c("name", "est")]
sig$`Table 2 CV%` <- round(sig$est * 100)
sig$Published <- unname(sd_names[sig$name])

# Six of the seven residual SDs reproduce the Table 2 column exactly. The
# seventh, the plasma salmeterol residual, does not: Data S1 `$SIGMA` gives
# 0.0461, whose square root is 21.5%, while Table 2 prints 22%. The
# assertion pins that to exactly one row differing by exactly one CV unit,
# so a genuine transcription error (which would move a row by far more)
# still fails the check. See Assumptions and deviations.
sig_diff <- sig$`Table 2 CV%` - sig$Published
stopifnot(sum(sig_diff != 0) == 1L,
          max(abs(sig_diff)) == 1,
          sig$name[sig_diff != 0] == "expSd")

dplyr::bind_rows(
  om |> dplyr::rename("Variance / SD" = est),
  sig |> dplyr::rename("Variance / SD" = est)
) |>
  dplyr::rename("Parameter" = name, "Derived CV%" = `Table 2 CV%`,
                "Table 2 CV%" = Published) |>
  knitr::kable(digits = 5,
               caption = "Every IIV variance and residual SD re-derived back to the Table 2 CV% column.")
Every IIV variance and residual SD re-derived back to the Table 2 CV% column.
Parameter Variance / SD Derived CV% Table 2 CV%
etalvc 0.02620 16 16
etalq 0.55400 86 86
etalvp 0.17100 43 43
etalcl 0.10300 33 33
etalkmet 0.15700 41 41
etalkurine 0.08100 29 29
etalurprod 0.42900 73 73
expSd 0.21471 21 22
expSd_ohsal 0.38079 38 38
expSdUrineSalLow 0.29103 29 29
expSdUrineSalMid 0.40743 41 41
expSdUrineSalHigh 0.56745 57 57
expSdUrineOhsalMid 0.37947 38 38
expSdUrineOhsalHigh 0.51381 51 51

Units

The analysis was run in molar units: doses were converted to nmol and concentrations to nmol/L (Methods 2.2), and the packaged model keeps that scale. The published tables report ng/mL, so the vignette converts with the molecular weights below. The salmeterol molecular weight is pinned by the worked examples in Data S1, which give AMT = 240.63 nmol for a 100 ug dose and AMT = 962.52 nmol for a 400 ug dose.

MW_SAL <- 415.57  # g/mol, salmeterol free base (C25H37NO4)
MW_OH  <- 431.57  # g/mol, alpha-hydroxysalmeterol (salmeterol + one oxygen)

nmol_per_ug <- 1 / (MW_SAL / 1000)
stopifnot(
  abs(100 * nmol_per_ug - 240.63) < 0.01,   # Data S1 scenario 2 worked example
  abs(400 * nmol_per_ug - 962.52) < 0.02    # Data S1 scenario 1 worked example
)
c(`100 ug (nmol)` = 100 * nmol_per_ug, `400 ug (nmol)` = 400 * nmol_per_ug)
#> 100 ug (nmol) 400 ug (nmol) 
#>      240.6333      962.5334

Structural checks

The paper reports two published quantities that the structural parameters determine exactly, with no simulation required.

Biexponential half-lives

The Discussion reports “average t1/2 for plasma salmeterol of 0.13 and 5 h, corresponding to the initial decline due to diffusion into the peripheral compartment and the onset of elimination, and a final decline due to elimination after reaching equilibrium”. Those are the two roots of the two-compartment characteristic equation at the typical parameter values.

th <- setNames(ui$iniDf$est, ui$iniDf$name)
V1 <- exp(th[["lvc"]]); V2 <- exp(th[["lvp"]])
Q  <- exp(th[["lq"]]);  CL <- exp(th[["lcl"]])
k12 <- Q / V1; k21 <- Q / V2; kel <- CL / V1
b <- k12 + k21 + kel
lambda <- c((b + sqrt(b^2 - 4 * k21 * kel)) / 2, (b - sqrt(b^2 - 4 * k21 * kel)) / 2)
t_half <- log(2) / lambda

data.frame(
  Phase = c("alpha (distribution)", "beta (terminal)"),
  Derived = round(t_half, 4),
  Published = c(0.13, 5)
) |>
  dplyr::rename("Half-life (h)" = Derived, "Thoueille 2026 Discussion (h)" = Published) |>
  knitr::kable(caption = "Biexponential half-lives derived from V1/F, V2/F, Q/F and CL_S/F.")
Biexponential half-lives derived from V1/F, V2/F, Q/F and CL_S/F.
Phase Half-life (h) Thoueille 2026 Discussion (h)
alpha (distribution) 0.1297 0.13
beta (terminal) 5.0054 5.00 

# The two sides are the same closed form evaluated on the same numbers, so
# this is an exact transcription check, not a stochastic comparison.
stopifnot(abs(t_half[1] - 0.13) < 0.005, abs(t_half[2] - 5) < 0.02)

Route fractions and urinary recovery

At typical values the model routes a fixed fraction of parent clearance to each elimination pathway. These closed forms follow from integrating the ODE system to infinity.

kmet_t   <- exp(th[["lkmet"]])
kurine_t <- exp(th[["lkurine"]])
kur_oh_t <- exp(th[["lkurine_ohsal"]])
CLM      <- exp(th[["lcl_ohsal"]])

frac_metab      <- kmet_t * V1 / CL                       # parent CL routed to metabolite
frac_urine_par  <- kurine_t * V1 / CL                      # dose recovered as parent in urine
frac_urine_meta <- kur_oh_t * kmet_t * V1^2 / (CL * CLM)   # dose recovered as metabolite in urine
frac_other      <- 1 - frac_metab - frac_urine_par         # residual k10 route

data.frame(
  Route = c("Systemic conversion to alpha-hydroxysalmeterol (k13)",
            "Direct urinary excretion of parent (k14)",
            "Residual non-specified elimination (k10)",
            "Nominal dose recovered as parent in urine",
            "Nominal dose recovered as metabolite in urine"),
  Fraction = round(100 * c(frac_metab, frac_urine_par, frac_other,
                           frac_urine_par, frac_urine_meta), 3)
) |>
  dplyr::rename("Percent" = Fraction) |>
  knitr::kable(caption = "Typical-value route fractions implied by the final parameter estimates.")
Typical-value route fractions implied by the final parameter estimates.
Route Percent
Systemic conversion to alpha-hydroxysalmeterol (k13) 69.326
Direct urinary excretion of parent (k14) 0.218
Residual non-specified elimination (k10) 30.456
Nominal dose recovered as parent in urine 0.218
Nominal dose recovered as metabolite in urine 1.951 

# The Introduction states parent salmeterol accounts for less than 5% of the
# dose excreted in urine, and the Discussion notes higher overall urine
# concentrations of the metabolite than of the parent.
stopifnot(
  frac_urine_par < 0.05,
  frac_urine_meta > frac_urine_par,
  frac_other > 0   # k10 is non-negative at the typical values
)

Note that frac_metab is 69%, not 100%: the abstract’s phrase “a complete parent conversion into alpha-hydroxysalmeterol” describes the modelling assumption that the conversion is irreversible, but the control stream retains a residual elimination route k10 = CL_S/V1 - k13 - k14 that carries the remaining 30% at the typical values. See Assumptions and deviations.

Virtual cohort and dosing regimens

Thoueille 2026 simulated five one-week regimens (Methods 2.2.3): the maximum authorised doses of 100 ug at 8 h and 16 h, 100 ug every 12 h (8 h and 20 h), and 200 ug once daily at 8 h; and the prohibited doses of 200 ug at 8 h and 16 h, and 200 ug every 12 h (8 h and 20 h). Bladder voiding was assumed every 4 h and prior to each inhalation, and urine volumes were generated from the constant urine-production rate. Predictions are reported for athletes / endurance-trained individuals and for healthy participants, on the USG-corrected branch.

A bladder voiding is an evid = 5 replacement record setting urine, urine_ohsal and urine_vol back to zero, which is exactly the “compartment closed and reopened” pattern the Data S1 data-formatting section describes.

DAYS       <- 7
VOID_EVERY <- 4
STATES     <- c("urine", "urine_ohsal", "urine_vol")

regimens <- tibble::tribble(
  ~regimen,                ~dose_times, ~dose_ug, ~status,
  "100 ug at 8/16 h",      c(8, 16),    100,      "Permitted",
  "100 ug at 8/20 h",      c(8, 20),    100,      "Permitted",
  "200 ug at 8 h",         c(8),        200,      "Permitted",
  "200 ug at 8/16 h",      c(8, 16),    200,      "Prohibited",
  "200 ug at 8/20 h",      c(8, 20),    200,      "Prohibited"
)

# Build one arm (regimen x population) as a self-contained event table.
# `id_offset` keeps subject IDs disjoint across arms: rxSolve treats id as
# the subject key, so colliding ids would silently merge subjects.
make_arm <- function(n, dose_times, dose_ug, athlete, obs_offsets,
                     obs_dvid = 1L, id_offset = 0L) {
  dose_t <- sort(as.vector(outer(dose_times, 24 * (0:(DAYS - 1)), "+")))
  last_dose <- max(dose_t)
  voids <- seq(0, 24 * DAYS, by = VOID_EVERY)

  per_subject <- dplyr::bind_rows(
    # bladder voiding: reset all three urine states
    tidyr::expand_grid(time = voids, cmt = STATES) |>
      dplyr::mutate(amt = 0, evid = 5L, dvid = NA_integer_),
    # inhaled doses, as an intravenous-like bolus into `central`
    tibble::tibble(time = dose_t, amt = dose_ug * nmol_per_ug,
                   evid = 1L, cmt = "central", dvid = NA_integer_),
    # observations: dvid alone (see the note below), one row per requested
    # endpoint so that each endpoint's residual-error draw is returned
    tidyr::expand_grid(time = last_dose + obs_offsets, dvid = obs_dvid) |>
      dplyr::mutate(amt = NA_real_, evid = 0L, cmt = NA_character_)
  ) |>
    # voiding precedes dosing precedes observation at any tied time
    dplyr::arrange(time, match(evid, c(5L, 1L, 0L)), dvid)

  tidyr::expand_grid(id = id_offset + seq_len(n), per_subject) |>
    dplyr::mutate(ATHLETE = athlete, USG_CORRECTED = 1, STUDY_SALM = 2,
                  last_dose = last_dose)
}

This model declares four endpoints (Cc, Cc_ohsal, urineSal, urineOhsal), so rxode2 builds a dvid-to-compartment map onto endpoint pseudo-compartments placed after the six ODE states. Observation rows therefore carry dvid and leave cmt as NA; every observable is returned as a column on every observation row regardless of which dvid the row names. useLinCmt = FALSE is passed to each solve because rxode2’s default ODE-to-linCmt auto-conversion corrupts the dvid mapping for multi-output models of this shape.

Typical-value replication of Tables 3 and S1

Table 3 (salmeterol and alpha-hydroxysalmeterol) and Table S1 (alpha-hydroxysalmeterol at 2, 3 and 4 h post-dose) report predicted urine concentration percentiles for each regimen and population. Every parameter enters the urinary concentration multiplicatively through a log-normal random effect, so the median of the simulated population is the typical-value prediction. Computing it with zeroRe() therefore reproduces the published median column deterministically, with no Monte Carlo noise and no dependence on the random seed.

“Directly after inhalation” is the limit of urine / urine_vol as time approaches the dose from above: the bladder has just been voided, so both the excreted amount and the accumulated volume start from zero and their ratio tends to kurine * A_central / urprod. A real sample cannot be taken at zero volume, so the vignette evaluates it 0.01 h (36 s) after the inhalation, which is the sense of the paper’s own gloss (“small volumes of urine immediately after inhalation resulting from the constant production of urine after voiding the bladder”).

OBS_OFFSETS <- c(0.01, 0.5, 2, 3, 3.999)  # h after the last dose
OBS_LABELS  <- c("Directly after inhalation", "30 min post-dose",
                 "2 h post-dose", "3 h post-dose", "4 h post-dose")

arms <- tidyr::expand_grid(regimens, athlete = c(1, 0)) |>
  dplyr::mutate(population = ifelse(athlete == 1,
                                    "Athletes/endurance-trained", "Healthy participants"),
                arm = dplyr::row_number())

typ_events <- dplyr::bind_rows(lapply(seq_len(nrow(arms)), function(i) {
  make_arm(n = 1, dose_times = arms$dose_times[[i]], dose_ug = arms$dose_ug[i],
           athlete = arms$athlete[i], obs_offsets = OBS_OFFSETS,
           id_offset = arms$arm[i] - 1L) |>
    dplyr::mutate(regimen = arms$regimen[i], population = arms$population[i])
}))
stopifnot(!anyDuplicated(unique(typ_events[, c("id", "time", "evid", "cmt")])))

typ <- rxode2::rxSolve(rxode2::zeroRe(mod), events = typ_events,
                       keep = c("regimen", "population", "last_dose"),
                       useLinCmt = FALSE, addDosing = FALSE) |>
  as.data.frame() |>
  dplyr::filter(!is.na(urineSal)) |>
  dplyr::mutate(
    offset  = round(time - last_dose, 3),
    timing  = OBS_LABELS[match(offset, round(OBS_OFFSETS, 3))],
    sal_ng  = urineSal   * MW_SAL / 1000,
    ohsal_ng = urineOhsal * MW_OH / 1000
  )
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etalurprod
#> 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: etalurprod
#> as a work-around try putting the mu-referenced expression on a simple line
#> ℹ omega/sigma items treated as zero: 'etalvc', 'etalq', 'etalvp', 'etalcl', 'etalkmet', 'etalkurine', 'etalurprod'
#> Warning: multi-subject simulation without without 'omega'

Table 3 – salmeterol 

pub_t3_sal <- tibble::tribble(
  ~population,                  ~regimen,           ~timing,                     ~published,
  "Athletes/endurance-trained", "100 ug at 8/16 h", "Directly after inhalation", 3.6,
  "Athletes/endurance-trained", "100 ug at 8/16 h", "30 min post-dose",          2.0,
  "Athletes/endurance-trained", "100 ug at 8/20 h", "Directly after inhalation", 3.6,
  "Athletes/endurance-trained", "100 ug at 8/20 h", "30 min post-dose",          1.9,
  "Athletes/endurance-trained", "200 ug at 8 h",    "Directly after inhalation", 6.9,
  "Athletes/endurance-trained", "200 ug at 8 h",    "30 min post-dose",          3.5,
  "Athletes/endurance-trained", "200 ug at 8/16 h", "Directly after inhalation", 7.3,
  "Athletes/endurance-trained", "200 ug at 8/16 h", "30 min post-dose",          3.9,
  "Athletes/endurance-trained", "200 ug at 8/20 h", "Directly after inhalation", 7.2,
  "Athletes/endurance-trained", "200 ug at 8/20 h", "30 min post-dose",          3.8,
  "Healthy participants",       "100 ug at 8/16 h", "Directly after inhalation", 1.3,
  "Healthy participants",       "100 ug at 8/16 h", "30 min post-dose",          0.8,
  "Healthy participants",       "100 ug at 8/20 h", "Directly after inhalation", 1.3,
  "Healthy participants",       "100 ug at 8/20 h", "30 min post-dose",          0.7,
  "Healthy participants",       "200 ug at 8 h",    "Directly after inhalation", 2.4,
  "Healthy participants",       "200 ug at 8 h",    "30 min post-dose",          1.3,
  "Healthy participants",       "200 ug at 8/16 h", "Directly after inhalation", 2.7,
  "Healthy participants",       "200 ug at 8/16 h", "30 min post-dose",          1.6,
  "Healthy participants",       "200 ug at 8/20 h", "Directly after inhalation", 2.6,
  "Healthy participants",       "200 ug at 8/20 h", "30 min post-dose",          1.5
)

cmp_sal <- pub_t3_sal |>
  dplyr::left_join(typ |> dplyr::select(population, regimen, timing, sal_ng),
                   by = c("population", "regimen", "timing")) |>
  dplyr::mutate(pct_diff = 100 * (sal_ng - published) / published)
stopifnot(nrow(cmp_sal) == 20L, !anyNA(cmp_sal$sal_ng))

cmp_sal |>
  dplyr::mutate(dplyr::across(c(sal_ng, pct_diff), \(x) round(x, 2))) |>
  dplyr::rename("Population" = population, "Regimen" = regimen, "Timing" = timing,
                "Thoueille 2026 median (ng/mL)" = published,
                "Model typical value (ng/mL)" = sal_ng, "% diff" = pct_diff) |>
  knitr::kable(caption = "Table 3 of Thoueille 2026, salmeterol urine concentration: published median vs packaged-model typical value.")
Table 3 of Thoueille 2026, salmeterol urine concentration: published median vs packaged-model typical value.
Population Regimen Timing Thoueille 2026 median (ng/mL) Model typical value (ng/mL) % diff
Athletes/endurance-trained 100 ug at 8/16 h Directly after inhalation 3.6 3.58 -0.52
Athletes/endurance-trained 100 ug at 8/16 h 30 min post-dose 2.0 1.94 -3.17
Athletes/endurance-trained 100 ug at 8/20 h Directly after inhalation 3.6 3.48 -3.29
Athletes/endurance-trained 100 ug at 8/20 h 30 min post-dose 1.9 1.84 -3.05
Athletes/endurance-trained 200 ug at 8 h Directly after inhalation 6.9 6.82 -1.16
Athletes/endurance-trained 200 ug at 8 h 30 min post-dose 3.5 3.55 1.38
Athletes/endurance-trained 200 ug at 8/16 h Directly after inhalation 7.3 7.16 -1.89
Athletes/endurance-trained 200 ug at 8/16 h 30 min post-dose 3.9 3.87 -0.69
Athletes/endurance-trained 200 ug at 8/20 h Directly after inhalation 7.2 6.96 -3.29
Athletes/endurance-trained 200 ug at 8/20 h 30 min post-dose 3.8 3.68 -3.05
Healthy participants 100 ug at 8/16 h Directly after inhalation 1.3 1.31 0.65
Healthy participants 100 ug at 8/16 h 30 min post-dose 0.8 0.77 -3.66
Healthy participants 100 ug at 8/20 h Directly after inhalation 1.3 1.26 -3.39
Healthy participants 100 ug at 8/20 h 30 min post-dose 0.7 0.72 2.86
Healthy participants 200 ug at 8 h Directly after inhalation 2.4 2.37 -1.25
Healthy participants 200 ug at 8 h 30 min post-dose 1.3 1.30 0.22
Healthy participants 200 ug at 8/16 h Directly after inhalation 2.7 2.62 -3.08
Healthy participants 200 ug at 8/16 h 30 min post-dose 1.6 1.54 -3.66
Healthy participants 200 ug at 8/20 h Directly after inhalation 2.6 2.51 -3.39
Healthy participants 200 ug at 8/20 h 30 min post-dose 1.5 1.44 -4.00 

# Structural: a mis-transcribed clearance, rate constant, dose or unit
# conversion moves the whole set by tens of percent. The two sides are
# deterministic, so no Monte Carlo tolerance is needed -- the residual
# spread below is the published values' one-decimal rounding plus the
# finite 0.01 h offset used for the "directly after" limit.
stopifnot(
  max(abs(cmp_sal$pct_diff)) < 8,
  median(abs(cmp_sal$pct_diff)) < 5
)

Table 3 – alpha-hydroxysalmeterol

The metabolite column of Table 3 spans values as small as 0.1 ng/mL, where the published one-decimal rounding is itself up to 50% of the value, so the comparison uses an absolute-plus-relative tolerance rather than a pure percentage.

pub_t3_oh <- tibble::tribble(
  ~population,                  ~regimen,           ~timing,                     ~published,
  "Athletes/endurance-trained", "100 ug at 8/16 h", "Directly after inhalation", 0.9,
  "Athletes/endurance-trained", "100 ug at 8/16 h", "30 min post-dose",          1.7,
  "Athletes/endurance-trained", "100 ug at 8/20 h", "Directly after inhalation", 0.4,
  "Athletes/endurance-trained", "100 ug at 8/20 h", "30 min post-dose",          1.2,
  "Athletes/endurance-trained", "200 ug at 8 h",    "Directly after inhalation", 0.1,
  "Athletes/endurance-trained", "200 ug at 8 h",    "30 min post-dose",          1.7,
  "Athletes/endurance-trained", "200 ug at 8/16 h", "Directly after inhalation", 1.7,
  "Athletes/endurance-trained", "200 ug at 8/16 h", "30 min post-dose",          3.3,
  "Athletes/endurance-trained", "200 ug at 8/20 h", "Directly after inhalation", 0.8,
  "Athletes/endurance-trained", "200 ug at 8/20 h", "30 min post-dose",          2.5,
  "Healthy participants",       "100 ug at 8/16 h", "Directly after inhalation", 1.6,
  "Healthy participants",       "100 ug at 8/16 h", "30 min post-dose",          2.5,
  "Healthy participants",       "100 ug at 8/20 h", "Directly after inhalation", 1.0,
  "Healthy participants",       "100 ug at 8/20 h", "30 min post-dose",          1.9,
  "Healthy participants",       "200 ug at 8 h",    "Directly after inhalation", 0.3,
  "Healthy participants",       "200 ug at 8 h",    "30 min post-dose",          2.1,
  "Healthy participants",       "200 ug at 8/16 h", "Directly after inhalation", 3.3,
  "Healthy participants",       "200 ug at 8/16 h", "30 min post-dose",          4.9,
  "Healthy participants",       "200 ug at 8/20 h", "Directly after inhalation", 2.1,
  "Healthy participants",       "200 ug at 8/20 h", "30 min post-dose",          3.8
)

cmp_oh <- pub_t3_oh |>
  dplyr::left_join(typ |> dplyr::select(population, regimen, timing, ohsal_ng),
                   by = c("population", "regimen", "timing")) |>
  dplyr::mutate(abs_diff = abs(ohsal_ng - published),
                tol = pmax(0.1, 0.15 * published))
stopifnot(nrow(cmp_oh) == 20L, !anyNA(cmp_oh$ohsal_ng))

cmp_oh |>
  dplyr::mutate(dplyr::across(c(ohsal_ng, abs_diff, tol), \(x) round(x, 2))) |>
  dplyr::rename("Population" = population, "Regimen" = regimen, "Timing" = timing,
                "Thoueille 2026 median (ng/mL)" = published,
                "Model typical value (ng/mL)" = ohsal_ng,
                "|difference|" = abs_diff, "Tolerance" = tol) |>
  knitr::kable(caption = "Table 3 of Thoueille 2026, alpha-hydroxysalmeterol urine concentration.")
Table 3 of Thoueille 2026, alpha-hydroxysalmeterol urine concentration.
Population Regimen Timing Thoueille 2026 median (ng/mL) Model typical value (ng/mL) |difference| Tolerance
Athletes/endurance-trained 100 ug at 8/16 h Directly after inhalation 0.9 0.93 0.03 0.14
Athletes/endurance-trained 100 ug at 8/16 h 30 min post-dose 1.7 1.67 0.03 0.26
Athletes/endurance-trained 100 ug at 8/20 h Directly after inhalation 0.4 0.44 0.04 0.10
Athletes/endurance-trained 100 ug at 8/20 h 30 min post-dose 1.2 1.20 0.00 0.18
Athletes/endurance-trained 200 ug at 8 h Directly after inhalation 0.1 0.11 0.01 0.10
Athletes/endurance-trained 200 ug at 8 h 30 min post-dose 1.7 1.67 0.03 0.26
Athletes/endurance-trained 200 ug at 8/16 h Directly after inhalation 1.7 1.85 0.15 0.26
Athletes/endurance-trained 200 ug at 8/16 h 30 min post-dose 3.3 3.33 0.03 0.49
Athletes/endurance-trained 200 ug at 8/20 h Directly after inhalation 0.8 0.88 0.08 0.12
Athletes/endurance-trained 200 ug at 8/20 h 30 min post-dose 2.5 2.40 0.10 0.38
Healthy participants 100 ug at 8/16 h Directly after inhalation 1.6 1.70 0.10 0.24
Healthy participants 100 ug at 8/16 h 30 min post-dose 2.5 2.46 0.04 0.38
Healthy participants 100 ug at 8/20 h Directly after inhalation 1.0 1.09 0.09 0.15
Healthy participants 100 ug at 8/20 h 30 min post-dose 1.9 1.86 0.04 0.28
Healthy participants 200 ug at 8 h Directly after inhalation 0.3 0.40 0.10 0.10
Healthy participants 200 ug at 8 h 30 min post-dose 2.1 2.01 0.09 0.32
Healthy participants 200 ug at 8/16 h Directly after inhalation 3.3 3.41 0.11 0.49
Healthy participants 200 ug at 8/16 h 30 min post-dose 4.9 4.92 0.02 0.74
Healthy participants 200 ug at 8/20 h Directly after inhalation 2.1 2.18 0.08 0.32
Healthy participants 200 ug at 8/20 h 30 min post-dose 3.8 3.73 0.07 0.57 

stopifnot(all(cmp_oh$abs_diff <= cmp_oh$tol))

Table S1 – alpha-hydroxysalmeterol at 2, 3 and 4 h post-dose

Table S1 samples the metabolite where its urinary concentration is largest and flattest, so the published values there are well clear of the rounding floor and support a percentage comparison. The 4 h sample is taken at 3.999 h because a bladder voiding occurs exactly at 4 h.

pub_s1 <- tibble::tribble(
  ~population,                  ~regimen,           ~timing,          ~published,
  "Athletes/endurance-trained", "100 ug at 8/16 h", "2 h post-dose",  2.2,
  "Athletes/endurance-trained", "100 ug at 8/16 h", "3 h post-dose",  2.3,
  "Athletes/endurance-trained", "100 ug at 8/16 h", "4 h post-dose",  2.3,
  "Athletes/endurance-trained", "100 ug at 8/20 h", "2 h post-dose",  1.9,
  "Athletes/endurance-trained", "100 ug at 8/20 h", "3 h post-dose",  2.0,
  "Athletes/endurance-trained", "100 ug at 8/20 h", "4 h post-dose",  2.0,
  "Athletes/endurance-trained", "200 ug at 8 h",    "2 h post-dose",  3.0,
  "Athletes/endurance-trained", "200 ug at 8 h",    "3 h post-dose",  3.3,
  "Athletes/endurance-trained", "200 ug at 8 h",    "4 h post-dose",  3.3,
  "Athletes/endurance-trained", "200 ug at 8/16 h", "2 h post-dose",  4.5,
  "Athletes/endurance-trained", "200 ug at 8/16 h", "3 h post-dose",  4.6,
  "Athletes/endurance-trained", "200 ug at 8/16 h", "4 h post-dose",  4.5,
  "Athletes/endurance-trained", "200 ug at 8/20 h", "2 h post-dose",  3.8,
  "Athletes/endurance-trained", "200 ug at 8/20 h", "3 h post-dose",  3.9,
  "Athletes/endurance-trained", "200 ug at 8/20 h", "4 h post-dose",  3.9,
  "Healthy participants",       "100 ug at 8/16 h", "2 h post-dose",  3.1,
  "Healthy participants",       "100 ug at 8/16 h", "3 h post-dose",  3.2,
  "Healthy participants",       "100 ug at 8/16 h", "4 h post-dose",  3.2,
  "Healthy participants",       "100 ug at 8/20 h", "2 h post-dose",  2.6,
  "Healthy participants",       "100 ug at 8/20 h", "3 h post-dose",  2.7,
  "Healthy participants",       "100 ug at 8/20 h", "4 h post-dose",  2.8,
  "Healthy participants",       "200 ug at 8 h",    "2 h post-dose",  3.7,
  "Healthy participants",       "200 ug at 8 h",    "3 h post-dose",  4.0,
  "Healthy participants",       "200 ug at 8 h",    "4 h post-dose",  4.2,
  "Healthy participants",       "200 ug at 8/16 h", "2 h post-dose",  6.3,
  "Healthy participants",       "200 ug at 8/16 h", "3 h post-dose",  6.5,
  "Healthy participants",       "200 ug at 8/16 h", "4 h post-dose",  6.5,
  "Healthy participants",       "200 ug at 8/20 h", "2 h post-dose",  5.3,
  "Healthy participants",       "200 ug at 8/20 h", "3 h post-dose",  5.6,
  "Healthy participants",       "200 ug at 8/20 h", "4 h post-dose",  5.6
)

cmp_s1 <- pub_s1 |>
  dplyr::left_join(typ |> dplyr::select(population, regimen, timing, ohsal_ng),
                   by = c("population", "regimen", "timing")) |>
  dplyr::mutate(pct_diff = 100 * (ohsal_ng - published) / published)
stopifnot(nrow(cmp_s1) == 30L, !anyNA(cmp_s1$ohsal_ng))

cmp_s1 |>
  dplyr::mutate(dplyr::across(c(ohsal_ng, pct_diff), \(x) round(x, 2))) |>
  dplyr::rename("Population" = population, "Regimen" = regimen, "Timing" = timing,
                "Thoueille 2026 median (ng/mL)" = published,
                "Model typical value (ng/mL)" = ohsal_ng, "% diff" = pct_diff) |>
  knitr::kable(caption = "Table S1 of Thoueille 2026: alpha-hydroxysalmeterol urine concentration 2 to 4 h post-dose.")
Table S1 of Thoueille 2026: alpha-hydroxysalmeterol urine concentration 2 to 4 h post-dose.
Population Regimen Timing Thoueille 2026 median (ng/mL) Model typical value (ng/mL) % diff
Athletes/endurance-trained 100 ug at 8/16 h 2 h post-dose 2.2 2.20 -0.06
Athletes/endurance-trained 100 ug at 8/16 h 3 h post-dose 2.3 2.27 -1.31
Athletes/endurance-trained 100 ug at 8/16 h 4 h post-dose 2.3 2.25 -2.01
Athletes/endurance-trained 100 ug at 8/20 h 2 h post-dose 1.9 1.79 -5.63
Athletes/endurance-trained 100 ug at 8/20 h 3 h post-dose 2.0 1.90 -5.04
Athletes/endurance-trained 100 ug at 8/20 h 4 h post-dose 2.0 1.91 -4.29
Athletes/endurance-trained 200 ug at 8 h 2 h post-dose 3.0 2.96 -1.35
Athletes/endurance-trained 200 ug at 8 h 3 h post-dose 3.3 3.23 -2.13
Athletes/endurance-trained 200 ug at 8 h 4 h post-dose 3.3 3.31 0.32
Athletes/endurance-trained 200 ug at 8/16 h 2 h post-dose 4.5 4.40 -2.28
Athletes/endurance-trained 200 ug at 8/16 h 3 h post-dose 4.6 4.54 -1.31
Athletes/endurance-trained 200 ug at 8/16 h 4 h post-dose 4.5 4.51 0.16
Athletes/endurance-trained 200 ug at 8/20 h 2 h post-dose 3.8 3.59 -5.63
Athletes/endurance-trained 200 ug at 8/20 h 3 h post-dose 3.9 3.80 -2.61
Athletes/endurance-trained 200 ug at 8/20 h 4 h post-dose 3.9 3.83 -1.83
Healthy participants 100 ug at 8/16 h 2 h post-dose 3.1 3.07 -0.86
Healthy participants 100 ug at 8/16 h 3 h post-dose 3.2 3.19 -0.36
Healthy participants 100 ug at 8/16 h 4 h post-dose 3.2 3.21 0.23
Healthy participants 100 ug at 8/20 h 2 h post-dose 2.6 2.52 -2.92
Healthy participants 100 ug at 8/20 h 3 h post-dose 2.7 2.67 -1.10
Healthy participants 100 ug at 8/20 h 4 h post-dose 2.8 2.72 -2.92
Healthy participants 200 ug at 8 h 2 h post-dose 3.7 3.49 -5.71
Healthy participants 200 ug at 8 h 3 h post-dose 4.0 3.88 -3.04
Healthy participants 200 ug at 8 h 4 h post-dose 4.2 4.06 -3.24
Healthy participants 200 ug at 8/16 h 2 h post-dose 6.3 6.15 -2.43
Healthy participants 200 ug at 8/16 h 3 h post-dose 6.5 6.38 -1.89
Healthy participants 200 ug at 8/16 h 4 h post-dose 6.5 6.41 -1.31
Healthy participants 200 ug at 8/20 h 2 h post-dose 5.3 5.05 -4.75
Healthy participants 200 ug at 8/20 h 3 h post-dose 5.6 5.34 -4.64
Healthy participants 200 ug at 8/20 h 4 h post-dose 5.6 5.44 -2.92 

stopifnot(
  max(abs(cmp_s1$pct_diff)) < 10,
  median(abs(cmp_s1$pct_diff)) < 5
)

Across all 70 published median values in Tables 3 and S1 the packaged model reproduces the paper to a median absolute deviation of 2.9%.

Stochastic simulation, Figure 3 and the Table 3 prediction intervals

Figure 3 and the percentile columns of Table 3 come from Monte Carlo simulations “accounting for both IIV and RUV” (Methods 2.2.3). The between-subject part alone is not enough: the urineSal / urineOhsal columns rxode2 returns are individual predictions, and the residual draw for each endpoint lives in the sim column of the observation row carrying that endpoint’s dvid. The event table below therefore emits one observation row per urine endpoint per time.

Residual magnitude depends on STUDY_SALM, which selects the Table 2 low / mid / high tier. The paper does not state which tier it simulated with. STUDY_SALM = 2 (Jessen 2021) is used here, and the percentile comparison below shows it is the choice that reproduces the published intervals – which is also the natural one, Jessen 2021 being the only pooled study contributing plasma metabolite data.

N_PER_ARM <- 200L
rxode2::rxSetSeed(20260825)
set.seed(20260825)

# Endpoint index (dvid) and the compartment number rxode2 reports it under.
endpoints <- ui$predDf[, c("var", "cmt")]
endpoints$dvid <- seq_len(nrow(endpoints))
urine_ep <- endpoints[endpoints$var %in% c("urineSal", "urineOhsal"), ]
stopifnot(nrow(urine_ep) == 2L)

# Sample the last inter-void window on a fine grid, offset from the void
# boundaries where urine_vol is exactly zero.
grid_offsets <- sort(unique(c(0.01, 0.5, seq(0.25, 3.75, by = 0.25))))

sto_events <- dplyr::bind_rows(lapply(seq_len(nrow(arms)), function(i) {
  make_arm(n = N_PER_ARM, dose_times = arms$dose_times[[i]], dose_ug = arms$dose_ug[i],
           athlete = arms$athlete[i], obs_offsets = grid_offsets,
           obs_dvid = urine_ep$dvid,
           id_offset = (arms$arm[i] - 1L) * N_PER_ARM) |>
    dplyr::mutate(regimen = arms$regimen[i], population = arms$population[i],
                  status = arms$status[i])
}))
stopifnot(!anyDuplicated(unique(sto_events[, c("id", "time", "evid", "cmt", "dvid")])))

sto <- rxode2::rxSolve(mod, events = sto_events,
                       keep = c("regimen", "population", "status", "last_dose"),
                       useLinCmt = FALSE, addDosing = FALSE) |>
  as.data.frame() |>
  dplyr::filter(!is.na(sim)) |>
  dplyr::mutate(
    offset  = round(time - last_dose, 3),
    analyte = ifelse(CMT == urine_ep$cmt[urine_ep$var == "urineSal"],
                     "Salmeterol", "alpha-hydroxysalmeterol"),
    # `sim` carries the endpoint's residual draw in nmol/L; convert with the
    # molecular weight of whichever analyte the row belongs to.
    conc_ng = sim * ifelse(analyte == "Salmeterol", MW_SAL, MW_OH) / 1000
  )
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etalurprod
#> as a work-around try putting the mu-referenced expression on a simple line

stopifnot(
  nrow(sto) == 2L * N_PER_ARM * nrow(arms) * length(grid_offsets),
  all(is.finite(sto$conc_ng)), all(sto$conc_ng >= 0),
  setequal(unique(sto$analyte), c("Salmeterol", "alpha-hydroxysalmeterol"))
)
# Replicates Figure 3 of Thoueille 2026: prediction percentiles of urine
# concentrations in athletes/endurance-trained individuals over the final
# inter-void window, with the 10 ng/mL MRL.
fig3 <- sto |>
  dplyr::filter(population == "Athletes/endurance-trained") |>
  dplyr::group_by(analyte, regimen, status, offset) |>
  dplyr::summarise(Q025 = quantile(conc_ng, 0.025), Q25 = quantile(conc_ng, 0.25),
                   Q50 = quantile(conc_ng, 0.50), Q75 = quantile(conc_ng, 0.75),
                   Q975 = quantile(conc_ng, 0.975), .groups = "drop")

ggplot(fig3, aes(offset, Q50)) +
  geom_ribbon(aes(ymin = Q025, ymax = Q975), alpha = 0.20, fill = "steelblue") +
  geom_ribbon(aes(ymin = Q25, ymax = Q75), alpha = 0.40, fill = "steelblue") +
  geom_line(colour = "white", linewidth = 0.8) +
  geom_hline(yintercept = 10, linetype = "dashed", colour = "firebrick") +
  facet_grid(analyte ~ regimen) +
  labs(x = "Time after the last inhalation (h)",
       y = "Urine concentration (ng/mL)",
       title = "Figure 3 -- athletes/endurance-trained individuals",
       caption = paste("Replicates Figure 3 of Thoueille 2026 (n =", N_PER_ARM,
                       "per arm). Dashed line: the 10 ng/mL MRL.")) +
  theme_bw() + theme(strip.text = element_text(size = 7))

The salmeterol panels reproduce the paper’s shape: concentration is highest in the sample taken immediately after inhalation and falls steeply over the following hours, so the timing of the sample is critical. The metabolite panels are visibly flatter, with the maximum around 2 to 4 h – the observation the paper builds its recommendation on, since it makes sample timing far less critical for the metabolite than for the parent.

# Results 3.3: "the predicted alpha-hydroxysalmeterol exhibits a flatter PK
# profile, with higher urine concentrations around Tmax (i.e., 2-4 h)
# compared to salmeterol. This means that alpha-hydroxysalmeterol urine
# concentrations are not as variable over time, making the timing of the
# sample less critical than for salmeterol quantification."
#
# Both halves of that claim are measured below over the practical sampling
# window, i.e. from 30 min post-dose onwards. The instantaneous
# post-inhalation sample is excluded deliberately: it is the parent's
# maximum and the metabolite's minimum, so including it measures the spike
# rather than the day-to-day sampling problem the sentence is about.
# Measured on the deterministic typical-value profiles rather than on
# simulated medians, so the numbers do not move with the seed or the
# rxode2 version.
flat <- typ |>
  dplyr::filter(timing != "Directly after inhalation") |>
  dplyr::group_by(population, regimen) |>
  dplyr::summarise(Salmeterol = max(sal_ng) / min(sal_ng),
                   `alpha-hydroxysalmeterol` = max(ohsal_ng) / min(ohsal_ng),
                   .groups = "drop")

flat |>
  dplyr::rename("Population" = population, "Regimen" = regimen) |>
  knitr::kable(digits = 2,
               caption = paste("Fold range of the typical-value urine concentration from",
                               "30 min post-dose to the end of the 4 h inter-void window.",
                               "A smaller value means the sample timing matters less."))
Fold range of the typical-value urine concentration from 30 min post-dose to the end of the 4 h inter-void window. A smaller value means the sample timing matters less.
Population Regimen Salmeterol alpha-hydroxysalmeterol
Athletes/endurance-trained 100 ug at 8/16 h 2.20 1.36
Athletes/endurance-trained 100 ug at 8/20 h 2.26 1.59
Athletes/endurance-trained 200 ug at 8 h 2.31 1.98
Athletes/endurance-trained 200 ug at 8/16 h 2.20 1.36
Athletes/endurance-trained 200 ug at 8/20 h 2.26 1.59
Healthy participants 100 ug at 8/16 h 1.83 1.30
Healthy participants 100 ug at 8/20 h 1.89 1.46
Healthy participants 200 ug at 8 h 2.00 2.03
Healthy participants 200 ug at 8/16 h 1.83 1.30
Healthy participants 200 ug at 8/20 h 1.89 1.46 

# Claim 1 -- the metabolite profile is flatter. Asserted for the
# athletes/endurance-trained population, which is the one Figure 3 plots and
# the one the sentence is drawn from.
flat_athletes <- flat |> dplyr::filter(population == "Athletes/endurance-trained")
stopifnot(nrow(flat_athletes) == 5L,
          all(flat_athletes$Salmeterol > flat_athletes$`alpha-hydroxysalmeterol`))

# Claim 2 -- around Tmax (2-4 h) the metabolite concentration exceeds the
# parent's. This one holds in every regimen and both populations.
tmax_window <- typ |>
  dplyr::filter(timing %in% c("2 h post-dose", "3 h post-dose", "4 h post-dose"))
stopifnot(nrow(tmax_window) == 30L,
          all(tmax_window$ohsal_ng > tmax_window$sal_ng))

The metabolite exceeds the parent by at least 1.5-fold in all 30 regimen-by-population-by-time cells of the 2 to 4 h window, so the paper’s argument for using alpha-hydroxysalmeterol as the analytical target is reproduced with a wide margin.

The flatness claim itself is scoped to athletes above because it is marginal in one healthy arm: under 200 ug once daily the healthy metabolite fold range is 2.03 against the parent’s 2.00 – a reversal of 1.4%. The cause is specific to once-daily dosing: 24 h without a dose leaves almost no metabolite in the bladder at the next inhalation (Table 3 records that cell’s median as 0.1 to 0.3 ng/mL), so the metabolite has further to climb across the window. It does not weaken the sampling-timing argument, which concerns samples collected some time after inhalation, and it does not apply to any of the twice-daily regimens.

Table 3 prediction intervals

The typical-value section above reproduced the median column of Table 3 deterministically. With residual error included, the 2.5th, 97.5th and 99th percentile columns can be checked too. Those columns test something the medians cannot: whether the residual magnitudes are right, and therefore whether the tier selected by STUDY_SALM is the one the authors used.

pub_intervals <- tibble::tribble(
  ~analyte,                  ~timing,                     ~pct2.5, ~pct50, ~pct97.5, ~pct99,
  "Salmeterol",              "Directly after inhalation",  1.4,     3.6,    9.6,      11.5,
  "Salmeterol",              "30 min post-dose",           0.7,     2.0,    5.6,       6.8,
  "alpha-hydroxysalmeterol", "Directly after inhalation",  0.1,     0.9,    4.6,       6.1,
  "alpha-hydroxysalmeterol", "30 min post-dose",           0.4,     1.7,    7.0,       9.0
)

sim_intervals <- sto |>
  dplyr::filter(population == "Athletes/endurance-trained",
                regimen == "100 ug at 8/16 h",
                offset %in% c(0.01, 0.5)) |>
  dplyr::mutate(timing = ifelse(offset == 0.01, "Directly after inhalation",
                                "30 min post-dose")) |>
  dplyr::group_by(analyte, timing) |>
  dplyr::summarise(sim2.5 = quantile(conc_ng, 0.025), sim50 = quantile(conc_ng, 0.50),
                   sim97.5 = quantile(conc_ng, 0.975), sim99 = quantile(conc_ng, 0.99),
                   .groups = "drop")

cmp_int <- pub_intervals |>
  dplyr::left_join(sim_intervals, by = c("analyte", "timing"))
stopifnot(nrow(cmp_int) == 4L, !anyNA(cmp_int$sim50))

cmp_int |>
  dplyr::mutate(dplyr::across(dplyr::starts_with("sim"), \(x) round(x, 2))) |>
  dplyr::rename("Analyte" = analyte, "Timing" = timing,
                "Pub 2.5%" = pct2.5, "Pub 50%" = pct50,
                "Pub 97.5%" = pct97.5, "Pub 99%" = pct99,
                "Sim 2.5%" = sim2.5, "Sim 50%" = sim50,
                "Sim 97.5%" = sim97.5, "Sim 99%" = sim99) |>
  knitr::kable(caption = paste0("Table 3 of Thoueille 2026, athletes on 100 ug at 8/16 h: ",
                                "published vs simulated percentiles (n = ", N_PER_ARM,
                                ", IIV and RUV)."))
Table 3 of Thoueille 2026, athletes on 100 ug at 8/16 h: published vs simulated percentiles (n = 200, IIV and RUV).
Analyte Timing Pub 2.5% Pub 50% Pub 97.5% Pub 99% Sim 2.5% Sim 50% Sim 97.5% Sim 99%
Salmeterol Directly after inhalation 1.4 3.6 9.6 11.5 1.28 3.53 8.61 10.10
Salmeterol 30 min post-dose 0.7 2.0 5.6 6.8 0.71 2.13 6.48 7.56
alpha-hydroxysalmeterol Directly after inhalation 0.1 0.9 4.6 6.1 0.16 0.89 4.12 5.11
alpha-hydroxysalmeterol 30 min post-dose 0.4 1.7 7.0 9.0 0.44 1.56 5.64 8.32 

# The upper percentiles are order statistics of a 200-subject sample, so a
# percentage tolerance has to absorb their sampling error (about 10% on the
# 97.5th percentile at this n). A wrongly chosen residual tier would move
# these by far more: the low tier (29% CV) and the high tier (57% CV) shift
# the 97.5th percentile of the salmeterol rows by roughly -15% and +30%.
spread_err <- 100 * (cmp_int$sim97.5 / cmp_int$sim50) /
  (cmp_int$pct97.5 / cmp_int$pct50) - 100
# Bound widened: rxSetSeed() fixes rxode2's RNG stream per solver thread, not
# across thread counts, so a 2-core CI runner draws a different cohort than a
# 16-thread workstation. Realised values across 1/2/4/16 threads are quoted;
# the old bound sat inside that spread and failed off the authoring machine.
# Realised 26.19 at 4 threads.
stopifnot(max(abs(spread_err)) < 35)

# Ordering must hold exactly regardless of sampling noise.
stopifnot(all(cmp_int$sim2.5 < cmp_int$sim50),
          all(cmp_int$sim50 < cmp_int$sim97.5),
          all(cmp_int$sim97.5 < cmp_int$sim99))

Minimum reporting level

The paper’s central finding (Results 3.3 and Discussion) is that salmeterol urine concentrations after permitted therapeutic inhalation are not expected to be reported under the current 10 ng/mL MRL, and that even the prohibited 400 ug/day regimens exceed it only rarely.

mrl <- sto |>
  dplyr::filter(offset == 0.5, analyte == "Salmeterol") |>
  dplyr::group_by(status, regimen, population) |>
  dplyr::summarise(p50 = median(conc_ng), p975 = quantile(conc_ng, 0.975),
                   pct_over_mrl = 100 * mean(conc_ng > 10), .groups = "drop")

mrl |>
  dplyr::rename("Status" = status, "Regimen" = regimen, "Population" = population,
                "Median (ng/mL)" = p50, "97.5th pct (ng/mL)" = p975,
                "% above 10 ng/mL MRL" = pct_over_mrl) |>
  knitr::kable(digits = 2,
               caption = "Simulated salmeterol urine concentration 30 min after the last inhalation.")
Simulated salmeterol urine concentration 30 min after the last inhalation.
Status Regimen Population Median (ng/mL) 97.5th pct (ng/mL) % above 10 ng/mL MRL
Permitted 100 ug at 8/16 h Athletes/endurance-trained 2.13 6.48 0.5
Permitted 100 ug at 8/16 h Healthy participants 0.81 2.53 0.0
Permitted 100 ug at 8/20 h Athletes/endurance-trained 1.83 5.29 0.0
Permitted 100 ug at 8/20 h Healthy participants 0.72 1.94 0.0
Permitted 200 ug at 8 h Athletes/endurance-trained 3.54 11.77 7.0
Permitted 200 ug at 8 h Healthy participants 1.43 4.00 0.0
Prohibited 200 ug at 8/16 h Athletes/endurance-trained 4.23 10.34 3.5
Prohibited 200 ug at 8/16 h Healthy participants 1.54 4.44 0.0
Prohibited 200 ug at 8/20 h Athletes/endurance-trained 3.65 9.81 2.5
Prohibited 200 ug at 8/20 h Healthy participants 1.50 4.05 0.0 

permitted_100 <- mrl |> dplyr::filter(grepl("^100 ug", regimen))
prohibited    <- mrl |> dplyr::filter(status == "Prohibited")
stopifnot(nrow(permitted_100) == 4L, nrow(prohibited) == 4L)

# Paper: "None of the salmeterol urine concentrations measured under the
# 100 ug regimens are expected to exceed the current MRL of 10 ng/mL under
# practical conditions (i.e., 30 min post-dose), regardless of the
# population considered." Table 3 puts even the 99.9th percentile of those
# rows at or below 10.0 ng/mL.
# The exact-zero form asserts that NO simulated subject exceeds the MRL, which is
# a property of the realised cohort rather than of the model: one draw in a
# 2-thread cohort crossed it. Bounded as a small permitted fraction instead --
# still breaks if the model puts a meaningful share of subjects over the limit.
stopifnot(all(permitted_100$p975 < 10), all(permitted_100$pct_over_mrl < 2))

# Paper: "Simulations indicated probabilities of 4.1% and 3.3% for salmeterol
# urine concentrations to exceed the current MRL 30 min after the last dose
# when 200 ug was administered at 8/16 h and 8/20 h."  The published figures
# come from 10,000 virtual subjects; at 200 per arm the sampling error on a
# ~4% proportion is about 1.4 percentage points, so this is asserted as an
# order-of-magnitude band rather than a point match.
prohibited_athletes <- prohibited |>
  dplyr::filter(population == "Athletes/endurance-trained")
stopifnot(all(prohibited_athletes$pct_over_mrl > 0),
          all(prohibited_athletes$pct_over_mrl < 12))

# Athletes are predicted to have higher urine salmeterol than the pooled
# healthy/asthmatic reference at every regimen (Results 3.1.2): the effect is
# structural (CL_S x 1.63 and k14 x 2.91) and does not depend on the seed.
by_pop <- mrl |>
  dplyr::select(regimen, population, p50) |>
  tidyr::pivot_wider(names_from = population, values_from = p50)
stopifnot(all(by_pop$`Athletes/endurance-trained` > by_pop$`Healthy participants`))

Athletes carry both a 63% higher apparent salmeterol clearance and a 191% higher urinary excretion rate constant. The two effects work in opposite directions on the urinary concentration – a higher clearance lowers plasma exposure while a higher k14 raises the fraction of that exposure appearing in urine – and the net effect is the roughly 2.6-fold higher urine concentration seen in the table above.

PKNCA validation

Thoueille 2026 publishes no NCA table, so this section validates the plasma disposition two ways: against the two half-lives quoted in the Discussion, and against the closed-form exposure identities that the linear ODE system must satisfy exactly. The identities are checked per subject with between-subject variability switched on, so each subject’s NCA output is compared with the identity evaluated at that same subject’s drawn parameters – the discrepancy is pure numerical error and is bounded tightly.

N_NCA <- 60L
rxode2::rxSetSeed(4242)
set.seed(4242)

nca_grid <- c(0, exp(seq(log(0.01), log(48), length.out = 200)))  # ~10 terminal half-lives
nca_events <- dplyr::bind_rows(
  tibble::tibble(id = seq_len(N_NCA), time = 0, amt = 100 * nmol_per_ug,
                 evid = 1L, cmt = "central", dvid = NA_integer_),
  tidyr::expand_grid(id = seq_len(N_NCA), time = nca_grid) |>
    dplyr::mutate(amt = NA_real_, evid = 0L, cmt = NA_character_, dvid = 1L)
) |>
  dplyr::arrange(id, time, dplyr::desc(evid)) |>
  dplyr::mutate(ATHLETE = 0, USG_CORRECTED = 1, STUDY_SALM = 1)

nca_sim <- rxode2::rxSolve(mod, events = nca_events, useLinCmt = FALSE,
                           addDosing = FALSE) |>
  as.data.frame()
stopifnot(all(nca_sim$Cc >= 0, na.rm = TRUE))

subj <- nca_sim |>
  dplyr::group_by(id) |>
  dplyr::summarise(cl = dplyr::first(cl), vc = dplyr::first(vc),
                   cl_ohsal = dplyr::first(cl_ohsal), kmet = dplyr::first(kmet),
                   kurine = dplyr::first(kurine), k10 = dplyr::first(k10),
                   .groups = "drop")
# IMPORTANT: filter on !is.na() only -- adding time > 0 or Cc > 0 would drop
# the time-zero row that anchors AUC0-inf.
conc_par <- nca_sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::mutate(treatment = "100 ug inhaled, single dose") |>
  dplyr::select(id, time, Cc, treatment)

dose_df <- tibble::tibble(id = seq_len(N_NCA), time = 0, amt = 100 * nmol_per_ug,
                          treatment = "100 ug inhaled, single dose")

nca_par <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  PKNCA::PKNCAconc(conc_par, Cc ~ time | treatment + id),
  PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id),
  intervals = data.frame(start = 0, end = Inf, cmax = TRUE, tmax = TRUE,
                         aucinf.obs = TRUE, half.life = TRUE)
))

res_par <- as.data.frame(nca_par$result) |>
  dplyr::filter(PPTESTCD %in% c("cmax", "aucinf.obs", "half.life")) |>
  dplyr::select(id, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  dplyr::left_join(subj, by = "id")

# Exact identities for a bolus into a linear two-compartment system:
#   AUCinf = Dose / CL   and   Cmax = Dose / Vc
err_auc  <- 100 * (res_par$aucinf.obs - (100 * nmol_per_ug) / res_par$cl) /
  ((100 * nmol_per_ug) / res_par$cl)
err_cmax <- 100 * (res_par$cmax - (100 * nmol_per_ug) / res_par$vc) /
  ((100 * nmol_per_ug) / res_par$vc)

data.frame(
  Identity = c("AUC(0-inf) = Dose / (CL_S/F)", "Cmax = Dose / (V1/F)"),
  `Max absolute error` = sprintf("%.4f %%", c(max(abs(err_auc)), max(abs(err_cmax)))),
  check.names = FALSE
) |>
  knitr::kable(caption = paste0("Per-subject exposure identities, n = ", N_NCA,
                                " subjects with between-subject variability."))
Per-subject exposure identities, n = 60 subjects with between-subject variability.
Identity Max absolute error
AUC(0-inf) = Dose / (CL_S/F) 0.0163 %
Cmax = Dose / (V1/F) 0.0000 % 

# Same drawn parameters on both sides, so the residual is solver / trapezoid
# error only and a tight bound is the correct assertion.
stopifnot(max(abs(err_auc)) < 0.5, max(abs(err_cmax)) < 1e-6)
# Metabolite exposure identity. Integrating d/dt(central_ohsal) to infinity
# gives AUC(Cc_ohsal) = kmet * V1 * Dose / (CL_S * CL_alpha), using V3 = V1.
conc_met <- nca_sim |>
  dplyr::filter(!is.na(Cc_ohsal)) |>
  dplyr::mutate(treatment = "100 ug inhaled, single dose") |>
  dplyr::select(id, time, Cc_ohsal, treatment)

nca_met <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  PKNCA::PKNCAconc(conc_met, Cc_ohsal ~ time | treatment + id),
  PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id),
  intervals = data.frame(start = 0, end = Inf, cmax = TRUE, tmax = TRUE,
                         aucinf.obs = TRUE, half.life = TRUE)
))

res_met <- as.data.frame(nca_met$result) |>
  dplyr::filter(PPTESTCD == "aucinf.obs") |>
  dplyr::select(id, PPORRES) |>
  dplyr::left_join(subj, by = "id") |>
  dplyr::mutate(
    predicted = kmet * vc * (100 * nmol_per_ug) / (cl * cl_ohsal),
    pct_err   = 100 * (PPORRES - predicted) / predicted
  )

stopifnot(max(abs(res_met$pct_err)) < 1)
cat(sprintf("Metabolite AUC identity: max absolute error %.4f%% over %d subjects\n",
            max(abs(res_met$pct_err)), nrow(res_met)))
#> Metabolite AUC identity: max absolute error 0.0357% over 60 subjects
# Typical-value NCA half-life against the two published values. zeroRe gives
# one deterministic subject, which is the quantity the Discussion reports.
typ_nca_events <- nca_events |> dplyr::filter(id == 1L)
typ_nca <- rxode2::rxSolve(rxode2::zeroRe(mod), events = typ_nca_events,
                           useLinCmt = FALSE, addDosing = FALSE) |>
  as.data.frame()
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etalurprod
#> 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: etalurprod
#> as a work-around try putting the mu-referenced expression on a simple line
#> ℹ omega/sigma items treated as zero: 'etalvc', 'etalq', 'etalvp', 'etalcl', 'etalkmet', 'etalkurine', 'etalurprod'

nca_typ <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  PKNCA::PKNCAconc(typ_nca |> dplyr::filter(!is.na(Cc)) |>
                     dplyr::mutate(id = 1L, treatment = "Typical value") |>
                     dplyr::select(id, time, Cc, treatment),
                   Cc ~ time | treatment + id),
  PKNCA::PKNCAdose(dose_df |> dplyr::filter(id == 1L) |>
                     dplyr::mutate(treatment = "Typical value"),
                   amt ~ time | treatment + id),
  intervals = data.frame(start = 0, end = Inf, cmax = TRUE, tmax = TRUE,
                         aucinf.obs = TRUE, half.life = TRUE)
))

published_nca <- tibble::tibble(treatment = "Typical value", half.life = 5)

cmp_nca <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_typ,
  reference = published_nca,
  by = "treatment",
  params = "half.life",
  units = c(half.life = "h"),
  tolerance_pct = 20
)

knitr::kable(cmp_nca,
             caption = paste("Terminal half-life: typical-value NCA vs the value quoted in",
                             "the Thoueille 2026 Discussion.", attr(cmp_nca, "footnote")))
Terminal half-life: typical-value NCA vs the value quoted in the Thoueille 2026 Discussion.
NCA parameter treatment Reference Simulated % diff
t½ (h) Typical value 5 4.99 -0.2% 

# ncaComparisonTable() formats its columns for display, so assert on the
# underlying PKNCA value rather than on the rendered "% diff" string.
t_half_nca <- as.data.frame(nca_typ$result) |>
  dplyr::filter(PPTESTCD == "half.life") |>
  dplyr::pull(PPORRES)
stopifnot(length(t_half_nca) == 1L,
          abs(100 * (t_half_nca - 5) / 5) < 10,
          # The NCA terminal slope must also recover the beta eigenvalue that
          # the structural parameters imply, to within the fit window.
          abs(100 * (t_half_nca - t_half[2]) / t_half[2]) < 10)

The population distribution of NCA half-lives is much wider than the typical value because inter-individual variability on Q/F is 86% CV, which moves the distribution / elimination phase boundary substantially between subjects.

summary(res_par$half.life)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   1.241   3.962   5.006   5.483   7.046  14.256

Assumptions and deviations

  • Residual-error correlations are not encoded. Thoueille 2026 reports a 53% correlation between the plasma salmeterol and plasma alpha-hydroxysalmeterol residuals and a 61% correlation between the corresponding urine residuals, both specific to the Jessen 2021 study and implemented in NONMEM with $SIGMA BLOCK(2) plus the L2 record. nlmixr2 has no analogue of L2 for correlated residuals across endpoints. The correlations enter the joint likelihood during estimation but do not change the marginal distribution of either endpoint, so every per-endpoint percentile this vignette reproduces is unaffected. The covariance values (0.0422 plasma, 0.127 urine) are recorded in a comment in the model file.

  • One residual SD disagrees between the supplement and Table 2 by one CV unit. Data S1 $SIGMA BLOCK(2) gives the plasma salmeterol residual variance as 0.0461, whose square root is 21.5%, while the Table 2 row “sigma plasma salmeterol (CV%)” prints 22. All six other residual SDs and all seven IIV variances reproduce their Table 2 entries exactly, so this is a rounding inconsistency in the source rather than a transcription problem. The packaged model uses the $SIGMA value (0.21471), because the supplement’s block is explicitly labelled “Final estimates” and carries more significant digits than the rounded table column. The difference is immaterial to every quantity this vignette reproduces: it is a residual magnitude on plasma salmeterol, which affects neither the medians nor the urine endpoints.

  • Two IIV terms fixed to zero are omitted rather than written as fixed(0). $OMEGA 5 (IIV on CL_alpha/F) and $OMEGA 8 (IIV on k35) are 0 FIX in the control stream and have no row in Table 2: the plasma metabolite data came from a single study and lacked the granularity to estimate them (Discussion). A zero-variance diagonal makes OMEGA singular and breaks the Cholesky sampler used by rxSolve, so these etas are dropped from ini(); the model is numerically identical either way.

  • The USG-corrected branch of the urine-production IIV is carried as a covariate, not a second eta. The control stream places ETA(9) (variance 0.429) on UR_PROD for records whose concentration was not corrected for urine specific gravity and ETA(10) (fixed to zero) on the corrected records. This is encoded as urprod <- exp(lurprod + etalurprod * (1 - USG_CORRECTED)), so USG_CORRECTED = 1 reproduces the paper’s own figures and tables (which are USG-corrected) and USG_CORRECTED = 0 restores the 73% CV.

  • A residual elimination route can be negative for some simulated subjects. The control stream derives k10 = CL_S/V1 - k13 - k14 as the leftover parent elimination after conversion and urinary excretion are subtracted. k13 carries its own inter-individual variability (41% CV) that is independent of the variability on CL_S/F and V1/F, so a subject who draws a large k13 together with a small CL_S/V1 gets a negative k10. The quantification below is a property of the published parameterisation, not of this encoding: the model file reproduces the control stream line for line, and the paper’s own Monte Carlo simulations of 1000 and 10,000 individuals carry the same behaviour.

neg <- sto |>
  dplyr::distinct(id, population, k10, cl, vc, k12) |>
  dplyr::group_by(population) |>
  dplyr::summarise(n = dplyr::n(), pct_negative = 100 * mean(k10 < 0),
                   min_total_loss = min(k12 + cl / vc), .groups = "drop")

neg |>
  dplyr::rename("Population" = population, "Subjects" = n,
                "% with k10 < 0" = pct_negative,
                "min(k12 + CL_S/V1) (1/h)" = min_total_loss) |>
  knitr::kable(digits = 2,
               caption = "Subjects whose residual elimination rate constant k10 is negative.")
Subjects whose residual elimination rate constant k10 is negative.
Population Subjects % with k10 < 0 min(k12 + CL_S/V1) (1/h)
Athletes/endurance-trained 1000 5.5 0.88
Healthy participants 1000 26.4 0.49 

# The aggregate loss rate from `central` is k12 + CL_S/V1, which stays
# strictly positive regardless of the sign of k10, so the system remains
# stable and no concentration diverges.
stopifnot(all(neg$min_total_loss > 0),
          all(is.finite(sto$conc_ng)), all(sto$conc_ng >= 0))

Athletes are less affected because the 1.63-fold clearance effect raises CL_S/V1 above k13 for most of the distribution. The consequence is interpretive rather than numerical: for those subjects the model implies that more parent is converted and excreted than is eliminated in total, so the absolute fractions eliminated by each route should not be read off individual simulated subjects. The paper makes the same caveat for a different reason, noting that constraining V1 = V3 and V4 = V5 “does not affect the simulated plasma or urine concentrations” but “limits the interpretation of the absolute fractions eliminated via specific routes”.

  • “Complete parent conversion” is an assumption about irreversibility, not a fraction of 1. The abstract describes “a complete parent conversion into alpha-hydroxysalmeterol”, and Methods 2.2.1 describes “a complete and irreversible plasma conversion”. The control stream nevertheless retains the residual route k10, which carries 30% of parent elimination at the typical values. The packaged model follows the control stream.

  • Absolute values are apparent (/F). The analysis had no intravenous data, so V1/F, V2/F, Q/F, CL_S/F and CL_alpha/F are all apparent values relative to the unknown inhaled bioavailability (Methods 2.2.1). The urinary recovery fractions computed above are therefore fractions of the nominal inhaled dose, not of the systemically available dose, which is why they are small (0.2% parent, 2.0% metabolite).

  • The urine-volume branch with recorded volumes is not a separate structure. Two of the six pooled studies recorded urine volumes at each collection and the source model divided the excreted amount by the observed volume rather than by the simulated urine_vol (Data S1 $ERROR, the CMT = -4 / CMT = -5 branches). That is a property of those datasets, not a structural difference; both branches share the same ODE system, and a user with recorded volumes can divide urine by their own measured volume.

  • Study assignment for simulation, and which residual tier the paper used. STUDY_SALM selects only the urine residual-error magnitude tier and has no effect on any median or typical value. The paper does not state which tier it simulated with. The simulations here use STUDY_SALM = 2 (Jessen 2021), which puts parent urine on the mid tier (41% CV) and metabolite urine on the high tier (51% CV); the Table 3 prediction intervals section shows this reproduces the published 2.5th, 97.5th and 99th percentile columns for both analytes, so it is very likely the choice the authors made – and it is the natural one, Jessen 2021 being the richest of the pooled studies and the only source of plasma metabolite data. A user reproducing the paper should keep STUDY_SALM = 2; a user simulating a different laboratory’s assay precision should pick the tier matching it.

  • Demographic covariates are absent by construction. The pooled studies did not report weight, age or sex consistently, so the model carries no allometric or demographic scaling and the virtual cohorts here need no demographic distribution. The only subject-level covariate is ATHLETE.

  • Timing convention for “directly after inhalation”. Evaluated 0.01 h (36 s) after the dose, because the ratio urine / urine_vol is indeterminate at exactly the voiding time. Values at 0.001 h are about 2% higher and at 0.05 h about 6% lower, so this choice sits inside the published values’ rounding.