Skip to contents

Atazanavir (Rekic 2011)

Source: vignettes/articles/Rekic_2011_atazanavir.Rmd
Rekic_2011_atazanavir.Rmd

Model and source

  • Citation: Rekic D, Clewe O, Roshammar D, Flamholc L, Sonnerborg A, Ormaasen V, Gisslen M, Abelo A, Ashton M. Bilirubin – a potential marker of drug exposure in atazanavir-based antiretroviral therapy. The AAPS Journal. 2011;13(4):598-605. doi:10.1208/s12248-011-9299-0. Absorption-phase parameters (ka, lag time) fixed from Colombo S, Buclin T, Cavassini M, et al. Population pharmacokinetics of atazanavir in patients with human immunodeficiency virus infection. Antimicrob Agents Chemother. 2006;50(11):3801-3808 (cited as ref 27).
  • Description: Population PK / PD model for atazanavir (boosted with ritonavir 100 mg QD) and its concentration-dependent effect on plasma bilirubin in adult antiretroviral-naive HIV-positive patients from the NORTHIV trial (Rekic 2011). Atazanavir disposition is described by a one-compartment model with first-order absorption and an absorption lag, fitted to log-transformed plasma atazanavir concentrations; ka and the lag time were fixed to the published values from the Colombo 2006 atazanavir popPK report (ref 27) because sparse absorption-phase sampling did not support their re-estimation, and CL/F and V/F were re-estimated under fixed allometric scaling on body weight centred at 70 kg (exponents 0.75 on CL/F and 1 on V/F, both fixed a priori). The bilirubin response is described by an indirect-response (turnover) model with concentration-dependent inhibition of the fractional turnover rate kout: dB/dt = kin - kout * (1 - Imax * Cc / (IC50 + Cc)) * B, with kin re-parameterised at steady state as kin = kout * Baseline. Inter-individual variability is supported only on V/F, CL/F (PK), and bilirubin baseline (PD); the paper notes that the data did not support IIV on the remaining PD parameters. PK residual variability is proportional; bilirubin residual variability is combined additive + proportional (the paper’s ‘slope-intercept’ model).
  • 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).

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 

mod <- rxode2::rxode(readModelDb("Rekic_2011_atazanavir"))
#> ℹ parameter labels from comments will be replaced by 'label()'
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:

mod_typical <- rxode2::zeroRe(mod, which = "omega")
sim_typical <- rxode2::rxSolve(mod_typical, events = events,
                               keep = c("treatment", "WT")) |>
  as.data.frame()
#> ℹ omega/sigma items treated as zero: 'etalvc', 'etalcl', 'etalrbase_bilirubin'
#> Warning: multi-subject simulation without without 'omega'

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.

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.

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.

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
}))
#> ℹ omega/sigma items treated as zero: 'etalvc', 'etalcl', 'etalrbase_bilirubin'
#> ℹ omega/sigma items treated as zero: 'etalvc', 'etalcl', 'etalrbase_bilirubin'
#> ℹ omega/sigma items treated as zero: 'etalvc', 'etalcl', 'etalrbase_bilirubin'

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.

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.")
Simulated steady-state NCA for atazanavir 300 mg QD over the final 24-hour interval.
Interval Start Interval End treatment N AUClast (h*umol/L) Cmax (umol/L) Cmin (umol/L) Tmax (h) Cav (umol/L) Ctrough (umol/L)
312 336 300 mg ATZ + 100 mg RTV QD 200 69.7 [43.2] 5.69 [32.8] 1.05 [166] 2.00 [2.00, 2.00] 2.90 [43.2] NC

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:

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."))
Side-by-side comparison of typical-individual exposure markers from Rekic 2011 Results vs. the simulated values.
Quantity Published Simulated
Cmax,ss (umol/L) ~4.5 5.62
Cmin,ss (umol/L) ~1.0 1.37
Cavg,ss (umol/L) 2.75 2.86
Bilirubin baseline (umol/L, typical) 7.69 (Table I, Baseline) 7.69
Bilirubin steady state (umol/L, typical) 34 (+/- 18.6) on therapy 35.77
Uninhibited bilirubin half-life (h, typical) 1.64 1.65
Inhibited bilirubin half-life at Cavg,ss (h, typical) 8.2 9.19

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).