---
title: "Aripiprazole (Koue 2007)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Aripiprazole (Koue 2007)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode2(readModelDb("Koue_2007_aripiprazole"))$reference`
* Description: `r rxode2::rxode2(readModelDb("Koue_2007_aripiprazole"))$description`
* Article: <https://doi.org/10.1248/bpb.30.2154>

# Population

Koue 2007 fit the popPK model to 68 healthy Japanese adult male volunteers
(age 20-32 years, mean 23.1; weight 48.2-75.2 kg, mean 60.7) enrolled across
three clinical-pharmacology trials (Pharmacokinetic Data section): a
single-dose trial (protocol 1; 6 mg ABILIFY tablet under fasting), a
multiple-dose trial (protocol 2; 3 mg once daily for 14 d, dosed 30 min
after breakfast), and an itraconazole-coadministration crossover trial
(protocol 3; 3 mg single dose alone, then re-dosed at 3 mg 7 days into a
21-day itraconazole 100 mg once-daily regimen, with a 5-week washout between
phases). Plasma aripiprazole was assayed by LC-MS/MS. The abstract reports
24 subjects in the single-dose trial; the Pharmacokinetic Data section
narrows that to 26 subjects, of whom 21 contributed analysable
concentration-time profiles to Figure 1.

CYP2D6 genotype was determined by PCR-RFLP and long-PCR for `*2`, `*4`,
`*5`, `*10`, `*14`, `*18`, `*36`, and `*41`; subjects were grouped into
three CYP2D6 strata (Table 1). Group 1 (extensive metabolizers, EM; n=14)
carried `*1/*1`, `*1/*2`, or `*2/*2`. Group 2 (intermediate metabolizers,
IM; n=39) carried `*1/*5`, `*1/*10`, `*2/*5`, `*2/*10`, or `*2/*41`. Group 3
(poor metabolizers, PM; n=15) carried `*5/*10`, `*10/*10`, or `*41/*41`.
CYP3A5 was not genotyped. Group 1 (EM) is the explicit reference category
in Eq. 4.

The same information is available programmatically via
`readModelDb("Koue_2007_aripiprazole")$population`.

# Source trace

Every parameter in the model file carries an inline source-location
comment. The table below collects the entries in one place. Koue 2007
Equations 1-7 specify the model structure explicitly; Table 2 reports the
final parameter point estimates and 95% confidence intervals.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `ltlag` (ALAG) | 0.805 h | Table 2 theta1 (Eq. 1) |
| `lka` (KA) | 2.65 1/h | Table 2 theta2 (Eq. 2) |
| `lvc` (V1/F per kg) | 3.84 L/kg | Table 2 theta3 (Eq. 3) |
| `lcl` (CL/F per kg, Group 1) | 0.0645 L/h/kg | Table 2 theta4 (Eq. 4) |
| `lvp` (V2/F per kg) | 1.54 L/kg | Table 2 theta8 (Eq. 5) |
| `lq` (Q/F per kg) | 0.168 L/h/kg | Table 2 theta9 (Eq. 6) |
| `e_2d6im_cl` (Group-2 additive shift on CL/F per kg) | -0.0135 L/h/kg | Table 2 theta5 (Eq. 4); paper reports as a positive decrement |
| `e_2d6pm_cl` (Group-3 additive shift on CL/F per kg) | -0.0293 L/h/kg | Table 2 theta6 (Eq. 4); paper reports as a positive decrement |
| `e_azole_cl` (Itraconazole additive shift on CL/F per kg) | -0.0181 L/h/kg | Table 2 theta7 (Eq. 4); paper reports as a positive decrement |
| IIV omega_ALAG (SD, log scale) | 0.110 | Table 2 omega_ALAG; variance = 0.110^2 = 0.0121 |
| IIV omega_KA (SD, log scale) | 0.910 | Table 2 omega_KA; variance = 0.910^2 = 0.8281 |
| IIV omega_V1/F (SD, log scale) | 0.381 | Table 2 omega_V1/F; variance = 0.381^2 = 0.1452 |
| IIV omega_CL/F (SD, log scale) | 0.397 | Table 2 omega_CL/F; variance = 0.397^2 = 0.1576 |
| IIV omega_V2/F (SD, log scale) | 0.365 | Table 2 omega_V2/F; variance = 0.365^2 = 0.1332 |
| IIV omega_Q/F (SD, log scale) | 0.262 | Table 2 omega_Q/F; variance = 0.262^2 = 0.0686 |
| Residual sigma (SD on log scale) | 0.166 | Table 2 sigma; Eq. 7 lognormal residual |
| CL/F covariate equation `CL/F = (theta4 - theta5*G2 - theta6*G3 - theta7*ITZ) * WT * exp(eta_CL/F)` | -- | Eq. 4; G2 = IM indicator, G3 = PM indicator, ITZ = azole coadministration indicator |
| V1/F covariate equation `V1/F = theta3 * WT * exp(eta_V1/F)` | -- | Eq. 3 |
| V2/F covariate equation `V2/F = theta8 * WT * exp(eta_V2/F)` | -- | Eq. 5 |
| Q/F covariate equation `Q/F = theta9 * WT * exp(eta_Q/F)` | -- | Eq. 6 |
| ALAG covariate equation `ALAG = theta1 * exp(eta_ALAG)` | -- | Eq. 1 |
| KA covariate equation `KA = theta2 * exp(eta_KA)` | -- | Eq. 2 |
| 2-cmt structure with first-order absorption and absorption lag | -- | Methods, NONMEM Analysis paragraph (ADVAN4 TRANS4) |
| Lognormal (exponential) residual error `C_ij = C*_ij * exp(eps_ij)`, eps ~ N(0, sigma^2) | -- | Eq. 7 |
| Concentration scaling `Cc (ng/mL) = central (mg) / Vc (L) * 1000` | -- | Figure y-axes; LC-MS/MS quantification in ng/mL |

# Virtual cohort

The Koue 2007 raw observations are not openly available. The cohort below
samples covariates compatible with the model's reference subject
(60.7 kg = population mean) and exposes each combination of the three
CYP2D6 phenotype groups and the two itraconazole-coadministration states
relevant to the paper's three dosing protocols.

```{r cohort}
set.seed(20071105)  # Koue 2007 publication month

scenarios <- tibble::tibble(
  scenario = c(
    "Group 1 EM (no ITZ)",
    "Group 2 IM (no ITZ)",
    "Group 3 PM (no ITZ)",
    "Group 1 EM + ITZ",
    "Group 2 IM + ITZ",
    "Group 3 PM + ITZ"
  ),
  WT           = 60.7,
  CYP2D6_EM    = c(1L, 0L, 0L, 1L, 0L, 0L),
  CYP2D6_PM    = c(0L, 0L, 1L, 0L, 0L, 1L),
  CONMED_AZOLE = c(0L, 0L, 0L, 1L, 1L, 1L)
) |>
  dplyr::mutate(id = dplyr::row_number(), .before = scenario)
```

# Simulation

Three event-table builders match the paper's three protocols: a single 6 mg
oral dose with sampling out to 168 h (protocol 1, Figure 1), a 3 mg once-daily
regimen for 14 d with sampling out to 600 h (protocol 2, Figure 2), and a
3 mg single dose with sampling out to 366 h with or without itraconazole
coadministration (protocol 3, Figure 3).

```{r build-events}
build_single_dose <- function(demo, dose_mg, tlast_h, obs_grid = NULL) {
  if (is.null(obs_grid)) {
    obs_grid <- sort(unique(c(seq(0, 24, by = 0.25),
                               seq(24, tlast_h, by = 1))))
  }
  doses <- demo |>
    dplyr::mutate(amt = dose_mg, evid = 1L, cmt = "depot", time = 0,
                  ii = NA_real_, addl = NA_integer_) |>
    dplyr::select(id, time, amt, evid, cmt, ii, addl, scenario,
                  WT, CYP2D6_EM, CYP2D6_PM, CONMED_AZOLE)
  obs <- demo |>
    dplyr::select(id, scenario, WT, CYP2D6_EM, CYP2D6_PM, CONMED_AZOLE) |>
    tidyr::crossing(time = obs_grid) |>
    dplyr::mutate(amt = NA_real_, evid = 0L, cmt = NA_character_,
                  ii = NA_real_, addl = NA_integer_)
  dplyr::bind_rows(doses, obs) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

build_multi_dose <- function(demo, dose_mg, n_doses, tau_h, tlast_h) {
  doses <- demo |>
    dplyr::mutate(amt = dose_mg, evid = 1L, cmt = "depot",
                  time = 0, ii = tau_h, addl = n_doses - 1L) |>
    dplyr::select(id, time, amt, evid, cmt, ii, addl, scenario,
                  WT, CYP2D6_EM, CYP2D6_PM, CONMED_AZOLE)
  final_dose_time <- (n_doses - 1L) * tau_h
  obs_grid <- sort(unique(c(
    seq(0, 24, by = 0.5),
    seq(24, final_dose_time, by = 6),
    seq(final_dose_time, final_dose_time + 24, by = 0.5),
    seq(final_dose_time + 24, tlast_h, by = 4)
  )))
  obs <- demo |>
    dplyr::select(id, scenario, WT, CYP2D6_EM, CYP2D6_PM, CONMED_AZOLE) |>
    tidyr::crossing(time = obs_grid) |>
    dplyr::mutate(amt = NA_real_, evid = 0L, cmt = NA_character_,
                  ii = NA_real_, addl = NA_integer_)
  dplyr::bind_rows(doses, obs) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

# Protocol 1: single 6 mg oral dose, 0-168 h, no ITZ -- restrict to the
# three no-ITZ scenarios.
events_p1 <- build_single_dose(
  scenarios |> dplyr::filter(CONMED_AZOLE == 0),
  dose_mg = 6, tlast_h = 168
)

# Protocol 2: 3 mg QD oral for 14 days, sampling out to 600 h, no ITZ --
# restrict to a single CYP2D6 group (the paper's Figure 2 pools n = 15
# subjects without stratifying by CYP2D6 group, but inspection of the
# protocol-2 row of Table 1 shows the subjects are predominantly Group 1
# and Group 2 -- the cohort below uses the population-mean weight at each
# phenotype to show the IIV across groups).
events_p2 <- build_multi_dose(
  scenarios |> dplyr::filter(CONMED_AZOLE == 0),
  dose_mg = 3, n_doses = 14L, tau_h = 24, tlast_h = 600
)

# Protocol 3: 3 mg single oral dose with vs without ITZ, 0-366 h, all six
# scenarios contribute (three CYP2D6 groups times two ITZ states).
events_p3 <- build_single_dose(
  scenarios,
  dose_mg = 3, tlast_h = 366
)

stopifnot(!anyDuplicated(unique(events_p1[, c("id", "time", "evid")])))
stopifnot(!anyDuplicated(unique(events_p2[, c("id", "time", "evid")])))
stopifnot(!anyDuplicated(unique(events_p3[, c("id", "time", "evid")])))
```

```{r simulate}
mod <- rxode2::rxode2(readModelDb("Koue_2007_aripiprazole"))

# Typical-value simulation (zero IIV) reproduces the deterministic
# population mean curves the paper overlays in Figures 1-3.
sim_typ_p1 <- rxode2::rxSolve(rxode2::zeroRe(mod), events = events_p1,
                              keep = c("scenario")) |>
  as.data.frame()
sim_typ_p2 <- rxode2::rxSolve(rxode2::zeroRe(mod), events = events_p2,
                              keep = c("scenario")) |>
  as.data.frame()
sim_typ_p3 <- rxode2::rxSolve(rxode2::zeroRe(mod), events = events_p3,
                              keep = c("scenario")) |>
  as.data.frame()
```

# Replicate published figures

```{r figure-1, fig.cap = "Replicates Figure 1 of Koue 2007. Typical-value plasma aripiprazole concentration vs. time following a single 6 mg oral dose, stratified by CYP2D6 phenotype group. The paper's Figure 1 plots observed mean concentrations as Group 1 (blue), Group 2 (green), and Group 3 (red); the solid curves below are the typical-value model predictions at the population-mean weight (60.7 kg)."}
sim_typ_p1 |>
  dplyr::filter(time > 0) |>
  ggplot2::ggplot(ggplot2::aes(time, Cc, color = scenario)) +
  ggplot2::geom_line(linewidth = 0.8) +
  ggplot2::scale_x_continuous(breaks = seq(0, 168, by = 24)) +
  ggplot2::labs(x = "Time (h)",
                y = "Aripiprazole plasma concentration (ng/mL)",
                color = "CYP2D6 phenotype",
                title = "Figure 1 -- single 6 mg oral dose, by CYP2D6 group")
```

```{r figure-2, fig.cap = "Replicates Figure 2 of Koue 2007. Typical-value plasma aripiprazole concentration vs. time for a 14-day course of 3 mg once-daily oral dosing followed by 12 days of follow-up sampling, stratified by CYP2D6 phenotype. The paper's Figure 2 pools all 15 subjects regardless of CYP2D6 phenotype; the three curves below show how the predicted mean profile shifts across the three phenotype strata, with Group 1 (EM) closest to the paper's published mean curve."}
sim_typ_p2 |>
  dplyr::filter(time > 0) |>
  ggplot2::ggplot(ggplot2::aes(time, Cc, color = scenario)) +
  ggplot2::geom_line(linewidth = 0.7) +
  ggplot2::scale_x_continuous(breaks = seq(0, 600, by = 120)) +
  ggplot2::labs(x = "Time (h)",
                y = "Aripiprazole plasma concentration (ng/mL)",
                color = "CYP2D6 phenotype",
                title = "Figure 2 -- 3 mg QD oral for 14 d, by CYP2D6 group")
```

```{r figure-3, fig.cap = "Replicates Figure 3 of Koue 2007. Typical-value plasma aripiprazole concentration vs. time following a single 3 mg oral dose with and without coadministered itraconazole (100 mg QD for 21 d, aripiprazole dosed on day 7 of itraconazole), shown on a log Y-axis. The paper's Figure 3 plots observed mean concentrations as aripiprazole alone (light blue) and aripiprazole with itraconazole (orange). The curves below stratify additionally by CYP2D6 phenotype to expose the joint effect of the two clearance covariates."}
sim_typ_p3 |>
  dplyr::filter(time > 0, Cc > 0) |>
  ggplot2::ggplot(ggplot2::aes(time, Cc, color = scenario,
                                linetype = factor(CONMED_AZOLE))) +
  ggplot2::geom_line(linewidth = 0.7) +
  ggplot2::scale_x_continuous(breaks = seq(0, 360, by = 48)) +
  ggplot2::scale_y_log10() +
  ggplot2::labs(x = "Time (h)",
                y = "Aripiprazole plasma concentration (ng/mL, log scale)",
                color = "CYP2D6 phenotype + ITZ",
                linetype = "Itraconazole (1 = yes)",
                title = "Figure 3 -- single 3 mg oral dose +/- itraconazole")
```

# PKNCA validation

PKNCA computes Cmax, Tmax, AUC0-inf, and apparent terminal half-life over
the protocol-1 single-dose 6 mg profile, stratified by CYP2D6 phenotype.
The grouping is on `scenario` so the per-phenotype values can be compared
directly against the half-lives that Koue 2007 reports in the Results and
Discussion section.

```{r pknca-p1}
sim_nca_p1 <- sim_typ_p1 |>
  dplyr::filter(!is.na(Cc), time > 0) |>
  dplyr::select(id, time, Cc, scenario)

dose_nca_p1 <- events_p1 |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, scenario)

conc_obj_p1 <- PKNCA::PKNCAconc(sim_nca_p1, Cc ~ time | scenario + id,
                                 concu = "ng/mL", timeu = "h")
dose_obj_p1 <- PKNCA::PKNCAdose(dose_nca_p1, amt ~ time | scenario + id,
                                 doseu = "mg")

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

nca_res_p1 <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  conc_obj_p1, dose_obj_p1, intervals = intervals_p1
))

nca_tbl_p1 <- as.data.frame(nca_res_p1$result) |>
  dplyr::select(scenario, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  dplyr::mutate(dplyr::across(where(is.numeric), \(x) signif(x, 3)))

knitr::kable(nca_tbl_p1,
             caption = paste0("Protocol-1 single-dose 6 mg PKNCA at the population-mean ",
                              "weight (60.7 kg), typical-value simulation. cmax in ng/mL; ",
                              "tmax and half.life in h; aucinf.obs in ng*h/mL."))
```

## Comparison against published half-lives

Koue 2007's Results and Discussion paragraph reports observed mean
terminal half-lives at 6 mg single dose of 44.2 h (Group 1), 47.7 h
(Group 2), and 76.3 h (Group 3), and at 3 mg single dose of 65.2 h
(aripiprazole alone) and 77.9 h (with itraconazole). These half-lives
were computed empirically by the authors from the observed concentration
data using post-hoc Bayesian individual estimates, NOT directly from the
population structural parameters; depending on how late the analytical
window extends, an empirical half-life can capture a mix of the alpha and
beta phases of a two-compartment disposition. The PKNCA half-lives below
are computed from the same typical-value structural model that produced
Figures 1-3; differences between the two values are documented under
Assumptions and deviations.

```{r pknca-vs-published}
published_p1 <- tibble::tibble(
  scenario   = c("Group 1 EM (no ITZ)", "Group 2 IM (no ITZ)", "Group 3 PM (no ITZ)"),
  hl_published = c(44.2, 47.7, 76.3)
)

comparison_p1 <- nca_tbl_p1 |>
  dplyr::select(scenario, half.life) |>
  dplyr::rename(hl_pknca = half.life) |>
  dplyr::left_join(published_p1, by = "scenario") |>
  dplyr::mutate(ratio = signif(hl_pknca / hl_published, 3))

knitr::kable(comparison_p1,
             caption = paste0("Apparent terminal half-life (h) -- PKNCA on the typical-value ",
                              "simulated protocol-1 profile vs. Koue 2007's reported observed ",
                              "mean half-lives from Results and Discussion."))
```

# Assumptions and deviations

* **Weight set to the population mean.** Every scenario above uses
  `WT = 60.7 kg`, the population-mean weight reported in the
  Pharmacokinetic Data section. The model's body-weight scaling is a
  fixed linear-per-kg multiplier on every clearance and volume
  parameter, so the relative-shape comparison against the paper's
  figures (which overlay observed individual data with a typical-value
  curve) is well-posed at any single weight; the choice of 60.7 kg
  matches the paper's population mean.
* **Single-dose protocol-1 cohort sized as published.** Figure 1 reports
  Group 1 n=5, Group 2 n=10, Group 3 n=11, which is fewer subjects than
  the genotype counts in Table 1 because not every Table-1 subject
  contributed to the protocol-1 concentration-time analysis. The
  vignette uses one virtual subject per phenotype because the simulation
  is typical-value (no IIV); reproducing the per-group n is unnecessary.
* **Multi-dose protocol-2 stratification.** The paper's Figure 2 pools
  all 15 protocol-2 subjects regardless of CYP2D6 group; the vignette's
  Figure 2 replicate plots all three phenotype curves so a downstream
  user can read off the CYP2D6 effect at steady state. The Group 1 (EM)
  curve is the closest match to the paper's published mean curve.
* **CYP2D6 reference category preserved.** Koue 2007's Eq. 4 takes
  Group 1 (EM) as the explicit reference, NOT IM as in the canonical
  three-level CYP2D6_PM + CYP2D6_EM register orientation. The model
  preserves the paper's EM reference by deriving the IM indicator
  inside `model()` as `1 - CYP2D6_EM - CYP2D6_PM`; the lcl parameter is
  therefore the EM-reference CL/F per kg (0.0645 L/h/kg).
* **omega values transcribed as SDs (not variances).** Koue 2007
  Table 2's omega column entries are the per-parameter standard
  deviations on the log scale (the text above Eq. 1-6 reads "variance
  of omega^2_ALAG" -- i.e., omega^2 is the variance and omega is the
  SD). nlmixr2's `~ value` syntax declares the variance, so each Table-2
  omega is squared in `ini()`.
* **Lognormal residual error.** Koue 2007 Eq. 7 writes
  `C_ij = C*_ij * exp(eps_ij)` with `eps ~ N(0, sigma^2)`. This is the
  lognormal / exponential residual structure; the model file uses
  `Cc ~ lnorm(expSd)` with `expSd = sigma = 0.166`.
* **PKNCA half-life vs paper's reported half-life.** The half-lives
  reported in Results and Discussion are empirical Bayesian post-hoc
  estimates from the observed data (e.g. log-linear regression on the
  terminal-phase observations of each subject's individual fitted
  profile), NOT direct functions of the typical-value structural
  parameters. PKNCA's typical-value half-life can therefore differ from
  the published Bayesian half-life even when the structural model is
  perfectly transcribed; that difference is a property of the
  observation grid and terminal-phase windowing the authors used to
  compute the published values and is not a transcription error.
* **Protocol-1 abstract n vs. Methods n.** The abstract reports 24
  contributing subjects in the single-dose trial; the Pharmacokinetic
  Data section reports 26 enrolled, of whom 21 are reflected in
  Figure 1's per-group ns (5 + 10 + 11). The model file's `population`
  metadata records 26 (the Pharmacokinetic Data section value, which is
  the authoritative source); the abstract value is noted in `notes`
  for traceability.
