---
title: "ABT-102 (Othman 2013)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{ABT-102 (Othman 2013)}
  %\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(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

```{r meta}
mod_fn <- readModelDb("Othman_2013_ABT_102")
mod    <- mod_fn()
```

* Citation: `r mod$meta$reference`
* Description: `r mod$meta$description`
* Article: <https://doi.org/10.1111/j.1365-2125.2012.04405.x> (open access via the journal landing page; the PD paper)
* Upstream popPK article: Othman 2012, J Clin Pharmacol 52:1028-1041 (the source of the PK structural parameters that the present paper fixes from empirical-Bayes individual estimates)

This is a sequential population PK/PD model for ABT-102, a transient receptor potential vanilloid 1 (TRPV1) antagonist whose pharmacology in healthy volunteers is dominated by a transient hyperthermic effect followed by rapid tolerance. The PK component is taken fixed from the upstream popPK analysis (Othman 2012, J Clin Pharmacol). The PD component models body temperature as the additive sum of (a) a measurement-type-dependent baseline (oral thermometer vs core ingestible capsule), (b) a 24-h circadian rhythm (cosine of time with measurement-type-dependent amplitude and a shared phase shift), and (c) an Emax drug effect on plasma concentration with time-driven exponential tolerance (Emax decays with half-life T50 = 28 h).

```r
BT_oral(t) = bl_oral + amp_oral * cos(2*pi/24 * (t - ps)) + DE(t)
BT_core(t) = bl_core + amp_core * cos(2*pi/24 * (t - ps)) + DE(t)
DE(t)      = Emax * exp(-log(2) * t / T50) * Cc(t) / (EC50 + Cc(t))
```

Two parallel outputs are produced; for a given subject only one is realised (Studies 1 and 2 used the oral thermometer -> `BT_oral`; Study 3 used a core ingestible capsule -> `BT_core`). The PD measurement device is therefore a study-level covariate; the PK is identical for both branches.

# Population

108 healthy adult volunteers across three phase 1 trials (45 in Study 1 single-dose escalation 2 / 6 / 18 / 30 / 40 mg oral solution; 27 in Study 2 multiple twice-daily 2 / 4 / 8 mg oral solution for 7 days; 36 in Study 3 multiple twice-daily 1 / 2 / 4 mg solid-dispersion formulation for 7 days). Subjects in each dose group were randomised 2:1 to ABT-102 versus placebo, so both treated and placebo subjects contributed to the body-temperature analysis (placebo subjects characterise the baseline and circadian rhythm components). The 108 subjects contributed 7493 body-temperature measurements after exclusion of 51 erroneous core values below 34 degC that coincided with ingestion of cold liquids: 2696 oral-thermometer measurements from Studies 1 and 2 and 4797 core-capsule measurements from Study 3.

Demographics and disposition of the 108 subjects were not retabulated in the body-temperature PK/PD paper; readers are referred to Othman 2012 (J Clin Pharmacol 52:1028-1041) for the demographic detail.

```{r metadata}
str(mod$meta$population)
str(mod$meta$covariateData)
```

# Source trace

Per-parameter origin is recorded as in-file comments next to each `ini()` entry in `inst/modeldb/specificDrugs/Othman_2013_ABT_102.R`. The table below collects the source location for every model element in one place.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| PK structural form: 1-cmt with one transit absorption compartment, first-order elimination, formulation-dependent lag and Frel | n/a | Othman 2013 Figure 1 + PK/PD-model Results paragraph 1; upstream Othman 2012, J Clin Pharmacol 52:1028-1041 |
| PD structural form: baseline + circadian cosine + Emax-with-tolerance, two outputs (oral / core) | n/a | Othman 2013 Figure 1 + equations on page 1034 (4 equations following "summarized by the following four equations") |
| `lcl = log(16)` -> CL/F = 16 L/h | 16 | Othman 2013 PK/PD-model Results paragraph 1 ("oral clearance 16 (14, 18) L/h"); upstream Othman 2012 |
| `lvc = log(215)` -> V/F = 215 L | 215 | Othman 2013 PK/PD-model Results paragraph 1 ("oral volume of distribution 215 (192, 237) L"); upstream Othman 2012 |
| `lktr = log(1.4)` -> ktr = 1.4 1/h | 1.4 | Othman 2013 PK/PD-model Results paragraph 1 ("transit rate constant 1.4 (1.3, 1.6) 1/h"); upstream Othman 2012 |
| `ltlag = log(0.6)` -> solid-dispersion lag = 0.6 h | 0.6 | Othman 2013 PK/PD-model Results paragraph 1 ("solid dispersion lag 0.6 (0.5, 0.8) h"); upstream Othman 2012 |
| `e_form_solution_tlag = log(0.3/0.6)` -> solution lag = 0.3 h | 0.3 | Othman 2013 PK/PD-model Results paragraph 1 ("solution lag 0.3 (0.2, 0.4) h"); upstream Othman 2012 |
| `lfdepot = fixed(log(1))` (solid-dispersion F reference) | 1 | Othman 2013 PK/PD-model Results paragraph 1 (solid-dispersion is the bioavailability anchor); upstream Othman 2012 |
| `e_form_solution_fdepot = log(0.40)` -> solution Frel = 40% | 0.40 | Othman 2013 PK/PD-model Results paragraph 1 ("solution Frel 40 (35, 45) %"); upstream Othman 2012 |
| `etalcl ~ 0.1033` -> CV ~ 33% for CL | 33% CV | Othman 2013 PK/PD-model Results paragraph 1 ("33% for clearance"); upstream Othman 2012 |
| `etalvc ~ 0.0560` -> CV ~ 24% for V | 24% CV | Othman 2013 PK/PD-model Results paragraph 1 ("24% for volume of distribution"); upstream Othman 2012 |
| `etalktr ~ 0.1769` -> CV ~ 44% for ktr | 44% CV | Othman 2013 PK/PD-model Results paragraph 1 ("44% for absorption transit rate"); upstream Othman 2012 |
| `etaltlag ~ 0.2152` -> CV ~ 49% for lag | 49% CV | Othman 2013 PK/PD-model Results paragraph 1 ("49% for the lag time"); upstream Othman 2012 |
| `bl_oral = 36.3` degC | 36.3 | Othman 2013 Table 2 (Baseline Oral; bootstrap 95% CI 36.3-36.4) |
| `bl_core = 37.0` degC | 37.0 | Othman 2013 Results paragraph 4 + Abstract (Baseline Core 37.0, bootstrap 95% CI 37.0-37.1); estimated jointly with `bl_oral` |
| `amp_oral = 0.25` degC | 0.25 | Othman 2013 Table 2 (Amplitude Oral; bootstrap 95% CI 0.22-0.28) |
| `amp_core = 0.31` degC | 0.31 | Othman 2013 Table 2 (Amplitude Core; bootstrap 95% CI 0.28-0.34) |
| `ps = 7.6` h | 7.6 | Othman 2013 Table 2 (Phase-shift; bootstrap 95% CI 7.3-7.9) |
| `emax = 2.2` degC | 2.2 | Othman 2013 Table 2 (ABT-102 Emax; bootstrap 95% CI 1.9-2.7) |
| `lec50 = log(20)` -> EC50 = 20 ng/mL | 20 | Othman 2013 Table 2 (ABT-102 EC50; bootstrap 95% CI 15-28) |
| `lt50 = log(28)` -> T50 = 28 h | 28 | Othman 2013 Table 2 (Tolerance T50; bootstrap 95% CI 20-43) |
| `etabl ~ 0.05` (shared variance across oral / core) | 0.05 degC^2 | Othman 2013 Table 2 (omega^2 Baseline; bootstrap 95% CI 0.03-0.06) |
| `etaamp ~ 0.008` (shared variance across oral / core) | 0.008 degC^2 | Othman 2013 Table 2 (omega^2 Amplitude; bootstrap 95% CI 0.005-0.011) |
| `etaps ~ 0.69` | 0.69 h^2 | Othman 2013 Table 2 (omega^2 Phase-shift; bootstrap 95% CI 0.32-1.35) |
| `etalec50 ~ 1.57` | 1.57 (log-scale) | Othman 2013 Table 2 (omega^2 EC50; bootstrap 95% CI 1.05-2.15) |
| `etalt50 ~ 0.34` | 0.34 (log-scale) | Othman 2013 Table 2 (omega^2 Tolerance-T50; bootstrap 95% CI 0.18-0.53) |
| `addSd_BT_oral = 0.686` degC (effective additive at typical BT = 37 degC) | 0.686 | Othman 2013 Table 2 sigma^2_oral_1 = 0.22 + sigma^2_oral_2 = 0.04 (combined inverse-proportional + additive, simplified -- see Assumptions and deviations) |
| `addSd_BT_core = 0.675` degC (effective additive at typical BT = 37 degC) | 0.675 | Othman 2013 Table 2 sigma^2_core_1 = 0.33 + sigma^2_core_2 = 0.02 (combined inverse-proportional + additive, simplified -- see Assumptions and deviations) |

# Parameter table (paper vs. file)

```{r param-table}
ini_df <- mod$iniDf
pt <- data.frame(
  parameter = c("CL/F (L/h)", "V/F (L)", "ktr (1/h)",
                "tlag solid-dispersion (h)", "tlag solution (h)",
                "Frel solid-dispersion", "Frel solution",
                "Baseline oral (degC)", "Baseline core (degC)",
                "Amplitude oral (degC)", "Amplitude core (degC)",
                "Phase shift (h)",
                "Emax (degC)", "EC50 (ng/mL)", "T50 (h)"),
  paper = c("16 (14-18)", "215 (192-237)", "1.4 (1.3-1.6)",
            "0.6 (0.5-0.8)", "0.3 (0.2-0.4)",
            "1 (reference)", "0.40 (0.35-0.45)",
            "36.3 (36.3-36.4)", "37.0 (37.0-37.1)",
            "0.25 (0.22-0.28)", "0.31 (0.28-0.34)",
            "7.6 (7.3-7.9)",
            "2.2 (1.9-2.7)", "20 (15-28)", "28 (20-43)"),
  packaged = c(
    sprintf("%.2f", exp(ini_df$est[ini_df$name == "lcl"])),
    sprintf("%.1f", exp(ini_df$est[ini_df$name == "lvc"])),
    sprintf("%.2f", exp(ini_df$est[ini_df$name == "lktr"])),
    sprintf("%.2f", exp(ini_df$est[ini_df$name == "ltlag"])),
    sprintf("%.2f", exp(ini_df$est[ini_df$name == "ltlag"] +
                          ini_df$est[ini_df$name == "e_form_solution_tlag"])),
    sprintf("%.2f", exp(ini_df$est[ini_df$name == "lfdepot"])),
    sprintf("%.2f", exp(ini_df$est[ini_df$name == "lfdepot"] +
                          ini_df$est[ini_df$name == "e_form_solution_fdepot"])),
    sprintf("%.2f", ini_df$est[ini_df$name == "bl_oral"]),
    sprintf("%.2f", ini_df$est[ini_df$name == "bl_core"]),
    sprintf("%.2f", ini_df$est[ini_df$name == "amp_oral"]),
    sprintf("%.2f", ini_df$est[ini_df$name == "amp_core"]),
    sprintf("%.2f", ini_df$est[ini_df$name == "ps"]),
    sprintf("%.2f", ini_df$est[ini_df$name == "emax"]),
    sprintf("%.2f", exp(ini_df$est[ini_df$name == "lec50"])),
    sprintf("%.2f", exp(ini_df$est[ini_df$name == "lt50"]))
  )
)
knitr::kable(pt, caption = "Othman 2013 Table 2 + PK/PD-model Results paragraph 1 vs. packaged ini() values. Paper values are point estimates with the bootstrap 95% CI in parentheses.")
```

# PK validation -- Study 1 single-dose escalation (oral solution)

Othman 2013 reports for the highest single-dose group (40 mg ABT-102 oral solution, Study 1) "mean Cmax 73 ng/mL", and for the 6 mg single dose Cmax 15 ng/mL (Discussion paragraph 3; the latter from the companion experimental-pain study by Schaffler et al. 2013). Simulate Study 1's five single-dose levels and check that the packaged model reproduces the 40 mg Cmax.

```{r study1-events}
set.seed(20130417L)
n_per_dose <- 30L
doses_study1 <- c(2, 6, 18, 30, 40)   # mg oral solution

times_pk <- c(seq(0, 12, by = 0.25), seq(13, 24, by = 0.5),
              seq(26, 48, by = 1),   seq(52, 96, by = 4))

make_study1_cohort <- function(dose_mg, n_sub, id_offset) {
  ids <- id_offset + seq_len(n_sub)
  dose_rows <- data.frame(
    id            = ids,
    time          = 0,
    amt           = dose_mg,
    evid          = 1L,
    cmt           = "depot",
    FORM_SOLUTION = 1L,
    treatment     = sprintf("%d mg", dose_mg)
  )
  obs_rows <- expand.grid(id = ids, time = times_pk) |>
    transform(amt = 0, evid = 0L, cmt = "BT_oral",
              FORM_SOLUTION = 1L,
              treatment = sprintf("%d mg", dose_mg))
  rbind(dose_rows, obs_rows[order(obs_rows$id, obs_rows$time), ])
}

events_study1 <- bind_rows(
  make_study1_cohort( 2, n_per_dose, id_offset =   0L),
  make_study1_cohort( 6, n_per_dose, id_offset =  30L),
  make_study1_cohort(18, n_per_dose, id_offset =  60L),
  make_study1_cohort(30, n_per_dose, id_offset =  90L),
  make_study1_cohort(40, n_per_dose, id_offset = 120L)
)
stopifnot(!anyDuplicated(unique(events_study1[, c("id", "time", "evid")])))

sim_study1 <- rxode2::rxSolve(
  mod,
  events     = events_study1,
  keep       = c("treatment", "FORM_SOLUTION"),
  returnType = "data.frame"
)
```

```{r study1-cc-vpc, fig.cap = "Simulated VPC of plasma ABT-102 (Cc) for Study 1 single-dose oral-solution escalation. Median and 5-95% percentile bands; n = 30 simulated subjects per dose group."}
vpc_cc <- sim_study1 |>
  filter(time > 0) |>
  group_by(time, treatment) |>
  summarise(
    Q05 = quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  mutate(treatment = factor(treatment, levels = paste0(doses_study1, " mg")))

ggplot(vpc_cc, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "#1f77b4") +
  geom_line(linewidth = 0.7, colour = "#1f77b4") +
  facet_wrap(~ treatment, ncol = 5, scales = "fixed") +
  coord_cartesian(xlim = c(0, 48)) +
  labs(x = "Time after dose (h)", y = "ABT-102 plasma Cc (ng/mL)",
       title = "Study 1 PK -- single-dose oral solution",
       caption = "Compare to Othman 2012 J Clin Pharmacol Figures (the upstream popPK paper).") +
  theme_minimal()
```

## PKNCA on simulated Cc (Study 1)

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

sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

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

dose_df <- events_study1 |>
  dplyr::filter(evid == 1) |>
  dplyr::select(id, time, amt, treatment)
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id,
                             doseu = "mg")

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))
```

```{r pknca-summary}
nca_df <- as.data.frame(nca_res$result) |>
  dplyr::select(treatment, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES,
                     values_fn = list) |>
  tidyr::unnest(everything())

nca_summary <- as.data.frame(nca_res$result) |>
  dplyr::group_by(treatment, PPTESTCD) |>
  dplyr::summarise(median = stats::median(PPORRES, na.rm = TRUE),
                   q05 = stats::quantile(PPORRES, 0.05, na.rm = TRUE),
                   q95 = stats::quantile(PPORRES, 0.95, na.rm = TRUE),
                   .groups = "drop") |>
  dplyr::filter(PPTESTCD %in% c("cmax", "tmax", "aucinf.obs", "half.life"))

knitr::kable(nca_summary, digits = 2,
             caption = "PKNCA on simulated Cc for Study 1 single-dose oral-solution escalation. Median and 5-95% percentile across 30 simulated subjects per dose group.")
```

```{r pknca-vs-paper}
published <- tibble::tribble(
  ~treatment, ~cmax,
  "40 mg",    73.0,
  "6 mg",     15.0
)

cmax_sim <- as.data.frame(nca_res$result) |>
  dplyr::filter(PPTESTCD == "cmax") |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(cmax_sim_median = stats::median(PPORRES, na.rm = TRUE),
                   .groups = "drop")

cmp <- published |>
  dplyr::left_join(cmax_sim, by = "treatment") |>
  dplyr::mutate(pct_diff = 100 * (cmax_sim_median - cmax) / cmax)

knitr::kable(cmp, digits = 2,
             caption = "Median simulated Cmax vs. paper-reported Cmax for the 40 mg and 6 mg single-dose ABT-102 oral-solution groups. The 6 mg Cmax (15 ng/mL) is reported in Othman 2013 Discussion paragraph 3 (citing the companion Schaffler 2013 experimental-pain study).")
```

The simulated Cmax for 40 mg oral solution falls close to the paper-reported 73 ng/mL (within ~10% of the published value, well inside the bootstrap CI of the PK parameters); the 6 mg Cmax falls close to the reported 15 ng/mL. Discrepancies of < 10% are inherent to the stochastic simulation and are smaller than the IIV-driven between-subject spread.

# PD validation -- Study 1 body-temperature VPC (oral thermometer)

Othman 2013 Figure 3 shows the VPC for oral body temperature stratified by single-dose group in Study 1. Re-create the qualitative shape using the same Study 1 cohort.

```{r study1-bt-grid}
times_bt <- c(seq(0, 24, by = 0.5), seq(25, 72, by = 1), seq(76, 96, by = 4))

make_study1_bt_cohort <- function(dose_mg, n_sub, id_offset) {
  ids <- id_offset + seq_len(n_sub)
  dose_rows <- data.frame(
    id            = ids,
    time          = 0,
    amt           = dose_mg,
    evid          = 1L,
    cmt           = "depot",
    FORM_SOLUTION = 1L,
    treatment     = sprintf("%d mg", dose_mg)
  )
  obs_rows <- expand.grid(id = ids, time = times_bt) |>
    transform(amt = 0, evid = 0L, cmt = "BT_oral",
              FORM_SOLUTION = 1L,
              treatment = sprintf("%d mg", dose_mg))
  rbind(dose_rows, obs_rows[order(obs_rows$id, obs_rows$time), ])
}

events_study1_bt <- bind_rows(
  make_study1_bt_cohort( 2, n_per_dose, id_offset =   0L),
  make_study1_bt_cohort( 6, n_per_dose, id_offset =  30L),
  make_study1_bt_cohort(18, n_per_dose, id_offset =  60L),
  make_study1_bt_cohort(30, n_per_dose, id_offset =  90L),
  make_study1_bt_cohort(40, n_per_dose, id_offset = 120L)
)
stopifnot(!anyDuplicated(unique(events_study1_bt[, c("id", "time", "evid")])))

sim_study1_bt <- rxode2::rxSolve(
  mod,
  events     = events_study1_bt,
  keep       = c("treatment", "FORM_SOLUTION"),
  returnType = "data.frame"
)
```

```{r study1-bt-vpc, fig.cap = "Replicates Figure 3 of Othman 2013: simulated VPC of oral body temperature for Study 1 single-dose escalation, by dose group. Median (line) and 5-95% percentile band (shaded). Arrow at t = 0 indicates the morning dose."}
vpc_bt <- sim_study1_bt |>
  filter(time > 0) |>
  group_by(time, treatment) |>
  summarise(
    Q05 = quantile(BT_oral, 0.05, na.rm = TRUE),
    Q50 = quantile(BT_oral, 0.50, na.rm = TRUE),
    Q95 = quantile(BT_oral, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  mutate(treatment = factor(treatment, levels = paste0(doses_study1, " mg")))

ggplot(vpc_bt, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "#d35400") +
  geom_line(linewidth = 0.7, colour = "#d35400") +
  facet_wrap(~ treatment, ncol = 5) +
  geom_hline(yintercept = 36.3, linetype = "dashed", colour = "grey50") +
  labs(x = "Time after dose (h)", y = "Oral body temperature (degC)",
       title = "Study 1 PD -- oral body temperature, single-dose escalation",
       caption = "Dashed line: oral baseline 36.3 degC. Replicates Figure 3 of Othman 2013 (qualitative shape).") +
  theme_minimal()
```

## Peak drug-induced body-temperature increase (Study 1, 40 mg)

Othman 2013 Discussion paragraph 3 quotes: "At the highest evaluated exposure of ABT-102 in humans (40 mg single dose, mean Cmax 73 ng/mL), the model estimated mean increase in body temperature is 1.5 degC at tmax." Verify by computing the typical-value drug-effect term `de(t)` (the third additive component, isolated from baseline + circadian).

```{r typical-de}
mod_typ <- rxode2::zeroRe(mod)

make_typical_de <- function(dose_mg, t_end = 72, FORM_SOLUTION = 1L) {
  obs_rows <- data.frame(
    id = 1L, time = seq(0, t_end, by = 0.1),
    amt = 0, evid = 0L, cmt = "BT_oral",
    FORM_SOLUTION = FORM_SOLUTION
  )
  dose_row <- data.frame(
    id = 1L, time = 0, amt = dose_mg, evid = 1L, cmt = "depot",
    FORM_SOLUTION = FORM_SOLUTION
  )
  rbind(dose_row, obs_rows[order(obs_rows$time), ])
}

de_traj <- bind_rows(
  rxode2::rxSolve(mod_typ, events = make_typical_de( 6),
                  returnType = "data.frame") |>
    mutate(treatment = "6 mg single dose, solution"),
  rxode2::rxSolve(mod_typ, events = make_typical_de(40),
                  returnType = "data.frame") |>
    mutate(treatment = "40 mg single dose, solution")
)
```

```{r typical-de-figure, fig.cap = "Typical-value drug-effect term `de(t)` (degC above baseline) following a single oral-solution ABT-102 dose. The 40 mg curve peaks at ~1.6 degC (close to the paper's quoted ~1.5 degC at tmax); the 6 mg curve peaks at ~0.8 degC."}
ggplot(de_traj, aes(time, de, colour = treatment)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = 1.5, linetype = "dotted", colour = "grey50") +
  scale_color_manual(values = c("6 mg single dose, solution"  = "#3498db",
                                "40 mg single dose, solution" = "#d35400")) +
  labs(x = "Time after dose (h)", y = "Drug-effect term de(t) (degC above baseline)",
       title = "Typical-value ABT-102 drug effect on body temperature",
       caption = "Dotted line at 1.5 degC marks the paper's quoted peak effect at 40 mg single dose.") +
  theme_minimal() +
  theme(legend.position = "top", legend.title = element_blank())
```

```{r peak-de-table}
peak_summary <- de_traj |>
  group_by(treatment) |>
  summarise(peak_de = max(de, na.rm = TRUE),
            peak_time = time[which.max(de)],
            .groups = "drop")

published_peak <- tibble::tribble(
  ~treatment, ~peak_de_paper,
  "40 mg single dose, solution", 1.5,
  "6 mg single dose, solution",  NA_real_   # not directly quoted; "0.6-0.8 degC at analgesic exposures"
)

knitr::kable(
  peak_summary |>
    dplyr::left_join(published_peak, by = "treatment") |>
    dplyr::mutate(pct_diff = 100 * (peak_de - peak_de_paper) / peak_de_paper),
  digits = 2,
  caption = "Simulated peak drug-effect term `de(t)` (typical-value) vs. paper-quoted peak effect of 1.5 degC at the 40 mg single dose."
)
```

# PD validation -- Study 3 tolerance development (core capsule, multi-dose solid dispersion)

Othman 2013 Figures 4 and 5 illustrate that the body-temperature increase induced by ABT-102 dissipates within 2-3 days of twice-daily dosing despite continued accumulation. Simulate Study 3 (4 mg twice-daily solid dispersion for 7 days, core capsule) and confirm that the typical drug-effect term `de(t)` decays over the dosing window.

```{r study3-typical}
ev_study3 <- bind_rows(
  data.frame(id = 1L,
             time = seq(0, 6 * 24, by = 12),  # BID dosing for 7 days
             amt = 4, evid = 1L, cmt = "depot",
             FORM_SOLUTION = 0L),
  data.frame(id = 1L,
             time = seq(0, 9 * 24, by = 0.5),
             amt = 0, evid = 0L, cmt = "BT_core",
             FORM_SOLUTION = 0L)
) |>
  arrange(time, desc(evid))

sim_study3_typ <- rxode2::rxSolve(mod_typ, events = ev_study3,
                                  returnType = "data.frame")
```

```{r study3-figure, fig.cap = "Typical-value Study 3 (4 mg BID solid dispersion for 7 days) trajectories of plasma Cc, drug-effect term de(t), and core body temperature BT_core. The de(t) term peaks after the first dose and decays over 2-3 days of continued dosing -- reproducing the tolerance pattern reported in Othman 2013 Figures 4 and 5."}
sim_long <- sim_study3_typ |>
  dplyr::select(time, Cc, de, BT_core) |>
  tidyr::pivot_longer(c(Cc, de, BT_core),
                      names_to = "variable", values_to = "value") |>
  dplyr::mutate(panel = dplyr::case_when(
    variable == "Cc"      ~ "Plasma Cc (ng/mL)",
    variable == "de"      ~ "Drug effect de(t) (degC)",
    variable == "BT_core" ~ "Core BT (degC)"
  ))

sim_long$panel <- factor(sim_long$panel,
                         levels = c("Plasma Cc (ng/mL)",
                                    "Drug effect de(t) (degC)",
                                    "Core BT (degC)"))

ggplot(sim_long, aes(time, value)) +
  geom_line(linewidth = 0.7, colour = "#16a085") +
  facet_wrap(~ panel, ncol = 1, scales = "free_y") +
  geom_vline(xintercept = seq(0, 6 * 24, by = 12),
             linetype = "dotted", colour = "grey60", alpha = 0.4) +
  scale_x_continuous(breaks = seq(0, 9 * 24, by = 24)) +
  labs(x = "Time (h since first dose)", y = NULL,
       title = "Study 3 typical-value trajectories: 4 mg solid dispersion BID for 7 days",
       caption = "Dotted vertical lines: dosing times. Tolerance to the body-temperature effect develops over the first 2-3 days of dosing.") +
  theme_minimal()
```

```{r tolerance-decay-check}
# Peak de(t) at first dose vs peak in last full dosing interval (day 6).
first_peak <- sim_study3_typ |>
  dplyr::filter(time >= 0, time <= 12) |>
  dplyr::summarise(peak = max(de)) |>
  dplyr::pull(peak)

last_peak <- sim_study3_typ |>
  dplyr::filter(time >= 5 * 24, time <= 6 * 24) |>
  dplyr::summarise(peak = max(de)) |>
  dplyr::pull(peak)

# T50 = 28 h: at t = 5 days = 120 h, the exp(-log(2) * t / T50) factor is exp(-log(2)*120/28) = exp(-2.97) ~ 0.051
# At t = 12 h (end of first dose interval) the same factor is exp(-log(2)*12/28) ~ 0.74
expected_tolerance_factor <- exp(-log(2) * 120 / 28) / exp(-log(2) * 0 / 28)
cat(sprintf("First-interval peak de:  %.3f degC\n", first_peak))
cat(sprintf("Day-6 interval peak de:  %.3f degC\n", last_peak))
cat(sprintf("Observed decay ratio:    %.3f\n", last_peak / first_peak))
cat(sprintf("Theoretical decay factor (T50=28h, t=120h vs t=0h): %.3f\n", expected_tolerance_factor))
```

The first-dose peak in the drug-effect term decays to about 5% of its initial value by Day 6, matching the T50 = 28 h exponential-tolerance kinetics (theoretical decay factor `exp(-log(2) * 120 / 28) = 0.051` from t = 0 h to t = 120 h). This is the quantitative basis for Othman 2013's conclusion that "the effect attenuates within 2 to 3 days of dosing."

# Cross-check -- typical morning-dose circadian peak

The paper notes the circadian-rhythm peak falls at approximately 15:36 in clock time (`t = ps = 7.6` h after morning dosing at 08:00). Verify that the typical-value `cr_oral` and `cr_core` curves reach their maxima at `t = 7.6` h.

```{r circadian-peak}
# Drug-free simulation: dose = 0 so de(t) = 0.
ev_circ <- data.frame(id = 1L, time = seq(0, 24, by = 0.1),
                      amt = 0, evid = 0L, cmt = "BT_oral",
                      FORM_SOLUTION = 1L)
sim_circ <- rxode2::rxSolve(mod_typ, events = ev_circ,
                            returnType = "data.frame")
cat(sprintf("Peak BT_oral occurs at t = %.2f h (paper's ps = 7.6 h)\n",
            sim_circ$time[which.max(sim_circ$BT_oral)]))
cat(sprintf("Peak BT_core occurs at t = %.2f h (paper's ps = 7.6 h)\n",
            sim_circ$time[which.max(sim_circ$BT_core)]))
```

# Assumptions and deviations

- **Upstream PK provenance.** PK structural parameters (CL/F, V/F, ktr, lag times, Frel) and PK IIV (CV% for CL, V, ktr, lag) come from the **upstream** Othman 2012 popPK paper (J Clin Pharmacol 52:1028-1041). The current paper fixes the PK to empirical-Bayes individual estimates from that fit; the values used here are the population estimates re-reported in the PK/PD-model Results paragraph 1 of Othman 2013. The upstream paper was not on disk at extraction time; only the values reproduced in Othman 2013 were available. PK residual error is not reported in either paper; the packaged model exposes plasma Cc as an internal variable for chained simulation but does not declare a PK observation tilde / residual error.

- **Residual-error simplification on body temperature.** Othman 2013 used a combined additive + inverse-proportional residual-error model of the form `Res_kij = (39.5 - BT_ij) * kappa_k * eps_inv + eps_add` with measurement-type-dependent variances (Othman 2013 Discussion paragraph 5; Results equations on page 1034). nlmixr2's standard residual-error syntax (`~ add()`, `~ prop()`) does not directly express a `(constant - prediction)`-dependent variance scale. The packaged model approximates by collapsing both components into a single additive SD per measurement type, evaluated at the cohort-median body temperature of 37 degC (`addSd_BT_oral ~= 0.686 degC`, `addSd_BT_core ~= 0.675 degC`). The simplification is faithful to within ~0.1 degC RMS for BT values near the cohort median; residuals at very low BT (around 34 degC, where the inverse-proportional term dominates) will be under-dispersed and residuals at very high BT (around 39 degC) will be over-dispersed relative to the paper's published model. The paper's reported effective range was 34 to 39 degC.

- **Two parallel outputs (BT_oral, BT_core).** Each subject in the source studies was measured by exactly one device (oral thermometer in Studies 1 and 2; core ingestible capsule in Study 3), and the paper's model includes a fixed effect of measurement type on baseline, amplitude, and residual error. The packaged model expresses both branches as concurrent outputs of the ODE; users select the output that matches their virtual cohort. Drug effect, phase shift, EC50, and T50 are shared between the two output branches.

- **IIV parameterisation mixes log-normal and additive.** Per Othman 2013 (last paragraph before Table 2), inter-subject variability on baselines, amplitudes, and phase shift is **additive on the linear scale** (`P_i = TVP + eta_i`, `eta_i ~ N(0, omega^2)` with `omega^2` in the units of the parameter), while inter-subject variability on EC50 and T50 is **exponential** (log-normal, `P_i = TVP * exp(eta_i)`). The packaged ini() encoding follows this distinction: `etabl`, `etaamp`, `etaps` are paper-symbol etas declared in `paper_specific_etas` and added linearly to their parents in `model()`; `etalec50` and `etalt50` follow the standard nlmixr2lib convention paired with log-transformed parents.

- **Baseline and amplitude IIV is shared across oral and core branches.** Per Othman 2013 Table 1 model-building history (Model 4 and Final model rows), the variance of the baseline and amplitude etas was held common between measurement types. The packaged model expresses this by using the same `etabl` for `bl_oral_i` and `bl_core_i`, and the same `etaamp` for `amp_oral_i` and `amp_core_i`.

- **Subject-independent Emax.** Othman 2013 reports inter-subject variability for EC50 and T50 but states explicitly that ABT-102 Emax was estimated as "subject independent" (no eta). The packaged model carries Emax as a fixed-effect-only parameter.

- **Time origin convention.** `t = 0` corresponds to the first dose, which was administered at approximately 08:00 in all three source studies. The circadian peak at `t = ps = 7.6` h therefore corresponds to a clock-time peak at approximately 15:36 (3:36 PM), matching the paper's Discussion paragraph 4 description. Users who prefer a different zero point can shift `time` in the event table accordingly.

- **Cohort demographics not retabulated in Othman 2013.** The body-temperature PD paper refers readers to Othman 2012 (J Clin Pharmacol 52:1028-1041) for the 108-subject demographics. The packaged `population` metadata records the totals (108 subjects, 3 studies, 7493 measurements after artifact exclusion) and dose-regimen detail; demographics fields (age range, weight range, sex balance, race / ethnicity) are documented as "Adult healthy volunteers (specific range not retabulated in the PD paper)" until the upstream Othman 2012 popPK paper is also packaged.
