Skip to contents

Benznidazole popPK and qPCR exposure-response (Assmus 2025)

Source: vignettes/articles/Assmus_2025_benznidazole.Rmd
Assmus_2025_benznidazole.Rmd

Model and source

This paper contributes two model files:

  • Assmus_2025_benznidazole – the population PK model (Assmus 2025 Table 2).
  • Assmus_2025_benznidazole_qpcr – the same PK layer extended with a cumulative AUC state and the paper’s beta binomial exposure-response model for post-treatment Trypanosoma cruzi qPCR positivity (Assmus 2025 Equation 2 and Table 3).
pk <- rxode2::rxode(readModelDb("Assmus_2025_benznidazole"))
#> ℹ parameter labels from comments will be replaced by 'label()'
pd <- rxode2::rxode(readModelDb("Assmus_2025_benznidazole_qpcr"))
#> ℹ parameter labels from comments will be replaced by 'label()'
  • Citation: Assmus F, Cruz C, Watson JA, White NJ, Adehin A, Hoglund RM, Blum de Oliveira B, Barreira F, Scandale I, Tarning J. Population pharmacokinetic-pharmacodynamic analysis of benznidazole monotherapy and combination therapy with fosravuconazole in chronic Chagas disease (BENDITA). PLoS Neglected Tropical Diseases. 2025;19(9):e0013522. doi:10.1371/journal.pntd.0013522.
  • Article: https://doi.org/10.1371/journal.pntd.0013522

Population

The BENDITA trial (NCT03378661) randomised 210 Bolivian adults (18-50 years, 50-80 kg) with chronic indeterminate Chagas disease to one of seven arms of 30: six active benznidazole regimens and a placebo arm. Benznidazole was given orally as Abarax tablets, split into a morning and an evening dose. The active arms were 150 mg twice daily for 8, 4 or 2 weeks; 150 mg once daily for 4 weeks alone or combined with fosravuconazole (E1224); and 300 mg once weekly (divided into two doses) for 8 weeks with fosravuconazole. Patients with signs of the chronic cardiac or digestive form of the disease were excluded.

Benznidazole was measured in dried blood spots (HPLC-MS/MS, LLOQ 50 ng/mL), so every disposition parameter in this model is a whole-blood apparent value. Sampling was sparse: predose plus nine protocol time points per patient, and the dataset was censored at 120 h after dose. The PK analysis used 175 of the 180 actively treated subjects (986 concentrations, 1.7 % below LLOQ) after excluding five subjects whose profiles were incompatible with their assigned arm. The exposure-response analysis used the modified intention-to-treat population, n = 201 including placebo. Baseline characteristics are in Assmus 2025 Table 1.

str(pk$population, max.level = 1)
#> List of 15
#>  $ species       : chr "human"
#>  $ n_subjects    : int 175
#>  $ n_studies     : int 1
#>  $ n_observations: int 986
#>  $ age_range     : chr "18-50 years"
#>  $ age_median    : chr "34 years"
#>  $ weight_range  : chr "50-80 kg"
#>  $ weight_median : chr "65 kg"
#>  $ sex_female_pct: num 68
#>  $ race_ethnicity: chr "All participants were Bolivian."
#>  $ disease_state : chr "Adults (18-50 years, 50-80 kg) with chronic indeterminate Chagas disease, confirmed by serological testing and "| __truncated__
#>  $ dose_range    : chr "Six active oral benznidazole arms (Abarax, Laboratorios ELEA): 150 mg twice daily for 8, 4 or 2 weeks; 150 mg o"| __truncated__
#>  $ regions       : chr "Bolivia (Cochabamba, Tarija and Sucre)."
#>  $ sampling      : chr "Sparse dried blood spot (DBS) sampling: predose and at follow-up visits on days 1-3 and weeks 2, 3, 4, 6 and 10"| __truncated__
#>  $ notes         : chr "BENDITA trial, ClinicalTrials.gov NCT03378661, conducted 2016-2018. 210 adults were randomised (30 per arm); th"| __truncated__

Source trace

Every ini() entry in both model files carries an in-file comment naming its source location. The table below collects them.

Equation / parameter Value Source location
Transit chain gut -> transit1 -> central, ktr = (1 + n) / MTT, n = 1 n/a Fig 2 (structural schematic)
One-compartment linear disposition n/a Results, Pharmacokinetics
lmtt (MTT) 0.75 h Table 2
lcl (CL/F) 1.30 L/h Table 2
lvc (V/F) 31.6 L Table 2
lfdepot (F) 1, fixed Table 2; Methods, popPK analysis (i)
e_wt_cl 0.75, fixed Methods, popPK analysis (i)
e_wt_vc 1.0, fixed Methods, popPK analysis (i)
Reference weight 65 kg Methods, popPK analysis (i); Table 1 median
e_conmed_fosravuconazole_cl +17.7 % Table 2
e_sexf_fdepot -12.9 % (men vs women) Table 2
etalmtt 61.6 % CV Table 2
etalcl 18.9 % CV Table 2
etalfdepot 10.2 % CV Table 2
expSd sqrt(0.076) Table 2 (variance of residual error)
d/dt(auc_central) <- Cc (AUCinf) n/a Methods, Definition of exposure metrics
logit(p) = b0 + b1 * x + b2 * Ct n/a Equation 2
e_auc_logit log(0.9995) Table 3, mITT incl. placebo, AUCinf column
e_ct_tcruzi_base_logit log(0.9021) Table 3, same column
b0_qpcr 3.496 Not published – digitised from Fig 8a; see Assumptions
Median baseline Ct 37.7 cycles Table 1 (PK/PD mITT column)
IC90 in DBS 7.61 mg/L (29.3 uM) Methods, Definition of exposure metrics
Historical target range in DBS 2.5-5 mg/L Methods, Definition of exposure metrics
Plasma-to-DBS scaling factor 0.84 S1 Table footnote a (after Bedor 2018)
Per-arm AUCinf / CMAX / TMAX / terminal t-half see below S1 Table (pooled column)

The validation reference values below come from S1 Table (Secondary pharmacokinetic parameter estimates, based on the final population pharmacokinetic model for benznidazole), a supplementary DOCX retrieved from the PLoS supporting-information endpoint and stored alongside the article PDF. It reports median (5th-95th percentile) AUCinf, CMAX, TMAX and terminal half-life for each of the six active arms, separately for women, men and pooled. The pooled column is used here.

Virtual cohort

Original participant data are held by IDDO under controlled access, so the figures below use virtual cohorts whose covariates approximate Assmus 2025 Table 1: body weight median 65 kg over the 50-80 kg enrolment window, 68 % female in the PK analysis set, and baseline qPCR cycle threshold median 37.7 cycles (range 30.2-40.0, ceiling 40).

set.seed(20250922)

N_PER_ARM <- 100L  # well under the 200/arm cap; medians are stable at this size

# Truncated normal via the inverse CDF (no rejection loop).
rtnorm <- function(n, mean, sd, lo, hi) {
  qnorm(runif(n, pnorm(lo, mean, sd), pnorm(hi, mean, sd)), mean, sd)
}

arms <- tibble::tribble(
  ~arm,                          ~ii, ~ndose, ~combo, ~split_day,
  "150 mg BID, 8 wk",             12,   112L,      0,      FALSE,
  "150 mg BID, 4 wk",             12,    56L,      0,      FALSE,
  "150 mg BID, 2 wk",             12,    28L,      0,      FALSE,
  "150 mg QD, 4 wk",              24,    28L,      0,      FALSE,
  "150 mg QD, 4 wk + E1224",      24,    28L,      0,      FALSE,
  "300 mg weekly, 8 wk + E1224", 168,     8L,      1,       TRUE
)
arms$combo[arms$arm == "150 mg QD, 4 wk + E1224"] <- 1
arms <- arms |> mutate(arm = factor(arm, levels = arm))

# The weekly arm gives 300 mg "divided in two doses" on one day per week
# (Methods, Treatment regimens), i.e. 150 mg morning + 150 mg evening. Every
# arm therefore administers 150 mg units.
dose_times <- function(a) {
  base <- a$ii * seq(0, a$ndose - 1)
  if (a$split_day) sort(c(base, base + 12)) else base
}

# Build a dosing + observation event table for one arm from a subject frame.
build_events <- function(subj, a, obs_times) {
  dt <- dose_times(a)
  bind_rows(
    subj |> tidyr::crossing(time = dt) |>
      mutate(amt = 150, evid = 1L, cmt = "depot"),
    subj |> tidyr::crossing(time = obs_times(max(dt))) |>
      mutate(amt = NA_real_, evid = 0L, cmt = "central")
  ) |>
    arrange(id, time, desc(evid))
}

# 1 h over the whole course plus a 10-half-life washout (AUCinf needs the tail),
# refined to 0.05 h across the interval following the final dose. By then every
# regimen has accumulated to steady state, so that interval's peak is the peak
# of the whole course, and 0.05 h resolves TMAX finely enough to compare against
# S1 Table's 1.74-1.99 h without grid quantisation. The window runs 10 h past
# the last dose, which covers TMAX even for a subject at the 95th percentile of
# mean transit time; it stays inside the 12 h twice-daily interval.
grid_full <- function(last) {
  sort(unique(c(seq(0, last + 240, by = 1),
                seq(last, last + 10, by = 0.05))))
}

offsets <- (seq_len(nrow(arms)) - 1L) * N_PER_ARM
events <- do.call(bind_rows, lapply(seq_len(nrow(arms)), function(i) {
  a <- arms[i, ]
  subj <- tibble(
    id             = offsets[i] + seq_len(N_PER_ARM),
    WT             = rtnorm(N_PER_ARM, 65, 7.5, 50, 80),
    SEXF           = as.numeric(seq_len(N_PER_ARM) > round(0.32 * N_PER_ARM)),
    CONMED_FOSRAVUCONAZOLE = a$combo,
    CT_TCRUZI_BASE = pmin(rtnorm(N_PER_ARM, 37.5, 2.0, 30.2, 40.0), 40),
    arm            = a$arm
  )
  build_events(subj, a, grid_full)
}))

# anyDuplicated() must see the raw rows: wrapping this in unique() would remove
# the duplicates first and the assertion could never go red.
stopifnot(!anyDuplicated(events[, c("id", "time", "evid")]))
stopifnot(dplyr::n_distinct(events$id) == N_PER_ARM * nrow(arms))

Simulation

A single solve of the exposure-response model supplies everything the sections below need: it contains the PK layer verbatim, plus the cumulative-AUC state and the qPCR-positivity probability.

# Every stochastic block seeds itself. base::set.seed() governs only the
# covariate draws above; the between-subject etas are drawn by rxode2's own
# RNG, which answers to rxSetSeed() alone.
rxode2::rxSetSeed(20250922)

sim <- rxode2::rxSolve(
  pd, events = events,
  keep = c("arm", "WT", "SEXF", "CONMED_FOSRAVUCONAZOLE", "CT_TCRUZI_BASE")
) |>
  as.data.frame() |>
  filter(!is.na(Cc))
nrow(sim)
#> [1] 738600

Structural checks

Two identities follow from the published parameters alone and are the sharpest tests available here.

Terminal half-life 

t_half <- function(cl_mult) log(2) * 31.6 / (1.30 * cl_mult)
tibble(
  Quantity  = c("Terminal half-life, monotherapy",
                "Terminal half-life, with fosravuconazole"),
  Model     = sprintf("%.1f h", c(t_half(1), t_half(1.177))),
  Published = c("~17 h", "~14 h")
) |>
  knitr::kable(caption = "Half-life implied by Table 2 vs Assmus 2025 Results, Pharmacokinetics.")
Half-life implied by Table 2 vs Assmus 2025 Results, Pharmacokinetics.
Quantity Model Published
Terminal half-life, monotherapy 16.8 h ~17 h
Terminal half-life, with fosravuconazole 14.3 h ~14 h

Total-course AUC equals dose / clearance

Absorption is complete and elimination is linear, so cumulative AUC to infinity over a whole course must equal total dose * F / CL for every subject. Running this on a typical-value cohort spanning both sexes, the full 50-80 kg weight range and all six arms turns it into a test of the covariate model as well: the expected value below is computed by hand from Assmus 2025 Table 2, not read back from the solve.

tv <- tidyr::crossing(arm_i = seq_len(nrow(arms)), SEXF = c(0, 1),
                      WT = c(50, 65, 80)) |>
  mutate(id = row_number())

tv_events <- do.call(bind_rows, lapply(seq_len(nrow(arms)), function(i) {
  a <- arms[i, ]
  subj <- tv |> filter(arm_i == i) |>
    transmute(id, WT, SEXF,
              CONMED_FOSRAVUCONAZOLE = a$combo,
              CT_TCRUZI_BASE = 37.7, arm = a$arm)
  # Only the washout endpoint is needed; doses drive the integration.
  build_events(subj, a, function(last) c(0, last + 240))
}))

tv_sim <- rxode2::rxSolve(rxode2::zeroRe(pd), events = tv_events,
                          omega = NA, sigma = NA,
                          keep = c("arm", "WT", "SEXF",
                                   "CONMED_FOSRAVUCONAZOLE")) |>
  as.data.frame()

tv_chk <- tv_sim |>
  group_by(id, arm, WT, SEXF, CONMED_FOSRAVUCONAZOLE) |>
  summarise(aucinf = max(auc_central), .groups = "drop") |>
  left_join(
    tv_events |> filter(evid == 1L) |>
      group_by(id) |> summarise(dose_total = sum(amt), .groups = "drop"),
    by = "id"
  ) |>
  mutate(
    cl_typ   = 1.30 * (WT / 65)^0.75 * (1 + 0.177 * CONMED_FOSRAVUCONAZOLE),
    f_typ    = 1 + (-0.129) * (1 - SEXF),
    expected = dose_total * f_typ / cl_typ,
    relerr   = aucinf / expected - 1
  )

cat(sprintf("max |AUCinf / (Dose * F / CL) - 1| over %d typical subjects = %.4f %%\n",
            nrow(tv_chk), 100 * max(abs(tv_chk$relerr))))
#> max |AUCinf / (Dose * F / CL) - 1| over 36 typical subjects = 0.0008 %
stopifnot(nrow(tv_chk) == nrow(arms) * 6, max(abs(tv_chk$relerr)) < 0.005)

The residual is numerical-integration error. Because the expectation is built from the published typical values, allometric exponents and both covariate coefficients, agreement to better than 0.5 % confirms the whole covariate model as well as the AUC bookkeeping state.

Replicate published figures

Figure 5 – simulated concentration-time profiles by arm

Assmus 2025 Fig 5 plots simulated median DBS concentrations for each arm against the in vitro IC90 (7.61 mg/L) and the historically accepted therapeutic range (2.5-5 mg/L in DBS). The paper’s headline reading of this figure is that only the twice-daily regimens keep concentrations above the IC90, while the lower 2.5 mg/L bound is exceeded for extended periods in every arm.

prof <- sim |>
  filter(time <= 24 * 7 * 10) |>
  group_by(arm, time) |>
  summarise(Q50 = median(Cc), .groups = "drop")

ggplot(prof, aes(time / 24, Q50)) +
  annotate("rect", xmin = -Inf, xmax = Inf, ymin = 2.5, ymax = 5,
           fill = "grey70", alpha = 0.4) +
  geom_hline(yintercept = 7.61, colour = "red") +
  geom_line(linewidth = 0.3) +
  facet_wrap(~arm, ncol = 2) +
  labs(x = "Time (days)", y = "Median benznidazole in DBS (mg/L)",
       caption = paste("Replicates Figure 5 of Assmus 2025. Red line = in vitro",
                       "IC90 (7.61 mg/L); grey band = 2.5-5 mg/L in DBS."))

The qualitative claim is directly testable. For each arm, the fraction of the on-treatment period during which the median profile sits above each threshold.

This fraction is only meaningful on an evenly-spaced time grid: any locally refined region would be over-weighted by mean(Q50 > threshold) relative to the rest of the course. The observation grid above is deliberately uniform at 1 h for the whole treatment period, with the fine 0.05 h sampling confined to the washout after the last dose, which this filter excludes.

last_dose <- events |> filter(evid == 1L) |>
  group_by(arm) |> summarise(last_dose = max(time), .groups = "drop")

thresholds <- prof |>
  left_join(last_dose, by = "arm") |>
  filter(time <= last_dose) |>
  group_by(arm) |>
  summarise(
    pct_above_25   = round(100 * mean(Q50 > 2.5)),
    pct_above_ic90 = round(100 * mean(Q50 > 7.61)),
    .groups = "drop"
  )

thresholds |>
  dplyr::rename("Arm" = arm,
                "% time > 2.5 mg/L" = pct_above_25,
                "% time > 7.61 mg/L (IC90)" = pct_above_ic90) |>
  knitr::kable(
    caption = paste("Median-profile time above target while on treatment.",
                    "Assmus 2025: only twice-daily dosing maintains levels above",
                    "the IC90; the 2.5 mg/L bound is exceeded in all arms."))
Median-profile time above target while on treatment. Assmus 2025: only twice-daily dosing maintains levels above the IC90; the 2.5 mg/L bound is exceeded in all arms.
Arm % time > 2.5 mg/L % time > 7.61 mg/L (IC90)
150 mg BID, 8 wk 100 89
150 mg BID, 4 wk 100 71
150 mg BID, 2 wk 100 89
150 mg QD, 4 wk 98 0
150 mg QD, 4 wk + E1224 82 0
300 mg weekly, 8 wk + E1224 21 0 

bid    <- grepl("BID", as.character(thresholds$arm))
weekly <- grepl("weekly", as.character(thresholds$arm))

stopifnot(
  # "only twice-daily dosing regimens maintained levels above the target
  # concentration". The separation is categorical, so it is asserted as such:
  # every non-twice-daily arm must spend EXACTLY zero time above the IC90.
  all(thresholds$pct_above_ic90[!bid] == 0),
  all(thresholds$pct_above_ic90[bid]   > 50),
  # "extended exposure above target concentrations was achieved in all
  # treatment arms" at the 2.5 mg/L bound. Every daily-dosing arm holds it
  # almost continuously; the weekly arm still clears it for about 1.5 days of
  # each 7-day cycle, which is the "extended exposure" the paper describes.
  all(thresholds$pct_above_25[!weekly] > 80),
  all(thresholds$pct_above_25[weekly]  > 15)
)

The separation is categorical rather than graded: no once-daily or weekly arm exceeds the IC90 at any point in its course, while every twice-daily arm spends most of its course above it. That is the paper’s claim reproduced as a binary, which is a far sharper test than any percentage comparison, and it is robust: the once-daily and weekly arms peak around 6.4-6.9 mg/L, comfortably short of the 7.61 mg/L IC90, so they cannot cross it.

The twice-daily percentage itself is deliberately asserted loosely, because it is a knife-edge statistic rather than a stable one. Over a 12 h interval the trough-to-peak ratio is exp(-kel * 12) = 0.61, and a median peak near 11 mg/L therefore puts the median trough near 6.6 mg/L – just below the 7.61 mg/L IC90. The median profile consequently crosses the threshold once per dosing interval, so the fraction of time above it lands around 80 % and shifts by several points with the virtual cohort’s weight draw. Reporting it is informative; pinning it to a tight tolerance would be pinning down sampling noise.

Against the lower, historically used 2.5 mg/L bound every arm qualifies, which is exactly the sensitivity to target choice that Assmus 2025 flags as the central uncertainty in interpreting these regimens: the in vivo target concentration is unknown.

Figure 4 – distribution of exposure metrics across arms 

by_id <- sim |>
  group_by(id, arm) |>
  summarise(cmax = max(Cc), aucinf = max(auc_central),
            p_qpcr = dplyr::last(pqpcr), .groups = "drop")

by_id |>
  select(arm, aucinf, cmax) |>
  pivot_longer(c(aucinf, cmax)) |>
  mutate(name = recode(name, aucinf = "AUCinf (mg*h/L)", cmax = "CMAX (mg/L)")) |>
  ggplot(aes(arm, value)) +
  geom_boxplot(outlier.size = 0.4) +
  facet_wrap(~name, scales = "free_x") +
  coord_flip() +
  labs(x = NULL, y = NULL,
       caption = "Replicates Figure 4a-b of Assmus 2025.")

Assmus 2025 reports that the standard 8-week arm has a 7.9-fold higher median AUCinf than the weekly arm, and that AUCinf is approximately dose-proportional across the 2-, 4- and 8-week monotherapy arms. S1 Table lets that first number be checked exactly: its pooled AUCinf medians are 11,336 and 1,428 mg*h/L, a ratio of 7.94.

These ratios are taken from the typical-value cohort rather than the stochastic one. Each arm draws its own covariates and etas, so a stochastic median carries a few percent of Monte-Carlo noise – enough to blur a claim of exact dose proportionality. On the typical-value ladder the ratios are noise-free, so the assertions below can be tight.

med <- tv_chk |> group_by(arm) |>
  summarise(auc = median(aucinf), .groups = "drop")
g <- function(a) med$auc[as.character(med$arm) == a]

ratios <- tibble(
  Comparison = c("8 wk BID / weekly + E1224",
                 "8 wk BID / 4 wk BID", "4 wk BID / 2 wk BID"),
  Simulated  = round(c(g("150 mg BID, 8 wk") / g("300 mg weekly, 8 wk + E1224"),
                       g("150 mg BID, 8 wk") / g("150 mg BID, 4 wk"),
                       g("150 mg BID, 4 wk") / g("150 mg BID, 2 wk")), 3),
  Published  = c("7.9 (S1 Table: 7.94)", "~2", "~2")
)
knitr::kable(ratios, caption = "Typical-value AUCinf ratios vs Assmus 2025 Results.")
Typical-value AUCinf ratios vs Assmus 2025 Results.
Comparison Simulated Published
8 wk BID / weekly + E1224 8.239 7.9 (S1 Table: 7.94)
8 wk BID / 4 wk BID 2.000 ~2
4 wk BID / 2 wk BID 2.000 ~2 

# Dose proportionality is exact by construction once absorption is complete and
# elimination is linear, so these two are pinned hard rather than loosely.
stopifnot(all(abs(ratios$Simulated[2:3] - 2) < 0.001))
# 8 wk vs weekly carries the fosravuconazole clearance effect as well as the
# 7-fold dose ratio: 7 * 1.177 = 8.24 structurally, against 7.9 observed.
stopifnot(abs(ratios$Simulated[1] - 7 * 1.177) < 0.01,
          abs(ratios$Simulated[1] - 7.94) < 0.5)

The 2-, 4- and 8-week monotherapy ratios are exactly 2.000, which is the signature of complete absorption and linear elimination. The 8-week-versus-weekly ratio is 8.24 structurally – the 7-fold total-dose ratio multiplied by the 17.7 % fosravuconazole clearance increase – against 7.94 in S1 Table. The remaining 4 % is adherence: the 8-week arm had the trial’s highest discontinuation rate (20 %) and lowest compliance (80 %), so its real-patient median AUCinf falls below the perfect-adherence value.

Sex effect

Assmus 2025 reports that AUCinf and CMAX were lower in men than in women, “but this effect was minor (< 20 %)”.

sex_eff <- by_id |>
  left_join(events |> distinct(id, SEXF), by = "id") |>
  group_by(SEXF) |>
  summarise(auc = median(aucinf), cmax = median(cmax), .groups = "drop")

sex_tab <- tibble(
  Metric = c("AUCinf", "CMAX"),
  `Men vs women (%)` = round(c(
    100 * (sex_eff$auc[sex_eff$SEXF == 0]  / sex_eff$auc[sex_eff$SEXF == 1]  - 1),
    100 * (sex_eff$cmax[sex_eff$SEXF == 0] / sex_eff$cmax[sex_eff$SEXF == 1] - 1)
  ), 1)
)
knitr::kable(sex_tab,
  caption = "Simulated sex difference (published: lower in men, magnitude < 20%).")
Simulated sex difference (published: lower in men, magnitude < 20%).
Metric Men vs women (%)
AUCinf -16.2
CMAX -9.2 
stopifnot(all(sex_tab[[2]] < 0), all(sex_tab[[2]] > -20))

PKNCA validation

S1 Table reports AUCinf, CMAX, TMAX and terminal half-life for every arm, so the comparison can be made arm by arm rather than against a pooled range. Two of those four quantities are whole-course descriptors and two are not, and matching each to the right simulation is what makes the comparison meaningful:

  • AUCinf and CMAX are cumulative / peak values over the entire treatment course, so they come from the full-regimen solve above.
  • TMAX is defined by S1 Table as “time after dose to reach the maximum concentration”. Because the maximum concentration occurs once the regimen has accumulated, this is a steady-state time-after-dose, not a single-dose one. It matters: the residual drug from previous doses is still declining when the new dose is absorbed, which pulls the peak earlier. A single 150 mg dose peaks at 2.3 h under this model, whereas the accumulated regimens peak 1.7-2.1 h after a dose – exactly the published range. Computing TMAX from a single dose would have manufactured a spurious ~25 % discrepancy.
  • Terminal half-life is regimen-independent, so it comes from a clean single-dose PKNCA run.
# Peak of the whole course, and how long after the most recent dose it occurred.
peak <- sim |>
  group_by(id, arm) |>
  slice_max(Cc, n = 1, with_ties = FALSE) |>
  ungroup() |>
  select(id, arm, tcmax = time)

wc <- peak |>
  left_join(events |> filter(evid == 1L) |> select(id, dtime = time),
            by = "id", relationship = "many-to-many") |>
  filter(dtime <= tcmax) |>
  group_by(id, arm) |>
  slice_max(dtime, n = 1, with_ties = FALSE) |>
  ungroup() |>
  mutate(tmax = tcmax - dtime) |>
  select(id, arm, tmax) |>
  left_join(by_id, by = c("id", "arm"))

stopifnot(nrow(wc) == N_PER_ARM * nrow(arms), all(wc$tmax > 0))

The single-dose run for half-life. Sub-LLOQ samples are dropped the way the paper’s own dataset was built (assay LLOQ 50 ng/mL, records censored at 120 h after dose), by truncating each profile after its last quantifiable sample – not by filtering on concentration, which would also delete the time-zero record and trigger PKNCA’s “AUC range starting before the first measurement” warning once per subject.

# Seeded independently of the multiple-dose block above, so that editing one
# cohort cannot silently shift the other's results.
set.seed(20250923)
rxode2::rxSetSeed(20250923)

LLOQ  <- 0.05   # mg/L = 50 ng/mL (Methods, Quantification of benznidazole)
# Single-dose cohort per arm, under the 200/arm cap. Sized at 150 rather than a
# token 50 because the arm-level MEDIANS are what the S1 Table comparison below
# asserts on, and MTT carries 61.6 % CV: at n = 50 the sampling error on a
# median TMAX is already several per cent, which would leave the 20 % tolerance
# hostage to which subjects happened to draw a long absorption.
N_SD  <- 150L
sd_times <- sort(unique(c(seq(0, 12, by = 0.1), seq(12, 120, by = 1))))

sd_events <- do.call(bind_rows, lapply(seq_len(nrow(arms)), function(i) {
  a <- arms[i, ]
  subj <- tibble(
    id   = 10000L + (i - 1L) * N_SD + seq_len(N_SD),
    WT   = rtnorm(N_SD, 65, 7.5, 50, 80),
    SEXF = as.numeric(seq_len(N_SD) > round(0.32 * N_SD)),
    CONMED_FOSRAVUCONAZOLE = a$combo,
    arm  = a$arm
  )
  bind_rows(
    subj |> mutate(time = 0, amt = 150, evid = 1L, cmt = "depot"),
    subj |> tidyr::crossing(time = sd_times) |>
      mutate(amt = NA_real_, evid = 0L, cmt = "central")
  ) |> arrange(id, time, desc(evid))
}))

sd_sim <- rxode2::rxSolve(pk, events = sd_events, keep = c("arm")) |>
  as.data.frame() |>
  filter(!is.na(Cc))
stopifnot(all(sd_sim$Cc >= 0),
          dplyr::n_distinct(sd_sim$id) == N_SD * nrow(arms))
sim_nca <- sd_sim |>
  group_by(id) |>
  filter(any(Cc >= LLOQ), time <= max(time[Cc >= LLOQ])) |>
  ungroup() |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, arm)

# Guarantee a time = 0 record per (id, arm); pre-dose Cc = 0 for oral dosing.
sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> dplyr::distinct(id, arm) |> dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, arm, time, .keep_all = TRUE) |>
  dplyr::arrange(id, arm, time)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | arm + id)

dose_df <- sd_events |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, arm)
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | arm + id)

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

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

Comparison against published NCA 

published <- tibble::tribble(
  ~arm,                          ~aucinf.obs, ~cmax, ~tmax, ~half.life,
  "150 mg BID, 8 wk",                  11336,  11.7,  1.80,       18.0,
  "150 mg BID, 4 wk",                   6275,  11.8,  1.77,       17.1,
  "150 mg BID, 2 wk",                   3058,  11.9,  1.74,       16.2,
  "150 mg QD, 4 wk",                    2837,  6.45,  1.99,       16.4,
  "150 mg QD, 4 wk + E1224",            2601,  6.17,  1.97,       14.2,
  "300 mg weekly, 8 wk + E1224",        1428,  6.92,  1.91,       13.8
) |> as.data.frame()

sim_long <- bind_rows(
  wc |> transmute(arm = as.character(arm), PPTESTCD = "aucinf.obs", PPORRES = aucinf),
  wc |> transmute(arm = as.character(arm), PPTESTCD = "cmax",       PPORRES = cmax),
  wc |> transmute(arm = as.character(arm), PPTESTCD = "tmax",       PPORRES = tmax),
  nca_res$result |>
    dplyr::filter(PPTESTCD == "half.life") |>
    transmute(arm = as.character(arm), PPTESTCD, PPORRES)
) |> as.data.frame()

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated     = sim_long,
  reference     = published,
  by            = "arm",
  units         = c(aucinf.obs = "mg*h/L", cmax = "mg/L",
                    tmax = "h", half.life = "h"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = paste("Simulated vs Assmus 2025 S1 Table (pooled column), by arm.",
                  "* differs by >20%.")
)
Simulated vs Assmus 2025 S1 Table (pooled column), by arm. * differs by >20%.
NCA parameter arm Reference Simulated % diff
Cmax (mg/L) 150 mg BID, 8 wk 11.7 11.3 -3.8%
Cmax (mg/L) 150 mg BID, 4 wk 11.8 10.7 -9.5%
Cmax (mg/L) 150 mg BID, 2 wk 11.9 11.5 -3.0%
Cmax (mg/L) 150 mg QD, 4 wk 6.45 6.99 +8.4%
Cmax (mg/L) 150 mg QD, 4 wk + E1224 6.17 6.01 -2.5%
Cmax (mg/L) 300 mg weekly, 8 wk + E1224 6.92 6.5 -6.1%
Tmax (h) 150 mg BID, 8 wk 1.8 1.7 -5.6%
Tmax (h) 150 mg BID, 4 wk 1.77 1.7 -4.0%
Tmax (h) 150 mg BID, 2 wk 1.74 1.85 +6.3%
Tmax (h) 150 mg QD, 4 wk 1.99 2.1 +5.5%
Tmax (h) 150 mg QD, 4 wk + E1224 1.97 2.02 +2.8%
Tmax (h) 300 mg weekly, 8 wk + E1224 1.91 2.12 +11.3%
AUC0-∞ (obs) (mg*h/L) 150 mg BID, 8 wk 11300 12800 +12.5%
AUC0-∞ (obs) (mg*h/L) 150 mg BID, 4 wk 6280 5860 -6.6%
AUC0-∞ (obs) (mg*h/L) 150 mg BID, 2 wk 3060 3330 +8.9%
AUC0-∞ (obs) (mg*h/L) 150 mg QD, 4 wk 2840 3250 +14.4%
AUC0-∞ (obs) (mg*h/L) 150 mg QD, 4 wk + E1224 2600 2610 +0.5%
AUC0-∞ (obs) (mg*h/L) 300 mg weekly, 8 wk + E1224 1430 1510 +5.5%
t½ (h) 150 mg BID, 8 wk 18 17.3 -3.9%
t½ (h) 150 mg BID, 4 wk 17.1 17.3 +0.9%
t½ (h) 150 mg BID, 2 wk 16.2 16.1 -0.8%
t½ (h) 150 mg QD, 4 wk 16.4 16.9 +3.0%
t½ (h) 150 mg QD, 4 wk + E1224 14.2 14.3 +0.9%
t½ (h) 300 mg weekly, 8 wk + E1224 13.8 14.6 +6.1%

Every one of the 24 comparisons is asserted, not merely displayed:

chk <- sim_long |>
  group_by(arm, PPTESTCD) |>
  summarise(sim = median(PPORRES), .groups = "drop") |>
  left_join(
    published |> pivot_longer(-arm, names_to = "PPTESTCD", values_to = "ref"),
    by = c("arm", "PPTESTCD")
  ) |>
  mutate(pct = 100 * (sim / ref - 1))

stopifnot(nrow(chk) == 24, !anyNA(chk$pct))
cat(sprintf("max |%% difference| across %d arm x parameter comparisons = %.1f %%\n",
            nrow(chk), max(abs(chk$pct))))
#> max |% difference| across 24 arm x parameter comparisons = 14.4 %
stopifnot(max(abs(chk$pct)) < 20)

The residual differences run in a consistent direction and have a single cause. Simulated AUCinf is a few per cent high in the arms with the poorest adherence – most visibly the 8-week arm (20 % discontinuation, 80 % compliance) and the 4-week 150 mg once-daily arm – because this simulation assumes every scheduled dose was taken, whereas S1 Table is derived from each patient’s actual dosing record. Simulated CMAX in the twice-daily arms runs a few per cent low for a different reason: the virtual cohort is a mixture in which 32 % of subjects are men carrying the 12.9 % lower relative bioavailability, and weight spans the whole 50-80 kg window. Nothing here was tuned; the worst single comparison is well inside the 20 % tolerance.

Exposure-response: qPCR positivity

Figure 8a – probability of a positive qPCR at a single follow-up visit 

b0    <- 3.496
b_auc <- log(0.9995)
b_ct  <- log(0.9021)

auc_grid <- tibble(aucinf = seq(0, 18000, length.out = 400)) |>
  mutate(p = plogis(b0 + b_auc * aucinf + b_ct * 37.7))

ggplot(auc_grid, aes(aucinf, p)) +
  geom_line() +
  geom_point(data = by_id, aes(aucinf, p_qpcr, colour = arm),
             size = 0.6, alpha = 0.6) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = "AUCinf (mg*h/L)",
       y = "Probability (qPCR+ at a single follow-up)",
       colour = NULL,
       caption = paste("Replicates Figure 8a of Assmus 2025. Curve drawn at the",
                       "median baseline Ct of 37.7; points are simulated subjects."))

The curve is reproduced by construction for the intercept, which was recovered from this very figure (see Assumptions). What the figure tests independently is where each arm’s simulated AUCinf distribution falls, which comes entirely from the PK layer:

by_id |>
  group_by(arm) |>
  summarise(auc = round(median(aucinf)), p = round(median(p_qpcr), 3),
            .groups = "drop") |>
  dplyr::rename("Arm" = arm, "Median AUCinf (mg*h/L)" = auc,
                "Median P(qPCR+ per visit)" = p) |>
  knitr::kable(caption = paste("Simulated exposure and predicted single-visit",
                               "qPCR positivity by arm."))
Simulated exposure and predicted single-visit qPCR positivity by arm.
Arm Median AUCinf (mg*h/L) Median P(qPCR+ per visit)
150 mg BID, 8 wk 12750 0.001
150 mg BID, 4 wk 5864 0.036
150 mg BID, 2 wk 3329 0.116
150 mg QD, 4 wk 3246 0.119
150 mg QD, 4 wk + E1224 2614 0.159
300 mg weekly, 8 wk + E1224 1506 0.254

Reading these against Assmus 2025 Fig 8a: the 8-week twice-daily arm sits at the right-hand end of the x-axis (10,000-15,000 mg*h/L) where the fitted probability has flattened near zero, while the weekly arm sits near 1,500-2,000 where the curve is still descending. That ordering, and a near-zero predicted probability for the longer arms, is what the paper shows.

Equation 3 – probability of at least one positive follow-up

Assmus 2025 Equation 3 converts the per-visit probability into the probability of at least one positive result across N follow-up visits, p* = 1 - (1 - p)^N. It is post-processing on p rather than part of the model, so it is applied here rather than in the model file.

by_id |>
  group_by(arm) |>
  summarise(p = median(p_qpcr), .groups = "drop") |>
  tidyr::crossing(N = c(1L, 3L, 5L, 7L)) |>
  mutate(p_star = round(1 - (1 - p)^N, 3)) |>
  select(-p) |>
  pivot_wider(names_from = N, values_from = p_star, names_prefix = "N = ") |>
  dplyr::rename("Arm" = arm) |>
  knitr::kable(caption = paste("P(at least one qPCR-positive follow-up) by",
                               "number of visits (Assmus 2025 Equation 3)."))
P(at least one qPCR-positive follow-up) by number of visits (Assmus 2025 Equation 3).
Arm N = 1 N = 3 N = 5 N = 7
150 mg BID, 8 wk 0.001 0.004 0.006 0.008
150 mg BID, 4 wk 0.036 0.105 0.168 0.227
150 mg BID, 2 wk 0.116 0.309 0.460 0.578
150 mg QD, 4 wk 0.119 0.315 0.468 0.587
150 mg QD, 4 wk + E1224 0.159 0.405 0.579 0.702
300 mg weekly, 8 wk + E1224 0.254 0.585 0.769 0.871

The number of follow-up visits varied by arm in BENDITA – shorter arms had more post-treatment visits – which is precisely why the paper adopted the proportion-based endpoint rather than a binary “any positive” endpoint.

What the paper concluded

The relationship encoded here is the placebo-inclusive fit. Assmus 2025 is explicit that it is driven almost entirely by the contrast between no drug and any drug: after excluding the placebo arm and one subject who took only four doses (n = 170), the AUCinf odds ratio moved to 0.9999 [0.9998, 1.0000], p = 0.282, and no exposure metric remained a significant predictor. Only the baseline Ct term survived, with an odds ratio between 0.85 and 0.88 across every model. This model should not be read as evidence of a dose-response relationship within the active regimens.

The per-arm predictions above make that concrete. The model assigns the weekly combination arm a per-visit positivity of about 0.26, rising to 0.87 over seven follow-up visits under Equation 3 – yet in the trial that arm behaved like the 4- and 8-week arms, with individual qPCR positivity at or below 20 %. The over-prediction is not a transcription error; it is the substantive finding. A relationship calibrated mostly on the placebo-versus-active contrast extrapolates poorly into the low-exposure end of the active range, which is why the paper concludes that time on treatment, rather than cumulative exposure alone, may drive the parasitological response, and why the once-weekly regimen performed better than its exposure would predict.

tibble::tribble(
  ~Predictor,                ~`OR, mITT incl. placebo (n=201)`, ~`OR, excl. placebo (n=170)`,
  "AUCinf (per mg*h/L)",     "0.9995 [0.9994, 0.9997] ***",     "0.9999 [0.9998, 1.0000]",
  "log AUCinf",              "0.236 [0.194, 0.288] ***",        "0.381 [0.117, 1.176]",
  "CMAX (per mg/L)",         "0.640 [0.581, 0.706] ***",        "0.966 [0.863, 1.073]",
  "Weeks of treatment",      "0.449 [0.366, 0.550] ***",        "0.854 [0.709, 1.009]",
  "Baseline Ct (per cycle)", "0.9021 [0.8024, 1.0142]",         "0.875 [0.765, 1.006]"
) |>
  knitr::kable(caption = paste("Assmus 2025 Table 3; the Ct row is taken from",
                               "the AUCinf column. *** p < 0.001. Only the",
                               "AUCinf row is encoded in the model file."))
Assmus 2025 Table 3; the Ct row is taken from the AUCinf column. *** p < 0.001. Only the AUCinf row is encoded in the model file.
Predictor OR, mITT incl. placebo (n=201) OR, excl. placebo (n=170)
AUCinf (per mg*h/L) 0.9995 [0.9994, 0.9997] *** 0.9999 [0.9998, 1.0000]
log AUCinf 0.236 [0.194, 0.288] *** 0.381 [0.117, 1.176]
CMAX (per mg/L) 0.640 [0.581, 0.706] *** 0.966 [0.863, 1.073]
Weeks of treatment 0.449 [0.366, 0.550] *** 0.854 [0.709, 1.009]
Baseline Ct (per cycle) 0.9021 [0.8024, 1.0142] 0.875 [0.765, 1.006]

Assumptions and deviations

  • b0_qpcr is figure-derived, not published. Assmus 2025 reports the two slopes of Equation 2 as odds ratios but never reports the intercept, without which no absolute probability can be computed. It was recovered by digitising the solid “Including placebo” fitted curve of Fig 8a (the caption calls it grey; it is the darker of the two fitted lines, and is distinct from the black placebo-arm observed points): the linear predictor at AUCinf = 0 and median baseline Ct is -0.388, i.e. p = 0.404 for an untreated patient of median parasite load. The digitisation validates itself – refitting the traced curve with a free slope returns an exposure odds ratio of 0.99955 against the published 0.9995, agreeing to four significant figures, which confirms the axis calibration. Residual uncertainty on the intercept is about +/- 0.05 logit units (p at zero exposure between 0.40 and 0.42 depending on the fitting window). Fig 8 does not state which Ct value its curves are drawn at; the median is assumed, following Methods (“utilizing the median Ct value at baseline”). Every other parameter in both model files comes from the paper’s tables.
  • Beta binomial overdispersion is not encoded. The paper fits the PD layer in glmmTMB with a beta binomial distribution parameterised per Morris 1997 but does not report the estimated overdispersion coefficient. Rather than invent one, the exposure-response layer is deterministic: pqpcr is a derived output column with no residual-error term. The PK layer’s residual error is unaffected.
  • Only AUCinf is encoded as the exposure driver. Assmus 2025 screens eight univariate exposure predictors against a Ct-only base model (Table 3: AUCinf, log AUCinf, CMAX, time above 3 and 6 mg/L, and three treatment-duration definitions – total days of dosing, weeks, and total duration in days). AUCinf is the paper’s leading metric, the subject of Fig 8a, and the only one that follows directly from integrating the PK layer. Representatives of the other families – CMAX (CMAX_BZN), time above target (T_ABOVE_TARGET) and weeks of treatment (DUR_BZN_WEEKS) – are recorded in covariatesDataExcluded with their published odds ratios but are not encoded as model states; the two remaining duration definitions are omitted as near-duplicates of the weeks metric. Time above target in particular is described by the paper itself as highly sensitive to a threshold that “remains unknown in vivo”.
  • CONMED_FOSRAVUCONAZOLE is an arm-level, time-fixed indicator, matching how the paper models it. In BENDITA it is confounded with benznidazole dose and schedule, since only the two combination arms received fosravuconazole.
  • Sex reference category. The register canonical SEXF uses 0 (male) as its reference, but Assmus 2025 reports typical values for a female subject. The published -12.9 % male effect is therefore applied as (1 + e_sexf_fdepot * (1 - SEXF)) so the Table 2 structural estimates stay verbatim – the same construction used by Bajaj_2017_nivolumab.R, Wada_2023_sparsentan.R and Birgersson_2019_artesunate.R.
  • Covariate distributions are assumed. Body weight is drawn from a normal truncated to the 50-80 kg enrolment window with a 65 kg median, sex is set to 68 % female (Table 1, PK analysis column), and baseline Ct is drawn from a normal truncated to 30.2-40.0 with the assay ceiling at 40. The paper reports medians and ranges only, not distributional shapes. Assmus 2025 Fig 5 fixes weight at exactly 65 kg for its own simulation; sampling weight here widens the percentiles slightly but leaves the medians essentially unchanged.
  • Perfect adherence is assumed. BENDITA had 80-97 % arm-level compliance, dose interruptions in up to 34 % of one arm, and five subjects excluded for possible treatment misallocation; none of that is modelled here. This is the main reason simulated AUCinf sits a few per cent above S1 Table in the lowest-compliance arms, and it is not a transcription error.
  • S1 Table is used as the per-arm NCA reference. It was not in the original source bundle; it was retrieved from the PLoS supporting-information endpoint during this extraction and is stored next to the article PDF as PMC12510642_S1_Table_secondary_pk.docx (along with S1 Text, S2 Text and S3 Table). Its values are themselves model-derived – “secondary pharmacokinetic parameter estimates, based on the final population pharmacokinetic model” – so this comparison tests whether this reimplementation reproduces the authors’ own model, which is exactly the question a library extraction needs to answer.
  • S3 Table (refined mITT, n = 186) is not extracted. The paper frames it as a sensitivity analysis, and the standing convention is to carry the authors’ final model rather than its robustness checks. Its odds ratios are close to the primary fit (AUCinf 1.000 [0.999, 1.000], Ct 0.902 [0.797, 1.020]).
  • auc_central is a bookkeeping state, declared via paper_specific_compartments, not a biological compartment.

Registry entries 

nlmixr2lib::modeldb |>
  dplyr::filter(grepl("^Assmus_2025", name)) |>
  dplyr::select(name, vignette) |>
  knitr::kable()
name vignette
Assmus_2025_benznidazole Assmus_2025_benznidazole
Assmus_2025_benznidazole_mouse Assmus_2025_benznidazole_mouse
Assmus_2025_benznidazole_qpcr Assmus_2025_benznidazole