---
title: "Amphenmulin PK/PD against Mycoplasma gallisepticum (Wang 2024)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Amphenmulin PK/PD against Mycoplasma gallisepticum (Wang 2024)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

Wang 2024 is a preclinical pharmacokinetic/pharmacodynamic integration study of
**amphenmulin**, a novel pleuromutilin derivative designed and synthesised in the
authors' own laboratory, against *Mycoplasma gallisepticum* strain S6. It is not a
marketed drug. The paper has three separable quantitative components, and this
package carries one model file for each.

* Article: <https://doi.org/10.1128/spectrum.03675-23>

```{r models}
mods <- c("Wang_2024_amphenmulin_chicken",
          "Wang_2024_amphenmulin_killrate",
          "Wang_2024_amphenmulin_pkpd_index")
for (m in mods) {
  cat("##", m, "\n")
  cat(rxode2::rxode2(readModelDb(m))$description, "\n\n")
}
```

| Model file | Paper component | Source |
|---|---|---|
| `Wang_2024_amphenmulin_chicken` | One-compartment intravenous PK in broiler chickens | Table 1, IV column |
| `Wang_2024_amphenmulin_killrate` | Sigmoid Emax kill rate driven by the in-vitro dynamic system PK | Table 2 (row 0-24 h) + Materials and Methods, In vitro dynamic model |
| `Wang_2024_amphenmulin_pkpd_index` | Inhibitory sigmoid Emax on the AUC24h/MIC index | Table 3, AUC24h/MIC column |

All three share the citation below.

```{r citation}
cat(rxode2::rxode2(readModelDb("Wang_2024_amphenmulin_chicken"))$reference)
```

# Population

The in-vivo component used **30 healthy adult yellow-feathered broiler chickens**
weighing 1.58 to 2.04 kg, randomised to three groups of 10 that each received a
single 20 mg/kg body-weight dose of amphenmulin by intravenous injection,
intramuscular injection, or oral gavage (Materials and Methods, In vivo
pharmacokinetics). Birds were caged, fed ad libitum, and received no antibiotics
or anticoccidial drugs. Plasma was assayed by HPLC-MS/MS with a lower limit of
quantification of 0.005 ug/mL. Only the intravenous group informs the packaged PK
model; see Assumptions and deviations for why.

The in-vitro components used *M. gallisepticum* standard strain S6 from the
National Center for Veterinary Culture Collection, Beijing, at a starting inoculum
of 10^7 CFU/mL. Amphenmulin MIC was 0.0039 ug/mL by broth microdilution and
0.0078 ug/mL by agar dilution, with MIC99 0.0077 ug/mL, MBC 0.0156 ug/mL and MPC
0.0500 ug/mL. Every in-vitro experiment was run in triplicate.

The same information is available programmatically from each model's `population`
metadata (`readModelDb("Wang_2024_amphenmulin_chicken")()$population`).

```{r constants}
MIC <- 0.0039  # ug/mL, Wang 2024 Results, Susceptibility assay (broth microdilution)
```

# Source trace

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

| Model | Equation / parameter | Value | Source location |
|---|---|---|---|
| chicken | `lcl` (CL) | 1.17 L/h/kg | Table 1, IV column |
| chicken | `lvc` (Vss) | 3.64 L/kg | Table 1, IV column |
| chicken | `addSd` | 0 | Not reported; held at 0 |
| chicken | `d/dt(central) <- -kel * central` | n/a | Discussion ("building a one-compartment elimination model is more straightforward") |
| killrate | `e0` (E0) | -0.0093 1/h | Table 2, row 0-24 h |
| killrate | `lemax` (Emax) | 0.1261 1/h | Table 2, row 0-24 h |
| killrate | `lec50` (EC50) | 0.0325 ug/mL | Table 2, row 0-24 h |
| killrate | `lhill` (N) | 0.9238 | Table 2, row 0-24 h |
| killrate | `lkel` | 0.3254 1/h | Materials and Methods, In vitro dynamic model; `log(2)/T1/2Ke` with T1/2Ke = 2.13 h from Table 1 |
| killrate | `lvc` | 300 mL | Materials and Methods, In vitro dynamic model (external compartment of the reaction chamber) |
| killrate | `log10_cfu0` | 7 | Materials and Methods, Static time-kill curves (10^7 CFU/mL) |
| killrate | `E = E0 + Emax * Ce^N / (EC50^N + Ce^N)` | n/a | Materials and Methods, Static time-kill curves (printed equation) |
| killrate | `C = C0 * exp(-k*t)` | n/a | Materials and Methods, In vitro dynamic model (printed equation) |
| pkpd_index | `e0` (E0) | -0.3845 log10 CFU/mL | Table 3, AUC24h/MIC column |
| pkpd_index | `emax` (Emax) | -2.4214 log10 CFU/mL | Table 3, AUC24h/MIC column |
| pkpd_index | `lec50` (EC50) | 1199.4720 h | Table 3, AUC24h/MIC column |
| pkpd_index | `lhill` (N) | 3.1997 | Table 3, AUC24h/MIC column |
| pkpd_index | `mic` | 0.0039 ug/mL | Results, Susceptibility assay |
| pkpd_index | `log10_cfu0` | 7 | Materials and Methods, In vitro dynamic model |
| pkpd_index | `E = E0 - (E0 - Emax) * Ce^N / (EC50^N + Ce^N)` | n/a | Materials and Methods, Integration and modeling of PK/PD (printed equation) |

Both printed equations were read from the article PDF. The two sigmoid forms are
genuinely different and are transcribed as printed: in the static time-kill
equation `Emax` is an **increment** above `E0`, whereas in the inhibitory equation
`Emax` is the **asymptote**. This matches the two distinct Phoenix built-in models
the authors named ("Sigmoid Emax model in Phoenix" and "inhibitory Sigmoid Emax
model"). See Assumptions and deviations.

# Part 1 - Intravenous pharmacokinetics in broiler chickens

The model is deterministic (Wang 2024 fitted no population model), so a single
typical-value profile is simulated on the paper's own intravenous sampling grid.

```{r chicken-sim}
mod_pk <- readModelDb("Wang_2024_amphenmulin_chicken")

# Wang 2024 Materials and Methods: IV sampling at 0, 5, 10, 15, 30, 45 min and
# 1, 2, 3, 4, 6, 8, 10, 12, 24 h.
t_iv <- c(0, 5, 10, 15, 30, 45) / 60
t_iv <- c(t_iv, 1, 2, 3, 4, 6, 8, 10, 12, 24)

ev_iv <- data.frame(
  id   = 1L,
  time = c(0, t_iv),
  evid = c(1L, rep(0L, length(t_iv))),
  amt  = c(20, rep(NA_real_, length(t_iv))),   # 20 mg/kg IV bolus
  cmt  = "central"
)

sim_iv <- rxode2::rxSolve(mod_pk, ev_iv, useLinCmt = FALSE,
                          returnType = "data.frame")
if (!"id" %in% names(sim_iv)) sim_iv$id <- 1L
sim_iv$treatment <- "IV 20 mg/kg"
head(sim_iv[, c("id", "time", "Cc", "treatment")])
```

```{r figure-3}
# Replicates the intravenous curve of Figure 3 of Wang 2024 (log concentration
# scale). Only the IV arm is packaged; see Assumptions and deviations.
sim_iv |>
  dplyr::filter(!is.na(Cc), Cc > 0) |>
  ggplot(aes(time, Cc)) +
  geom_line() +
  geom_point() +
  scale_y_log10() +
  labs(x = "Time (h)", y = "Amphenmulin plasma concentration (ug/mL)",
       title = "Figure 3 (IV arm) - amphenmulin 20 mg/kg IV in broiler chickens",
       caption = "Replicates the intravenous curve of Figure 3 of Wang 2024.")
```

## PKNCA validation against Table 1

```{r pknca-chicken}
sim_nca <- sim_iv |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

# IV bolus: the time-zero row is the post-dose concentration C0 = Dose/Vss, which
# the simulation already produces, so no defensive Cc = 0 row is added here (that
# would be correct for an extravascular model, not for an IV bolus).
stopifnot(any(sim_nca$time == 0), sim_nca$Cc[sim_nca$time == 0] > 0)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id,
                             concu = "ug/mL", timeu = "h")

dose_df <- ev_iv |>
  dplyr::filter(evid == 1) |>
  dplyr::select(id, time, amt) |>
  dplyr::mutate(treatment = "IV 20 mg/kg")

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id,
                             doseu = "mg/kg", route = "intravascular",
                             duration = 0)

intervals_iv <- data.frame(
  start = 0, end = Inf,
  auclast = TRUE, aucinf.obs = TRUE, half.life = TRUE,
  cl.obs = TRUE, mrt.obs = TRUE, vss.obs = TRUE
)

nca_iv <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj, dose_obj,
                                         intervals = intervals_iv))
```

```{r nca-compare-chicken}
# Wang 2024 Table 1, IV column (mean +/- SD over 10 birds). Cmax and Tmax are
# reported as "--" for the IV arm and are therefore not compared.
published_iv <- tibble::tribble(
  ~treatment,    ~auclast, ~aucinf.obs, ~half.life, ~cl.obs, ~mrt.obs, ~vss.obs,
  "IV 20 mg/kg", 18.76,    18.89,       2.13,       1.17,    3.13,     3.64
)

cmp_iv <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_iv,
  reference = published_iv,
  by        = "treatment",
  units     = c(auclast = "h*ug/mL", aucinf.obs = "h*ug/mL", half.life = "h",
                cl.obs = "L/h/kg", mrt.obs = "h", vss.obs = "L/kg"),
  tolerance_pct = 20
)

knitr::kable(
  cmp_iv,
  caption = "Simulated vs. published intravenous NCA (Wang 2024 Table 1). * differs from reference by >20%.",
  align = c("l", "l", "r", "r", "r")
)
```

No row exceeds the 20% tolerance. CL and Vss are recovered exactly because they
are the model's primary parameters; the half-life agrees to 1.2% and the MRT to
0.6%, which is the substantive check that Vss is the right single-compartment
volume. The AUC values sit about 9% low. That gap is arithmetic, not an encoding
error: Table 1 reports the mean of the per-bird CL and the mean of the per-bird
AUC separately, and `mean(Dose/AUC_i)` is not `Dose/mean(AUC_i)`. Simulating from
the reported mean CL of 1.17 L/h/kg necessarily returns `20/1.17 = 17.1`
h\*ug/mL rather than the separately averaged 18.89. The gap is well inside the
reported SD of 6.03.

Note that `ncaComparisonTable()` labels the clearance and volume rows "CL/F" and
"Vss/F"; for the intravenous route F is 1 by definition, so these are CL and Vss.

# Part 2 - Static time-kill curves and the kill-rate model

Wang 2024 Table 2 reports the sigmoid Emax kill-rate fit over seven time windows.
The packaged model carries the 0-24 h row, which the paper identifies as the
optimal fit (R = 0.9936) and plots in Figure 4b. All seven are tabulated here.

```{r table-2}
table2 <- tibble::tribble(
  ~window,  ~Emax,   ~EC50,   ~E0,      ~N,      ~R,
  "0-24",   0.1261,  0.0325,  -0.0093,  0.9238,  0.9936,
  "0-36",   0.1817,  0.0914,  -0.0359,  0.4264,  0.9856,
  "0-48",   0.3277,  0.9135,  -0.0361,  0.3760,  0.9914,
  "6-24",   0.0805,  0.0082,  -0.0373,  2.7101,  0.9748,
  "6-36",   0.0821,  0.0053,  -0.0550,  2.9068,  0.9649,
  "6-48",   0.1743,  0.0526,  -0.0524,  0.4175,  0.9820,
  "12-48",  0.2880,  0.5245,  -0.0615,  0.2919,  0.9839
)

table2 |>
  dplyr::mutate(`E0 + Emax` = E0 + Emax) |>
  dplyr::rename("Time window (h)" = window, "Emax (1/h)" = Emax,
                "EC50 (ug/mL)" = EC50, "E0 (1/h)" = E0, "Hill N" = N,
                "Saturating kill rate (1/h)" = `E0 + Emax`) |>
  knitr::kable(digits = 4,
               caption = "Wang 2024 Table 2. The packaged model uses the 0-24 h row. The final column is E0 + Emax, the kill rate the printed equation approaches at saturating concentration.")
```

## The ODE encoding reproduces the fitted kill rate exactly

The model integrates `d/dt(bact) = -ln(10) * E * bact`, so `log10(bact)` must fall
at exactly `E` log10 CFU/mL per hour at any constant concentration. Setting the
elimination rate constant to zero recovers the static experiment.

```{r killrate-static}
mod_kr <- rxode2::rxode2(readModelDb("Wang_2024_amphenmulin_killrate"))
mod_static <- mod_kr |> rxode2::ini(lkel = log(1e-10))

kr_events <- function(conc, times) {
  data.frame(
    id   = 1L,
    time = c(0, times),
    evid = c(1L, rep(0L, length(times))),
    amt  = c(conc * 300, rep(NA_real_, length(times))),  # vc = 300 mL
    cmt  = "central",
    dvid = c(NA, rep(1L, length(times)))
  )
}

# Closed form of the Table 2 row 0-24 h equation.
kill_closed <- function(ce, E0 = -0.0093, Emax = 0.1261,
                        EC50 = 0.0325, N = 0.9238) {
  E0 + Emax * ce^N / (EC50^N + ce^N)
}

mult <- c(0, 1, 2, 4, 6, 8, 16, 32, 64)
static_check <- lapply(mult, function(k) {
  s <- rxode2::rxSolve(mod_static, kr_events(k * MIC, c(0, 24)),
                       useLinCmt = FALSE, returnType = "data.frame")
  tibble::tibble(
    `MIC multiple` = k,
    `Concentration (ug/mL)` = k * MIC,
    `ODE kill rate (1/h)` = -(s$logCfu[2] - s$logCfu[1]) / 24,
    `Table 2 closed form (1/h)` = kill_closed(k * MIC)
  )
}) |> dplyr::bind_rows() |>
  dplyr::mutate(`Absolute difference` = abs(`ODE kill rate (1/h)` -
                                              `Table 2 closed form (1/h)`))

knitr::kable(static_check, digits = c(0, 4, 6, 6, 10),
             caption = "The ODE encoding reproduces the Table 2 (0-24 h) sigmoid Emax kill rate to solver tolerance at every tested concentration, including the drug-free control.")
```

The agreement is to about 1e-7, i.e. solver tolerance, which confirms the
`-ln(10)` conversion between the log10-scale kill rate the paper reports and the
linear-density ODE state.

```{r figure-4b}
# Replicates Figure 4b of Wang 2024: fitted kill rate vs amphenmulin
# concentration over the 0-24 h window.
curve_df <- tibble::tibble(conc = 10^seq(log10(MIC / 2), log10(64 * MIC),
                                         length.out = 200)) |>
  dplyr::mutate(kill = kill_closed(conc))

ggplot(curve_df, aes(conc, kill)) +
  geom_line() +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_vline(xintercept = 0.0325, linetype = "dashed", colour = "grey40") +
  annotate("text", x = 0.0325, y = -0.008, hjust = -0.1, size = 3,
           label = "EC50 = 0.0325 ug/mL") +
  scale_x_log10() +
  labs(x = "Amphenmulin concentration (ug/mL, log scale)",
       y = "Kill rate (log10 CFU/mL per h)",
       title = "Figure 4b - sigmoid Emax kill rate, 0-24 h window",
       caption = "Replicates Figure 4b of Wang 2024. Negative kill rate = net mycoplasma growth.")
```

## Check against the observed 48 h reductions

The Results section reports the observed 48 h reductions from Figure 4a. Those are
compared below against the **0-48 h** parameter row, which is the row fitted to
that window.

```{r kill-48h}
kill_48 <- function(ce) kill_closed(ce, E0 = -0.0361, Emax = 0.3277,
                                    EC50 = 0.9135, N = 0.3760)

tibble::tibble(
  `MIC multiple` = c(1, 64),
  `Concentration (ug/mL)` = c(1, 64) * MIC,
  `Observed 48 h reduction (log10 CFU/mL)` = c(0.44, 3.85),
  `Model 48 h reduction (log10 CFU/mL)` = kill_48(c(1, 64) * MIC) * 48
) |>
  knitr::kable(digits = 4,
               caption = "Wang 2024 Results (Analysis of static time-kill curves) reports 0.44 log10 CFU/mL at 1 MIC and 3.85 log10 CFU/mL at 64 MIC over 48 h; between 2 and 32 MIC the reduction was 0.57 to 2.70 log10 CFU/mL. The model row fitted to the 0-48 h window is shown for comparison.")
```

The 64 MIC endpoint is reproduced within about 10%. The 1 MIC endpoint is not: the
0-48 h fit has a very shallow Hill slope (N = 0.376) and an EC50 of 0.9135 ug/mL,
far above the tested range, so it under-predicts the low-concentration end. This
is a property of the published fit, not of the encoding - the closed form above is
Table 2 row 0-48 h evaluated directly. It is one reason the paper identifies the
0-24 h window as optimal.

## The in-vitro dynamic system

Coupling the Table 2 kill rate to the apparatus PK reproduces Figure 5.

```{r figure-5}
regimens <- c(0.1, 0.2, 0.5, 0.8, 1.0, 1.5)  # target peak concentration, ug/mL
obs_t <- sort(unique(c(seq(0, 72, by = 0.25), 0, 24, 48, 72)))

dyn_events <- function(conc, id) {
  dose <- data.frame(id = id, time = c(0, 24, 48), evid = 1L,
                     amt = conc * 300, cmt = "central", dvid = NA_integer_)
  obs <- data.frame(id = id, time = obs_t, evid = 0L,
                    amt = NA_real_, cmt = "central", dvid = 1L)
  dplyr::bind_rows(dose, obs) |> dplyr::arrange(time, dplyr::desc(evid))
}

ev_dyn <- dplyr::bind_rows(
  lapply(seq_along(regimens), function(i) {
    e <- dyn_events(regimens[i], id = i)
    e$regimen <- sprintf("%.1f ug/mL", regimens[i])
    e
  })
)
stopifnot(!anyDuplicated(unique(ev_dyn[, c("id", "time", "evid")])))

sim_dyn <- rxode2::rxSolve(mod_kr, ev_dyn, useLinCmt = FALSE,
                           keep = "regimen", returnType = "data.frame") |>
  dplyr::filter(!is.na(logCfu))
```

```{r figure-5a}
# Replicates Figure 5a of Wang 2024: amphenmulin concentration-time curves in the
# dynamic model at each dosing regimen.
sim_dyn |>
  ggplot(aes(time, Cc, colour = regimen)) +
  geom_line() +
  scale_y_log10() +
  scale_x_continuous(breaks = seq(0, 72, by = 12)) +
  labs(x = "Time (h)", y = "Amphenmulin concentration (ug/mL, log scale)",
       colour = "Peak concentration",
       title = "Figure 5a - concentration-time curves in the dynamic model",
       caption = "Replicates Figure 5a of Wang 2024.")
```

```{r figure-5b}
# Replicates Figure 5b of Wang 2024: dynamic time-kill curves.
sim_dyn |>
  ggplot(aes(time, logCfu, colour = regimen)) +
  geom_line() +
  geom_hline(yintercept = 7 - 3, linetype = "dashed", colour = "grey40") +
  annotate("text", x = 2, y = 7 - 3 + 0.15, hjust = 0, size = 3,
           label = "bactericidal threshold (3 log10 reduction)") +
  scale_x_continuous(breaks = seq(0, 72, by = 12)) +
  labs(x = "Time (h)", y = "Mycoplasma density (log10 CFU/mL)",
       colour = "Peak concentration",
       title = "Figure 5b - dynamic time-kill curves",
       caption = "Replicates Figure 5b of Wang 2024.")
```

### An independent check of the composition

The kill-rate parameters were fitted to **static** concentrations and the dynamic
apparatus was a **separate** experiment; Wang 2024 never fitted the two jointly.
Composing them is therefore a falsifiable prediction, and the paper's own
narrative conclusions about the dynamic experiment are the test.

```{r dynamic-check}
net_change <- sim_dyn |>
  dplyr::group_by(regimen) |>
  dplyr::summarise(`Net change over 0-72 h (log10 CFU/mL)` =
                     dplyr::last(logCfu) - dplyr::first(logCfu),
                   .groups = "drop") |>
  dplyr::mutate(`Bactericidal (>= 3 log10 reduction)` =
                  ifelse(`Net change over 0-72 h (log10 CFU/mL)` <= -3, "yes", "no"))

knitr::kable(net_change, digits = 2,
             caption = "Predicted net change in mycoplasma density over the three-day course, from the static kill-rate parameters composed with the apparatus PK.")
```

Wang 2024 Results (In vitro pharmacokinetics and the effects on *M.
gallisepticum*) states two things about this experiment: that at 0.1 to 0.5 ug/mL
the counts fell by **at most 2.82 log10 CFU/mL** over 0-72 h, and that **at 0.8
ug/mL or above amphenmulin was bactericidal** (a reduction of at least 3 log10
CFU/mL). The composed model reproduces both: the largest reduction across the 0.1
to 0.5 ug/mL regimens is about 2.6 log10 CFU/mL, and the bactericidal threshold is
crossed at exactly 0.8 ug/mL. Neither statement was used to build the model, so
this is genuine external corroboration of the composition rather than a
self-consistency check.

## PKNCA on the dynamic-system exposures

```{r pknca-dyn}
sim_dyn_nca <- sim_dyn |>
  dplyr::filter(!is.na(Cc), regimen %in% c("0.1 ug/mL", "1.5 ug/mL")) |>
  dplyr::select(id, time, Cc, regimen)

conc_dyn <- PKNCA::PKNCAconc(sim_dyn_nca, Cc ~ time | regimen + id,
                             concu = "ug/mL", timeu = "h")

dose_dyn <- ev_dyn |>
  dplyr::filter(evid == 1, regimen %in% c("0.1 ug/mL", "1.5 ug/mL")) |>
  dplyr::select(id, time, amt, regimen)

dose_dyn_obj <- PKNCA::PKNCAdose(dose_dyn, amt ~ time | regimen + id,
                                 doseu = "ug", route = "intravascular",
                                 duration = 0)

nca_dyn <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  conc_dyn, dose_dyn_obj,
  intervals = data.frame(start = 0, end = 24, cmax = TRUE, auclast = TRUE)
))
```

```{r nca-compare-dyn}
# Wang 2024 Results (In vitro pharmacokinetics and the effects on M.
# gallisepticum) reports the 0-24 h AUC across the six regimens as 0.34 to 6.19
# h*ug/mL and the measured first-interval Cmax as 0.12 to 1.52 ug/mL. The range
# endpoints correspond to the lowest and highest regimen.
published_dyn <- tibble::tribble(
  ~regimen,     ~cmax, ~auclast,
  "0.1 ug/mL",  0.12,  0.34,
  "1.5 ug/mL",  1.52,  6.19
)

cmp_dyn <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_dyn,
  reference = published_dyn,
  by        = "regimen",
  units     = c(cmax = "ug/mL", auclast = "h*ug/mL"),
  tolerance_pct = 20
)

knitr::kable(
  cmp_dyn,
  caption = "Simulated vs. published in-vitro dynamic-system exposures over 0-24 h. * differs from reference by >20%.",
  align = c("l", "l", "r", "r", "r")
)
```

Both Cmax values and the AUC at the lowest regimen agree within 20%. The AUC at
the highest regimen is starred. The cause is visible in the published numbers
themselves: a first-order system with a single elimination rate constant must give
the same AUC/Cmax ratio at every dose level, but Wang 2024's reported values give
`0.34/0.12 = 2.8` h at the lowest regimen and `6.19/1.52 = 4.1` h at the highest,
a 44% spread that no linear model can reproduce. The packaged model uses the
apparatus design constant the paper specifies (`k = log(2)/2.13 h`, giving
`1/k = 3.07` h, between the two), and the parameters were **not** tuned toward
either endpoint.

# Part 3 - The PK/PD index model

Wang 2024 Table 3 fits the inhibitory sigmoid Emax model against two PK/PD
indices. AUC24h/MIC correlates better (R = 0.9657 vs 0.8995) and is the
parameterisation packaged here.

```{r table-3}
table3 <- tibble::tribble(
  ~index,        ~Emax,    ~EC50,      ~E0,      ~N,      ~R,      ~target,
  "AUC24h/MIC",  -2.4214,  1199.4720,  -0.3845,  3.1997,  0.9657,  904.05,
  "Cmax/MIC",    -3.3571,  441.4603,   -0.3575,  1.3842,  0.8995,  190.11
)

table3 |>
  dplyr::rename("PK/PD index" = index, "Emax (log10 CFU/mL)" = Emax,
                "EC50" = EC50, "E0 (log10 CFU/mL)" = E0, "Hill N" = N,
                "Published target for 1 log10 reduction" = target) |>
  knitr::kable(digits = 4, caption = "Wang 2024 Table 3 and the published PK/PD targets.")
```

## Back-solving the published targets

The paper states that a 1 log10 CFU/mL reduction requires AUC24h/MIC = 904.05 h or
Cmax/MIC = 190.11. Substituting Table 3 into the printed equation is an
independent confirmation of both the equation form and the sign convention.

```{r target-backsolve}
inhib_emax <- function(ce, E0, Emax, EC50, N) {
  E0 - (E0 - Emax) * ce^N / (EC50^N + ce^N)
}

backsolve <- table3 |>
  dplyr::rowwise() |>
  dplyr::mutate(
    `E at the published target (log10 CFU/mL)` =
      inhib_emax(target, E0, Emax, EC50, N),
    `Index giving exactly -1 log10` =
      stats::uniroot(function(x) inhib_emax(x, E0, Emax, EC50, N) + 1,
                     c(1e-3, 1e6))$root
  ) |>
  dplyr::ungroup() |>
  dplyr::mutate(`% difference from published target` =
                  100 * (`Index giving exactly -1 log10` - target) / target) |>
  dplyr::select(index, target, `E at the published target (log10 CFU/mL)`,
                `Index giving exactly -1 log10`,
                `% difference from published target`) |>
  dplyr::rename("PK/PD index" = index, "Published target" = target)

knitr::kable(backsolve, digits = c(0, 2, 4, 2, 1),
             caption = "The printed inhibitory sigmoid Emax equation evaluated at Wang 2024's published PK/PD targets returns approximately -1 log10 CFU/mL for both indices, confirming the equation form and the signed-effect convention.")
```

Both indices return close to -1 log10 CFU/mL, with the residual discrepancy
attributable to rounding of the four-significant-figure parameters. The AUC24h/MIC
target reproduces to 2.1%; the Cmax/MIC target to 9.2%. The larger gap on
Cmax/MIC is expected given its shallower Hill slope (N = 1.38 versus 3.20) and
correspondingly higher sensitivity to rounding in EC50.

Critically, evaluating the **alternative** sign convention (treating `Emax` as an
increment rather than the asymptote, as in the static time-kill equation) returns
-1.08 and -1.16 respectively, i.e. worse for both indices. The two equations are
therefore correctly transcribed as different forms.

```{r figure-6}
# Replicates Figure 6 of Wang 2024: anti-mycoplasma effect vs PK/PD index.
fig6 <- dplyr::bind_rows(lapply(seq_len(nrow(table3)), function(i) {
  r <- table3[i, ]
  ce <- 10^seq(0, log10(5 * r$target), length.out = 150)
  tibble::tibble(index = r$index, ce = ce, target = r$target,
                 E = inhib_emax(ce, r$E0, r$Emax, r$EC50, r$N))
}))

ggplot(fig6, aes(ce, E)) +
  geom_line() +
  geom_hline(yintercept = -1, linetype = "dotted") +
  geom_vline(aes(xintercept = target), linetype = "dashed", colour = "grey40") +
  facet_wrap(~index, scales = "free_x") +
  scale_x_log10() +
  labs(x = "PK/PD index (log scale)",
       y = "Change in mycoplasma count over 24 h (log10 CFU/mL)",
       title = "Figure 6 - inhibitory sigmoid Emax fits",
       caption = paste("Replicates Figure 6 of Wang 2024. Dashed line: published",
                       "target for a 1 log10 reduction (dotted horizontal line)."))
```

## The packaged ODE reproduces the per-interval change

```{r pkpd-ode}
mod_idx <- readModelDb("Wang_2024_amphenmulin_pkpd_index")

idx_check <- lapply(c(0, 300, 904.05, 2000), function(aucmic) {
  ev <- data.frame(id = 1L, time = c(0, 24), evid = 0L,
                   AUC_AMPH = aucmic * MIC)
  s <- rxode2::rxSolve(mod_idx, ev, useLinCmt = FALSE,
                       returnType = "data.frame")
  tibble::tibble(
    `AUC24h/MIC (h)` = aucmic,
    `AUC_AMPH covariate (h*ug/mL)` = aucmic * MIC,
    `ODE change over 24 h (log10 CFU/mL)` = s$Cc[2] - s$Cc[1],
    `Table 3 closed form (log10 CFU/mL)` =
      inhib_emax(aucmic, -0.3845, -2.4214, 1199.4720, 3.1997)
  )
}) |> dplyr::bind_rows()

knitr::kable(idx_check, digits = c(2, 4, 6, 6),
             caption = "The packaged ODE reproduces the Table 3 per-interval change exactly at each 24 h boundary, including the untreated control (AUC24h/MIC = 0, which returns E0).")
```

# Assumptions and deviations

* **Only the intravenous route is packaged in the chicken PK model.** Wang 2024
  Table 1 also reports oral and intramuscular NCA, but neither can be represented
  by a one-compartment model built on the intravenous parameters:
  * The published oral Cmax of 0.73 ug/mL is **2.26 times** `F * Dose / Vss =
    0.0588 * 20 / 3.64 = 0.323` ug/mL, which is the absolute ceiling on the
    concentration a one-compartment model can reach after an oral dose at F =
    5.88%, attained only in the limit of instantaneous absorption with no
    elimination. Equivalently, the published oral Cmax/AUC ratio of 0.658 1/h
    exceeds `kel = 0.321` 1/h, which is the maximum a one-compartment model can
    produce. Encoding the oral route with the paper's own `ka = 3.15` 1/h
    (Discussion, estimated by the method of residuals) yields Cmax = 0.25 ug/mL
    and Tmax = 0.81 h against the published 0.73 ug/mL and 0.25 h. The likeliest
    explanation is a distribution phase, so that the true central volume is
    smaller than Vss, but Wang 2024 reports no Vc or Vz from which a
    two-compartment model could be built.
  * The intramuscular route is feasible in magnitude (`F * Dose / Vss = 2.87`
    ug/mL against a published Cmax of 1.93 ug/mL) but Wang 2024 reports no
    intramuscular absorption rate constant, and the reported Tmax of 0.38 h is
    incompatible with the reported terminal rate constant of 0.13 1/h: flip-flop
    absorption at 0.13 1/h would place Tmax at 4.7 h. Inventing a `ka` was not an
    option.
* **Vss is used as the one-compartment volume** in the chicken model. This is
  corroborated by two independent identities: `Vss/CL = 3.11` h against the
  reported MRT of 3.13 h, and `log(2) * Vss/CL = 2.16` h against the reported
  T1/2Ke of 2.13 h.
* **CL and AUC in Table 1 are not exactly reciprocal.** `Dose/CL = 17.1` h\*ug/mL
  against the separately reported mean AUC0-inf of 18.89 h\*ug/mL. Both are means
  of per-bird estimates, and `mean(Dose/AUC_i)` exceeds `Dose/mean(AUC_i)`. CL is
  taken as the model primary; the 9.5% gap is inside the reported SD.
* **Text versus printed equation for the static kill model.** Wang 2024's Results
  text describes Emax as "the maximum kill rate", but the printed equation
  `E = E0 + Emax * Ce^N / (EC50^N + Ce^N)` makes Emax an increment above E0, so
  the saturating kill rate is `E0 + Emax`. The printed equation is used, both
  because it is the more specific statement and because it matches the Phoenix
  built-in Sigmoid Emax model the authors named. For the packaged 0-24 h row the
  difference is 0.1168 versus 0.1261 1/h.
* **The two sigmoid equations use different sign and parameterisation
  conventions**, transcribed as printed. In the static time-kill model E is a
  **kill rate** with positive meaning killing, and E0 is negative because the
  untreated control grows. In the PK/PD-index model E is a **signed change** in
  log10 CFU/mL with negative meaning killing, and Emax is the asymptote rather
  than an increment. The back-solve of the published targets above confirms this
  reading numerically for the second equation.
* **The kill-rate model composes two separately fitted components.** Wang 2024
  fitted the Table 2 kill rates to static concentrations, and ran the dynamic
  apparatus as a separate experiment; the two were never fitted jointly. The
  packaged `Wang_2024_amphenmulin_killrate` model couples them, which is the
  paper's stated PK/PD integration written as an ODE system. Setting the
  elimination rate constant to zero recovers the static experiment exactly, as
  demonstrated above. The composition is corroborated, not assumed: it reproduces
  the paper's independent narrative conclusions about the dynamic experiment
  (a maximum of 2.82 log10 CFU/mL reduction at 0.1 to 0.5 ug/mL, and a
  bactericidal effect from 0.8 ug/mL upward).
* **Only the 0-24 h kill-rate window is packaged.** Wang 2024 fitted seven
  overlapping windows; the paper identifies 0-24 h as optimal and plots it in
  Figure 4b. All seven parameter sets are tabulated above so any of them can be
  substituted with `rxode2::ini()`.
* **Only the AUC24h/MIC index is packaged** in the PK/PD-index model, following
  the paper's own finding that it is the better correlate (R = 0.9657 versus
  0.8995). The Cmax/MIC parameterisation is tabulated and plotted above.
* **The apparatus elimination rate is a design constant, not a fitted value.**
  It is set from the chicken intravenous half-life of 2.13 h. The measured
  in-vitro AUC and Cmax values reported in the Results are mutually inconsistent
  with any single first-order rate constant (see the starred row above); nothing
  was tuned to accommodate them.
* **No between-subject variability and no residual error are encoded.** Wang 2024
  reports mean and SD of individual NCA estimates for the chicken data and only
  the correlation coefficient R for each sigmoid fit; it publishes no population
  model, no OMEGA, and no residual SD. Every residual SD is therefore held at
  zero and the models are intended for deterministic typical-value simulation.
  The Table 1 SDs are reproduced in the `population` metadata so a user can add
  variability deliberately.
* **The kill-rate model has no carrying-capacity term**, so the drug-free control
  grows without bound. Over the 48 to 72 h horizon of the source experiments this
  is immaterial (about 0.22 log10 CFU/mL gained per 24 h), but the model should
  not be extrapolated to long drug-free periods.
* **`useLinCmt = FALSE` is passed to every `rxSolve()` call.** The kill-rate model
  is multi-output (`Cc` and `logCfu`), and rxode2's automatic ODE-to-linCmt
  conversion corrupts the dvid-to-cmt mapping for such models.
* **No non-paper-derived parameter values.** Every value in all three model files
  comes from the article text or its tables; nothing was digitised from a figure,
  supplied by correspondence, or carried from an upstream model.
