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

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

# Population and biological context

Coproporphyrin I (CPI) is a heme-biosynthesis byproduct produced predominantly in erythrocytes, with ALA-synthase considered the rate-limiting step of its formation (Barnett 2018 Introduction). CPI is a selective in-vitro substrate of OATP1B1 and OATP1B3 with no involvement of OATP2B1 or renal uptake transporters in the available data. Its plasma concentrations are reported to be stable between occasions and have low interindividual variability (< 25% CV) in subjects with SLCO1B1 c.521 T>C wildtype, making it a candidate endogenous biomarker for OATP1B-mediated drug-drug interactions.

Barnett 2018 developed the semi-mechanistic CPI model jointly with the rifampicin and rosuvastatin popPK models in 12 healthy adult male volunteers (Lai et al. 2016 cohort, SLCO1B1 wildtype). The cohort underwent three study occasions (rifampicin-only on OCC1, rosuvastatin-only on OCC2, combined rifampicin + rosuvastatin on OCC3) with a 7-day washout between occasions; baseline CPI plasma and urine samples were collected on every occasion in addition to post-dose monitoring. The 420 CPI plasma samples (OCC1 = 144, OCC2 = 144, OCC3 = 132) plus 102 CPI urine samples were fit simultaneously in NONMEM under FOCE. Identifiability of the structural model was confirmed in DAISY.

The same context is available programmatically via `readModelDb("Barnett_2018_coproporphyrin_I")$population`.

# Source trace

| Equation / parameter | Value         | Source location                                       |
|----------------------|---------------|-------------------------------------------------------|
| `lksyn`              | log(12.7)     | Barnett 2018 Table 1, CPI row 'k syn (nM/h)' = 12.7 (SE 6%) -- but see Errata below on the 'nM/h' label |
| `lclb`               | log(12.3)     | Barnett 2018 Table 1, CPI row 'CL b,CPI (L/h)' = 12.3 |
| `lclr`               | log(1.64)     | Barnett 2018 Table 1, CPI row 'CL R,CPI (L/h)' = 1.64 |
| `lvc` (baseline)     | log(6.59)     | Barnett 2018 Table 1, CPI row 'V CPI (L)' = 6.59      |
| `e_rif_vc`           | -0.4841       | Barnett 2018 Table 1, CPI row 'V CPI (L) (RIF)' = 3.4; derived as 3.4 / 6.59 - 1 |
| `lki`                | log(1.15)     | Barnett 2018 Table 1, CPI row 'Ki CPI (uM) c' = 1.15 (footnote c: unbound Ki = 0.13 uM after RIF fu = 0.11) |
| `propSd` (plasma)    | 0.139         | Barnett 2018 Table 1, CPI row 'r prop (%) - plasma' = 13.9 |
| `addSd` (plasma)     | 0.001 nmol/L FIXED | Barnett 2018 Table 1, CPI row 'r add (nM) - plasma' = 0.001 FIXED |
| `propSd_Ucpi`        | 0.342         | Barnett 2018 Table 1, CPI row 'r prop (%) - urine' = 34.2 |
| `addSd_Ucpi`         | 2.69 nmol/L   | Barnett 2018 Table 1, CPI row 'r add (nM) - urine' = 2.69 |
| ODE form             | n/a           | Barnett 2018 Eq. 3 (baseline) + Eq. 4 (RIF inhibition) + Eq. 5 (urine accumulation) |
| Steady-state baseline (analytic) | `ksyn / (CLb + CLr)` | Derived from Eq. 3 dC/dt = 0; equals 12.7 / (12.3 + 1.64) = 0.911 nmol/L (= nM) |

## Units of every ODE term (dimensional analysis)

| Term in `d/dt(central) = ksyn - (CLb_eff + CLr) * Ccpi`              | Units                              |
|---------------------------------------------------------------------|------------------------------------|
| `central` (state)                                                   | nmol                               |
| `ksyn`                                                              | nmol / h                           |
| `CLb_eff`, `CLr`                                                    | L / h                              |
| `Ccpi = central / vc`                                               | nmol / L                           |
| `(CLb_eff + CLr) * Ccpi`                                            | L / h * nmol / L = nmol / h        |
| `d/dt(central)`                                                     | nmol / h ok                        |
| `d/dt(urine) = CLr * Ccpi`                                          | nmol / h ok                        |

# Steady-state check (no RIF, deterministic typical-value)

This verifies the analytic baseline `Css = ksyn / (CLb + CLr)` and the absence of spurious drift in the ODE solver.

```{r ss-check}
mod <- readModelDb("Barnett_2018_coproporphyrin_I")
mod_typical <- rxode2::zeroRe(mod)

make_cpi_events <- function(t_end = 120, dt = 2, conmed_rif = 0, crif = 0, occ = 1) {
  data.frame(
    id   = 1L,
    time = seq(0, t_end, by = dt),
    evid = 0L,
    amt  = 0,
    cmt  = "Cc",
    OCC  = occ,
    CONMED_RIF = conmed_rif,
    CP_RIF_UM  = crif
  )
}

ss_sim <- rxode2::rxSolve(mod_typical, events = make_cpi_events(t_end = 200))
cat("CPI typical-value baseline (no RIF):\n")
cat("  Ccpi(t=0)  :", round(ss_sim$Cc[1], 4), "nmol/L\n")
cat("  Ccpi(t=200):", round(tail(ss_sim$Cc, 1), 4), "nmol/L\n")
cat("  Drift over 200 h:", signif(diff(range(ss_sim$Cc)), 3), "nmol/L\n")
cat("  Analytic Css (ksyn / (CLb + CLr) = 12.7 / 13.94):",
    round(12.7 / (12.3 + 1.64), 4), "nmol/L\n")
stopifnot(diff(range(ss_sim$Cc)) < 1e-6)
```

The simulator holds the analytic baseline indefinitely within numerical tolerance.

# Perturbation-recovery (no RIF, displaced initial condition)

Displacing the central state to 0.5x and 2x the steady-state amount should give a monotone recovery to the baseline.

```{r perturb}
ev <- make_cpi_events(t_end = 100, dt = 0.5)
ic_amount_baseline <- 0.911 * 6.59  # Css * Vcpi at baseline

sim_low  <- rxode2::rxSolve(mod_typical, events = ev,
                            inits = c(central = 0.5 * ic_amount_baseline, urine = 0))
sim_high <- rxode2::rxSolve(mod_typical, events = ev,
                            inits = c(central = 2.0 * ic_amount_baseline, urine = 0))

cat("Perturbation recovery toward Css = 0.911 nmol/L:\n")
cat("  Start  0.5x:", round(sim_low$Cc[1], 4),
    "  End:", round(tail(sim_low$Cc, 1), 4), "nmol/L\n")
cat("  Start  2.0x:", round(sim_high$Cc[1], 4),
    "  End:", round(tail(sim_high$Cc, 1), 4), "nmol/L\n")
```

```{r perturb-plot}
recovery <- dplyr::bind_rows(
  sim_low  |> as.data.frame() |> dplyr::mutate(start = "0.5 x Css"),
  sim_high |> as.data.frame() |> dplyr::mutate(start = "2.0 x Css")
)
ggplot(recovery, aes(time, Cc, colour = start)) +
  geom_line(linewidth = 1) +
  geom_hline(yintercept = 0.911, linetype = "dashed", alpha = 0.6) +
  labs(x = "Time (h)", y = "CPI plasma Ccpi (nmol/L)",
       title = "Perturbation recovery toward steady-state baseline",
       colour = "Initial state",
       caption = "Both trajectories monotonically approach Css = 0.911 nmol/L.") +
  theme_minimal()
```

# Rifampicin perturbation (reproduces Figure 3 upper-panel pattern)

Barnett 2018 Figure 3 shows the VPC of CPI plasma data superimposed under the rifampicin phase (CPI (RIF), upper-right panel). The clinical observation is that a single 600 mg oral rifampicin dose raises CPI plasma 2-5x over the next ~8 h, with recovery toward baseline by 24 h. Reproduce the typical-value trajectory by simulating the rifampicin model first and feeding its central-compartment concentration as `CP_RIF_UM` on the CPI event table.

```{r rif-perturb}
# Step 1: simulate the rifampicin typical-value PK over 24 h
mod_rif <- rxode2::zeroRe(readModelDb("Barnett_2018_rifampicin"))
ev_rif <- rxode2::et(amt = 600, cmt = "depot", time = 0) |>
  rxode2::et(time = seq(0.05, 24, by = 0.05), cmt = "Cc") |>
  as.data.frame()
ev_rif$OCC <- 1
sim_rif <- rxode2::rxSolve(mod_rif, events = ev_rif) |> as.data.frame()
crif_fn <- approxfun(sim_rif$time, sim_rif$Cc, rule = 2, yleft = 0)

# Step 2: simulate CPI under RIF perturbation; the perturbation window
# covers t = 0 to 24 h with CONMED_RIF = 1 throughout
ev_cpi_rif <- data.frame(
  id = 1L,
  time = seq(0, 24, by = 0.1),
  evid = 0L,
  amt = 0,
  cmt = "Cc",
  OCC = 1,
  CONMED_RIF = 1
)
ev_cpi_rif$CP_RIF_UM <- crif_fn(ev_cpi_rif$time)
sim_cpi_rif <- rxode2::rxSolve(mod_typical, events = ev_cpi_rif) |> as.data.frame()

ev_cpi_base <- ev_cpi_rif
ev_cpi_base$CONMED_RIF <- 0
ev_cpi_base$CP_RIF_UM <- 0
sim_cpi_base <- rxode2::rxSolve(mod_typical, events = ev_cpi_base) |> as.data.frame()

cpi_compare <- dplyr::bind_rows(
  sim_cpi_base |> dplyr::transmute(time, Ccpi = Cc, phase = "baseline"),
  sim_cpi_rif  |> dplyr::transmute(time, Ccpi = Cc, phase = "rifampicin 600 mg PO")
)
ggplot(cpi_compare, aes(time, Ccpi, colour = phase)) +
  geom_line(linewidth = 1) +
  labs(x = "Time (h)", y = "CPI plasma Ccpi (nmol/L)",
       title = "Typical-value CPI trajectory under a single 600 mg RIF dose",
       colour = NULL,
       caption = "Compare to the CPI (RIF) panel of Barnett 2018 Figure 3 (upper right).") +
  theme_minimal()

auc_baseline <- with(sim_cpi_base, sum(diff(time) * (head(Cc, -1) + tail(Cc, -1)) / 2))
auc_rif      <- with(sim_cpi_rif,  sum(diff(time) * (head(Cc, -1) + tail(Cc, -1)) / 2))
aucr <- auc_rif / auc_baseline
cat("Typical-value CPI AUC ratio (0-24 h) under RIF vs baseline:", round(aucr, 2), "\n")
cat("Barnett 2018 Power-calculations text reports mean simulated CPI AUCR for the\n")
cat("I/Ki = 1 (rifampicin) scenario as 3.46. Sim AUCR =", round(aucr, 2),
    " (typical-value, no IIV).\n")
```

# Mass-balance check at the analytic baseline

At steady state with no RIF, the production rate (`ksyn = 12.7 nmol/h`) must exactly balance the elimination rate (`(CLb + CLr) * Css = (12.3 + 1.64) * 0.911 = 12.7 nmol/h`):

```{r mass-balance}
ksyn       <- 12.7
clb        <- 12.3
clr        <- 1.64
css        <- ksyn / (clb + clr)
elim_rate  <- (clb + clr) * css
cat("Production rate  :", ksyn, "nmol/h\n")
cat("Elimination rate :", round(elim_rate, 6), "nmol/h\n")
stopifnot(abs(ksyn - elim_rate) < 1e-9)

# 24-hour cumulative urine excretion = CLr * Css * 24
ucpi_24h_analytic <- clr * css * 24
ucpi_24h_sim <- tail(sim_cpi_base$urine, 1)
cat("Analytic 24-h cumulative urinary CPI:", round(ucpi_24h_analytic, 3), "nmol\n")
cat("Simulated 24-h cumulative urinary CPI:", round(ucpi_24h_sim, 3),     "nmol\n")
```

# Comparison against published values

Barnett 2018 reports the following CPI characteristics that the typical-value model should approximately match:

| Quantity                                  | Barnett 2018 reported value      | Simulated typical-value |
|-------------------------------------------|----------------------------------|-------------------------|
| Baseline plasma CPI (Css)                 | ~0.5-2 nmol/L typical CPI range  | `r round(12.7 / (12.3 + 1.64), 3)` nmol/L |
| Mean CPI AUCR under RIF (I/Ki = 1)         | 3.46 (Power calculations text)   | `r round(auc_rif / auc_baseline, 2)`             |
| Biliary fraction of CPI elimination        | ~88% (Table 1 narrative)         | `r round(100 * 12.3 / (12.3 + 1.64), 1)`%       |

# Assumptions and deviations

* **`ksyn` units.** Barnett 2018 Table 1 reports `k syn` with label `nM/h`, but dimensional analysis of the steady-state equation `Css = ksyn / (CLb + CLr)` only resolves if `ksyn` has units of `nmol/h` (amount/time), not `nmol/L/h` (concentration/time). The model encodes `ksyn` in `nmol/h`. With `Css = 12.7 / 13.94 = 0.911 nmol/L`, the value matches the expected CPI baseline range in healthy SLCO1B1-wildtype adults and reproduces the paper's reported 88% biliary-excretion fraction.
* **`Vcpi` covariate effect as instantaneous binary shift.** Barnett 2018 reports two values for the CPI volume of distribution (Vcpi = 6.59 L baseline, 3.4 L during the rifampicin phase) and encodes this as a binary covariate at the period level. The model file encodes the effect multiplicatively as `vc <- exp(lvc + etalvc) * (1 + e_rif_vc * CONMED_RIF)` with `e_rif_vc = -0.4841`. The shift is structurally non-physiological (an empirical capture of the OATP1B-inhibition effect on apparent CPI distribution) and is applied instantaneously when CONMED_RIF turns on; in clinical simulations this produces a discontinuous change in Ccpi at the moment CONMED_RIF switches from 0 to 1. The initial condition `central(0) <- Css_baseline * vc` deliberately uses the CONMED_RIF-modified `vc` so that `Ccpi(0) = Css_baseline` (0.911 nmol/L) regardless of the CONMED_RIF value at t = 0 -- the simulation always starts at the pre-perturbation baseline and the subsequent dynamics evolve under the chosen perturbation.
* **`ksyn` IIV is very poorly estimated.** Barnett 2018 Table 1 reports IIV on `ksyn` as 2.9% CV with RSE 318.7% (essentially un-identified). The model encodes the point estimate but the resulting `etalksyn ~ 0.000841` variance is effectively zero. Users running stochastic simulations should treat the ksyn IIV component as nominal.
* **No demographics data.** As with the sibling rifampicin and rosuvastatin models, Barnett 2018 does not tabulate per-subject body-weight / age / region data for the n = 12 cohort.
* **`urine` and `serum` as compartment names.** The model uses `urine` (cumulative renal excretion amount) as a compartment name. The nlmixr2lib canonical compartment register (`R/conventions.R`) does not include `urine`, so `checkModelConventions()` will flag this as a non-canonical compartment name. Precedent: `inst/modeldb/endogenous/Aksenov_2018_uricAcid.R` uses `serum` + `urine` for the same purpose and is the existing example of an endogenous turnover model with simultaneous plasma + cumulative-urine outputs.
