---
title: "Phenytoin (Yukawa 1990)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Phenytoin (Yukawa 1990)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Yukawa_1990_phenytoin"))$reference`
* Description: `r rxode2::rxode(readModelDb("Yukawa_1990_phenytoin"))$description`
* Article: <https://doi.org/10.1248/cpb.38.1973>

Yukawa, Higuchi, and Aoyama (Department of Hospital Pharmacy, Kyushu University
Hospital) re-analysed routine therapeutic-drug-monitoring records from 334
Japanese epileptic outpatients to update the Michaelis-Menten population PK
parameters of phenytoin (PHT). The analysis covers 756 paired daily-dose / steady-state
serum concentration records collected during chronic oral PHT maintenance therapy.
Three candidate bioavailability sub-models for the powder formulation were fit
(Methods, page 1974; constant Model 1, exponential Model 2, linear Model 3);
Model 2 was preferred on predictive performance (Table IV: lowest mean absolute
prediction error, MAE = 32.47 mg/d) and is reproduced here.

# Population

The cohort included 334 outpatients (170 male, 164 female; 49.1 percent female)
aged 0.6 to 71.1 years (mean 24.3, SD 14.1) and weighing 9.0 to 115.0 kg
(mean 49.1, SD 15.5). Mean daily PHT dose was 225.8 mg/d (SD 73.1) and mean
steady-state PHT serum concentration was 9.78 ug/mL (SD 7.77; therapeutic
range 10-20 ug/mL). 101 patients were on PHT monotherapy and 233 were on PHT
combined with one or more of phenobarbital, carbamazepine, valproate, primidone,
clonazepam, sultiame, ethotoin, ethosuximide, acetazolamide, or diazepam
(Yukawa 1990 Table I). All patients had normal renal and hepatic function;
concurrent therapy was not altered during the analysis window. Steady-state
serum concentrations were drawn 2-5 hours post-dose at least 30 days after
any dose change. Source: Yukawa 1990 Tables I and II and Methods Data-Sources
(page 1973).

The same information is available programmatically via the model's
`population` metadata.

```{r population-meta}
mod <- readModelDb("Yukawa_1990_phenytoin")
str(rxode2::rxode(mod)$population)
```

# Source trace

The per-parameter origin is recorded as an in-file comment next to each
`ini()` entry in `inst/modeldb/specificDrugs/Yukawa_1990_phenytoin.R`. The
table below collects the same information in one place for review.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lvmax` (Vm, mg/d) at reference covariates | 325 | Table III Model 2 (SEM 9.18) |
| `lkm` (Km, mg/L) at reference covariates | 2.41 | Table III Model 2 (SEM 0.25); 1 ug/mL = 1 mg/L |
| `lvc` (Vc, L; not from this paper) | 36 (fixed) | Phenytoin literature 0.6 L/kg x 60 kg standard; required only for ODE-based simulation |
| `lka` (ka, 1/day; not from this paper) | 36 (fixed; = 1.5 1/h x 24) | Phenytoin tablet literature 1.5 1/h; required only for ODE-based simulation |
| `e_wt_vmax` (power exponent of WT/60 on Vmax) | 0.737 (fixed) | Table III Model 2 (SEM 0.036); Methods Eq. 2 |
| `e_child_km` (multiplicative factor on Km when CHILD = 1) | 0.752 (fixed) | Table III Model 2 (SEM 0.073); Methods Eq. 3 |
| `e_conmed_aed_vmax` (multiplicative factor on Vmax when CONMED_AED = 1) | 1.08 (fixed) | Table III Model 2 (SEM 0.039); Methods Eq. 2 |
| `e_conmed_aed_km` (multiplicative factor on Km when CONMED_AED = 1) | 1.32 (fixed) | Table III Model 2 (SEM 0.165); Methods Eq. 3 |
| `e_form_powder_f` (theta in F_powder = 1 - exp(-theta/D)) | 9.92 (fixed) | Table III Model 2 (SEM 0.52); Methods Eq. 4 |
| `etalvmax` (omega^2, IIV Vmax) | 0.0372 | Table III Model 2: 19.3 percent CV (SEM 2.2); omega^2 = 0.193^2 |
| `etalkm` (omega^2, IIV Km) | 0.398 | Table III Model 2: 63.1 percent CV (SEM 5.1); omega^2 = 0.631^2 |
| `propSd` (proportional residual SD on Cc) | 0.103 | Table III Model 2: theta_E = 10.3 percent (SEM 0.5); paper applies the residual to predicted dose R (Eq. 6), here mapped to Cc |
| ODE `d/dt(central) = ka * depot - vmax * Cc / (km + Cc)` | n/a | Methods Eq. 1 (page 1974) re-cast in time-resolved form |
| `f(depot) = 1 - FORM_POWDER * (1 - 1 + exp(-9.92/DOSE_PHT_MGKGD))` | n/a | Methods Eq. 4 (page 1974), Model 2 branch |

# Virtual cohort

Original observed data are not publicly available. The validation below uses
a grid of typical patient covariate combinations spanning the analysis
covariate space (weight, age stratum, co-medication, formulation, daily dose).
For each cell of the grid we simulate steady-state PHT pharmacokinetics under
a TID dosing regimen and compare the resulting average steady-state
concentration to the algebraic Michaelis-Menten prediction of Yukawa 1990
Methods Eq. 1.

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

# Grid of typical patients. Daily dose is encoded as mg/kg/d so each weight
# stratum receives a clinically appropriate amount and the resulting effective
# input rate R*F stays well below Vmax_ij (which itself scales as
# (WT/60)^0.737). The 3, 4, 5, 6 mg/kg/d levels cover the therapeutic window
# for both paediatric and adult phenytoin maintenance dosing.
grid <- expand.grid(
  WT             = c(15, 30, 60, 80),       # kg; paediatric to adult
  CHILD          = c(0L, 1L),                # 0 = adult, 1 = age < 15 yr
  CONMED_AED     = c(0L, 1L),                # 0 = PHT alone, 1 = on at least one co-AED
  FORM_POWDER    = c(0L, 1L),                # 0 = tablet (F = 1), 1 = powder
  dose_per_kg    = c(2.5, 3.5, 4.5, 5),      # mg/kg/d (kept below saturation)
  stringsAsFactors = FALSE
)
grid$daily_dose_mg <- grid$WT * grid$dose_per_kg
grid$id <- seq_len(nrow(grid))

# Phenytoin is split TID at 0, 8, and 16 hours every 24 hours.
tau   <- 24       # dosing interval (hours of the day)
n_day <- 30       # number of days simulated
sim_end_h <- n_day * 24

# DOSE_PHT_MGKGD is the patient's own total daily dose / WT (mg/kg/d)
grid$DOSE_PHT_MGKGD <- grid$daily_dose_mg / grid$WT

# Build the event table. Time is in DAYS to match the model's units$time = "day".
events <- grid %>%
  rowwise() %>%
  do({
    g <- .
    dose_amt <- g$daily_dose_mg / 3        # mg per TID dose
    dose_times_h <- seq(0, sim_end_h - 1, by = tau / 3)  # 0, 8, 16, 24, ...
    # Observation grid: every 2 hours over the final 24 hours (steady state)
    obs_times_h  <- seq(sim_end_h - 24, sim_end_h, by = 1)
    out <- bind_rows(
      tibble(
        time = dose_times_h / 24, evid = 1L, amt = dose_amt, cmt = "depot",
        Cc = NA_real_
      ),
      tibble(
        time = obs_times_h / 24, evid = 0L, amt = 0, cmt = NA_character_,
        Cc = NA_real_
      )
    )
    out$id           <- g$id
    out$WT           <- g$WT
    out$CHILD        <- g$CHILD
    out$CONMED_AED   <- g$CONMED_AED
    out$FORM_POWDER  <- g$FORM_POWDER
    out$DOSE_PHT_MGKGD <- g$DOSE_PHT_MGKGD
    out$daily_dose_mg <- g$daily_dose_mg
    out
  }) %>%
  ungroup() %>%
  arrange(id, time, desc(evid))

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

# Simulation

Simulate typical-value PK (between-subject random effects zeroed) so the
steady-state concentrations match the model's typical-value prediction
exactly and can be compared directly to the algebraic Michaelis-Menten
solution.

```{r simulate}
mod <- readModelDb("Yukawa_1990_phenytoin")
mod_typical <- rxode2::zeroRe(mod)

sim <- rxode2::rxSolve(
  mod_typical,
  events = as.data.frame(events),
  keep   = c("WT", "CHILD", "CONMED_AED", "FORM_POWDER",
             "DOSE_PHT_MGKGD", "daily_dose_mg")
) |>
  as.data.frame()

head(sim)
```

# Replicate published figures

## Figure 2 -- relative bioavailability of phenytoin powder vs daily dose

Yukawa 1990 Figure 2 shows the powder relative bioavailability F as a function
of daily dose mg/kg/d for the three competing models. Model 2 (the preferred
fit) gives `F = 1 - exp(-9.92 / D)`. Replicating the curve algebraically:

```{r figure-2-powder-F}
fig2 <- tibble(
  D_mgkgd = seq(0.1, 12, by = 0.05),
  F_model1 = 0.840,
  F_model2 = 1 - exp(-9.92 / D_mgkgd),
  F_model3 = pmax(0, 1 - 0.0325 * D_mgkgd)
) |>
  pivot_longer(starts_with("F_"), names_to = "model", values_to = "F")

ggplot(fig2, aes(D_mgkgd, F, colour = model)) +
  geom_line(linewidth = 0.8) +
  scale_colour_manual(
    values = c(F_model1 = "grey50", F_model2 = "black", F_model3 = "steelblue"),
    labels = c(F_model1 = "Model 1 (constant 0.840)",
               F_model2 = "Model 2 (1 - exp(-9.92/D))",
               F_model3 = "Model 3 (1 - 0.0325 * D)")
  ) +
  geom_hline(yintercept = 1, linetype = "dotted") +
  labs(x = "Daily dose D (mg/kg/d)", y = "Relative bioavailability F (powder)",
       colour = "Model",
       title = "Figure 2 - Powder relative bioavailability vs daily dose",
       caption = "Replicates Figure 2 of Yukawa 1990. Model 2 is the preferred fit.")
```

The Model 2 curve is implemented inside the packaged model via
`fdepot = (1 - FORM_POWDER) + FORM_POWDER * (1 - exp(-e_form_powder_f / DOSE_PHT_MGKGD))`,
with `e_form_powder_f = 9.92` taken directly from Table III. At low daily
doses (< 2 mg/kg/d) F approaches 1, matching the Discussion claim that
"the percent of relative bioavailability of PHT powder was nearly 100 percent
in less than 2 mg/kg of daily dose" (page 1976).

# PKNCA validation

Compute steady-state NCA over the final 24-h dosing interval. Because the
model is steady-state by design (the paper fits a regression of daily dose
against steady-state concentration; see Methods Eq. 1), the natural NCA
target is the average concentration over a dosing interval at steady state
(`cav.tau`), which the algebraic Michaelis-Menten relation predicts directly:
at steady state `R * F = Vmax * Cav / (Km + Cav)`, so
`Cav = R * F * Km / (Vmax - R * F)`.

```{r pknca}
# rxSolve output rows are observation rows by default (one per requested time);
# restrict further to the steady-state final-day window so PKNCA computes
# steady-state metrics over a single dosing interval.
sim_nca <- sim |>
  filter(!is.na(Cc), time >= 29, time <= 30) |>
  mutate(scenario = paste0(
    "WT=", WT, "kg ",
    ifelse(CHILD == 1, "child ", "adult "),
    ifelse(CONMED_AED == 1, "+coAED ", "alone "),
    ifelse(FORM_POWDER == 1, "powder ", "tablet "),
    daily_dose_mg, "mg/d")) |>
  select(id, time, Cc, scenario) |>
  distinct(id, time, .keep_all = TRUE)

dose_df <- events |>
  filter(evid == 1) |>
  mutate(scenario = paste0(
    "WT=", WT, "kg ",
    ifelse(CHILD == 1, "child ", "adult "),
    ifelse(CONMED_AED == 1, "+coAED ", "alone "),
    ifelse(FORM_POWDER == 1, "powder ", "tablet "),
    daily_dose_mg, "mg/d")) |>
  select(id, time, amt, scenario)

conc_obj <- PKNCA::PKNCAconc(
  as.data.frame(sim_nca), Cc ~ time | scenario + id,
  concu = "mg/L", timeu = "day"
)
dose_obj <- PKNCA::PKNCAdose(
  as.data.frame(dose_df), amt ~ time | scenario + id,
  doseu = "mg"
)

# Steady-state interval: the final 24 hours = 1 day of the simulation.
ss_start <- 30 - 1
ss_end   <- 30
intervals <- data.frame(
  start    = ss_start,
  end      = ss_end,
  cmax     = TRUE,
  tmax     = TRUE,
  cmin     = TRUE,
  cav      = TRUE,
  auclast  = TRUE
)

nca_data <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- suppressWarnings(PKNCA::pk.nca(nca_data))
nca_tbl  <- as.data.frame(nca_res$result)
head(nca_tbl[, c("scenario", "id", "PPTESTCD", "PPORRES")], 12)
```

## Comparison against the published algebraic relationship

The Yukawa 1990 model is a steady-state population-PK model: Methods Eq. 1
states `R * F = Vmax * Cpss / (Km + Cpss)`, with all individual covariate
adjustments (allometric WT on Vmax; co-anticonvulsants on Vmax and Km;
pediatric-status factor on Km; powder-formulation dose-dependent F) baked
into Vmax_ij, Km_ij, and F_ij. Solving for steady-state concentration:

`Cpss = R * F * Km / (Vmax - R * F)`

For each grid cell we compute the algebraic Cpss directly and compare against
the simulated steady-state average concentration extracted by PKNCA.

```{r compare-algebraic}
# Algebraic prediction per Methods Eq. 1 (Model 2)
algebraic <- grid |>
  mutate(
    vmax_ij  = 325 * (WT / 60)^0.737 * 1.08^CONMED_AED,
    km_ij    = 2.41 * 0.752^CHILD * 1.32^CONMED_AED,
    F_powder = 1 - exp(-9.92 / DOSE_PHT_MGKGD),
    F_ij     = (1 - FORM_POWDER) + FORM_POWDER * F_powder,
    R_eff    = daily_dose_mg * F_ij,
    Cpss_algebraic = ifelse(R_eff < vmax_ij,
                            R_eff * km_ij / (vmax_ij - R_eff),
                            NA_real_)
  )

cav_tbl <- nca_tbl |>
  filter(PPTESTCD == "cav") |>
  group_by(scenario, id) |>
  summarise(Cav_sim = first(PPORRES), .groups = "drop")

comparison <- algebraic |>
  left_join(cav_tbl, by = "id") |>
  mutate(rel_diff_pct = 100 * (Cav_sim - Cpss_algebraic) / Cpss_algebraic)

# Show a representative subset (60 kg adult patients)
comparison |>
  filter(WT == 60, FORM_POWDER == 0) |>
  select(WT, CHILD, CONMED_AED, daily_dose_mg, Cpss_algebraic, Cav_sim,
         rel_diff_pct) |>
  knitr::kable(digits = 3,
               caption = "60 kg adult tablet patients - algebraic vs simulated steady-state Cav. CHILD = 0 throughout (60 kg implies adult).")
```

```{r compare-summary}
cat(sprintf(
  "Median |relative difference|: %.2f%%\n",
  median(abs(comparison$rel_diff_pct), na.rm = TRUE)
))
cat(sprintf(
  "95th percentile |relative difference|: %.2f%%\n",
  quantile(abs(comparison$rel_diff_pct), 0.95, na.rm = TRUE)
))
```

For typical-value (zero random-effects) simulations, the algebraic and
ODE-simulated steady-state averages agree to within a few percent across the
covariate grid. The remaining residual reflects (a) the finite number of
TID doses simulated before the steady-state interval, and (b) the temporal
average over the dosing interval being computed by PKNCA's trapezoidal rule
on a discrete sampling grid rather than analytically. Both are simulation
artefacts rather than model-structural disagreements with Yukawa 1990 Eq. 1.

## Steady-state Cc time course over the final dosing interval

```{r ss-time-course}
sim |>
  filter(time >= 29, time <= 30,
         WT == 60, FORM_POWDER == 0, CHILD == 0) |>
  mutate(scenario = paste0(
    ifelse(CONMED_AED == 1, "+coAED ", "alone "),
    daily_dose_mg, " mg/d"),
    time_h = (time - 29) * 24) |>
  ggplot(aes(time_h, Cc, colour = scenario)) +
  geom_line(linewidth = 0.7) +
  geom_hline(yintercept = 10, linetype = "dotted", colour = "grey40") +
  geom_hline(yintercept = 20, linetype = "dotted", colour = "grey40") +
  labs(x = "Time within dosing interval (h)",
       y = "PHT serum concentration (mg/L)",
       colour = "Scenario",
       title = "Steady-state Cc over the final 24 h",
       caption = "60 kg adult tablet patients on TID dosing. Dotted lines = therapeutic range 10-20 ug/mL.")
```

# Assumptions and deviations

- **Vc and ka are not from Yukawa 1990.** The paper's published model is a
  steady-state regression of daily dose against steady-state concentration
  (Methods Eq. 1) and does not require either parameter; both are needed
  only as ODE structural constants for time-resolved simulation. The model
  file fixes `Vc = 36 L` (= 0.6 L/kg x 60 kg standard, a widely cited
  phenytoin total-drug Vd from product-label and pharmacology references)
  and `ka = 1.5 1/h = 36 1/day` (a typical phenytoin tablet absorption
  rate from the literature, matching the Hennig 2015 phenytoin model file
  in this registry). Steady-state Cpss is invariant to either choice; both
  values affect only the time taken to reach steady state and the
  intra-interval Cmax/Cmin oscillation around the steady-state mean.
- **Residual error mapping from R-space to Cc-space.** Yukawa 1990 Methods
  Eq. 6 specifies a 10.3 percent CV proportional residual on the predicted
  daily dose R, not on Cpss. Forward simulation in nlmixr2 carries the
  residual on the observation Cc. At steady state R is monotonic in Cpss
  but the elasticity Km / (Km + Cpss) means a 10.3 percent CV on R
  corresponds to a larger CV on Cpss in the saturated regime; the model
  file applies `propSd = 0.103` directly to Cc as a re-parameterization,
  preserving the magnitude reported in the source while shifting the noise
  axis. This is a structural simplification of the published statistical
  model.
- **IIV variance from CV percent.** Yukawa 1990 Table III reports CV
  percent for Vm and Km IIV. The model file uses the squared-fractional-CV
  approximation `omega^2 ~ (CV/100)^2` (giving 0.0372 for 19.3 percent CV
  and 0.398 for 63.1 percent CV) consistent with the typical NONMEM
  exponential-IIV-on-positive-parameter encoding and with the registry's
  Hennig 2015 phenytoin convention.
- **CONMED_AED orientation flipped from source.** The paper's `CO` indicator
  is 1 when the patient is on PHT alone (Methods Eqs. 2-3); the canonical
  `CONMED_AED` indicator inverts this so 0 is the monotherapy reference,
  matching the rest of the `CONMED_*` family in the registry. The same
  inversion applies to `CHILD` (paper's `AGE` indicator = 1 for adults;
  canonical `CHILD = 1 - AGE_indicator`) and `FORM_POWDER` (paper's `BA`
  indicator = 1 for tablets; canonical `FORM_POWDER = 1 - BA_indicator`).
  The numerical coefficients (theta_AGE = 0.752, theta_coVm = 1.08,
  theta_coKm = 1.32, theta_BA2 = 9.92) are unchanged because the
  multiplicative-factor formulation absorbs the orientation flip cleanly.
- **`DOSE_PHT_MGKGD` is a per-record regressor.** The Yukawa 1990 powder
  bioavailability formula F_powder = 1 - exp(-9.92 / D) requires the
  patient's own current daily PHT dose normalised by body weight. Users
  must supply this column on every dose record (compute as total daily
  dose mg/d divided by current WT in kg). For tablet records
  (FORM_POWDER = 0) the value is multiplied by 0 in the F expression and
  has no effect, but a non-NA placeholder is still required so the
  expression evaluates.
