---
title: "Triflusal (Park 2014)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Triflusal (Park 2014)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r readModelDb("Park_2014_triflusal")()$meta$reference`
* Description: `r readModelDb("Park_2014_triflusal")()$meta$description`
* Article: <https://doi.org/10.1186/2050-6511-15-75>

# Population

Park 2014 enrolled 38 healthy Korean adult male volunteers in a randomized,
open-label, multiple-dose, two-period, two-treatment, comparative crossover
study at Kyungpook National University Hospital Clinical Trial Center
(KNUH CTC) between August and September 2008. Four subjects dropped out;
data from the remaining 34 subjects (mean age 24.1 +/- 1.7 years, mean
height 176.1 +/- 4.9 cm, mean weight 70.8 +/- 9.0 kg) were used to construct
the population PK and PD models. Volunteers were required to be within 20%
of their ideal body weight per the Devine-style formula
`IBW (kg) = (height_cm - 100) * 0.9` and to pass routine physical and
laboratory screening. Baseline demographics are reported in
Park 2014 Table 1.

The dosing regimen during each period was a single 900 mg oral loading
dose of triflusal on Day 1 followed by a 600 mg/day maintenance dose
(two 300 mg capsules once daily) from Day 2 through Day 9. Plasma HTB was
sampled at 0, 24, 48, 96, 144, 168, 192, 192.5, 193, 194, 196, 199, 202,
and 216 h after the first dose; platelet aggregation was assessed at the
same times excluding the dense post-Day-9 sub-sampling. Population PK was
fit to 476 HTB plasma concentrations and population PD to 340 platelet-
aggregation observations.

The structured `population` metadata is available programmatically via
`readModelDb("Park_2014_triflusal")()$meta$population`.

# Source trace

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| Structural PK form (`d/dt(depot) = -ka * depot`; `d/dt(central) = ka * depot - kel * central`) | n/a | Methods: Population PK modeling -- "ADVAN2 TRANS2" (1-cmt with first-order absorption); the paper's `k_f` plays the role of the canonical `ka` |
| `kf` (= `lka` canonical) | 0.341 1/h | Park 2014 Table 2 theta_3 (RSE 15.1%) |
| `CL/F` (= `lcl` canonical) | 0.200 L/h at 71.65 kg | Park 2014 Table 2 theta_1 (RSE 2.7%) |
| `V/F` (= `lvc` canonical) | 8.300 L at 71.65 kg | Park 2014 Table 2 theta_2 (RSE 2.7%) |
| `(WT/71.65)^0.845` weight effect on CL/F | exponent 0.845 (estimated) | Park 2014 Table 2 theta_4 (RSE 17.4%) |
| `(WT/71.65)^1` weight effect on V/F | exponent fixed at 1 | Park 2014 Table 2 V/F equation `V_d/F = theta_2 * (weight/71.65)` (no separate exponent reported) |
| `EC50` (= `lec50` canonical) | 84.9 ug/mL | Park 2014 Table 2 theta_5 (RSE 4.0%) |
| Hill exponent `gamma` (= `lhill` canonical) | 19.2 | Park 2014 Table 2 theta_6 (RSE 22.4%) |
| `P(IPA) = Cc^gamma / (EC50^gamma + Cc^gamma)` | n/a | Methods: Population PK-PD model development, binary probability model equation |
| BSV on CL/F | 14.9% CV | Park 2014 Table 2 omega_1^2 (RSE 21.0%; shrinkage 0.6%) |
| BSV on V/F  |  9.5% CV | Park 2014 Table 2 omega_2^2 (RSE 57.8%; shrinkage 40.2%) |
| BSV on k_f  | 88.0% CV | Park 2014 Table 2 omega_3^2 (RSE 26.5%; shrinkage 10.5%) |
| BSV on EC50 | 21.8% CV | Park 2014 Table 2 omega_4^2 (RSE 28.4%; shrinkage 8.3%) |
| BSV on gamma | 0 (fixed) | Park 2014 Table 2 omega_5^2 (NE) |
| Proportional residual error | sigma^2 = 0.098 | Park 2014 Table 2 sigma^2 (RSE 6.5%; shrinkage 11.5%); propSd = sqrt(0.098) |

# Virtual cohort

The published trial data are not available. We assemble a virtual cohort
of 200 healthy male volunteers with weights drawn from a truncated normal
distribution that matches the reported 70.8 +/- 9.0 kg (53.3 - 89.7 kg
range) of Park 2014 Table 1.

```{r cohort}
set.seed(20140523L)

n_subjects <- 200L

cohort <- tibble::tibble(
  id = seq_len(n_subjects),
  WT = pmin(pmax(rnorm(n_subjects, mean = 70.8, sd = 9.0), 53.3), 89.7)
)

# Dosing: 900 mg loading at t=0, then 600 mg/day from Day 2 to Day 9.
dose_times <- c(0, seq(from = 24, to = 192, by = 24))
dose_amts  <- c(900, rep(600, length(dose_times) - 1L))

dose_events <- expand.grid(
  id   = cohort$id,
  idx  = seq_along(dose_times)
) |>
  dplyr::mutate(
    time = dose_times[idx],
    amt  = dose_amts[idx],
    evid = 1L,
    cmt  = "depot"
  ) |>
  dplyr::select(id, time, amt, evid, cmt) |>
  dplyr::left_join(cohort, by = "id")

# Observation times: a dense grid + the paper's specific PK sampling points
obs_times <- sort(unique(c(
  seq(0, 216, by = 0.5),
  c(0, 24, 48, 96, 144, 168, 192, 192.5, 193, 194, 196, 199, 202, 216)
)))

obs_events <- expand.grid(
  id   = cohort$id,
  time = obs_times
) |>
  dplyr::mutate(
    amt  = NA_real_,
    evid = 0L,
    cmt  = NA_character_
  ) |>
  dplyr::left_join(cohort, by = "id")

events <- dplyr::bind_rows(dose_events, obs_events) |>
  dplyr::arrange(id, time, dplyr::desc(evid))

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- readModelDb("Park_2014_triflusal")()

sim <- rxode2::rxSolve(mod, events = events, keep = c("WT")) |>
  as.data.frame()
```

For deterministic typical-value replication (e.g., reproducing the paper's
Css,min and accumulation factor exactly), the random effects are zeroed:

```{r simulate-typical}
mod_typical <- mod |> rxode2::zeroRe()

events_typical <- events |> dplyr::filter(id == 1L)
events_typical$WT <- 71.65

sim_typical <- rxode2::rxSolve(mod_typical, events = events_typical) |>
  as.data.frame()
```

# Replicate published figures

## Figure 2(a) -- HTB concentration over time

Park 2014 Figure 2(a) shows individual and median HTB plasma concentration
profiles across the multi-dose regimen. The simulated VPC below uses the
200-subject virtual cohort.

```{r figure-2a}
vpc_df <- sim |>
  dplyr::filter(time %in% c(0, 24, 48, 96, 144, 168, 192, 192.5, 193, 194,
                            196, 199, 202, 216)) |>
  dplyr::group_by(time) |>
  dplyr::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"
  )

ggplot(vpc_df, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "steelblue") +
  geom_line(colour = "steelblue", linewidth = 0.8) +
  geom_point(colour = "steelblue") +
  labs(
    x = "Time after first dose (h)",
    y = "HTB plasma concentration (ug/mL)",
    title = "Figure 2(a) -- simulated HTB concentration over time",
    subtitle = "Median and 5th/95th percentile across 200 virtual subjects",
    caption = "Replicates Figure 2(a) of Park 2014. Sampling times match the paper."
  ) +
  theme_minimal()
```

## Figure 5(b) -- predicted concentration vs probability of IPA

Park 2014 Figure 5(b) plots the predicted HTB concentration vs probability
of IPA curve, with circles marking the model-predicted IPA probabilities
at each individual PD sampling timepoint. The Hill exponent gamma = 19.2
is steep (quantal-like response).

```{r figure-5b}
conc_grid <- seq(0, 200, by = 0.5)
ec50 <- 84.9
hill <- 19.2
prob_curve <- conc_grid^hill / (ec50^hill + conc_grid^hill)

ggplot(data.frame(Cc = conc_grid, prob = prob_curve),
       aes(Cc, prob)) +
  geom_line(linewidth = 0.8, colour = "firebrick") +
  geom_vline(xintercept = ec50, linetype = "dashed", colour = "grey50") +
  geom_hline(yintercept = 0.5, linetype = "dashed", colour = "grey50") +
  annotate("text", x = ec50 + 2, y = 0.05,
           label = paste0("EC50 = ", ec50, " ug/mL"), hjust = 0) +
  labs(
    x = "HTB plasma concentration (ug/mL)",
    y = "Probability of IPA",
    title = "Figure 5(b) -- predicted concentration vs probability of IPA",
    caption = "Replicates Figure 5(b) of Park 2014."
  ) +
  theme_minimal()
```

# PKNCA validation

The paper reports two derived steady-state quantities under "Population
PK modeling" / "Discussion": Css,min = 103.5 ug/mL and accumulation
factor = 2.26 ("estimated using the values of final population PK
parameters"). PKNCA is used here on the typical-value simulation across
the steady-state interval (the Day-8 to Day-9 window, the last full
dosing interval before the dense intensive sampling) to extract Cmin,
Cmax, AUC0-tau, and Cavg for comparison.

```{r pknca}
tau <- 24

# Use the typical-value simulation for the published-value comparison.
# rxSolve may drop the `id` column for single-subject simulations; add it
# back explicitly before passing to PKNCA.
sim_nca <- sim_typical |>
  dplyr::filter(!is.na(Cc), time >= 168, time <= 192) |>
  dplyr::mutate(id = 1L, treatment = "typical") |>
  dplyr::select(id, time, Cc, treatment)

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

dose_df <- dose_events |>
  dplyr::filter(id == 1L) |>
  dplyr::select(id, time, amt) |>
  dplyr::mutate(treatment = "typical")

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

intervals <- data.frame(
  start    = 168,
  end      = 192,
  cmin     = TRUE,
  cmax     = TRUE,
  tmax     = TRUE,
  auclast  = TRUE,
  cav      = TRUE
)

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

knitr::kable(
  as.data.frame(nca_res$result),
  digits = 3,
  caption = paste0(
    "PKNCA NCA over the Day-8 to Day-9 steady-state interval ",
    "(168 - 192 h post-first-dose) on the typical-value simulation."
  )
)
```

## Comparison against published Css,min and accumulation factor

```{r compare}
sim_cmin_ss <- min(sim_typical$Cc[sim_typical$time >= 168 &
                                  sim_typical$time <= 192])
# Accumulation factor (AUC basis) per closed-form for 1-cmt first-order:
# F_acc = 1 / (1 - exp(-kel * tau))
kel_typ <- 0.200 / 8.300
acc_factor_closed <- 1 / (1 - exp(-kel_typ * tau))

dplyr::tibble(
  Quantity = c("Css,min (ug/mL)", "AUC accumulation factor"),
  Published = c(103.5, 2.26),
  Simulated = round(c(sim_cmin_ss, acc_factor_closed), 2)
) |>
  knitr::kable(caption = "Comparison against Park 2014 published values.")
```

# Binary probability PD -- post-processing

The PD endpoint in Park 2014 is the binary indicator IPA (1 = inhibition,
0 = non-inhibition), dichotomised at platelet aggregation < 74%. The model
exposes `prob_ipa` as a deterministic transform of `Cc`; binary draws
are generated downstream from `prob_ipa`. The chunk below reproduces
Park 2014 Table 3 (observed vs predicted IPA counts at the PD sampling
times) by averaging the model-predicted `prob_ipa` across the virtual
cohort.

```{r table-3}
table3_times <- c(0, 24, 48, 96, 144, 168, 192, 196, 202, 216)
n_sim <- length(unique(sim$id))

table3 <- sim |>
  dplyr::filter(time %in% table3_times) |>
  dplyr::group_by(time) |>
  dplyr::summarise(
    mean_prob_ipa  = mean(prob_ipa, na.rm = TRUE),
    expected_count = round(mean_prob_ipa * n_sim),
    expected_pct   = round(100 * mean_prob_ipa),
    .groups        = "drop"
  )

park2014_table3 <- dplyr::tibble(
  time = table3_times,
  paper_observed_pct  = c(3, 6, 32, 85, 74, 79, 88, 100, 100, 82),
  paper_predicted_pct = c(0, 6, 32, 74, 82, 88, 88, 100, 100, 88)
)

dplyr::inner_join(table3, park2014_table3, by = "time") |>
  dplyr::transmute(
    `Time (h)` = time,
    `Mean prob_ipa` = round(mean_prob_ipa, 3),
    `Simulated expected IPA (%)` = expected_pct,
    `Paper predicted IPA (%)`    = paper_predicted_pct,
    `Paper observed IPA (%)`     = paper_observed_pct
  ) |>
  knitr::kable(caption = paste0(
    "Park 2014 Table 3 -- observed vs predicted IPA counts. ",
    "Simulated expected IPA (%) is 100 * mean(prob_ipa) over the ",
    n_sim, "-subject virtual cohort."
  ))
```

For a Bernoulli realization of the binary outcome:

```{r bernoulli, eval = FALSE}
# Per-subject binary draws on the PD sampling grid
ipa_draws <- sim |>
  dplyr::filter(time %in% table3_times) |>
  dplyr::mutate(ipa = rbinom(dplyr::n(), size = 1, prob = prob_ipa))
```

# Assumptions and deviations

- The PD endpoint is binary (IPA = 1 if platelet aggregation < 74% else 0).
  The packaged model exposes the deterministic probability `prob_ipa` as
  a derived quantity in `model()`; the model file declares only `Cc` as
  the formal observation variable (with proportional residual error per
  Park 2014 Table 2). Binary realizations are drawn downstream via
  `rbinom(n, 1, prob_ipa)` -- see the chunk above. This decision keeps
  the model file aligned with the nlmixr2 convention that observation
  variables carry standard residual-error families (`add`, `prop`,
  `lnorm`); the paper's binary likelihood is encoded as a documented
  post-processing step rather than as an in-model `dbinom(prob, 1)` call.
- The cohort body-weight distribution is sampled from a truncated normal
  approximation of the reported mean +/- SD because the trial-level
  individual demographics are not publicly available.
- Race / ethnicity is set to 100% Asian to match the reported Korean
  cohort; no explicit race-effect covariate is in the model.
- The V/F covariate equation `V/F = theta_2 * (WT/71.65)` is encoded with
  `e_wt_vc <- fixed(1)`. The paper reports the equation without a
  separately estimated exponent, consistent with the authors' Discussion
  statement that "V/F was directly proportional to weight (exponent = 1)".
- The Figure 1 caption lists PK sampling days (2, 3, 5, 7, 8, 9) but the
  Methods text "from Day 2 to 9" is the authoritative dosing schedule and
  is consistent with the paper's reported steady-state Css,min = 103.5
  ug/mL and AUC accumulation factor 2.26. The vignette and the model
  metadata reflect daily dosing across Day 2 - Day 9 inclusive.
- The paper's `kf` (first-order HTB formation rate constant) is encoded
  as the canonical `lka`. The math is identical to a 1-compartment
  first-order absorption model (NONMEM ADVAN2 TRANS2); the depot
  compartment carries triflusal and converts to HTB in central at rate
  `ka == kf`. F (bioavailability) absorbs the unknown fraction of
  triflusal converted to HTB, so CL/F and V/F are apparent HTB-disposition
  parameters.
