Skip to contents

Danoprevir HCV PK/VK (Canini 2015)

Source: vignettes/articles/Canini_2015_danoprevir.Rmd
Canini_2015_danoprevir.Rmd

Model and source

  • Citation: Canini L, Chatterjee A, Guedj J, Lemenuel-Diot A, Brennan B, Smith PF, Perelson AS. (2015). A pharmacokinetic/viral kinetic model to evaluate the treatment effectiveness of danoprevir against chronic HCV. Antiviral Therapy 20(5):469-477. doi:10.3851/IMP2879.
  • Description: Combined pharmacokinetic / viral-kinetic (PK/VK) model for the HCV NS3/4A protease inhibitor danoprevir in adults with chronic hepatitis C (Canini 2015). PK is a two-compartment disposition (central = plasma, peripheral1 = tissue including liver) with zero-order absorption of duration Tk0 into central after a post-dose lag Tlag, followed by first-order elimination. The viral-kinetic layer follows the standard Neumann 1998 model with target cells assumed to remain constant (only the productively-infected cell pool infected and the free-virion pool virus are carried as ODE states). Drug effectiveness e(t) = C1 / (EC50 + C1) is a Hill-Emax function with Emax fixed at 1 and Hill coefficient fixed at 1 (the Canini 2015 model-selection winner over h = 2 and h = 3 by AICc/BIC); e(t) blocks virion production so d/dt(virus) = (1 - e) * p * infected - c * virus. Because Canini 2015 does not estimate the virion production rate p or the infection rate beta separately, the packaged model fixes p = 1 as a mathematical scaling (the identifiability constraint p * beta_T = d * c at the pre-treatment steady state fully determines beta_T = d * c and I(0) = c * V0). All 10 estimated typical values from Canini 2015 Table 2 (Tlag, Tk0, V1, ke, k12, k21, V0, c, d, EC50) and all 10 IIVs are carried; cohort-specific V0 offsets (Table 2 cohorts 2-5) and the screened HCV genotype covariate are documented in vignette Errata but not applied as covariates in the packaged model (see notes in population and in the vignette Assumptions section). PK is in per-hour units to match the paper’s Table 2 reporting; the viral-kinetic /day-scale rate constants c and d are converted inline to /h so the integrated system runs on a single hour-scale time axis.
  • Article: https://doi.org/10.3851/IMP2879

Canini et al 2015 fit an integrated PK plus viral-kinetic (PK/VK) model to 14 days of oral danoprevir monotherapy in adults with chronic HCV (Phase 1 single ascending dose, 40 randomised patients, 32 active drug across 5 cohorts). PK is a two-compartment disposition (central plasma + peripheral tissue including the liver) with zero-order absorption of duration Tk0 into central after a post-dose lag Tlag, and first-order elimination. The viral-kinetic layer is the standard Neumann 1998 target-cell model with target cells assumed constant, so only the productively-infected cell pool infected and the free-virion pool virus are carried as ODE states. Drug effectiveness follows a Hill-Emax function e(t) = C1 / (EC50 + C1) with Emax fixed at 1 and Hill exponent fixed at 1 (best of h = 1, 2, 3 by AICc and BIC).

Population

The packaged model is fit to the 32 active-drug patients (8 per cohort allocated 8:2 active:placebo, with 5-6 non-responders receiving active drug in cohort 5) from a randomised, double-blind, placebo-controlled, 5-cohort Phase 1 single-ascending-dose study of oral danoprevir monotherapy for 14 days (Canini 2015 Table 1). The cohort structure is:

Cohort Dose regimen Patient type Median initial VL (log10 IU/mL) Median 14-day decline (log10 IU/mL)
1 100 mg twice daily Treatment-naive 5.8 2.0
2 100 mg three times daily Treatment-naive 6.2 2.7
3 200 mg twice daily Treatment-naive 6.3 2.3
4 200 mg three times daily Treatment-naive 6.4 3.9
5 300 mg twice daily PEG-IFN / RBV non-responder 6.5 2.7

Genotype distribution across the 40 patients was 30 percent GT1a, 55 percent GT1b, and 15 percent GT1 with undetermined subtype (Canini 2015 Methods, Patients paragraph). HCV genotype was screened as a covariate but not retained in the final PK/VK model. Age, weight, and sex are not reported in the trimmed source; the Canini 2015 paper cites the underlying Phase 1 design (reference [9]) for baseline demographics.

The same information is available programmatically via the model’s population metadata (rxode2::rxode(readModelDb("Canini_2015_danoprevir"))$meta$population).

Source trace

The per-parameter origin is recorded as an in-file comment next to each ini() entry in inst/modeldb/specificDrugs/Canini_2015_danoprevir.R. The table below collects the source location for every model equation and every parameter.

Equation / parameter Value Source location
Two-compartment PK ODEs n/a Canini 2015 Methods, PK model paragraph + Figure 1A; Table 2 heading.
Zero-order absorption over Tk0 after lag Tlag n/a Canini 2015 Methods, PK model paragraph (‘after a lag-time, Tlag, danoprevir is absorbed into the first or central compartment … following a zero order absorption law, with rate constant k0 = D/Tk0’).
Neumann viral-kinetic ODE system with T constant n/a Canini 2015 Methods, VK model paragraph (‘the kinetics of viral decline were assumed to follow the standard model developed by Neumann et al., with target cells assumed to remain constant’).
Steady-state initial conditions I(0) = c V0 / p, V(0) = V0 n/a Canini 2015 Methods (standard Neumann pre-treatment steady-state derivation; p = 1 fixed as a scaling because Canini 2015 does not identify p separately, see Assumptions section).
Emax effectiveness e = C1^h / (EC50^h + C1^h), Emax = 1 n/a Canini 2015 Methods, Danoprevir pharmacodynamics paragraph (‘drug effectiveness e(t) varies as a function of C1, according to the Emax model shown in Equation 1, where the maximum effectiveness is assumed to be 1’).
Hill coefficient h = 1 (FIXED) 1 Canini 2015 Results, PK/VK model paragraph (‘Fitting the PK/VK model for h=1 (AICc=6097, BIC=6098), h=2 (AICc=6123, BIC=6154) and h=3 (AICc=6154, BIC=6182), showed that h=1 provided the best model’).
Tlag = 0.50 h (IIV 63.5% CV) 0.50 Canini 2015 Table 2.
Tk0 = 0.58 h (IIV 122% CV) 0.58 Canini 2015 Table 2.
V1 = 2100 L (V1/F1 effective; IIV 90% CV) 2100 Canini 2015 Table 2.
ke = 1.02 /h (IIV 30.3% CV) 1.02 Canini 2015 Table 2.
k12 = 0.25 /h (IIV 79.4% CV) 0.25 Canini 2015 Table 2.
k21 = 0.81 /h (IIV 89.2% CV) 0.81 Canini 2015 Table 2.
V0 = 5.4e5 IU/mL (cohort-1 reference; IIV 118% CV) 540000 Canini 2015 Table 2 (cohort-1 row).
c = 5.28 /day (IIV 38.4% CV) 5.28 Canini 2015 Table 2 (converted to /h in model() by dividing by 24).
d = 0.18 /day (IIV 55.9% CV) 0.18 Canini 2015 Table 2 (converted to /h in model() by dividing by 24).
EC50 = 0.0082 ng/mL (IIV 152% CV) 0.0082 Canini 2015 Table 2.
PK additive residual SD a = 0.25 ng/mL 0.25 Canini 2015 Results (‘combined error model … additive error term (a = 0.25 +- 0.012 ng/ml)’).
PK proportional residual SD b = 0.61 0.61 Canini 2015 Results (‘proportional error term (b = 0.61 +- 0.02)’).
VK additive residual SD on log10 = 0.29 0.29 Canini 2015 Results (‘additive error model was found to best describe the residual error for the VK data (a = 0.29 +- 0.0094 log10 IU/ml)’).

Errata

Three points from the source paper are worth calling out because the packaged model does not encode them literally.

  • Abstract vs Methods contradiction on the effectiveness driver. The paper’s abstract states that the effectiveness varies “with the predicted danoprevir concentration inside the SECOND compartment”, but the Methods “Danoprevir pharmacodynamics” paragraph explicitly writes “e(t) varies as a function of C1, according to the Emax model” – i.e. the FIRST-compartment plasma concentration. Canini 2015 does not report V2 anywhere, so a second-compartment concentration cannot be evaluated in ng/mL against the tabulated EC50 = 0.0082 ng/mL; likewise the paper’s Discussion compares EC50 with plasma Cmin / Cmax (‘EC50 estimate, which seems small compared to Cmax, is only 2.6 to 7.7 x higher than the Cmin of the cohort with the highest and lowest dose in treatment-naive patients’). The packaged model follows the Methods statement and drives effectiveness with C1 = central / V1. The abstract wording appears to be a transcription slip.

  • Cohort-specific typical V0 (Table 2 cohorts 2-5). Canini 2015 Table 2 reports different typical baseline viral loads by cohort: 5.4e5 (cohort 1, reference); 1.9e6 (cohort 2, P = 0.08); 1.8e6 (cohort 3, P = 0.048); 2.6e6 (cohort 4, P = 0.012); 2.7e6 (cohort 5, P = 0.0092). Two of the four P-values are above the 0.05 threshold, the sample size is ~8 per cohort, and the between-cohort spread (0.5 log10) is small relative to the reported 118 percent IIV on V0 (~2.4 log10 within +/- 2 SD). The packaged model uses the cohort-1 reference typical value (V0 = 5.4e5) with the paper’s IIV of 118 percent CV; simulated between-subject viral-load ranges therefore span the observed between-cohort values without encoding cohort as a covariate.

  • Neumann parameters p and beta not identifiable from the reported data. Canini 2015 estimates only c, d, V0, EC50 in the VK layer. The classical Neumann model at pre-treatment steady state constrains p x beta_T = d x c but not p and beta_T separately, so one degree of freedom must be fixed by convention. The packaged model follows the Neumann-family convention of fixing p = 1 as a mathematical scaling; the derived quantities are then beta_T = d x c and I(0) = c x V0. p = 1 is not a biological rate; users interested in reading infected on a biological hepatocyte-count scale would need to rescale I and beta_T by their preferred p.

Virtual cohort

The Canini 2015 individual-subject data are not publicly available. This vignette simulates a virtual cohort matching the 5 dosing regimens of Canini 2015 Table 1 with 20 subjects per cohort (100 total) drawn from the model’s between-subject variability distributions. This is well below the 200-per-arm cap for validation vignettes and renders quickly.

set.seed(20260725)

n_per_cohort <- 20L
dose_regimens <- tibble::tribble(
  ~cohort, ~dose_mg, ~n_per_day, ~interval_h,        ~patient_type,
       1L,   100.0,        2L,        12.0,     "treatment-naive",
       2L,   100.0,        3L,         8.0,     "treatment-naive",
       3L,   200.0,        2L,        12.0,     "treatment-naive",
       4L,   200.0,        3L,         8.0,     "treatment-naive",
       5L,   300.0,        2L,        12.0,     "non-responder"
)

sim_dur_h <- 14 * 24  # 14-day monotherapy window (Canini 2015 Methods)

# Helper: build one cohort's event table (dose rows + observation grid).
# id_offset disjoint per cohort so ids never collide when bind_rows'ed.
make_cohort <- function(cohort_id, dose_mg, interval_h, n_subj, id_offset,
                        patient_type, dur_h) {
  ids <- id_offset + seq_len(n_subj)

  dose_times <- seq(0, dur_h - 0.1, by = interval_h)
  # Two observation grids: dense on the PK observable Cc during the first
  # dosing interval to capture Cmax; sparser on the VK observable Vlog10
  # across the 14-day treatment window. rxode2 requires cmt = <observable
  # name> for multi-output models (Cc from the `central` state, Vlog10 from
  # the `virus` state); this is the pattern that the pattern-2 warning in
  # known-vignette-failure-patterns.md permits when the algebraic
  # observables are declared AFTER all ODE states in the model body (as
  # they are here  -  the auto-injected cmt(Cc), cmt(Vlog10) slots sit at
  # positions 5, 6 with no downstream ODE state to renumber).
  obs_times_pk <- sort(unique(c(
    seq(0,   interval_h,     by = 0.25),
    seq(interval_h,   24,    by = 1),
    seq(24,           dur_h, by = 2)
  )))
  obs_times_vk <- sort(unique(c(
    seq(0,  24,    by = 2),
    seq(24, dur_h, by = 8)
  )))

  dose_rows <- expand.grid(id = ids, time = dose_times) |>
    dplyr::mutate(
      evid = 1L,
      amt  = dose_mg,
      cmt  = "central",
      rate = -2,             # rxode2: use dur(central) = d1 for zero-order abs
      cohort = cohort_id,
      patient_type = patient_type,
      dose_mg_per_dose = dose_mg,
      interval_h = interval_h
    )
  obs_pk_rows <- expand.grid(id = ids, time = obs_times_pk) |>
    dplyr::mutate(
      evid = 0L,
      amt  = 0,
      cmt  = "Cc",
      rate = 0,
      cohort = cohort_id,
      patient_type = patient_type,
      dose_mg_per_dose = dose_mg,
      interval_h = interval_h
    )
  obs_vk_rows <- expand.grid(id = ids, time = obs_times_vk) |>
    dplyr::mutate(
      evid = 0L,
      amt  = 0,
      cmt  = "Vlog10",
      rate = 0,
      cohort = cohort_id,
      patient_type = patient_type,
      dose_mg_per_dose = dose_mg,
      interval_h = interval_h
    )
  dplyr::bind_rows(dose_rows, obs_pk_rows, obs_vk_rows) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

events <- dplyr::bind_rows(purrr::pmap_dfr(
  list(
    cohort_id    = dose_regimens$cohort,
    dose_mg      = dose_regimens$dose_mg,
    interval_h   = dose_regimens$interval_h,
    patient_type = dose_regimens$patient_type,
    id_offset    = (dose_regimens$cohort - 1L) * n_per_cohort
  ),
  make_cohort,
  n_subj = n_per_cohort,
  dur_h  = sim_dur_h
))

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid", "cmt")])))
nrow(events)
#> [1] 29780

Simulation 

mod <- readModelDb("Canini_2015_danoprevir")

sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("cohort", "patient_type", "dose_mg_per_dose", "interval_h"),
  # useLinCmt = FALSE guards against the multi-output ODE->linCmt auto-
  # conversion that can corrupt dvid mapping for models with multiple
  # algebraic observables backed by different ODE states (see
  # known-vignette-failure-patterns.md pattern 5b).
  useLinCmt = FALSE
) |> as.data.frame()
#> ℹ parameter labels from comments will be replaced by 'label()'

Replicate Figure 2 - Viral decline during 14 days of danoprevir

Canini 2015 Figure 2 shows individual-patient viral-load trajectories over 14 days of danoprevir monotherapy for each cohort, with the median 14-day decline reported in Table 1 (2.0 log10 IU/mL for cohort 1 up to 3.9 log10 IU/mL for cohort 4). The VPC below replicates the aggregate cohort-level trajectories.

sim_vpc <- sim |>
  dplyr::group_by(cohort, dose_mg_per_dose, interval_h, time) |>
  dplyr::summarise(
    Vlog10_Q05 = quantile(Vlog10, 0.05, na.rm = TRUE),
    Vlog10_Q50 = quantile(Vlog10, 0.50, na.rm = TRUE),
    Vlog10_Q95 = quantile(Vlog10, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  dplyr::mutate(
    day = time / 24,
    cohort_label = paste0(
      "Cohort ", cohort, ": ", dose_mg_per_dose, " mg every ",
      interval_h, " h"
    )
  )

ggplot(sim_vpc, aes(day, Vlog10_Q50)) +
  geom_ribbon(aes(ymin = Vlog10_Q05, ymax = Vlog10_Q95), alpha = 0.25) +
  geom_line(linewidth = 0.7) +
  facet_wrap(~cohort_label, ncol = 3) +
  labs(
    x = "Time (days)",
    y = expression(log[10]~"(HCV RNA IU/mL)"),
    title = "Figure 2 replication - 14-day viral decline per cohort",
    caption = paste0(
      "Median (line) and 5th-95th percentile (band) of simulated ",
      n_per_cohort, " subjects per cohort. Replicates the aggregate ",
      "trajectories underlying Canini 2015 Figure 2."
    )
  )

Median 14-day viral decline vs Canini 2015 Table 1 

baseline_median <- sim |>
  dplyr::filter(time == 0) |>
  dplyr::group_by(cohort) |>
  dplyr::summarise(Vlog10_baseline = median(Vlog10, na.rm = TRUE),
                   .groups = "drop")

nadir_median <- sim |>
  dplyr::filter(time >= 24) |>
  dplyr::group_by(id, cohort) |>
  dplyr::summarise(Vlog10_end = Vlog10[which.max(time)], .groups = "drop") |>
  dplyr::group_by(cohort) |>
  dplyr::summarise(Vlog10_end_median = median(Vlog10_end, na.rm = TRUE),
                   .groups = "drop")

published_decline <- tibble::tribble(
  ~cohort, ~published_decline_log10,
       1L,                       2.0,
       2L,                       2.7,
       3L,                       2.3,
       4L,                       3.9,
       5L,                       2.7
)

decline_tbl <- baseline_median |>
  dplyr::left_join(nadir_median, by = "cohort") |>
  dplyr::left_join(published_decline, by = "cohort") |>
  dplyr::mutate(
    simulated_decline_log10 = Vlog10_baseline - Vlog10_end_median,
    pct_diff = 100 * (simulated_decline_log10 - published_decline_log10) /
      published_decline_log10
  ) |>
  dplyr::select(cohort, Vlog10_baseline, Vlog10_end_median,
                simulated_decline_log10, published_decline_log10, pct_diff) |>
  dplyr::rename(
    "Cohort"                                = cohort,
    "Sim baseline (log10 IU/mL)"            = Vlog10_baseline,
    "Sim day-14 median (log10 IU/mL)"       = Vlog10_end_median,
    "Sim 14-day decline (log10 IU/mL)"      = simulated_decline_log10,
    "Published 14-day decline (log10 IU/mL)" = published_decline_log10,
    "Percent difference"                    = pct_diff
  )

knitr::kable(decline_tbl, digits = 2,
             caption = paste0(
               "Simulated vs published (Canini 2015 Table 1) median 14-day ",
               "log10 viral decline. Because the packaged model uses a single ",
               "typical V0 = 5.4e5 (cohort-1 reference), simulated baselines ",
               "will be similar across cohorts and the reported vs simulated ",
               "declines compare the DYNAMIC (percent-inhibition) effect of ",
               "the dose regimen, not the between-cohort baseline shift."))
Simulated vs published (Canini 2015 Table 1) median 14-day log10 viral decline. Because the packaged model uses a single typical V0 = 5.4e5 (cohort-1 reference), simulated baselines will be similar across cohorts and the reported vs simulated declines compare the DYNAMIC (percent-inhibition) effect of the dose regimen, not the between-cohort baseline shift.
Cohort Sim baseline (log10 IU/mL) Sim day-14 median (log10 IU/mL) Sim 14-day decline (log10 IU/mL) Published 14-day decline (log10 IU/mL) Percent difference
1 5.69 3.48 2.21 2.0 10.34
2 5.81 2.17 3.64 2.7 34.77
3 5.69 2.89 2.79 2.3 21.41
4 5.80 2.65 3.15 3.9 -19.18
5 5.74 3.13 2.61 2.7 -3.31

Effectiveness comparison against Canini 2015 Table 3

Canini 2015 Table 3 reports the effectiveness averaged over the first 2 days (e2) and over the full 14-day dosing period (e14) for each cohort. The averages are computed here on the same simulated cohorts. The Canini 2015 point-estimate typical PK profile (obtained by zeroing the random effects with rxode2::zeroRe()) is used so the comparison is the model’s typical-value effectiveness, matching Table 3’s tabulated typical values.

mod_typical <- rxode2::zeroRe(mod)
#> ℹ parameter labels from comments will be replaced by 'label()'

# Use one subject per cohort under typical parameters. Reuse the cohort
# helper with n_subj = 1 and disjoint id offsets so we don't need to
# rebuild the whole event structure.
events_typical <- dplyr::bind_rows(purrr::pmap_dfr(
  list(
    cohort_id    = dose_regimens$cohort,
    dose_mg      = dose_regimens$dose_mg,
    interval_h   = dose_regimens$interval_h,
    patient_type = dose_regimens$patient_type,
    id_offset    = dose_regimens$cohort - 1L
  ),
  make_cohort,
  n_subj = 1L,
  dur_h  = sim_dur_h
))

sim_typical <- rxode2::rxSolve(
  mod_typical,
  events    = events_typical,
  keep      = c("cohort", "dose_mg_per_dose", "interval_h"),
  useLinCmt = FALSE
) |> as.data.frame()
#> ℹ omega/sigma items treated as zero: 'etaltlag', 'etald1', 'etalvc', 'etalke', 'etalk12', 'etalk21', 'etalc', 'etald', 'etalrbase', 'etalec50'
#> Warning: multi-subject simulation without without 'omega'

# e2 average over first 2 days; e14 average over 14 days.
eff_summary <- sim_typical |>
  dplyr::group_by(cohort) |>
  dplyr::summarise(
    e2_sim  = mean(e_t[time <= 48],            na.rm = TRUE),
    e14_sim = mean(e_t[time <= sim_dur_h],     na.rm = TRUE),
    .groups = "drop"
  )

published_eff <- tibble::tribble(
  ~cohort, ~e2_pub, ~e14_pub,
       1L,   0.946,    0.953,
       2L,   0.981,    0.990,
       3L,   0.966,    0.974,
       4L,   0.985,    0.994,
       5L,   0.974,    0.982
)

eff_tbl <- eff_summary |>
  dplyr::left_join(published_eff, by = "cohort") |>
  dplyr::mutate(
    e2_pct_diff  = 100 * (e2_sim  - e2_pub)  / e2_pub,
    e14_pct_diff = 100 * (e14_sim - e14_pub) / e14_pub
  ) |>
  dplyr::rename(
    "Cohort"                     = cohort,
    "Sim e2 (0-2 day mean)"      = e2_sim,
    "Sim e14 (0-14 day mean)"    = e14_sim,
    "Published e2 (Table 3)"     = e2_pub,
    "Published e14 (Table 3)"    = e14_pub,
    "e2 percent difference"      = e2_pct_diff,
    "e14 percent difference"     = e14_pct_diff
  )

knitr::kable(eff_tbl, digits = 3,
             caption = paste0(
               "Simulated vs published (Canini 2015 Table 3 PK/VK model rows) ",
               "typical-value average danoprevir effectiveness. e2 is averaged ",
               "over 0-48 h; e14 is averaged over 0-14 days."))
Simulated vs published (Canini 2015 Table 3 PK/VK model rows) typical-value average danoprevir effectiveness. e2 is averaged over 0-48 h; e14 is averaged over 0-14 days.
Cohort Sim e2 (0-2 day mean) Sim e14 (0-14 day mean) Published e2 (Table 3) Published e14 (Table 3) e2 percent difference e14 percent difference
1 0.912 0.932 0.946 0.953 -3.584 -2.237
2 0.939 0.972 0.981 0.990 -4.274 -1.831
3 0.932 0.956 0.966 0.974 -3.514 -1.826
4 0.944 0.978 0.985 0.994 -4.212 -1.601
5 0.939 0.965 0.974 0.982 -3.558 -1.700

PKNCA validation

The paper reports Cmax between 2.67 and 589 ng/mL across all cohorts. Since per-cohort NCA summaries are not tabulated, the PKNCA check below validates that (a) simulated Cmax and Tmax are order-of-magnitude consistent with the paper’s range, (b) Cmin exceeds EC50 for each cohort (matching the paper’s Discussion claim that ‘EC50 estimate … is only 2.6 to 7.7 x higher than the Cmin of the cohort with the highest and lowest dose in treatment-naive patients’ – i.e. EC50 > Cmin, so Cmin < EC50).

# PKNCA on the first dosing interval (single-dose Cmax and Tmax per subject).
sim_nca <- sim |>
  dplyr::filter(!is.na(Cc), time <= 24) |>
  dplyr::select(id, time, Cc, cohort)

# Guarantee a time-zero row per (id, cohort) at Cc = 0 (extravascular).
sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> 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(sim_nca, Cc ~ time | cohort + id)

# First-dose window per subject: dose at time 0, interval defined per cohort.
dose_df <- events |>
  dplyr::filter(evid == 1, time == 0) |>
  dplyr::select(id, time, amt, cohort)
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | cohort + id)

# Compute Cmax + Tmax over the first-dose window (0 to end of first interval,
# max = 12 h for BID, 8 h for TID cohorts). Using 12 h captures both.
intervals <- data.frame(
  start = 0,
  end   = 12,
  cmax  = TRUE,
  tmax  = TRUE,
  auclast = TRUE
)

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

nca_summary <- as.data.frame(nca_res$result) |>
  dplyr::group_by(cohort, PPTESTCD) |>
  dplyr::summarise(median = median(PPORRES, na.rm = TRUE), .groups = "drop") |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = median) |>
  dplyr::rename(
    "Cohort"        = cohort,
    "Cmax (ng/mL)"  = cmax,
    "Tmax (h)"      = tmax,
    "AUC0-12 (ng*h/mL)" = auclast
  )

knitr::kable(nca_summary, digits = 3,
             caption = paste0(
               "Simulated first-dose PKNCA summary (median across ",
               n_per_cohort, " subjects). Canini 2015 reports Cmax across ",
               "cohorts between 2.67 and 589 ng/mL; per-cohort NCA values are ",
               "not tabulated so this is a range consistency check."))
Simulated first-dose PKNCA summary (median across 20 subjects). Canini 2015 reports Cmax across cohorts between 2.67 and 589 ng/mL; per-cohort NCA values are not tabulated so this is a range consistency check.
Cohort AUC0-12 (ng*h/mL) Cmax (ng/mL) Tmax (h)
1 50.972 28.178 1.500
2 136.463 32.130 1.750
3 64.291 53.851 1.375
4 140.644 54.500 1.625
5 153.918 124.839 1.125

Cmin vs EC50 spot-check

The Canini 2015 Discussion reports minimal effectiveness values of 0.706 (cohort 1) and 0.877 (cohort 4). For the Hill-Emax function with Emax = 1 and h = 1, e_min = Cmin / (EC50 + Cmin) inverts to Cmin / EC50 = e_min / (1 - e_min) = 2.40 (cohort 1) and 7.13 (cohort 4). The Discussion sentence “EC50 … is 2.6 to 7.7 x higher than the Cmin” reverses the direction of the ratio implied by the tabulated minimal-effectiveness values (see Errata / Assumptions); the ratios below are reported as Cmin / EC50 so the sim vs paper comparison uses the direction that is internally consistent with the paper’s 0.706 / 0.877 numbers.

# Cmin at end of first dosing interval per subject.
cmin_by_cohort <- sim |>
  dplyr::group_by(id, cohort, interval_h) |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::slice(which.min(abs(time - interval_h[1]))) |>
  dplyr::ungroup() |>
  dplyr::group_by(cohort) |>
  dplyr::summarise(
    Cmin_median_ngmL = median(Cc, na.rm = TRUE),
    .groups = "drop"
  ) |>
  dplyr::mutate(
    EC50_ngmL         = 0.0082,
    ratio_Cmin_EC50   = Cmin_median_ngmL / EC50_ngmL,
    emin_at_Cmin      = ratio_Cmin_EC50 / (1 + ratio_Cmin_EC50)
  ) |>
  dplyr::rename(
    "Cohort"              = cohort,
    "Median Cmin (ng/mL)" = Cmin_median_ngmL,
    "EC50 (ng/mL)"        = EC50_ngmL,
    "Ratio Cmin / EC50"   = ratio_Cmin_EC50,
    "e at Cmin"           = emin_at_Cmin
  )

knitr::kable(cmin_by_cohort, digits = 4,
             caption = paste0(
               "Median simulated Cmin at end of first dosing interval per ",
               "cohort. Canini 2015 Discussion implies (via the 0.706 / ",
               "0.877 minimal-effectiveness values for cohorts 1 / 4) a ",
               "Cmin / EC50 ratio of 2.40 (cohort 1) to 7.13 (cohort 4). ",
               "The simulated cohort-1 ratio should be near 2.40; cohort-4 ",
               "ratio in the packaged model is larger than the paper's ",
               "implied value because the paper's reported Cmin appears to ",
               "have been computed from a different sampling window than ",
               "our end-of-first-interval slice (see Assumptions)."))
Median simulated Cmin at end of first dosing interval per cohort. Canini 2015 Discussion implies (via the 0.706 / 0.877 minimal-effectiveness values for cohorts 1 / 4) a Cmin / EC50 ratio of 2.40 (cohort 1) to 7.13 (cohort 4). The simulated cohort-1 ratio should be near 2.40; cohort-4 ratio in the packaged model is larger than the paper’s implied value because the paper’s reported Cmin appears to have been computed from a different sampling window than our end-of-first-interval slice (see Assumptions).
Cohort Median Cmin (ng/mL) EC50 (ng/mL) Ratio Cmin / EC50 e at Cmin
1 0.0311 0.0082 3.7956 0.7915
2 0.3161 0.0082 38.5494 0.9747
3 0.0653 0.0082 7.9579 0.8884
4 0.5244 0.0082 63.9489 0.9846
5 0.0804 0.0082 9.8009 0.9074

Assumptions and deviations

  • Effectiveness driver = C1 (plasma), not C2 (tissue). The abstract wording contradicts the Methods statement. The packaged model follows the explicit Methods statement (‘e(t) varies as a function of C1’). See the Errata section above for the discriminating evidence (V2 unreported; Discussion compares EC50 to plasma Cmin).

  • Cohort-specific V0 subsumed by IIV. The packaged model uses the cohort-1 reference typical V0 = 5.4e5 IU/mL. The paper’s cohort-specific V0 typical values (1.9e6 to 2.7e6 for cohorts 2-5) are between 0.5 and 0.7 log10 IU/mL above the reference and are well within the paper’s 118 percent IIV on V0. This choice keeps the packaged model portable to simulation contexts outside the original 5-cohort study design; a user who wants to reproduce Canini 2015 cohort-specific baseline VLs can either post-scale the simulated baseline (multiply Vlog10 at t = 0 by the cohort-specific ratio in log10 space) or bias the etalrbase distribution.

  • HCV genotype screened but not retained. Canini 2015 tested HCV genotype (GT1A vs GT1B vs undetermined GT1) as a covariate on the PK/VK parameters (Methods, ‘HCV genotype was tested as a covariate in the model to study its effect on the PK/VK parameters’) but no genotype effect appears in Table 2, so the covariate was not retained. It is captured in covariatesDataExcluded for provenance and is not referenced in model().

  • Neumann parameters p, beta fixed by convention (p = 1). Because Canini 2015 estimates only c, d, V0, EC50 in the VK layer, p and beta are related by the pre-treatment steady-state constraint p * beta_T = d * c and are individually unidentifiable. The packaged model fixes p = 1 as a mathematical scaling and derives beta_T = d * c, I(0) = c * V0. This does NOT change the observable behaviour (Vlog10 trajectories) but shifts the interpretation of the infected state away from a biological hepatocyte-count scale. Users who want a biological-scale infected interpretation can rescale infected and beta_t by their preferred p (e.g. from Neumann 1998).

  • Age / weight / sex not encoded. Baseline demographics beyond the cohort structure are not reported in the trimmed source (cited to the original Phase 1 report, reference [9]); no demographic covariates are used in the final PK/VK model per Canini 2015 Table 2.

  • Discussion Cmin ratio wording. The Canini 2015 Discussion states that “EC50 … is 2.6 to 7.7 x higher than the Cmin”, but the same paragraph reports minimal effectiveness values of 0.706 (cohort 1) and 0.877 (cohort 4). For the Hill-Emax e = C / (EC50 + C) with h = 1 and Emax = 1, 0.706 implies Cmin / EC50 = 2.40 (i.e., Cmin is 2.4 x HIGHER than EC50) and 0.877 implies Cmin / EC50 = 7.13. We treat the “higher than” direction in the Discussion sentence as a transcription slip and report Cmin / EC50 in the spot-check table above.