---
title: "Nusinersen (Biliouris 2018)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Nusinersen (Biliouris 2018)}
  %\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

* Citation: `r rxode2::rxode(readModelDb("Biliouris_2018_nusinersen"))$reference`
* Description: `r rxode2::rxode(readModelDb("Biliouris_2018_nusinersen"))$description`
* Article: <https://doi.org/10.1002/psp4.12323>

The on-disk source is the paper's simulation NONMEM control stream
(`$PROBLEM SMA MONKEY`, `$SIMULATION (12345678) ONLYSIM
SUBPROBLEM=1000`). All `$THETA`, `$OMEGA`, and `$SIGMA` entries carry
the `FIX` flag, so every value in the packaged `ini()` block is wrapped
in `fixed()`. See "Assumptions and deviations" below for the paper-text
sections not on disk with the control stream.

# Population

The Biliouris 2018 model was originally fitted to cynomolgus-monkey
plasma and CNS-tissue data (the `$PROBLEM` tag "SMA MONKEY" identifies
the fit dataset). The on-disk simulation control stream then re-uses
the monkey-fit parameter estimates for a paediatric-human extrapolation:
volumes are linearly weight-scaled from a 2.8 kg reference, rate
constants use a slight `(WT/2.8)^(-0.08)` scaling, and the CSF
physiological volume `V_CSF` is overridden by a stepwise age function
that matches human paediatric CSF anatomy (120, 130, 135, 140 mL for
age classes < 0.25, 0.25 to < 0.5, 0.5 to < 1, 1 to < 2 years; adult
reference 150 mL). Full baseline demographics for the fitted monkey
cohort are not on disk with the supplementary simulation stream; refer
to Biliouris 2018 main text Methods and Table 1 (published in CPT
Pharmacometrics Syst Pharmacol) for `n_subjects`, age, and sex
distributions.

The population metadata packaged with the model is available
programmatically:

```{r population}
mod_fun <- readModelDb("Biliouris_2018_nusinersen")
str(mod_fun()$population, max.level = 1)
```

# Source trace

Every value below is traced to the paper's simulation NONMEM control
stream (on disk as `PMID_30043511.pdf`, five pages containing the
`$PROBLEM SMA MONKEY` model verbatim). The full `$THETA(k) FIX` /
`$OMEGA k FIX` label is repeated as the in-file comment on each `ini()`
entry.

| Structural element | Value | Source location |
|--------------------|-------|-----------------|
| `V1 = 150 mL` (CSF adult) | `lv_csf = fixed(log(150))` | `$THETA(1) 150 FIX ; (1. V1) mL` |
| `V2 = 937 mL` (plasma) | `lv_plasma = fixed(log(937))` | `$THETA(2) 937 FIX ; (2. V2) mL` |
| `V3 = 1.91 mL` (cervical SC) | `lv_cerv = fixed(log(1.91))` | `$THETA(3) 1.91 FIX ; (3. V3) mL` |
| `V4 = 53.8 mL` (brain) | `lv_brain = fixed(log(53.8))` | `$THETA(4) 53.8 FIX ; (4. V4) mL` |
| `V6 = 1.08 mL` (lumbar SC) | `lv_lumb = fixed(log(1.08))` | `$THETA(20) 1.08 FIX ; (20. V6) mL` |
| `V8 = 1.52 mL` (thoracic SC) | `lv_thor = fixed(log(1.52))` | `$THETA(25) 1.52 FIX ; (25. V8) mL` |
| `V9 = 2.11 mL` (pons) | `lv_pons = fixed(log(2.11))` | `$THETA(28) 2.11 FIX ; (28. V9) mL` |
| `K12 = 0.0891 1/h` (CSF -> plasma) | `lk_csf_plasma` | `$THETA(9) 0.0891 FIX ; (9. K12)` |
| `K20 = 0.206 1/h` (plasma elim) | `lkel` | `$THETA(10) 0.206 FIX ; (10. K20)` |
| `K13/K31, K14/K41, K25/K52, K16/K61, K47/K74, K18/K81, K19/K91` | all `fixed(log(...))` | `$THETA(5..8, 11..12, 21..24, 26..27, 29..30)` |
| `V1(AGE)` overrides: 120/130/135/140 mL for AGE < 0.25 / < 0.5 / < 1 / < 2 y | in `model()` `ifelse` chain | `$PK IF (AGE.LT.0.25) TVV1 = 120` etc. |
| `Volume WT scaling`: `V = TVV * (WT/2.8)` for V2..V9 (not V1) | `wt_ratio = WT/2.8` in `model()` | `$PK TVV_i = THETA(_i) * (WT/2.8)`; V1's `(WT/2.8)` is commented out |
| `Rate WT scaling`: `K = TVK * (WT/2.8)^(-0.08)` | `wt_kfac = wt_ratio^e_wt_k` | `$PK TVK_ij = THETA(_) * 1/((WT/2.8)**0.08)` |
| Twelve exponential IIVs on K20, K12, V1, V2, K13, K14, K41, K16, K18, K81, K19, K91 | `~ fixed(<var>)` diagonal | `$OMEGA` variances 0.421, 1.3, 0.854, 0.832, 0.453, 13.6, 0.153, 0.124, 0.28, 0.129, 0.654, 0.289 (all FIX) |
| Proportional residual SD per tissue | `propSd_Ccsf` = 0.714, `propSd` = 0.412, `propSd_Ccerv` = 0.101, `propSd_Cbrain` = 1.77, `propSd_Clumb` = 0.312, `propSd_Cthor` = 0.0239, `propSd_Cpons` = 0.0103 | `$THETA(13..19) FIX ; ERR1..ERR9`; `$SIGMA 1 FIX` combined with `W = SQRT(THETA^2 * IPRED^2)` makes THETA the SD-fraction |

The nine ODEs replicate the paper's `$DES` block one-for-one; the
in-line comments in `inst/modeldb/specificDrugs/Biliouris_2018_nusinersen.R`
give the paper-to-model mapping compartment-by-compartment.

# Virtual cohort

The paper does not provide observed data. The figures below use two
paediatric-focused virtual populations for a single 12 mg intrathecal
lumbar-puncture bolus (the fixed dose used in the ISIS 396443 clinical
programme).

**Cohort 1 (typical-value trajectories)**: one representative
five-year-old (15 kg, AGE 5 y). All random effects zeroed.
`V_CSF = 150 mL` (adult regime).

**Cohort 2 (pediatric VPC)**: 100 virtual paediatric subjects with
weight uniformly in 5 to 20 kg and age uniformly in 0.1 to 2 y. The
age range exercises the `V_CSF` step function (each of the four
paediatric bins is sampled).

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

obs_grid <- c(0, 0.5, 1, 2, 4, 8, 12, 24, 36, 48, 72, 96, 120, 144, 168)

# Cohort 1 -- single-subject typical trajectories
cohort1_events <- {
  dose <- data.frame(id = 1L, time = 0, amt = 12, cmt = "csf",
                     evid = 1L, dvid = NA_integer_,
                     WT = 15, AGE = 5, cohort = "typical (15 kg, 5 y)",
                     stringsAsFactors = FALSE)
  obs <- expand.grid(id = 1L, time = obs_grid,
                     dvid = 1:7, stringsAsFactors = FALSE)
  obs$amt   <- NA_real_
  obs$cmt   <- c("Ccsf", "Cc", "Ccerv", "Cbrain",
                 "Clumb", "Cthor", "Cpons")[obs$dvid]
  obs$evid  <- 0L
  obs$WT    <- 15
  obs$AGE   <- 5
  obs$cohort <- "typical (15 kg, 5 y)"
  rbind(dose, obs)
}

# Cohort 2 -- 100-subject paediatric VPC
n_ped <- 100L
cohort2_events <- {
  ped_wt  <- runif(n_ped, 5, 20)
  ped_age <- runif(n_ped, 0.1, 2)
  dose_rows <- data.frame(
    id    = seq_len(n_ped) + 100L,
    time  = 0,
    amt   = 12,
    cmt   = "csf",
    evid  = 1L,
    dvid  = NA_integer_,
    WT    = ped_wt,
    AGE   = ped_age,
    cohort = "paediatric VPC (5-20 kg, 0.1-2 y)",
    stringsAsFactors = FALSE
  )
  obs_rows <- do.call(rbind, lapply(seq_len(n_ped), function(i) {
    o <- expand.grid(id = 100L + i, time = obs_grid,
                     dvid = c(1L, 2L),  # Ccsf and Cc only, to keep the run light
                     stringsAsFactors = FALSE)
    o$amt   <- NA_real_
    o$cmt   <- c("Ccsf", "Cc")[o$dvid]
    o$evid  <- 0L
    o$WT    <- ped_wt[i]
    o$AGE   <- ped_age[i]
    o$cohort <- "paediatric VPC (5-20 kg, 0.1-2 y)"
    o
  }))
  rbind(dose_rows, obs_rows)
}

events <- rbind(cohort1_events, cohort2_events)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid", "dvid")])))
```

# Simulation

```{r simulate}
mod <- readModelDb("Biliouris_2018_nusinersen")
sim <- rxode2::rxSolve(mod, events = events,
                       keep = c("WT", "AGE", "cohort", "dvid"))
sim_df <- as.data.frame(sim)
```

Typical-value trajectories (Cohort 1) use the zero-random-effect
package model to replicate the paper's semi-mechanistic prediction. VPC
percentiles (Cohort 2) apply the paper's fixed OMEGA variances.

```{r simulate-typical}
mod_typ <- rxode2::zeroRe(mod)
sim_typ_df <- as.data.frame(
  rxode2::rxSolve(mod_typ, events = cohort1_events,
                  keep = c("WT", "AGE", "cohort"))
)
```

# Typical-value time course across sampled tissues

The seven `~ prop()` residual-error clauses correspond to the seven
sampled matrices in the paper's `$ERROR` block: CSF, plasma, cervical
spinal cord, brain, lumbar spinal cord, thoracic spinal cord, and
pons. Each panel below traces its algebraic observable (`Ccsf`, `Cc`,
`Ccerv`, `Cbrain`, `Clumb`, `Cthor`, `Cpons`) over 7 days after a
single 12 mg intrathecal bolus in the typical 15 kg / 5 y subject.

```{r typical-trajectories}
tissue_labels <- c(
  Ccsf   = "CSF",
  Cc     = "Plasma",
  Ccerv  = "Cervical spinal cord",
  Cbrain = "Brain",
  Clumb  = "Lumbar spinal cord",
  Cthor  = "Thoracic spinal cord",
  Cpons  = "Pons"
)

sim_typ_long <- sim_typ_df |>
  dplyr::select(time, dplyr::all_of(names(tissue_labels))) |>
  tidyr::pivot_longer(-time, names_to = "obs", values_to = "conc") |>
  dplyr::mutate(tissue = factor(tissue_labels[obs],
                               levels = unname(tissue_labels)))

ggplot(sim_typ_long, aes(time, pmax(conc, 1e-3))) +
  geom_line() +
  facet_wrap(~tissue, scales = "free_y") +
  scale_y_log10() +
  labs(x = "Time after 12 mg IT bolus (h)",
       y = "Concentration (ng/mL)",
       title = "Typical-value concentrations across the seven sampled matrices",
       caption = paste0("Semi-mechanistic 9-compartment nusinersen model ",
                        "(Biliouris 2018); zero random effects; ",
                        "typical 15 kg, 5 y subject."))
```

# Virtual paediatric VPC on CSF and plasma

CSF (dose site) and plasma (systemic) are the two matrices where the
paper's model is most directly informative for clinical dosing.

```{r vpc-csfplasma}
vpc_df <- sim_df |>
  dplyr::filter(cohort == "paediatric VPC (5-20 kg, 0.1-2 y)") |>
  dplyr::select(id, time, WT, AGE, dvid, Ccsf, Cc) |>
  dplyr::mutate(matrix = c("CSF (ng/mL)", "Plasma (ng/mL)")[dvid],
                conc   = ifelse(dvid == 1L, Ccsf, Cc)) |>
  dplyr::filter(!is.na(conc), conc > 0)

vpc_summary <- vpc_df |>
  dplyr::group_by(matrix, time) |>
  dplyr::summarise(
    Q05 = quantile(conc, 0.05, na.rm = TRUE),
    Q50 = quantile(conc, 0.50, na.rm = TRUE),
    Q95 = quantile(conc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(vpc_summary, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line() +
  facet_wrap(~matrix, scales = "free_y") +
  scale_y_log10() +
  labs(x = "Time after 12 mg IT bolus (h)",
       y = "Concentration (ng/mL)",
       title = "Paediatric VPC (n = 100, 5-20 kg, 0.1-2 y): CSF + plasma",
       caption = "5th, 50th, 95th percentiles of virtual subjects.")
```

# PKNCA on typical-value plasma exposure

The paper does not tabulate NCA parameters for the semi-mechanistic
extrapolation, so the block below serves as a sanity summary rather
than a validation-table comparison: it confirms that the packaged
model reproduces plausible plasma Cmax, Tmax, terminal half-life, and
AUC after a single 12 mg IT dose in the typical paediatric subject.

```{r pknca}
# rxSolve drops the id column for single-subject sims; add it back so
# PKNCA gets a per-subject key.
if (!"id" %in% names(sim_typ_df)) sim_typ_df$id <- 1L

pk_conc <- sim_typ_df |>
  dplyr::select(id = id, time = time, Cc = Cc, cohort = cohort) |>
  dplyr::filter(!is.na(Cc))

# Time-zero guarantee -- PKNCA needs a time = 0 row per (id, cohort).
pk_conc <- dplyr::bind_rows(
  pk_conc,
  pk_conc |> dplyr::distinct(id, cohort) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, cohort, time, .keep_all = TRUE) |>
  dplyr::arrange(id, cohort, time)

conc_obj <- PKNCA::PKNCAconc(pk_conc, Cc ~ time | cohort + id)

dose_df <- cohort1_events |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, cohort)
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | cohort + id)

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

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

nca_res$result |>
  dplyr::select(cohort, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  dplyr::rename(
    "Cohort"                       = cohort,
    "Cmax (ng/mL)"                 = cmax,
    "Tmax (h)"                     = tmax,
    "AUC[0,inf] observed (ng.h/mL)" = aucinf.obs,
    "t1/2 (h)"                     = half.life
  ) |>
  knitr::kable(digits = 3,
               caption = "Plasma NCA on the typical-value replication (n = 1).")
```

# Assumptions and deviations

- **On-disk source is the simulation control stream, not the full paper text.**
  The five-page PDF on file contains only the `$PROBLEM SMA MONKEY`
  simulation model (`$SIMULATION (12345678) ONLYSIM SUBPROBLEM=1000`)
  with every `$THETA`, `$OMEGA`, and `$SIGMA` value carrying the
  `FIX` flag. All packaged values are traced verbatim to the on-disk
  block. Baseline demographics for the fitted monkey cohort, the
  paper's narrative Methods, and any published NCA comparison values
  live in Biliouris 2018 main text (CPT PSP 2018;7(9):581-592) and
  are not reproduced here.
- **`V_CSF` piecewise age function is a paediatric-human anatomical
  overlay.** The paper's `$PK` block replaces the fitted V1 = 150 mL
  with anatomically scaled paediatric values via `IF (AGE.LT.0.25)
  TVV1 = 120` etc. The `WT` scaling term `(WT/2.8)` is explicitly
  commented out for V1 in the control stream, so `V_CSF` is set by
  the age function alone.
- **Reference weight 2.8 kg (newborn / juvenile-scale)** is coded
  verbatim from the `$PK` block's `(WT/2.8)` factor. Applied linearly
  to V2..V9 and as `(WT/2.8)^(-0.08)` to every rate constant, per the
  simulation-file description ("scale exp -0.08").
- **Twelve exponential IIVs, twelve diagonal `$OMEGA` variances.** The
  paper's `$PK` block explicitly comments out ETA on K31, K25, K52,
  V3, V4, V6, K61, K47, K74, V8, V9, so those individual parameters
  are typical-value only. The `etalk_csf_brain ~ fixed(13.6)` variance
  is unusually large (log-scale variance 13.6 implies a CV so large
  that stochastic VPCs on K14 will explode) and is retained verbatim
  from the paper's `$OMEGA` block; users running stochastic VPCs
  should be aware.
- **Residual error interpretation.** The paper's `$ERROR` block is
  `W = SQRT(THETA(k)^2 * IPRED^2)` and `Y = IPRED + W * ERR(1)` with
  `$SIGMA 1 FIX`, so `W = |THETA(k)| * |IPRED|` and `THETA(k)` is the
  proportional-SD fraction on the linear (untransformed) scale. That
  is what `propSd_<output>` records.
- **Bioavailability scaling.** The dose is in mg but volumes are in
  mL, so `f(csf) <- 1e6` converts each dose to ng in the state and
  the observations come out directly in ng/mL. No paper-side conversion
  is implied.
- **Compartment naming.** Cervical, lumbar, and thoracic spinal-cord
  compartments and the pons compartment are declared as
  `paper_specific_compartments` (they do not have canonical entries in
  the `brain_<region>` / `spinal_cord_<segment>` families as of this
  writing). Plasma is mapped to the canonical `central`, matching the
  Luu 2017 nusinersen model.

# Related models

* `Luu_2017_nusinersen` -- simpler 4-compartment nusinersen popPK in
  paediatric SMA patients, plasma + CSF observations only (Luu 2017,
  J Clin Pharmacol; `modellib("Luu_2017_nusinersen")`).
