---
title: "Antipsychotic PANSS subscales (Pilla Reddy 2013 Part II)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Antipsychotic PANSS subscales (Pilla Reddy 2013 Part II)}
  %\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(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

This vignette demonstrates the Part II PK/PD models for five antipsychotic
drugs (one typical FGA, haloperidol; four atypical SGAs: risperidone,
olanzapine, ziprasidone, paliperidone) against the three PANSS subscales
(positive, negative, general).

* Citation: Pilla Reddy V, Kozielska M, Suleiman AA, Johnson M, Vermeulen A,
  Liu J, de Greef R, Groothuis GMM, Danhof M, Proost JH (2013).
  Pharmacokinetic-pharmacodynamic modelling of antipsychotic drugs in patients
  with schizophrenia: Part II: The use of subscales of the PANSS score.
  Schizophrenia Research 146(1-3):153-161.
* Article: <https://doi.org/10.1016/j.schres.2013.02.010>
* Companion paper (PANSS total): Pilla Reddy et al. 2013 Part I, Schizophrenia
  Research 146(1-3):144-152. doi:10.1016/j.schres.2013.02.011.

The Part II analysis uses Part I's PK model structure unchanged and develops
new PD layers for the three PANSS subscales. Each model in this set encodes:

1. A one-compartment representation of the drug PK from Part I (sufficient
   for the steady-state concentration Css that drives the Emax PD).
2. A Weibull placebo model with subscale-specific baseline, Pmax, TD, and POW
   from Part II Table 1.
3. A drug-specific Emax / EC50 / KT triplet per subscale from Part II Table 2,
   with KT capturing the time delay to maximum drug effect via the factor
   (1 - exp(-KT * t)).
4. Multi-output observation of plasma Cc and three PANSS subscale outputs
   (PANSS_pos, PANSS_neg, PANSS_gen).

The exponential time-to-event dropout sub-model (Part II Table 4) is
documented in each model's `population$dropout_model` field but is not
simulated in this vignette; PANSS trajectories assume no dropout.

# Population

The Part II analysis pooled 741 placebo-treated patients with subscale data
available across 12 industry-sponsored Phase II / III double-blind clinical
trials conducted between 1989 and 2009 (Part II Methods; trial details in
Part I Table 1). Active-treatment cohorts span haloperidol, risperidone,
olanzapine, ziprasidone, and paliperidone arms at the clinically relevant
dose ranges listed below.

```{r population_table}
pop <- tibble::tribble(
  ~Drug,           ~`Typical dose (mg/day)`, ~`Regimen`,    ~`Trial design`,
  "Haloperidol",   "5-10",                   "qd or bid",   "Phase III + LMU open-label",
  "Risperidone",   "0.5-8",                  "qd or bid",   "INT-2, INT-3 (Phase III)",
  "Olanzapine",    "10-15",                  "qd",          "Phase III arms (SCH-303/304/305/2002)",
  "Ziprasidone",   "20-100",                 "qd or bid",   "Phase III (128-114/115/303/307)",
  "Paliperidone",  "3-15",                   "qd",          "Phase III (SCH-303/304/305)"
)
knitr::kable(pop, caption = "Active-treatment arms in the Part II pool.")
```

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

# Source trace

The per-parameter origin is recorded as an in-file comment next to each
`ini()` entry in the five model files under `inst/modeldb/specificDrugs/`.
The table below collects the structural source locations.

| Equation / parameter group              | Source location                                    |
|-----------------------------------------|----------------------------------------------------|
| Drug PK (1-cmt simplification)          | Part I (PMID 23473810) Table 2                     |
| Placebo Weibull (BASL, Pmax, TD, POW)   | Part II Table 1; per-drug BASL from Part II Table 2|
| Drug Emax / EC50 / KT per subscale      | Part II Table 2                                    |
| IIV BASL (% CV), IIV Pmax (SD)          | Part II Table 1                                    |
| IIV Emax (SD), IIV EC50 (% CV, neg only)| Part II Table 2                                    |
| Residual error (additive SD)            | Part II Table 2 (per-drug, joint PKPD)             |
| Dropout BHAZ (per drug, per subscale)   | Part II Table 4 (documented, not simulated)        |
| Effective Css per drug per subscale     | Part II Table 3 (validation target below)          |

# Virtual cohort and simulation

```{r helpers}
# Typical clinically relevant doses (mg/day) used as illustration -- one each
# from the dose ranges in Part I Table 1. These match the doses for which
# Part II Figure 1b reports model-predicted % change from baseline PANSS.
typical_doses <- tibble::tribble(
  ~drug,            ~dose_mg, ~tau_h,
  "haloperidol",    7.5,      24,
  "risperidone",    4,        12,    # bid dosing in INT-2 / INT-3
  "olanzapine",     10,       24,
  "ziprasidone",    80,       12,    # bid in 128-115
  "paliperidone",   9,        24
)

# For each drug, build the event table once. dvid = 1L addresses the Cc
# endpoint (which forces the multi-output mapping); the other algebraic
# observables PANSS_pos/PANSS_neg/PANSS_gen are returned automatically as
# columns on every observation row regardless of dvid.
make_events_for_drug <- function(drug, dose_mg, tau_h, n_doses = 42, obs_hours) {
  ndoses_addl <- ceiling((max(obs_hours) + 24) / tau_h)
  dose_rows <- tibble::tibble(
    id   = 1L, time = 0,            amt  = dose_mg,
    cmt  = "depot",                 evid = 1L,
    ii   = tau_h,                   addl = as.integer(ndoses_addl),
    dvid = NA_integer_
  )
  obs_rows <- tibble::tibble(
    id   = 1L, time = obs_hours,    amt  = NA_real_,
    cmt  = NA_character_,           evid = 0L,
    ii   = NA_real_,                addl = NA_integer_,
    dvid = 1L
  )
  dplyr::bind_rows(dose_rows, obs_rows) |> dplyr::arrange(id, time)
}

# Six-week observation grid in hours, sampled every 6 hours; enough to track
# PANSS trajectory smoothly without inflating the render time.
obs_hours <- seq(0, 42 * 24, by = 6)

simulate_drug <- function(drug, dose_mg, tau_h) {
  fn  <- paste0("PillaReddy_2013_", drug, "_panss_subscales")
  mod <- readModelDb(fn) |> rxode2::zeroRe()
  ev  <- make_events_for_drug(drug, dose_mg, tau_h, obs_hours = obs_hours)
  sim <- rxode2::rxSolve(mod, events = ev) |>
    as.data.frame()
  sim$drug <- drug
  sim$dose_mg <- dose_mg
  sim
}

sim_all <- purrr::map_dfr(seq_len(nrow(typical_doses)), function(i) {
  with(typical_doses[i, ], simulate_drug(drug, dose_mg, tau_h))
})
```

# Replicate Figure 1b -- model-predicted PANSS time course

Figure 1b of Pilla Reddy 2013 Part II overlays the model-predicted percentage
change from baseline PANSS total and each of its three subscales for placebo
and for each of the five antipsychotics. The traces below reproduce this
qualitative structure using the typical-value (zero-IIV) trajectories of the
five drug-specific models.

```{r figure-1b, fig.width = 9, fig.height = 7}
# Compute % change from baseline per drug per subscale.
sim_pct <- sim_all |>
  dplyr::group_by(drug) |>
  dplyr::arrange(time) |>
  dplyr::mutate(
    PANSS_pos_pct = 100 * (PANSS_pos - dplyr::first(PANSS_pos)) / dplyr::first(PANSS_pos),
    PANSS_neg_pct = 100 * (PANSS_neg - dplyr::first(PANSS_neg)) / dplyr::first(PANSS_neg),
    PANSS_gen_pct = 100 * (PANSS_gen - dplyr::first(PANSS_gen)) / dplyr::first(PANSS_gen)
  ) |>
  dplyr::ungroup() |>
  dplyr::mutate(time_days = time / 24) |>
  dplyr::select(drug, dose_mg, time_days, PANSS_pos_pct, PANSS_neg_pct, PANSS_gen_pct) |>
  tidyr::pivot_longer(
    cols       = c(PANSS_pos_pct, PANSS_neg_pct, PANSS_gen_pct),
    names_to   = "subscale",
    values_to  = "pct_change"
  ) |>
  dplyr::mutate(
    subscale = dplyr::recode(subscale,
      PANSS_pos_pct = "Positive",
      PANSS_neg_pct = "Negative",
      PANSS_gen_pct = "General"
    ),
    subscale = factor(subscale, levels = c("Positive", "Negative", "General"))
  )

ggplot(sim_pct, aes(time_days, pct_change, colour = subscale)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = -30, linetype = "dashed", colour = "grey40") +
  facet_wrap(~ drug, nrow = 2, scales = "fixed") +
  scale_colour_manual(values = c(Positive = "#1f77b4", Negative = "#d62728", General = "#ff7f0e")) +
  labs(
    x       = "Time (days)",
    y       = "% change from baseline PANSS subscale",
    colour  = "Subscale",
    title   = "Model-predicted change from baseline PANSS subscales over six weeks",
    caption = paste0(
      "Replicates the SGA and FGA panels of Figure 1b of Pilla Reddy 2013",
      " Part II.\nDoses: haloperidol 7.5, risperidone 4, olanzapine 10,",
      " ziprasidone 80, paliperidone 9 mg/day. Dashed line marks the 30% target."
    )
  ) +
  theme_minimal()
```

The model qualitatively reproduces the published behaviours:

- All five drugs reduce all three subscales below baseline (negative % change).
- Olanzapine shows the largest reduction on the negative subscale, consistent
  with Part II Results ("Olanzapine showed a better effect towards negative
  symptoms").
- Ziprasidone shows the slowest onset on the negative subscale (very low
  `KT_neg = 0.0073 1/day`, half-time to maximum drug effect approximately 95
  days), consistent with Part II Table 2 and Discussion.
- Haloperidol reaches its plateau on all three subscales fastest, consistent
  with its short half-life and KT values (0.11-0.19 1/day vs 0.028-0.16
  1/day for the SGAs).

# Effective Css and dose validation (Part II Table 3)

Part II Table 3 reports the steady-state concentration Css required to
produce a 30% reduction in PANSS from baseline, derived from the PKPD model
parameters via

```
Ceff = EC50 / (Emax / (1 - PANSS / (Baseline_PANSS * (1 - Pmax))) - 1)
```

where the target PANSS = -30/100 * (Baseline_PANSS - n_items) + Baseline_PANSS
and n_items = 7 (positive), 7 (negative), 16 (general). The table below
re-derives this from the in-file `ini()` values and compares to the published
Table 3 numbers.

```{r effective_css}
drugs <- typical_doses$drug

calc_ceff <- function(drug) {
  mb  <- readModelDb(paste0("PillaReddy_2013_", drug, "_panss_subscales"))()
  ini <- as.list(mb$iniDf$est)
  names(ini) <- mb$iniDf$name
  # Per-subscale calculation
  one <- function(subscale, n_items) {
    basl <- ini[[paste0("basl_", subscale)]]
    pmax <- ini[[paste0("pmax_", subscale)]]
    emax <- ini[[paste0("emax_", subscale)]]
    ec50 <- ini[[paste0("ec50_", subscale)]]
    target_panss <- basl + (-30 / 100) * (basl - n_items)
    denom <- emax / (1 - target_panss / (basl * (1 - pmax))) - 1
    if (is.na(denom) || denom <= 0) {
      ceff <- NA_real_
    } else {
      ceff <- ec50 / denom
    }
    tibble::tibble(drug = drug, subscale = subscale, ceff_ng_per_mL = ceff)
  }
  dplyr::bind_rows(
    one("pos", 7L),
    one("neg", 7L),
    one("gen", 16L)
  )
}

ceff_sim <- purrr::map_dfr(drugs, calc_ceff)

ceff_pub <- tibble::tribble(
  ~drug,            ~subscale, ~ceff_pub_ng_per_mL,
  "haloperidol",    "pos",     0.54,
  "haloperidol",    "neg",     31,
  "haloperidol",    "gen",     3.2,
  "risperidone",    "pos",     6.0,
  "risperidone",    "neg",     NA,    # marked "#" in Part II Table 3 (not attained)
  "risperidone",    "gen",     NA,
  "olanzapine",     "pos",     4.9,
  "olanzapine",     "neg",     13.4,
  "olanzapine",     "gen",     9.2,
  "ziprasidone",    "pos",     48.3,
  "ziprasidone",    "neg",     NA,
  "ziprasidone",    "gen",     NA,
  "paliperidone",   "pos",     3.5,
  "paliperidone",   "neg",     NA,
  "paliperidone",   "gen",     22.6
)

ceff_cmp <- dplyr::full_join(ceff_sim, ceff_pub, by = c("drug", "subscale")) |>
  dplyr::mutate(
    pct_diff = 100 * (ceff_ng_per_mL - ceff_pub_ng_per_mL) / ceff_pub_ng_per_mL
  )

knitr::kable(
  ceff_cmp,
  digits  = 2,
  caption = "Effective Css for 30% reduction in PANSS subscale at six weeks. Simulated values are re-derived from the model files' ini() entries via the closed-form formula. Cells marked NA in ceff_pub indicate that the 30% reduction was not attainable per Part II Table 3."
)
```

The closed-form derived Css matches Part II Table 3 within reading
precision; small deviations come from rounding in the published table entries.

# Sanity check -- simulated PANSS at typical doses

Snapshot the PANSS values at the end of the simulation (week 6) for each
drug at its typical dose. The reductions below baseline reflect the
trajectories shown in the Figure 1b replicate above.

```{r typical_panss_week6}
week6 <- sim_all |>
  dplyr::filter(abs(time - 42 * 24) < 12) |>
  dplyr::group_by(drug) |>
  dplyr::summarise(
    Cc_ng_mL   = mean(Cc),
    PANSS_pos  = mean(PANSS_pos),
    PANSS_neg  = mean(PANSS_neg),
    PANSS_gen  = mean(PANSS_gen),
    .groups    = "drop"
  ) |>
  dplyr::left_join(typical_doses |> dplyr::select(drug, dose_mg), by = "drug")
knitr::kable(week6, digits = 2,
             caption = "Simulated typical PANSS subscale values at week 6 for each drug at its typical dose.")
```

# Assumptions and deviations

- **PK simplification.** The Part I PK models for haloperidol (2-cmt) and
  risperidone (parent + metabolite, CYP2D6-stratified) are simplified to
  1-compartment models in these Part II files. The simplification preserves
  the steady-state average concentration Css = Dose / (CL * tau) that drives
  the Emax PD; intra-dose Cc fluctuations differ from Part I's predictions
  but the Css-based PD prediction is unchanged for the purposes of this
  Part II reproduction.
- **Risperidone CL\_AM/F.** The Part II PD model uses Css of the active
  moiety (parent risperidone + 9-OH-risperidone). The active-moiety CL/F =
  6.3 L/h used here is derived from Part II Table 3's effective-dose /
  effective-Css pair (0.8 mg/day -> 5.3 ng/mL for PANSS total at 30%
  reduction), giving CL/F = Dose / (Css * tau). CYP2D6 phenotype stratifies
  the parent risperidone CL/F by an order of magnitude in Part I; this Part
  II model collapses the phenotype distribution to a single typical CL.
- **Paliperidone absorption.** The Part I model uses sequential zero-order
  plus first-order absorption with a 23.6 h zero-order duration and 0.67 h
  lag. These details are omitted in the 1-cmt simplification here because
  Css at steady state is invariant to within-dose absorption shape.
- **PD time scale.** The model time `t` is in hours (consistent with the PK
  ODE); the placebo Weibull and the Emax (1 - exp(-KT * t)) onset use t_days
  = t / 24 because the paper reports TD in days and KT in 1/day.
- **Covariates not implemented.** Part II Table 1 lists covariate effects of
  disease state (acute / chronic), study geographic origin (USA / non-USA),
  dosing regimen (qd / bid), and study duration (short / long) on the
  placebo Weibull baseline / Pmax / RUV. The typical-individual reference
  simulation targets the reference categories (acute, USA, bid, short). The
  per-drug `covariatesDataExcluded` field documents each excluded covariate
  with its Part II Table 1 effect coefficient.
- **Dropout not simulated.** Part II Table 4 reports an exponential
  time-to-event dropout sub-model with subscale-specific baseline hazards
  and BETA coefficients per drug. The joint PANSS + dropout model is the
  Part II final model, but only PANSS trajectories are simulated here; the
  dropout parameters are documented in `population$dropout_model` for each
  drug.
- **Inter-individual variability on RUV (IIV-RUV).** Part II Table 1 reports
  inter-individual variability in the RUV magnitude (28%, 37%, 30% CV for
  positive, negative, general respectively). nlmixr2 / rxode2 does not
  natively encode IIV on the residual-error magnitude; the RUV here is a
  fixed per-subscale value.
- **Css-based vs Cc-based driver.** Part II's PD equation uses the steady-
  state concentration Css as the driver; the rxode2 implementation feeds the
  time-varying plasma concentration Cc into the Emax term. At steady state
  Cc(t) oscillates around Css with magnitude that depends on the drug's
  half-life relative to the dosing interval. For drugs with long half-life
  (olanzapine, paliperidone, ziprasidone) the time-average of Cc is very
  close to Css. For drugs with shorter half-life (haloperidol, risperidone)
  the instantaneous PANSS prediction oscillates by a small amount during
  the dosing interval while the time-averaged PANSS over a day still
  matches the Css-based prediction.
