---
title: "Tilmicosin PK/PD against Pasteurella multocida (Chen 2023)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Tilmicosin PK/PD against Pasteurella multocida (Chen 2023)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5)
library(nlmixr2lib)
library(rxode2)
library(dplyr)
library(ggplot2)
```

# Model and source

* Citation: Chen Y, Ji X, Zhang S, Wang W, Zhang H, Ding H. (2023).
  Pharmacokinetic/pharmacodynamic integration of tilmicosin against
  *Pasteurella multocida* in a piglet tissue cage model.
  Frontiers in Veterinary Science 10:1260990.
* Article: <https://doi.org/10.3389/fvets.2023.1260990>

This vignette validates the packaged model `Chen_2023_tilmicosin` against the
published Chen 2023 results.

Chen 2023 implanted silicone tissue cages subcutaneously in piglets, infected
them with *P. multocida* serovar D:7, dosed tilmicosin orally at 30, 40, 50 or
60 mg/kg once daily for three days, and related the drug exposure inside the
cage to the change in bacterial count. The PK/PD integration is a **sigmoidal
Emax model on a PK/PD index** (Section 2.6):

$$E = E_0 + \frac{(E_{max} - E_0) \times C_e^N}{EC_{50}^N + C_e^N}$$

where `E` is the change in bacterial count (log10 CFU/mL) over a 24 h dosing
interval, `E0` is the corresponding change in the untreated control, and `Ce`
is a PK/PD index. Chen 2023 evaluated three candidate indices and found
`AUC24h/MIC` to correlate best with effect (R^2 = 0.92, versus 0.90 for
`Cmax/MIC` and 0.83 for `%T>MIC`); Table 2 reports parameters for the
`AUC24h/MIC` fit only, which is what the packaged model encodes.

Two features of this source shape the packaged model:

* **There is no PK model to package.** Chen 2023 analysed the tissue-cage-fluid
  concentrations non-compartmentally in WinNonlin 6.1 (Section 2.5) and reports
  only NCA summaries (Table 1). Drug exposure therefore enters the packaged
  model as an externally supplied per-interval covariate `AUC_TILM`, exactly as
  the paper's own PK/PD integration consumed it.
* **`E` is a positive reduction magnitude.** In the printed equation `E` rises
  from `E0` at zero exposure to `Emax` at saturating exposure, and both Table 2
  values are positive. The sign convention is confirmed numerically below: the
  equation reproduces the paper's own reported 1/3-log and 3/4-log thresholds.

# Population

Ten healthy crossbred piglets (Duroc x Landrace x Yorkshire), 25-30 kg, were
randomised to four treatment groups and one infected-but-untreated control
group, with two piglets and four tissue cages per group and one male and one
female per group (Sections 2.2 and 2.4). Tissue cages were food-grade silicone
tubes (65 mm long, 13 mm inner / 18 mm outer diameter, twelve holes at each
end) implanted subcutaneously on both sides of the neck and allowed four weeks
to recover before use.

Each sterile cage was inoculated with 1 mL of a ~1.5 x 10^8 CFU/mL
*P. multocida* suspension; after 48 h incubation, only cages reaching
~10^7 CFU/mL entered the experiment. The challenge strain was C44-15
(serovar D:7) from the China Institute of Veterinary Drug Control. The
tilmicosin MIC was **0.25 ug/mL in both tissue-cage fluid and tryptic soy
broth** by CLSI (2013) microdilution (Section 3.1) -- the paper notes that
measuring the MIC in the actual biological matrix matters, and here the two
media agreed.

Bacteria were counted at 24 h after each administration and at 48 and 72 h
after the final dose; the untreated control stayed around 10^7 CFU/mL
throughout (Section 3.3).

```{r population-table}
pop <- rxode2::rxode(readModelDb("Chen_2023_tilmicosin"))$population
data.frame(
  field = c("species", "n_subjects", "weight_range", "organism", "dose_range"),
  value = c(pop$species, pop$n_subjects, pop$weight_range,
            pop$organism, pop$dose_range),
  stringsAsFactors = FALSE
) |>
  knitr::kable(caption = "Selected population metadata from the packaged model.")
```

# Source trace

Every `ini()` parameter in the model file carries an in-file comment pointing
to its source location. The table below collects them in one place.

| Equation / parameter | Value | Source |
|---|---|---|
| Sigmoidal Emax PK/PD-index equation | `E = E0 + (Emax - E0) * Ce^N / (EC50^N + Ce^N)` | Section 2.6 |
| PK/PD index `Ce` | `AUC_TILM / mic` (i.e. AUC24h/MIC) | Section 2.6; best index per Section 3.4 (R^2 = 0.92) |
| `Emax` (log10 CFU/mL) | 1.09 | Table 2 |
| `E0` (log10 CFU/mL) | 0.003 | Table 2 |
| `EC50` (h) | 26.66 | Table 2 |
| `N` (Hill slope, unitless) | 2.69 | Table 2 |
| `mic` (ug/mL) | 0.25 (FIXED, measured) | Section 3.1 |
| `log10_cfu0` (log10 CFU/mL) | 7 (FIXED, design) | Sections 2.2 and 3.3 |
| Bacterial density ODE | `d/dt(bact) = -ln(10) * (E / 24) * bact` | encoding of the Section 2.6 per-interval model (see below) |
| Residual SD on log10 CFU/mL | FIXED 0 | not reported (only R^2 given, Section 3.4) |
| `AUC_TILM` covariate values | Table 1 (dosing days 1-3) + Section 3.2 (24-48 h after final dose) | Table 1, Section 3.2 |

The paper fitted only the **per-interval** reduction and counted bacteria only
at 24 h interval boundaries. The packaged ODE spreads each interval's reduction
uniformly across that interval, so `log10(bact)` falls by exactly `E` over each
24 h interval and the trajectory matches the paper's model at every time point
the paper actually observed. This is verified numerically in the encoding
checks below.

# Study exposures (model input)

Chen 2023 Table 1 gives the NCA `AUC0-24h` for each dose group on each of the
three dosing days. Section 3.2 additionally gives the AUC over the 24-48 h
window after the final dose. Together these are the four consecutive 24 h
exposures the model needs to span the paper's full observation period through
96 h.

```{r exposures}
MIC <- 0.25  # ug/mL; Chen 2023 Section 3.1

exposures <- tibble::tribble(
  ~dose_mg_kg, ~interval, ~auc_ug_h_mL, ~source,
  30L, 1L,  3.57, "Table 1",
  30L, 2L,  4.87, "Table 1",
  30L, 3L,  6.03, "Table 1",
  30L, 4L,  6.31, "Section 3.2",
  40L, 1L,  5.40, "Table 1",
  40L, 2L,  6.92, "Table 1",
  40L, 3L,  8.61, "Table 1",
  40L, 4L,  8.35, "Section 3.2",
  50L, 1L,  7.49, "Table 1",
  50L, 2L, 10.47, "Table 1",
  50L, 3L, 12.30, "Table 1",
  50L, 4L, 11.69, "Section 3.2",
  60L, 1L,  8.96, "Table 1",
  60L, 2L, 11.33, "Table 1",
  60L, 3L, 13.65, "Table 1",
  60L, 4L, 12.87, "Section 3.2"
) |>
  mutate(aucmic_h = auc_ug_h_mL / MIC)

exposures |>
  mutate(across(c(auc_ug_h_mL, aucmic_h), ~ round(.x, 2))) |>
  rename("Dose (mg/kg)" = dose_mg_kg, "24 h interval" = interval,
         "AUC0-24h (ug*h/mL)" = auc_ug_h_mL, "Source" = source,
         "AUC24h/MIC (h)" = aucmic_h) |>
  knitr::kable(caption = "Reproduces Table 1 (AUC column) of Chen 2023, plus the 24-48 h post-final-dose exposures from Section 3.2. MIC = 0.25 ug/mL.")
```

# Setup

```{r model-load}
mod <- rxode2::rxode(readModelDb("Chen_2023_tilmicosin"))

# Build an observation-only event table for one arm. The model has no dosing
# compartment: tilmicosin exposure enters solely through the piecewise-constant
# AUC_TILM covariate (one value per 24 h interval), so every record is an
# observation (evid = 0).
build_arm_events <- function(auc_by_interval, t_end = 96, by_h = 0.5) {
  times <- seq(0, t_end, by = by_h)
  # Interval index 1..n for [0,24), [24,48), ...; the closing time point of the
  # final interval stays inside that interval.
  idx <- pmin(floor(times / 24) + 1L, length(auc_by_interval))
  data.frame(
    id       = 1L,
    time     = times,
    AUC_TILM = auc_by_interval[idx],
    amt      = 0,
    evid     = 0
  )
}

simulate_arm <- function(auc_by_interval, label) {
  ev  <- build_arm_events(auc_by_interval)
  sim <- as.data.frame(rxode2::rxSolve(mod, events = ev, keep = "AUC_TILM"))
  sim$arm <- label
  sim
}

arm_auc <- split(exposures$auc_ug_h_mL, exposures$dose_mg_kg)
arm_auc <- arm_auc[order(as.integer(names(arm_auc)))]
```

# Replicate Figure 3: effect versus AUC24h/MIC

Chen 2023 Figure 3 plots the sigmoidal Emax relationship between each PK/PD
index and the antibacterial effect. The panel for `AUC24h/MIC` (the best index,
R^2 = 0.92) is reproduced below from the Table 2 point estimates, with the
paper's three reported effect thresholds marked and the 16 studied exposures
from the table above rug-plotted along the x-axis.

```{r figure-3, fig.height = 4.6}
pars <- as.list(rxode2::rxode(readModelDb("Chen_2023_tilmicosin"))$theta)
e0   <- pars$e0
emax <- exp(pars$lemax)
ec50 <- exp(pars$lec50)
hill <- exp(pars$lhill)

# Chen 2023 Section 2.6 equation, evaluated directly from the packaged
# parameters (same algebra the model() block uses).
effect_of_index <- function(ce) e0 + (emax - e0) * ce^hill / (ec50^hill + ce^hill)

curve_df <- data.frame(aucmic_h = seq(0, 80, length.out = 401)) |>
  mutate(effect = effect_of_index(aucmic_h))

thresholds <- data.frame(
  label    = c("1/3-log", "1/2-log", "3/4-log"),
  effect   = c(1 / 3, 1 / 2, 3 / 4)
) |>
  mutate(aucmic_h = vapply(effect, function(tg) {
    stats::uniroot(function(ce) effect_of_index(ce) - tg,
                   interval = c(1e-3, 500), tol = 1e-10)$root
  }, numeric(1)))

ggplot(curve_df, aes(aucmic_h, effect)) +
  geom_line(linewidth = 0.8, colour = "steelblue") +
  geom_hline(yintercept = emax, linetype = "dotted", colour = "grey50") +
  geom_segment(data = thresholds,
               aes(x = aucmic_h, xend = aucmic_h, y = 0, yend = effect),
               linetype = "dashed", colour = "firebrick") +
  geom_point(data = thresholds, aes(aucmic_h, effect),
             colour = "firebrick", size = 2) +
  geom_text(data = thresholds,
            aes(aucmic_h, effect, label = sprintf("%s (%.2f h)", label, aucmic_h)),
            hjust = -0.08, vjust = 1.6, size = 3, colour = "firebrick") +
  geom_rug(data = exposures, aes(aucmic_h, y = NULL), sides = "b", alpha = 0.6) +
  labs(x = "AUC24h/MIC (h)",
       y = expression("E = "*log[10]~"CFU/mL reduction per 24 h interval"),
       title = "Replicates Figure 3 (AUC24h/MIC panel) of Chen 2023",
       caption = paste("Dotted line = Emax (1.09). Rug = the 16 studied 24 h exposures.",
                       "Thresholds back-solved from the Table 2 parameters.")) +
  coord_cartesian(xlim = c(0, 80), ylim = c(0, 1.15))
```

# Replicate Figure 2: time-kill curves

Figure 2 of Chen 2023 shows the bacterial count in tissue-cage fluid over the
three-day dosing course and the follow-up period. The packaged model reproduces
the trajectory by integrating the per-interval reduction; the untreated control
arm is simulated by setting `AUC_TILM = 0`, which collapses the sigmoid term
and leaves the near-zero control change `E0 = 0.003`.

```{r figure-2, fig.height = 5}
sims <- bind_rows(
  simulate_arm(rep(0, 4), "Control (untreated)"),
  do.call(rbind, lapply(names(arm_auc), function(d) {
    simulate_arm(arm_auc[[d]], sprintf("%s mg/kg", d))
  }))
)
sims$arm <- factor(sims$arm,
                   levels = c("Control (untreated)",
                              sprintf("%s mg/kg", names(arm_auc))))

ggplot(sims, aes(time, Cc, colour = arm)) +
  geom_vline(xintercept = c(0, 24, 48), linetype = "dotted", colour = "grey65") +
  geom_hline(yintercept = 4, linetype = "dashed", colour = "grey40") +
  geom_line(linewidth = 0.8) +
  scale_x_continuous(breaks = seq(0, 96, by = 12)) +
  scale_y_continuous(breaks = seq(3, 8, by = 1)) +
  labs(x = "Time (h)",
       y = expression(log[10]~"CFU/mL"),
       colour = NULL,
       title = "Replicates Figure 2 of Chen 2023 (typical-value)",
       caption = paste("Dotted verticals = the three once-daily doses.",
                       "Dashed horizontal = 3-log10 reduction from the 10^7 CFU/mL baseline.")) +
  theme(legend.position = "bottom")
```

The simulated curves reproduce the qualitative findings of Section 3.3 and the
Discussion: the 30 mg/kg arm falls well short of a 3-log10 reduction, the 50
and 60 mg/kg arms cross the bactericidal 3-log10 threshold, and the 50 and 60
mg/kg curves are nearly superimposed -- the paper's observation that the two
highest doses were not significantly different because the maximum effect had
already been reached at 50 mg/kg.

Because the exposures over the 24-48 h window after the final dose remain high
(AUC24h/MIC of 46.8 h and 51.5 h at 50 and 60 mg/kg), the curves keep falling
after dosing stops -- the slow-elimination effect the Discussion highlights.

# Validation: cumulative bacterial reduction

Section 3.3 reports the total reduction in bacterial count for each dose group.
The Discussion states these accumulate "within 96 h after administration", i.e.
across four consecutive 24 h intervals: the three dosing intervals plus the
24-48 h window after the final dose. The comparison below uses exactly that
window.

```{r validate-reduction}
reported_total <- c(`30` = 1.48, `40` = 2.82, `50` = 3.39, `60` = 3.52)
reported_sd    <- c(`30` = 0.13, `40` = 0.10, `50` = 0.11, `60` = 0.15)

reduction_tbl <- do.call(rbind, lapply(names(arm_auc), function(d) {
  s <- simulate_arm(arm_auc[[d]], d)
  data.frame(
    dose      = as.integer(d),
    simulated = s$Cc[s$time == 0] - s$Cc[s$time == 96],
    reported  = reported_total[[d]],
    rep_sd    = reported_sd[[d]],
    stringsAsFactors = FALSE
  )
})) |>
  mutate(
    pct_diff = 100 * (simulated - reported) / reported,
    flag     = ifelse(abs(pct_diff) > 20, "*", "")
  )

reduction_tbl |>
  mutate(simulated = round(simulated, 3),
         reported  = sprintf("%.2f +/- %.2f", reported, rep_sd),
         pct_diff  = sprintf("%+.1f%%", pct_diff)) |>
  select(-rep_sd) |>
  rename("Dose (mg/kg)" = dose,
         "Simulated reduction at 96 h" = simulated,
         "Chen 2023 Section 3.3" = reported,
         "Difference" = pct_diff,
         ">20% flag" = flag) |>
  knitr::kable(caption = "Cumulative log10 CFU/mL reduction over four 24 h intervals (0-96 h) versus the published totals.")
```

Three of the four dose groups agree closely -- 30 mg/kg (1.484 vs 1.48) and
60 mg/kg (3.516 vs 3.52) reproduce the published totals to within 0.01 log10,
and 50 mg/kg agrees to 0.11 log10 (3.2%). No group exceeds the 20% flag
threshold.

The 40 mg/kg group is the largest discrepancy: the model predicts 2.40 versus
the reported 2.82 log10 (-14.8%). This is a property of the published
parameters and exposures, not of the packaging -- the Table 1 exposures for
40 mg/kg, pushed through the Table 2 sigmoid, cannot produce 2.82 log10. Note
that the reported totals are not monotone in a way the exposures explain: the
30-to-40 mg/kg step in reported reduction (1.48 to 2.82, +1.34) is far larger
than the 50-to-60 mg/kg step (3.39 to 3.52, +0.13), while the corresponding
exposure increments are comparable. No parameter was tuned to close this gap.

# Validation: PK/PD index thresholds

Table 2 also reports the `AUC24h/MIC` magnitude required for a 1/3-log,
1/2-log and 3/4-log reduction per 24 h interval. Back-solving those thresholds
from the Table 2 parameters is a direct internal-consistency check on the
encoded equation.

```{r validate-thresholds}
threshold_tbl <- thresholds |>
  mutate(
    table2   = c(19.64, 33.45, 35.64),
    abstract = c(19.65, 23.86, 35.77),
    e_at_table2 = effect_of_index(table2)
  ) |>
  transmute(
    Reduction               = label,
    `Model-implied AUC24h/MIC (h)` = round(aucmic_h, 2),
    `Table 2 (h)`           = table2,
    `Abstract (h)`          = abstract,
    `E at the Table 2 value` = round(e_at_table2, 3),
    `Target E`              = round(effect, 3)
  )

knitr::kable(threshold_tbl,
             caption = "Thresholds back-solved from the packaged Table 2 parameters versus the values printed in Chen 2023.")
```

The 1/3-log and 3/4-log thresholds are reproduced: the model implies 19.59 h
and 35.72 h against printed values of 19.64 h and 35.64 h (Table 2) or 19.65 h
and 35.77 h (Abstract) -- agreement to 0.3% or better, which is what confirms
both the equation's printed form and the reduction-magnitude sign convention.

**The 1/2-log threshold is internally inconsistent in the source.** Substituting
Table 2's 33.45 h into Table 2's own equation and parameters gives
E = 0.707 log10 CFU/mL, i.e. very nearly a 3/4-log reduction rather than a
1/2-log one; the Abstract's 23.86 h gives E = 0.466. The value consistent with
the published parameters is 25.01 h. The packaged model encodes the four fitted
parameters (which are self-consistent), so it necessarily disagrees with the
printed 1/2-log threshold. This is recorded under Errata below; no parameter was
adjusted.

# Validation: encoding checks

PKNCA is not applicable to this model: there is no PK component and no
concentration observable to run non-compartmental analysis on. The paper's own
PK analysis was non-compartmental and its outputs (`AUC0-24h`, `Cmax`,
`T>MIC`) are *inputs* to the packaged model, reproduced verbatim in the
exposures table above rather than re-derived. In place of NCA, the checks below
verify that the ODE encoding faithfully implements the published algebraic
model.

```{r encoding-checks}
# Check 1 -- interval-boundary exactness. The ODE must drop log10(bact) by
# exactly the algebraic E over each 24 h interval.
check1 <- do.call(rbind, lapply(names(arm_auc), function(d) {
  s   <- simulate_arm(arm_auc[[d]], d)
  b   <- s$Cc[s$time %in% c(0, 24, 48, 72, 96)]
  data.frame(dose = as.integer(d), interval = 1:4,
             ode_drop = -diff(b),
             algebraic = effect_of_index(arm_auc[[d]] / MIC))
})) |>
  mutate(abs_err = abs(ode_drop - algebraic))

# Check 2 -- untreated control. AUC_TILM = 0 must give a per-interval change
# of exactly E0 and leave the count essentially at the 10^7 CFU/mL baseline.
ctl <- simulate_arm(rep(0, 4), "control")
check2 <- c(per_interval_change = -diff(ctl$Cc[ctl$time %in% c(0, 24)]),
            e0_parameter = e0,
            log10_cfu_at_96h = ctl$Cc[ctl$time == 96])

# Check 3 -- saturation. As AUC_TILM grows without bound the per-interval
# reduction must approach Emax.
check3 <- c(reduction_at_huge_exposure =
              -diff(simulate_arm(rep(1e6, 4), "sat")$Cc[c(1, 49)]),
            emax_parameter = emax)

cat(sprintf("Check 1  max |ODE drop - algebraic E| over 16 intervals: %.3e\n",
            max(check1$abs_err)))
cat(sprintf("Check 2  control per-interval change %.6f vs E0 %.6f; log10 CFU at 96 h = %.4f\n",
            check2[["per_interval_change"]], check2[["e0_parameter"]],
            check2[["log10_cfu_at_96h"]]))
cat(sprintf("Check 3  reduction at saturating exposure %.6f vs Emax %.6f\n",
            check3[["reduction_at_huge_exposure"]], check3[["emax_parameter"]]))
```

```{r encoding-checks-table}
check1 |>
  mutate(across(c(ode_drop, algebraic), ~ round(.x, 6)),
         abs_err = signif(abs_err, 3)) |>
  rename("Dose (mg/kg)" = dose, "24 h interval" = interval,
         "ODE log10 drop" = ode_drop, "Algebraic E" = algebraic,
         "|difference|" = abs_err) |>
  knitr::kable(caption = "Check 1: the integrated ODE reproduces the algebraic per-interval reduction at every 24 h boundary, to solver tolerance.")
```

All three checks pass: the integrated trajectory matches the algebraic
per-interval model to solver tolerance, the untreated control reproduces `E0`
and stays at the 10^7 CFU/mL baseline, and saturating exposure drives the
per-interval reduction to `Emax`.

**Dimensional analysis.** `AUC_TILM` is in ug*h/mL and `mic` in ug/mL, so
`aucmic = AUC_TILM / mic` has units of h -- matching the `EC50` unit printed in
Table 2 ("EC50 (h)"), which is itself a useful confirmation that the index is
`AUC24h/MIC` and not a concentration ratio. `kill_log10` is in log10 CFU/mL per
24 h interval, and dividing by 24 h before multiplying by `ln(10) * bact` gives
`d/dt(bact)` in CFU/mL/h.

# Assumptions and deviations

* **No PK model.** Chen 2023 published no compartmental PK model for tilmicosin
  in tissue-cage fluid (Section 2.5 states PK parameters were computed with
  non-compartmental methods in WinNonlin 6.1). Drug exposure therefore enters
  the packaged model only through the `AUC_TILM` covariate. Simulating a
  regimen not in Table 1 requires an external source of per-interval AUC, for
  example a published tilmicosin popPK or PBPK model in swine; none was
  substituted here.
* **Within-interval uniformity of the kill rate.** The published model is
  defined per 24 h interval, not as an instantaneous rate. The packaged ODE
  spreads each interval's reduction uniformly across the interval. This is an
  encoding choice, not a published result; it is exact at every 24 h boundary
  (the only times Chen 2023 counted bacteria -- Section 2.4), so no observed
  data point is affected, but the shape *within* an interval should not be
  interpreted as a paper finding. Check 1 above verifies the boundary
  exactness.
* **Only the AUC24h/MIC sub-model is packaged.** Chen 2023 fitted the sigmoid
  to all three indices but reports parameter values (Table 2) for
  `AUC24h/MIC` only; for `Cmax/MIC` and `%T>MIC` just the R^2 values (0.90 and
  0.83) are given. The other two sub-models are therefore not packaged, because
  their `E0`, `Emax`, `EC50` and `N` are unpublished. No values were invented.
* **No between-subject variability.** Chen 2023 fitted the PK/PD relationship by
  sigmoid regression across pooled tissue cages and reported no random-effect
  variances, so the packaged model has no `eta` parameters and is a
  typical-value model. The Table 1 standard deviations describe variability in
  the *exposure inputs*, not IIV on the PD parameters, and are not encoded as
  IIV.
* **Residual SD fixed at zero.** Only the coefficient of determination
  (R^2 = 0.92) is reported for the fit; no residual standard deviation is
  given. `addSd <- fixed(0)` makes the model deterministic for typical-value
  simulation rather than encoding an invented error magnitude.
* **Baseline count taken as exactly 10^7 CFU/mL.** Section 2.2 selected cages
  with "bacterial counts around 10^7 CFU/mL" and Section 3.3 reports the
  untreated control "remained around 10^7 CFU/mL". The packaged
  `log10_cfu0 <- fixed(7)` adopts that nominal value; the paper gives no more
  precise mean baseline.
* **MIC treated as a fixed model parameter.** The MIC (0.25 ug/mL) is a
  measured property of strain C44-15 rather than a fitted parameter, and is
  wrapped in `fixed()`. It is exposed as a parameter rather than hard-coded so
  the model can be applied to isolates of different susceptibility, which is
  the point of a PK/PD-index model.
* **Cumulative-reduction window inferred as 96 h.** Section 3.3 reports total
  reductions without restating the window. The Discussion resolves it ("the
  number of bacteria could be reduced by 3 log10CFU/mL within 96 h after
  administration"), and the four-interval sum reproduces the 30 and 60 mg/kg
  totals to within 0.01 log10, which corroborates the 96 h reading.

# Errata and internal inconsistencies in the source

No published erratum or corrigendum for this article was located. The following
are internal inconsistencies within the article as published; none was resolved
by adjusting a parameter.

* **The 1/2-log threshold is inconsistent with the paper's own parameters.**
  Table 2 gives 33.45 h and the Abstract gives 23.86 h for the `AUC24h/MIC`
  needed for a 1/2-log10 CFU/mL reduction. Substituted back into the Section
  2.6 equation with the Table 2 parameters these give E = 0.707 and E = 0.466
  respectively; the self-consistent value is **25.01 h**. Table 2's 33.45 h is
  the more serious of the two, corresponding to nearly a 3/4-log reduction and
  so almost duplicating the adjacent 3/4-log row (35.64 h). The 1/3-log and
  3/4-log thresholds in both locations are consistent with the parameters to
  within 0.3%.
* **Abstract and Results disagree on all three thresholds.** The Abstract
  reports 19.65 h, 23.86 h and 35.77 h; Section 3.4 and Table 2 report 19.64 h,
  33.45 h and 35.64 h. For the 1/3-log and 3/4-log values the difference is
  negligible (and the Abstract's 35.77 h is marginally closer to the
  model-implied 35.72 h than Table 2's 35.64 h). Only the 1/2-log value differs
  materially, as above.
* **The 40 mg/kg cumulative reduction is not reproducible from the published
  exposures and parameters** (2.40 predicted versus 2.82 reported; -14.8%).
  The other three dose groups reproduce to within 3.2%. See the cumulative
  reduction section above.
* **Section 2.6 defines `E` via `(Emax - E0)`, so `Emax` is the effect at
  saturation, not the increment over baseline.** This is stated only implicitly
  in the paper's parameter descriptions ("Emax is the maximum antibacterial
  effect during the dosing interval"). The reading is confirmed by the
  reproduction of the 1/3-log and 3/4-log thresholds above; the alternative
  reading (`E = E0 - Emax * ...`, treating `Emax` as an increment) does not
  reproduce them.
