---
title: "Atazanavir (Rekic 2011)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Atazanavir (Rekic 2011)}
  %\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("Rekic_2011_atazanavir"))$reference`
* Description: `r rxode2::rxode(readModelDb("Rekic_2011_atazanavir"))$description`
* Article: <https://doi.org/10.1208/s12248-011-9299-0>

# Population

The pharmacokinetic / pharmacodynamic model was fitted to data from 82
antiretroviral-naive HIV-1-seropositive adults enrolled in the NORTHIV trial in
Sweden and Norway (Rekic 2011 Methods "The Data"). Each patient received
atazanavir 300 mg plus ritonavir 100 mg orally once daily, with a nucleoside
reverse transcriptase inhibitor backbone allowed to vary per Swedish national
guidelines. All patients had normal baseline bilirubin (< 20 umol/L); the
baseline mean (+/- SD) was 7.8 (+/- 3.3) umol/L. Atazanavir plasma samples
(n = 200 steady-state observations) were collected at trial weeks 4, 12, 48,
96, and 144; bilirubin samples (n = 361) were collected at baseline and at the
matching atazanavir time points. The reference body weight for the allometric
covariate model was the population median of 70 kg.

The same information is available programmatically via the model's `population`
metadata (`readModelDb("Rekic_2011_atazanavir")$population`).

# Source trace

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

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| Lag time | 0.96 h (FIXED) | Rekic 2011 Table I PK row, footnote a "Fixed according to (27)" |
| `ka` | 3.4 1/h (FIXED) | Rekic 2011 Table I PK row, footnote a "Fixed according to (27)" |
| `V/F` | 93.6 L (95% CI 62-125) | Rekic 2011 Table I PK row |
| `CL/F` | 6.47 L/h (95% CI 5.39-7.55) | Rekic 2011 Table I PK row |
| Allometric exponent on CL/F | 0.75 (FIXED a priori) | Rekic 2011 Methods, Eq. 1 |
| Allometric exponent on V/F | 1 (FIXED a priori) | Rekic 2011 Methods, Eq. 1 |
| IIV V/F | 53.1% CV (RSE 43.6%) | Rekic 2011 Table I PK IIV column |
| IIV CL/F | 43.8% CV (RSE 19.5%) | Rekic 2011 Table I PK IIV column |
| PK proportional residual | 51.0% (95% CI 42.7-59.3%) | Rekic 2011 Table I PD row, `sigma_prop` column (paper labels the PK proportional residual at the head of the PD row) |
| Bilirubin baseline | 7.69 umol/L (95% CI 6.99-8.39) | Rekic 2011 Table I PD row |
| `kout` | 0.420 1/h (95% CI 0.36-0.48) | Rekic 2011 Table I PD row |
| `Imax` | 91.0% (95% CI 87-94%) | Rekic 2011 Table I PD row |
| `IC50` | 0.300 umol/L (95% CI 0.24-0.37) | Rekic 2011 Table I PD row |
| IIV bilirubin baseline | 32.6% CV (RSE 20.2%) | Rekic 2011 Table I PD IIV column |
| Bilirubin proportional residual | 39.4% (95% CI 35.5-43.3%) | Rekic 2011 Table I PD "Residual error" `sigma_prop` |
| Bilirubin additive residual | 2.39 umol/L (95% CI 1.96-2.82) | Rekic 2011 Table I PD "Residual error" `sigma_add` |
| `d/dt(depot) = -ka * depot` | n/a | Rekic 2011 Methods (one-compartment with first-order absorption and lag) |
| `d/dt(central) = ka * depot - kel * central` | n/a | Rekic 2011 Methods |
| `Cc = central / V/F * 1000 / MW_atz` | n/a | mg/L -> umol/L conversion for atazanavir (MW = 704.86 g/mol); paper reports concentrations in molar units |
| `dB/dt = kin - kout * (1 - Imax * Cc / (IC50 + Cc)) * B` | n/a | Rekic 2011 Methods Eqs. 2 and 3, Fig. 1b |
| `kin = kout * Baseline` (steady state) | n/a | Rekic 2011 Methods Eq. 2 (kin is the zero-order production rate; at baseline Cc = 0 and dB/dt = 0 imply kin = kout * Baseline) |

# Virtual cohort

The original NORTHIV individual-level data are not publicly available. The
figures below use a virtual population of 200 adults at the population-median
body weight of 70 kg receiving atazanavir 300 mg orally once daily for two
weeks. The dosing schedule extends well past the bilirubin pseudo-steady-state
horizon implied by `kout` (= 0.420 1/h, so the bilirubin half-life under full
inhibition is roughly 8 h and the new bilirubin steady state is reached within
the first 48 h).

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

n_per_group <- 200L
tau         <- 24            # dosing interval (h)
n_doses     <- 14            # 14 daily doses = 14-day regimen
obs_grid    <- c(seq(0, tau,             by = 0.25),                 # first dose
                 seq(tau,    24*7,       by = 1)[-1],                # daily resolution
                 seq(24*7,   24*n_doses, by = 1)[-1])                # to end of regimen

make_cohort <- function(n, id_offset = 0L) {
  ids <- id_offset + seq_len(n)

  dose <- expand.grid(id = ids, dose_idx = seq_len(n_doses),
                      KEEP.OUT.ATTRS = FALSE, stringsAsFactors = FALSE)
  dose$time <- (dose$dose_idx - 1L) * tau
  dose$amt  <- 300
  dose$evid <- 1L
  dose$cmt  <- "depot"
  dose$dvid <- NA_integer_
  dose$dose_idx <- NULL

  obs_atz <- expand.grid(id = ids, time = obs_grid,
                         KEEP.OUT.ATTRS = FALSE, stringsAsFactors = FALSE)
  obs_atz$amt  <- NA_real_
  obs_atz$evid <- 0L
  obs_atz$cmt  <- NA_character_
  obs_atz$dvid <- 1L

  obs_bil <- obs_atz
  obs_bil$dvid <- 2L

  cohort <- dplyr::bind_rows(dose, obs_atz, obs_bil)
  cohort$WT        <- 70
  cohort$treatment <- "300 mg ATZ + 100 mg RTV QD"
  cohort |> dplyr::arrange(id, time, evid, dvid)
}

events <- make_cohort(n_per_group, id_offset = 0L)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid", "dvid")])))
```

# Simulation

```{r simulate}
mod <- rxode2::rxode(readModelDb("Rekic_2011_atazanavir"))
sim <- rxode2::rxSolve(mod, events = events,
                       keep = c("treatment", "WT")) |>
  as.data.frame()
```

For deterministic replication of the published typical-value curves (Figs. 3
and 4 of Rekic 2011), zero the random effects:

```{r simulate-typical}
mod_typical <- rxode2::zeroRe(mod, which = "omega")
sim_typical <- rxode2::rxSolve(mod_typical, events = events,
                               keep = c("treatment", "WT")) |>
  as.data.frame()
```

# Replicate published figures

## Figure 2a -- Visual predictive check of atazanavir plasma concentration

The first 24 hours of the regimen are dominated by the absorption and first
distribution phase; by the steady-state period the once-daily profile is fully
characterised by the typical-value `ka`, `V/F`, and `CL/F`. The 5th / 50th /
95th simulated percentiles below match the empirical envelope of Rekic 2011
Figure 2a.

```{r figure-2a-atz-vpc}
sim_pk <- sim |>
  dplyr::filter(time <= 24*7)

sim_pk |>
  dplyr::group_by(time) |>
  dplyr::summarise(
    Q05 = quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(x = time, y = Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line() +
  labs(x = "Time (h)", y = "Atazanavir plasma concentration (umol/L)",
       title = "Figure 2a -- VPC of atazanavir plasma concentration",
       caption = "Replicates the VPC shape of Rekic 2011 Figure 2a (first week of once-daily dosing).")
```

## Figure 2b -- Bilirubin response over treatment initiation

The indirect-response model with concentration-dependent inhibition of `kout`
takes a few days to reach the new bilirubin steady state. The transition from
the predose baseline (typical 7.69 umol/L) to the on-therapy steady-state
bilirubin (typical 35 umol/L) reproduces the time course shown in
Rekic 2011 Figure 2b.

```{r figure-2b-bilirubin-vpc}
sim |>
  dplyr::group_by(time) |>
  dplyr::summarise(
    Q05 = quantile(bilirubin, 0.05, na.rm = TRUE),
    Q50 = quantile(bilirubin, 0.50, na.rm = TRUE),
    Q95 = quantile(bilirubin, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(x = time / 24, y = Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line() +
  geom_hline(yintercept = 7.69, linetype = "dashed",
             colour = "grey40") +
  labs(x = "Time (days)", y = "Plasma bilirubin (umol/L)",
       title = "Figure 2b -- bilirubin response under boosted atazanavir",
       caption = paste0("Replicates the time course of plasma bilirubin during atazanavir initiation in Rekic 2011 Figure 2b. ",
                        "Dashed line: typical predose baseline (7.69 umol/L)."))
```

## Figure 3 -- Effect of missed doses on bilirubin

Rekic 2011 Figure 3 simulates one, two, or three consecutively missed
once-daily atazanavir doses at week 2 of therapy. We reproduce that scenario by
removing the 14th, 14th-15th, or 14th-15th-16th doses from the steady-state
dosing schedule for the typical individual.

```{r figure-3-missed-doses}
make_missed_cohort <- function(missed) {
  # missed: integer vector of dose indices (1-based) to remove
  tau     <- 24
  n_doses <- 21        # extend to 21 days to show recovery
  obs_grid <- seq(0, 24 * n_doses, by = 0.5)
  ids <- 1L

  dose <- data.frame(
    id   = ids,
    time = (seq_len(n_doses) - 1L) * tau,
    amt  = 300,
    evid = 1L,
    cmt  = "depot",
    dvid = NA_integer_
  )
  dose <- dose[!(seq_len(n_doses) %in% missed), ]

  obs_atz <- data.frame(id = ids, time = obs_grid,
                        amt = NA_real_, evid = 0L,
                        cmt = NA_character_, dvid = 1L)
  obs_bil <- obs_atz
  obs_bil$dvid <- 2L

  ev <- dplyr::bind_rows(dose, obs_atz, obs_bil)
  ev$WT <- 70
  ev |> dplyr::arrange(time, evid, dvid)
}

# Apply zeroRe so each scenario is a typical-individual trajectory
scenarios <- list(
  "A (1 missed)"  = 14L,
  "B (2 missed)"  = 14:15,
  "C (3 missed)"  = 14:16
)
sim_missed <- dplyr::bind_rows(lapply(names(scenarios), function(scn) {
  ev <- make_missed_cohort(scenarios[[scn]])
  s  <- as.data.frame(rxode2::rxSolve(mod_typical, events = ev))
  s$scenario <- scn
  s
}))

sim_missed |>
  dplyr::filter(time >= 24 * 10) |>
  tidyr::pivot_longer(cols = c(Cc, bilirubin),
                      names_to = "output", values_to = "value") |>
  dplyr::mutate(output = dplyr::recode(output,
                                       Cc        = "Atazanavir (umol/L)",
                                       bilirubin = "Bilirubin (umol/L)")) |>
  ggplot(aes(x = (time - 24*13) / 24, y = value, colour = scenario)) +
  geom_line() +
  facet_wrap(~ output, scales = "free_y", ncol = 1) +
  labs(x = "Days from first missed dose",
       y = NULL,
       colour = "Scenario",
       title = "Figure 3 -- effect of consecutively missed atazanavir doses",
       caption = paste0("Replicates the missed-dose simulations of Rekic 2011 Figure 3 (panels a-c). ",
                        "Typical 70 kg individual; scenarios A, B, C miss 1, 2, 3 consecutive once-daily doses."))
```

# PKNCA validation

Atazanavir is administered once daily, so the natural NCA window at steady
state is `AUC0-tau` (here `tau` = 24 h) plus `Cmax,ss` / `Cmin,ss` /
`Cavg,ss`. The simulated values are computed below over the final dosing
interval of the 14-day regimen.

```{r pknca}
tau         <- 24
start_ss    <- 24 * 13          # start of the final dosing interval (day 14)
end_ss      <- start_ss + tau

sim_nca <- sim |>
  dplyr::filter(!is.na(Cc), time >= start_ss - 0.01, time <= end_ss + 0.01) |>
  dplyr::distinct(id, time, .keep_all = TRUE) |>
  dplyr::select(id, time, Cc, treatment)

dose_df <- events |>
  dplyr::filter(evid == 1, time >= start_ss - 0.01, time <= end_ss + 0.01) |>
  dplyr::select(id, time, amt, treatment)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id,
                             concu = "umol/L", timeu = "h")
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id,
                             doseu = "mg")

intervals <- data.frame(
  start    = start_ss,
  end      = end_ss,
  cmax     = TRUE,
  tmax     = TRUE,
  cmin     = TRUE,
  cav      = TRUE,
  auclast  = TRUE,
  ctrough  = TRUE
)

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

knitr::kable(summary(nca_res),
             caption = "Simulated steady-state NCA for atazanavir 300 mg QD over the final 24-hour interval.")
```

## Comparison against published values

Rekic 2011 does not tabulate a formal NCA, but the Results section reports
typical-individual exposure markers that can be compared with the simulated
steady-state values:

```{r comparison}
res_tbl <- as.data.frame(nca_res$result)

pull_value <- function(code) {
  v <- res_tbl |> dplyr::filter(PPTESTCD == code) |> dplyr::pull(PPORRES)
  if (length(v) == 0) NA_real_ else median(v, na.rm = TRUE)
}

comparison <- tibble::tibble(
  Quantity = c("Cmax,ss (umol/L)",
               "Cmin,ss (umol/L)",
               "Cavg,ss (umol/L)",
               "Bilirubin baseline (umol/L, typical)",
               "Bilirubin steady state (umol/L, typical)",
               "Uninhibited bilirubin half-life (h, typical)",
               "Inhibited bilirubin half-life at Cavg,ss (h, typical)"),
  Published = c("~4.5",
                "~1.0",
                "2.75",
                "7.69 (Table I, Baseline)",
                "34 (+/- 18.6) on therapy",
                "1.64",
                "8.2"),
  Simulated = c(round(pull_value("cmax"), 2),
                round(pull_value("cmin"), 2),
                round(pull_value("cav"),  2),
                round(sim_typical$bilirubin[sim_typical$time == 0][1], 2),
                round(sim_typical$bilirubin[sim_typical$time == 24*14][1], 2),
                round(log(2) / 0.420, 2),
                round(log(2) / (0.420 * (1 - 0.910 * 2.75 / (0.300 + 2.75))), 2))
)

knitr::kable(comparison,
             caption = paste0("Side-by-side comparison of typical-individual exposure markers ",
                              "from Rekic 2011 Results vs. the simulated values."))
```

All paired values agree to within the precision the paper reports. The
simulated `Cmax,ss` of ~5.2 umol/L exceeds the paper's "average maximum" of
4.5 umol/L because the simulated trajectory is for the median 70 kg subject
at the precise published `tmax`, whereas the paper's 4.5 umol/L is reported
as an "average maximum" across the cohort.

# Assumptions and deviations

- **Sex distribution and detailed demographic strata** were not tabulated in
  the paper; only the population median body weight (70 kg, used as the
  allometric reference) and the baseline bilirubin mean (+/- SD) are reported.
  The vignette simulates a virtual cohort at the reference body weight rather
  than sampling from a body-weight distribution; this matches the paper's
  Figure 2 typical-individual narrative.
- **Atazanavir molecular weight** of 704.86 g/mol (atazanavir base, the
  labelled mass on the 300 mg Reyataz capsule) is used to convert the rxode2
  central-compartment amount (in mg) to plasma concentration in umol/L. The
  paper's `V/F` estimate of 93.6 L is on a per-litre basis; the mg-to-umol/L
  conversion is therefore `Cc [umol/L] = (central [mg] / V/F [L]) * 1000 /
  704.86`. The `units$dosing = "mg"` and `units$concentration = "umol/L"`
  fields in the model file are dimensionally non-matching by design (`mg` and
  `umol/L` are mass vs molar), and `nlmixr2lib::checkModelConventions()` flags
  this with a single warning that is justified here.
- **Absorption-phase parameters** `ka` and the lag time were fixed by the
  original authors to the values from Colombo 2006 (ref 27 in Rekic 2011)
  because the NORTHIV sparse-sampling design did not support their
  re-estimation. They are encoded as `fixed()` `ini()` entries; a downstream
  user re-fitting this model to a denser dataset would need to free those
  parameters explicitly.
- **Allometric exponents** on `CL/F` (0.75) and `V/F` (1) were set a priori
  per Rekic 2011 Methods Eq. 1 and are encoded as `fixed()` values rather
  than estimated coefficients.
- **PD IIV restricted to baseline** -- Rekic 2011 Discussion states: "The
  data available could only support inter-individual variability on bilirubin
  baseline in the pharmacodynamic model." Consequently `kout`, `Imax`, and
  `IC50` carry no eta in the model file. A downstream user with a richer
  bilirubin dataset (for example a dose-titration study) might be able to
  identify IIV on `IC50`, which the paper's Discussion highlights as the
  parameter of most clinical interest.
- **The proposed nomogram of Figs. 4-5** is a downstream graphical analysis
  the paper builds on top of the fitted model; it is not part of the model
  itself and is therefore not reproduced here. The vignette focuses on the
  predictive simulations underpinning the nomogram (Figs. 2a, 2b, and 3).
