---
title: "Guhong injection PK/PD in myocardial ischemia/reperfusion rats (Chen 2024)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Guhong injection PK/PD in myocardial ischemia/reperfusion rats (Chen 2024)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

Guhong injection (GHI) is a Chinese herbal compound preparation of
N-acetyl-L-glutamine plus an aqueous safflower (*Carthamus tinctorius* L.)
extract, used clinically for cardio- and cerebrovascular disease. Chen 2024
gave a single tail-vein injection of GHI at 2.5, 5 or 10 mL/kg to male
Sprague-Dawley rats subjected to myocardial ischemia/reperfusion (MI/R) injury,
quantified eight constituents in plasma by LC-MS/MS, and fitted a
two-compartment intravenous disposition model to each constituent in each dose
group. Four myocardial-injury biomarkers were measured in parallel and linked
to the plasma concentrations through direct-effect sigmoid-Emax models.

The paper contributes **eight** independent model files - one per quantified
constituent, because the authors fitted each constituent separately:

```{r models}
guhong_models <- c(
  "Chen_2024_nAcetylglutamine_rat",
  "Chen_2024_hydroxysafflorYellowA_rat",
  "Chen_2024_chlorogenicAcid_rat",
  "Chen_2024_pCoumaricAcid_rat",
  "Chen_2024_rutin_rat",
  "Chen_2024_hyperoside_rat",
  "Chen_2024_kaempferol3ORutinoside_rat",
  "Chen_2024_kaempferol3OGlucoside_rat"
)
```

* Citation: `r readModelDb("Chen_2024_chlorogenicAcid_rat")()$reference`
* Article: [J Food Drug Anal 2024;32(1):79-102](https://doi.org/10.38212/2224-6614.3492)

The `description` of one member, for orientation:

```{r description}
cat(strwrap(readModelDb("Chen_2024_chlorogenicAcid_rat")()$description, 76), sep = "\n")
```

# Population

Adult male Sprague-Dawley rats weighing 320-340 g were randomised to six groups
(Chen 2024 Sections 2.2-2.3): sham, MI/R model, three GHI dose groups (GHI-L
2.5 mL/kg, GHI-M 5 mL/kg, GHI-H 10 mL/kg) and a verapamil hydrochloride
0.95 mg/kg positive control. MI/R was induced by 30 min of left-anterior-
descending coronary artery ligation followed by 1 h of reperfusion. Saline
(10 mL/kg) was given to the sham and model groups. Plasma was drawn from the
retro-orbital venous plexus at 2, 5, 10, 20, 40, 60, 90, 120, 240 and 360 min.
The pharmacokinetic parameters (Tables 2-9) come from n = 6 rats per GHI dose
group; the biomarker series (Figure 3) from n = 3 per group. Only male rats
were studied, so no sex covariate applies, and no age or race/ethnicity data
are reported.

The same information is available programmatically:

```{r population}
str(readModelDb("Chen_2024_chlorogenicAcid_rat")()$population)
```

# Analytes, content in GHI, and administered dose

Each constituent's dose is the GHI volume dose multiplied by its content in the
preparation (Chen 2024 Results 3.2). This is what the `DOSE_GHI_MLKG` covariate
encodes, and it is what turns one injection into eight simultaneous doses.

```{r analytes}
analytes <- tibble::tribble(
  ~model,                                 ~analyte,                    ~abbrev, ~content, ~content_units, ~dose_units, ~conc_units, ~pk_table,
  "Chen_2024_nAcetylglutamine_rat",       "N-acetyl-L-glutamine",      "NAG",   30.0,     "mg/mL",        "mg/kg",     "ug/mL",     7L,
  "Chen_2024_hydroxysafflorYellowA_rat",  "hydroxysafflor yellow A",   "HSYA",   1.0,     "mg/mL",        "mg/kg",     "ug/mL",     3L,
  "Chen_2024_chlorogenicAcid_rat",        "chlorogenic acid",          "CGA",   12.9,     "ug/mL",        "ug/kg",     "ng/mL",     2L,
  "Chen_2024_pCoumaricAcid_rat",          "p-coumaric acid",           "pCA",   78.3,     "ug/mL",        "ug/kg",     "ng/mL",     8L,
  "Chen_2024_rutin_rat",                  "rutin",                     "RT",    11.7,     "ug/mL",        "ug/kg",     "ng/mL",     9L,
  "Chen_2024_hyperoside_rat",             "hyperoside",                "HYP",    3.4,     "ug/mL",        "ug/kg",     "ng/mL",     4L,
  "Chen_2024_kaempferol3ORutinoside_rat", "kaempferol-3-O-rutinoside", "K-3-R", 65.7,     "ug/mL",        "ug/kg",     "ng/mL",     6L,
  "Chen_2024_kaempferol3OGlucoside_rat",  "kaempferol-3-O-glucoside",  "K-3-G",  8.7,     "ug/mL",        "ug/kg",     "ng/mL",     5L
)

dose_groups <- tibble::tibble(
  group = factor(c("GHI-L", "GHI-M", "GHI-H"), levels = c("GHI-L", "GHI-M", "GHI-H")),
  ghi_dose = c(2.5, 5, 10)
)

dosing <- tidyr::expand_grid(analytes, dose_groups) |>
  dplyr::mutate(amt = content * ghi_dose)

dosing |>
  dplyr::select(analyte, abbrev, content, content_units, group, ghi_dose, amt, dose_units) |>
  dplyr::rename(
    "Analyte" = analyte, "Abbrev." = abbrev,
    "Content in GHI" = content, "Content units" = content_units,
    "Dose group" = group, "GHI dose (mL/kg)" = ghi_dose,
    "Analyte dose" = amt, "Dose units" = dose_units
  ) |>
  knitr::kable(caption = "Analyte dose per GHI dose group (Chen 2024 Results 3.2 and Section 2.3).")
```

# Source trace

Every `ini()` entry carries an in-file comment naming its source cell. The table
below collects the provenance by kind; the per-value comments in
`inst/modeldb/specificDrugs/Chen_2024_*_rat.R` are the authoritative record.

| Equation / parameter | Source location |
|----------------------|-----------------|
| `d/dt(central)`, `d/dt(peripheral1)` (two-compartment IV disposition) | Chen 2024 Section 2.9 ("we first investigated the PK parameters with a two-compartment model") and formulas 3-12 |
| `lvc_ghil` / `lvc_ghim` / `lvc_ghih` (V1, L/kg) | Tables 2-9, row `V1`, columns GHI-L / GHI-M / GHI-H |
| `lvp_ghil` / `lvp_ghim` / `lvp_ghih` (V2, L/kg) | Tables 2-9, row `V2` |
| `lcl_ghil` / `lcl_ghim` / `lcl_ghih` (CL1, L/h/kg) | Tables 2-9, row `CL1` |
| `lq_ghil` / `lq_ghim` / `lq_ghih` (Q, L/h/kg) | Tables 2-9, row `Q` (formula 12: `Q = k12 * V1`) |
| `kel <- cl / vc`, `k12 <- q / vc`, `k21 <- q / vp` | Chen 2024 formulas 5, 6 and 12 (`CL = k * V`, `Q = k12 * V1`) |
| `E = Emax * C^gamma / (ED50^gamma + C^gamma)` | Chen 2024 Section 3.6, PK/PD equation form |
| `lemax_ckmb_*`, `lec50_ckmb_*`, `lhill_ckmb_*` | Table 12 (delta CK-MB), row for the analyte within the dose-group block |
| `lemax_ima_*`, `lec50_ima_*`, `lhill_ima_*` | Table 13 (delta IMA) |
| `lemax_ctni_*`, `lec50_ctni_*`, `lhill_ctni_*` | Table 14 (delta cTn I) |
| `lemax_hbdh_*`, `lec50_hbdh_*`, `lhill_hbdh_*` | Table 15 (delta alpha-HBDH) |
| Which analyte/biomarker/dose combinations were fitted at all | Table 10 (PLSR equations) via Section 3.6: only analytes with a negative PLSR coefficient were carried into Tables 12-15 |
| Analyte content in GHI (used to convert the GHI volume dose) | Chen 2024 Results 3.2 |
| `DOSE_GHI_MLKG` levels 2.5 / 5 / 10 mL/kg | Chen 2024 Section 2.3 |
| Residual error, between-subject variability | Not reported anywhere in Chen 2024; fixed at zero (see Assumptions) |

# Virtual cohort

Chen 2024 reports no between-subject variability and no residual-error model,
so every packaged model is deterministic. A "cohort" is therefore one typical
rat per analyte per dose group - 24 profiles in total. Simulating more subjects
would reproduce the same curve 200 times and add nothing.

```{r cohort}
set.seed(20240315)

obs_times_min <- c(2, 5, 10, 20, 40, 60, 90, 120, 240, 360)
# Dense grid for NCA / plotting, augmented with the paper's nominal sampling times.
grid_h <- sort(unique(c(
  obs_times_min / 60,
  seq(0, 6, length.out = 241)
)))

make_profile <- function(row, id) {
  ev <- rxode2::et(amt = row$amt, cmt = "central") |>
    rxode2::et(grid_h, cmt = "Cc") |>
    as.data.frame()
  ev$id <- id
  ev$DOSE_GHI_MLKG <- row$ghi_dose
  ev$group <- as.character(row$group)
  ev$model <- row$model
  ev$abbrev <- row$abbrev
  ev
}

events <- dplyr::bind_rows(
  lapply(seq_len(nrow(dosing)), function(i) make_profile(dosing[i, ], id = i))
)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
nrow(events)
```

# Simulation

Each model is solved over its own three dose-group profiles. `keep =` carries
the grouping columns through `rxSolve()` so no post-hoc join is needed.

```{r simulate}
simulate_one <- function(model_name) {
  mod <- readModelDb(model_name)
  ev <- events[events$model == model_name, , drop = FALSE]
  out <- rxode2::rxSolve(
    mod, ev,
    keep = c("group", "abbrev"),
    returnType = "data.frame"
  )
  # rxSolve must return every profile it was given.
  stopifnot(length(unique(out$id)) == length(unique(ev$id)))
  out$model <- model_name
  out
}

sims <- lapply(guhong_models, simulate_one)
names(sims) <- guhong_models

sim <- dplyr::bind_rows(lapply(sims, function(x) {
  x[, c("id", "time", "Cc", "group", "abbrev", "model")]
})) |>
  dplyr::mutate(group = factor(group, levels = c("GHI-L", "GHI-M", "GHI-H")))

dplyr::glimpse(sim)
```

## Structural check: C(0) equals dose / V1

For an intravenous bolus into a two-compartment system the initial
concentration is exactly `dose / V1`. This is a closed-form identity, so it is
checked to numerical tolerance rather than eyeballed.

```{r c0-check}
published_v1 <- tibble::tribble(
  ~abbrev, ~group,   ~V1,
  "CGA",   "GHI-L",  0.672,   "CGA",   "GHI-M", 0.496,  "CGA",   "GHI-H", 0.389,
  "HSYA",  "GHI-L",  0.218,   "HSYA",  "GHI-M", 0.126,  "HSYA",  "GHI-H", 0.498,
  "HYP",   "GHI-L",  0.493,   "HYP",   "GHI-M", 0.340,  "HYP",   "GHI-H", 0.309,
  "K-3-G", "GHI-L",  0.540,   "K-3-G", "GHI-M", 0.590,  "K-3-G", "GHI-H", 0.475,
  "K-3-R", "GHI-L",  0.546,   "K-3-R", "GHI-M", 0.603,  "K-3-R", "GHI-H", 0.490,
  "NAG",   "GHI-L",  0.223,   "NAG",   "GHI-M", 0.024,  "NAG",   "GHI-H", 0.339,
  "pCA",   "GHI-L",  0.462,   "pCA",   "GHI-M", 0.677,  "pCA",   "GHI-H", 0.723,
  "RT",    "GHI-L",  0.379,   "RT",    "GHI-M", 0.466,  "RT",    "GHI-H", 0.239
)

c0_check <- sim |>
  dplyr::filter(time == 0) |>
  dplyr::select(abbrev, group, Cc) |>
  dplyr::mutate(group = as.character(group)) |>
  dplyr::left_join(dosing |>
                     dplyr::mutate(group = as.character(group)) |>
                     dplyr::select(abbrev, group, amt),
                   by = c("abbrev", "group")) |>
  dplyr::left_join(published_v1, by = c("abbrev", "group")) |>
  dplyr::mutate(expected = amt / V1, rel_err = abs(Cc - expected) / expected)

stopifnot(max(c0_check$rel_err) < 1e-8)
sprintf("C(0) = dose/V1 for all %d profiles (max relative error %.1e)",
        nrow(c0_check), max(c0_check$rel_err))
```

# Replicate published figures

## Figure 2 - plasma concentration-time curves

```{r figure-2, fig.height = 7}
# Replicates Figure 2 of Chen 2024: plasma concentration-time curves of the
# eight analytes in the three GHI dose groups.
sim |>
  dplyr::filter(time > 0, time <= 6) |>
  dplyr::left_join(analytes |> dplyr::select(abbrev, analyte, conc_units),
                   by = "abbrev") |>
  dplyr::mutate(panel = paste0(analyte, " (", conc_units, ")")) |>
  ggplot(aes(time * 60, Cc, colour = group)) +
  geom_line(linewidth = 0.6) +
  facet_wrap(~panel, scales = "free_y", ncol = 2) +
  scale_x_continuous(breaks = c(0, 60, 120, 240, 360)) +
  labs(x = "Time (min)", y = "Plasma concentration",
       colour = "Dose group",
       title = "Figure 2 - plasma concentration-time curves",
       caption = "Replicates Figure 2 of Chen 2024 (typical-value profiles; the paper plots observed mean +/- SD).") +
  theme(legend.position = "bottom")
```

The simulated curves are the model's typical-value predictions, whereas Chen
2024 Figure 2 plots the observed group means. The rank order of analytes
(N-acetyl-L-glutamine highest, in ug/mL; hyperoside lowest, in ng/mL), the
dose ordering within each panel, and the biphasic shape all reproduce. The
simulated concentration at the first sampling time (2 min) runs roughly 2-3x
below the plotted observed mean for most analytes; see Assumptions below.

## Figure 4 - concentration-effect relationships

```{r figure-4, fig.height = 6}
# Replicates the concentration-effect relationship underlying Figure 4 of
# Chen 2024, for the four analytes the paper identified as the
# pharmacodynamic substance basis of GHI.
pd_models <- c("Chen_2024_nAcetylglutamine_rat",
               "Chen_2024_hydroxysafflorYellowA_rat",
               "Chen_2024_chlorogenicAcid_rat",
               "Chen_2024_pCoumaricAcid_rat")

pd_long <- dplyr::bind_rows(lapply(pd_models, function(mn) {
  x <- sims[[mn]]
  keep <- intersect(c("dCKMB", "dIMA", "dCTNI", "dHBDH"), names(x))
  out <- x[, c("time", "Cc", "group", "abbrev", keep), drop = FALSE]
  tidyr::pivot_longer(out, cols = dplyr::all_of(keep),
                      names_to = "biomarker", values_to = "effect")
})) |>
  dplyr::mutate(
    group = factor(group, levels = c("GHI-L", "GHI-M", "GHI-H")),
    biomarker = factor(
      biomarker,
      levels = c("dCKMB", "dIMA", "dCTNI", "dHBDH"),
      labels = c("delta CK-MB (ng/mL)", "delta IMA (U/mL)",
                 "delta cTn I (ng/mL)", "delta alpha-HBDH (ng/mL)")
    )
  )

pd_long |>
  dplyr::filter(time > 0, time <= 6, effect > 0) |>
  ggplot(aes(Cc, effect, colour = group)) +
  geom_path(linewidth = 0.6) +
  facet_grid(biomarker ~ abbrev, scales = "free") +
  scale_x_log10() +
  labs(x = "Plasma concentration (ng/mL for CGA and pCA; ug/mL for HSYA and NAG)",
       y = "Effect: GHI group minus MI/R model group",
       colour = "Dose group",
       title = "Concentration-effect curves for the four PLSR-selected analytes",
       caption = paste("Replicates the PK/PD relationship of Chen 2024 Figure 4 and Tables 12-15.",
                       "Blank panels are combinations for which no equation was fitted.")) +
  theme(legend.position = "bottom")
```

## Effect magnitudes versus the observed biomarker differences

Chen 2024 Figure 3 plots each biomarker for the model group and the three GHI
groups. The observed reduction that the GHI groups achieve at late times sets
the scale that `Emax` must be able to reach. The check below confirms that the
maximum simulated effect over the 6 h profile never exceeds the fitted `Emax`,
and reports it against the biomarker's observed range.

```{r effect-magnitude}
emax_published <- tibble::tribble(
  ~biomarker,               ~observed_model_group, ~observed_units,
  "delta CK-MB (ng/mL)",    32,                    "ng/mL",
  "delta IMA (U/mL)",       115,                   "U/mL",
  "delta cTn I (ng/mL)",    18,                    "ng/mL",
  "delta alpha-HBDH (ng/mL)", 160,                 "ng/mL"
)

pd_long |>
  dplyr::filter(time > 0) |>
  dplyr::group_by(biomarker, abbrev, group) |>
  dplyr::summarise(max_effect = max(effect), .groups = "drop") |>
  dplyr::filter(max_effect > 0) |>
  dplyr::group_by(biomarker) |>
  dplyr::summarise(
    n_curves = dplyr::n(),
    min_max_effect = min(max_effect),
    max_max_effect = max(max_effect),
    .groups = "drop"
  ) |>
  dplyr::left_join(emax_published, by = "biomarker") |>
  dplyr::rename(
    "Biomarker"                       = biomarker,
    "Fitted curves"                   = n_curves,
    "Smallest peak effect"            = min_max_effect,
    "Largest peak effect"             = max_max_effect,
    "MI/R model group level (Fig. 3)" = observed_model_group,
    "Units"                           = observed_units
  ) |>
  knitr::kable(
    digits = 1,
    caption = paste("Simulated peak effect per biomarker against the MI/R model-group level",
                    "read from Chen 2024 Figure 3. Every peak effect is a plausible fraction",
                    "of the model-group level, confirming the Emax scale and units.")
  )
```

# PKNCA validation

Non-compartmental analysis of the simulated profiles, grouped by analyte and
dose group, compared against the published half-life and AUC.

```{r pknca}
# The observation grid starts at time 0, so every profile already carries a
# time-zero record; for an intravenous bolus its value is dose/V1, not zero.
# Filter on !is.na() only - a `time > 0` or `Cc > 0` filter would drop it and
# trigger PKNCA's "AUC range starting before the first measurement" warning.
sim_nca <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::mutate(treatment = paste(abbrev, group, sep = " / ")) |>
  dplyr::select(id, time, Cc, treatment) |>
  dplyr::arrange(id, treatment, time)

stopifnot(all(tapply(sim_nca$time, sim_nca$id, min) == 0))

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

dose_df <- events |>
  dplyr::filter(evid == 1) |>
  dplyr::mutate(treatment = paste(abbrev, group, sep = " / ")) |>
  dplyr::select(id, time, amt, treatment)

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id)

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

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

## Comparison against published NCA

Chen 2024 reports `t1/2beta` and `AUC(0-inf)` for every analyte in every dose
group (Tables 2-9). The terminal half-life is the direct test of the fitted
disposition parameters.

```{r nca-compare}
published_nca <- tibble::tribble(
  ~treatment,       ~half.life, ~aucinf.obs,
  "CGA / GHI-L",    0.919,      102.146,
  "CGA / GHI-M",    0.829,      277.105,
  "CGA / GHI-H",    0.520,      421.197,
  "HSYA / GHI-L",   0.732,       40.679,
  "HSYA / GHI-M",   1.386,      115.716,
  "HSYA / GHI-H",   0.972,      190.363,
  "HYP / GHI-L",    0.669,        8.252,
  "HYP / GHI-M",    0.312,       19.295,
  "HYP / GHI-H",    0.313,       34.504,
  "K-3-G / GHI-L",  0.674,       25.109,
  "K-3-G / GHI-M",  0.691,       35.448,
  "K-3-G / GHI-H",  0.385,      122.652,
  "K-3-R / GHI-L",  0.563,      172.917,
  "K-3-R / GHI-M",  0.621,      236.738,
  "K-3-R / GHI-H",  0.348,      819.682,
  "NAG / GHI-L",    0.814,      553.382,
  "NAG / GHI-M",    0.462,     1060.065,
  "NAG / GHI-H",    0.709,     1079.308,
  "pCA / GHI-L",    0.487,      400.705,
  "pCA / GHI-M",    0.392,      969.468,
  "pCA / GHI-H",    0.321,     1579.385,
  "RT / GHI-L",     0.541,       36.886,
  "RT / GHI-M",     0.492,      129.725,
  "RT / GHI-H",     0.348,      217.984
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_res,
  reference = published_nca,
  by        = "treatment",
  params    = c("half.life", "aucinf.obs"),
  units     = c(half.life = "h", aucinf.obs = "conc*h"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = paste("Simulated vs. published NCA for all 24 analyte/dose-group profiles.",
                  "* differs from the published value by more than 20%.")
)
```

### Half-life: agreement

The simulated terminal half-life reproduces the published `t1/2beta` closely
for every profile. Quantitatively:

```{r half-life-agreement}
hl <- as.data.frame(nca_res) |>
  dplyr::filter(PPTESTCD == "half.life") |>
  dplyr::select(treatment, simulated = PPORRES) |>
  dplyr::left_join(published_nca |> dplyr::select(treatment, published = half.life),
                   by = "treatment") |>
  dplyr::mutate(pct_diff = 100 * (simulated - published) / published)

sprintf(
  "Terminal half-life: median absolute difference %.1f%%, maximum %.1f%% (n = %d profiles)",
  stats::median(abs(hl$pct_diff)), max(abs(hl$pct_diff)), nrow(hl)
)

hl |>
  dplyr::arrange(dplyr::desc(abs(pct_diff))) |>
  dplyr::rename("Analyte / dose group" = treatment,
                "Simulated t1/2 (h)"   = simulated,
                "Published t1/2beta (h)" = published,
                "Difference (%)"       = pct_diff) |>
  knitr::kable(digits = c(0, 3, 3, 1),
               caption = "Simulated terminal half-life vs. Chen 2024 Tables 2-9.")
```

The half-life is a pure function of `V1`, `V2`, `CL1` and `Q`, so this
agreement confirms that the four disposition parameters transcribed from
Tables 2-9 are mutually consistent and that the packaged parameterisation
(`kel = CL1/V1`, `k12 = Q/V1`, `k21 = Q/V2`) is the one the authors fitted.

### AUC: a published-value inconsistency, not a model error

Every `aucinf.obs` row is starred. The simulated AUC is, by construction,
`dose / CL1` - the closed-form result for an intravenous bolus - and PKNCA
recovers it to within a couple of percent (the residual is the trapezoidal
and terminal-extrapolation error of a numerical NCA on a finite grid):

```{r auc-identity}
auc <- as.data.frame(nca_res) |>
  dplyr::filter(PPTESTCD == "aucinf.obs") |>
  dplyr::select(treatment, simulated = PPORRES)

published_cl1 <- tibble::tribble(
  ~abbrev, ~group,  ~CL1,
  "CGA",   "GHI-L", 2.669,   "CGA",   "GHI-M", 1.939,  "CGA",   "GHI-H", 2.554,
  "HSYA",  "GHI-L", 1.220,   "HSYA",  "GHI-M", 0.831,  "HSYA",  "GHI-H", 1.130,
  "HYP",   "GHI-L", 5.336,   "HYP",   "GHI-M", 4.595,  "HYP",   "GHI-H", 4.915,
  "K-3-G", "GHI-L", 5.255,   "K-3-G", "GHI-M", 4.123,  "K-3-G", "GHI-H", 4.4558,
  "K-3-R", "GHI-L", 5.158,   "K-3-R", "GHI-M", 4.134,  "K-3-R", "GHI-H", 4.481,
  "NAG",   "GHI-L", 1.577,   "NAG",   "GHI-M", 1.133,  "NAG",   "GHI-H", 1.633,
  "pCA",   "GHI-L", 3.317,   "pCA",   "GHI-M", 2.808,  "pCA",   "GHI-H", 3.469,
  "RT",    "GHI-L", 3.795,   "RT",    "GHI-M", 2.872,  "RT",    "GHI-H", 2.813
)

auc_check <- dosing |>
  dplyr::mutate(group = as.character(group)) |>
  dplyr::left_join(published_cl1, by = c("abbrev", "group")) |>
  dplyr::mutate(treatment = paste(abbrev, group, sep = " / "),
                dose_over_cl1 = amt / CL1) |>
  dplyr::left_join(auc, by = "treatment") |>
  dplyr::left_join(published_nca |> dplyr::select(treatment, published = aucinf.obs),
                   by = "treatment") |>
  dplyr::mutate(
    sim_vs_identity_pct = 100 * (simulated - dose_over_cl1) / dose_over_cl1,
    published_over_sim  = published / simulated
  )

sprintf(
  "Simulated AUCinf equals dose/CL1 to within %.2f%% for all %d profiles.",
  max(abs(auc_check$sim_vs_identity_pct)), nrow(auc_check)
)

auc_check |>
  dplyr::select(treatment, amt, CL1, dose_over_cl1, simulated, published, published_over_sim) |>
  dplyr::rename(
    "Analyte / dose group"      = treatment,
    "Dose"                      = amt,
    "Published CL1 (L/h/kg)"    = CL1,
    "dose / CL1"                = dose_over_cl1,
    "Simulated AUCinf"          = simulated,
    "Published AUC(0-inf)"      = published,
    "Published / simulated"     = published_over_sim
  ) |>
  knitr::kable(digits = c(0, 2, 3, 1, 1, 1, 1),
               caption = paste("The simulated AUC matches the closed-form dose/CL1 identity to within",
                               "numerical NCA error. The published AUC column does not: it exceeds",
                               "dose/CL1 by a factor that varies by analyte.",
                               "See Assumptions and deviations."))
```

The ratio in the last column is roughly 5-9x for the six analytes reported in
ng/mL and roughly 11-20x for N-acetyl-L-glutamine and hydroxysafflor yellow A -
the two analytes whose plasma concentrations Chen 2024 plots in ug/mL (Figure
2B and 2F) while labelling the AUC rows of Tables 3 and 7 as `ug/L*h`. A
1000-fold unit mislabelling accounts for most of the gap for those two, but
after correcting for it a residual factor of order 10 remains for all eight.
The published `AUC(0-inf)` column therefore cannot be reconciled with the
published `CL1` column, and the packaged model follows `CL1` (which, unlike
AUC, is internally consistent with `V1`, `V2`, `Q`, `t1/2alpha` and
`t1/2beta`). This is a defect in the source table, not in the extraction; no
parameter was tuned.

# Assumptions and deviations

- **Parameterisation.** Chen 2024 Tables 2-9 report both a macro-parameter set
  (`V1`, `V2`, `CL1`, `Q`) and a micro-constant set (`K10`, `K12`, `K21`). They
  are not exactly mutually consistent, because each is the arithmetic mean of
  six individually fitted rats and the mapping between them is non-linear
  (`mean(K10) * mean(V1) != mean(CL1)`). The packaged models use the
  macro-parameter set, which reproduces the published `t1/2alpha` and
  `t1/2beta` markedly better: for chlorogenic acid GHI-L, `V1/V2/CL1/Q` gives
  `t1/2alpha = 0.034 h` and `t1/2beta = 0.917 h` against published values of
  0.034 h and 0.919 h, whereas `V1/K10/K12/K21` gives 0.032 h and 0.903 h.
  For hydroxysafflor yellow A GHI-L the micro-constant route gives
  `t1/2beta = 0.40 h` against a published 0.732 h, while the macro route gives
  0.70 h.

- **`AUC(0-t)` and `AUC(0-inf)` are not reproducible from the published
  parameters.** See the PKNCA section above. The AUC columns of Tables 2-9
  exceed `dose / CL1` by an analyte-dependent factor and were not used to
  derive any packaged parameter.

- **Concentration units of two analytes are mislabelled in the tables.**
  Chen 2024 Figure 2B and 2F plot hydroxysafflor yellow A and
  N-acetyl-L-glutamine in ug/mL, consistent with their content in GHI (1 and
  30 mg/mL) and with the 50- and 200-fold sample dilutions described in Section
  2.6.2, but the AUC rows of Tables 3 and 7 are labelled `ug/L*h` like the
  other six analytes. The packaged models use ug/mL for these two and ng/mL for
  the other six; the `ED50` values in Tables 12-15 are consistent with that
  reading (for example the N-acetyl-L-glutamine `ED50` values of 23.6-278 fall
  inside its ug/mL concentration range, and would be far below the measurable
  range if read as ng/mL).

- **No between-subject variability and no residual error were reported.** Chen
  2024 fitted each rat individually in DAS v3.2.6 and tabulates the mean +/- SD
  of the resulting parameter estimates; it publishes no random-effects model,
  no `OMEGA`, and no residual-error term. Every parameter is therefore
  `fixed()` at the published mean and every residual SD is `fixed(0)`, per the
  library convention for deterministic typical-value models. The published SDs
  are estimate-to-estimate spread across six rats, which conflates true
  between-animal variability with per-animal estimation error, so they were not
  reinterpreted as an `omega`.

- **Dose-group parameters are selected, not interpolated.** The three dose
  groups were fitted independently; Chen 2024 estimated no dose-covariate
  function. The models therefore switch parameter sets on `DOSE_GHI_MLKG` at
  the midpoints of the studied levels (3.75 and 7.5 mL/kg). A `DOSE_GHI_MLKG`
  outside 2.5-10 mL/kg extrapolates the nearest fitted group and is not
  supported by the source data.

- **Unfitted analyte/biomarker/dose combinations return an effect of zero.**
  Chen 2024 carried an analyte into the PK/PD tables only where its PLSR
  regression coefficient against that biomarker in that dose group was negative
  (Section 2.9 and Table 10). For the remaining combinations the paper gives no
  `Emax`, `ED50` or `gamma`, so the packaged models contribute nothing for them
  rather than extrapolating a neighbouring dose group's curve. Hyperoside had
  no negative coefficient anywhere and therefore carries the PK layer only.
  Blank panels in the Figure 4 replication above are exactly these cases.

- **Several PK/PD parameters sit at estimation boundaries.** `gamma` is
  reported as 0.01 for eight of the 45 fitted curves and `ED50` as 0.01 for
  three of them - values that flatten the sigmoid to a near-constant effect
  and that look like the lower bound of the DAS search space rather than
  informative estimates. They are packaged as published; users fitting
  derivative models should treat those curves as unidentified.

- **The simulated 2-minute concentration runs below the observed group mean.**
  The typical-value profiles reach roughly a third to two thirds of the
  concentration plotted at the first sampling time in Chen 2024 Figure 2, while
  `C(0) = dose / V1` matches the extrapolated intercept well (for p-coumaric
  acid GHI-L, `C(0) = 424 ng/mL` against an observed 2-minute mean of about
  440 ng/mL). The published `t1/2alpha` values of 0.005-0.045 h place the
  entire distribution phase before the first 2-minute sample, so the fitted
  alpha phase decays faster than the observed data can resolve. This is a
  property of the published fit, not of the extraction.

- **The PLSR layer is not encoded.** The paper's headline method - partial
  least squares regression of the eight analyte concentrations against each
  biomarker (Table 10) - is a data-analysis step used to decide which
  analyte/biomarker pairs to fit, not a dynamical model. Its coefficients are
  recorded in the model descriptions and in the source-trace table above, and
  they determine which sigmoid-Emax parameters exist, but there is nothing in
  it to express as an ODE.

- **No parameter value came from outside the paper.** Every `ini()` entry
  traces to Chen 2024 Tables 2-9 (disposition) or Tables 12-15 (PK/PD); the
  dose conversion uses the analyte contents in Results 3.2 and the GHI dose
  levels in Section 2.3. Figure 2 and Figure 3 axis labels were read to settle
  concentration and biomarker units, and Figure 3 was read to sanity-check the
  `Emax` scale, but no numeric parameter was digitised from a figure.
