---
title: "PF-04878691 TLR7 agonist for chronic hepatitis C (Jones 2011)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{PF-04878691 TLR7 agonist for chronic hepatitis C (Jones 2011)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                     fig.width = 7, fig.height = 4)
library(nlmixr2lib)
library(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
set.seed(20260613L)
```

# Model and source

Jones HM, Chan PLS, van der Graaf PH, Webster R. Use of modelling and
simulation techniques to support decision making on the progression of
PF-04878691, a TLR7 agonist being developed for hepatitis C.
*Br J Clin Pharmacol.* 2012;73(1):77-92.

* Article: <https://doi.org/10.1111/j.1365-2125.2011.04047.x>
* ClinicalTrials.gov: NCT00810758 (PF-04878691 multiple-dose escalation, healthy volunteers).

This paper develops four sequentially-fit population PK/PD models from
two clinical data sets and uses the chain to predict the antiviral
efficacy of PF-04878691 (a toll-like-receptor-7 / TLR7 agonist) in
chronic hepatitis C (HCV) patients. The PK and biomarker fits were
performed on a Phase 1 multiple-dose escalation study in 24 healthy
adult volunteers (PF-04878691, 3 / 6 / 9 mg orally twice weekly for
2 weeks); the OAS-viral-load relationship was fit on a previously
published Phase 1 study of CPG-10101, a TLR9 agonist, in 39 chronic-HCV
patients (Jones 2011 reference [16] = McHutchison 2007). The four
nlmixr2lib model files are:

| File                                            | Layer                                         | Source table                |
|-------------------------------------------------|-----------------------------------------------|-----------------------------|
| `Jones_2011_PF04878691`                         | PK only (2-compartment + time-varying CL)     | Table 1                     |
| `Jones_2011_PF04878691_oas`                     | PK + OAS gene-expression fold change          | Tables 1, 2                 |
| `Jones_2011_PF04878691_lymphocyte`              | PK + absolute lymphocyte count                | Tables 1, 3                 |
| `Jones_2011_PF04878691_viralLoad`               | PK + OAS + HCV viral RNA                      | Tables 1, 2, 4              |

The viral-load file combines the PF-04878691 PK + OAS sub-models with
the CPG-10101-fit OAS-viral-load relationship; this chain mirrors the
paper's Figure 10 simulation that drove the no-progress decision.

```{r load-models}
mod_pk    <- readModelDb("Jones_2011_PF04878691")
mod_oas   <- readModelDb("Jones_2011_PF04878691_oas")
mod_lymph <- readModelDb("Jones_2011_PF04878691_lymphocyte")
mod_vl    <- readModelDb("Jones_2011_PF04878691_viralLoad")
```

# Population

PK / OAS / lymphocyte data: 24 healthy adult volunteers (Methods,
"Clinical TLR7 study data"), age 21-55 (median 34) years, body weight
57-97 (median 79) kg, 92% male. Six active subjects + two placebo per
3 / 6 / 9 mg cohort; placebo data are not used in any model. Two SAEs
in the 9 mg cohort prematurely terminated the study with four subjects
withdrawing during active treatment after two doses.

OAS-viral-load fit (Methods, "Clinical TLR9 study data"): 39 chronic
HCV patients of 60 randomised, dosed subcutaneously with CPG-10101
0.25 / 1 / 4 / 10 / 20 mg twice weekly or 0.5 / 0.75 mg/kg once weekly,
for 4 weeks (Jones 2011 reference [16]).

The same information is available programmatically:

```{r population}
rxode2::rxode(mod_pk)$population
```

# Source trace

Per-parameter origins are recorded as in-file comments next to each
`ini()` entry; the table below collects them in one place.

| Equation / parameter | Value | Source |
|----------------------|-------|--------|
| Two-compartment + time-varying CL: `CL(t) = CL_SS + CL_TIME * exp(-kdeg * t)` | n/a | Methods, "Population PK model"; Table 1 |
| `cl_ss` (= CLF, L/h/kg)                       | 1.7    | Table 1, TH1 (CLF)                       |
| `cl_time` initial offset (= CL0 - CLF)        | 1.8    | Derived from Table 1 (CL0 - CLF)         |
| `kdeg` (= DEG, 1/h)                           | 0.24   | Table 1, TH7                             |
| `vc` (Vc, L/kg)                               | 3.3    | Table 1, TH2                             |
| `lq` (Q, L/h/kg)                              | 0.74   | Table 1, TH4                             |
| `lvp` (Vp, L/kg)                              | 21     | Table 1, TH5                             |
| `lka` (ka, 1/h)                               | 0.078  | Table 1, TH3                             |
| IIV(`lcl`) (= IIV CLF)                        | 0.067  | Table 1, OM1                             |
| IIV(`lka`)                                    | 0.19   | Table 1, OM3                             |
| PK `propSd`                                   | 0.046  | Table 1, SIG1                            |
| `dOAS/dt = kin * (1 + slope * Cc^gamma) - kout * OAS`, with `kin = rbase * kout` | n/a | Methods, "Population PK-OAS and PK-lymphocyte models" |
| OAS `lkout` (1/h)                             | 0.034  | Table 2, TH1                             |
| OAS `lslope` (per (ng/mL)^gamma)              | 3.5    | Table 2, TH2                             |
| OAS `lrbase` (fold change)                    | 0.96   | Table 2, TH3                             |
| OAS `lgamma`                                  | 1.6    | Table 2, TH4                             |
| OAS IIVs (`kout`, `rbase`)                    | 1.7, 0.18 | Table 2, OM1 / OM3                    |
| OAS `propSd`                                  | 0.19   | Table 2, SIG1                            |
| `dLYMPH/dt = kin - kout * (1 + slope * Cc^gamma) * LYMPH`, with `kin = rbase * kout` | n/a | Methods, "Population PK-OAS and PK-lymphocyte models" |
| Lymph `lkout` (1/h)                           | 0.044  | Table 3, TH1                             |
| Lymph `lslope` (per (ng/mL)^gamma)            | 0.44   | Table 3, TH2                             |
| Lymph `lrbase` (cells/uL)                     | 1890   | Table 3, TH3 (paper-printed unit "pg/mL" is a typo; see Errata) |
| Lymph `lgamma`                                | 2.2    | Table 3, TH4                             |
| Lymph IIVs (`kout`, `slope`, `rbase`)         | 0.19, 0.20, 0.051 | Table 3, OM1 / OM2 / OM3      |
| Lymph `propSd`                                | 0.021  | Table 3, SIG1                            |
| `vload = BASE + Imax * oas_fc_above^gamma / (VO50^gamma + oas_fc_above^gamma)`, `oas_fc_above = oas / rbase_oas - 1` | n/a | Methods, "Population OAS-viral load model" |
| VL `lbase_vl` (log10 copies/mL)               | 7.3    | Table 4, TH1                             |
| VL `limax_vl` (|Imax|, log10 copies/mL; signed inside model() as `-imax_abs`) | 2.7 | Table 4, TH2 (= -2.7) |
| VL `lvo50_vl` (fold change above baseline)    | 3.6    | Table 4, TH3                             |
| VL `lgamma_vl`                                | 0.68   | Table 4, TH4                             |
| VL IIVs (Imax, VO50)                          | 0.29, 0.25 | Table 4, OM2 / OM3 (see Errata)      |
| VL `addSd_vload` (log10 scale)                | 0.41   | Table 4, TH5                             |

# Virtual cohort

The Jones 2011 PF-04878691 cohort was a small Phase 1 trial; we
recreate a comparable virtual cohort for the PK / OAS / lymphocyte
simulations.

```{r cohort-pk}
make_phase1_cohort <- function(n_per_dose, dose_mg, wt_median = 79,
                               wt_sd = 10, id_offset = 0L,
                               obs_times = unique(c(seq(0, 24, by = 1),
                                                     seq(24, 264, by = 6),
                                                     seq(264, 312, by = 12)))) {
  ids <- id_offset + seq_len(n_per_dose)
  tibble::tibble(id = ids,
                 dose_mg = dose_mg,
                 WT = pmax(50, rnorm(n_per_dose, wt_median, wt_sd)))
}

build_events <- function(cohort, dose_times = c(0, 72, 168, 240),
                         obs_times = unique(c(seq(0, 24, by = 1),
                                              seq(30, 312, by = 6)))) {
  doses <- cohort |>
    tidyr::expand_grid(time = dose_times) |>
    dplyr::transmute(id, time, evid = 1L, amt = dose_mg,
                     cmt = "depot", WT, dose_mg, treatment)
  obs <- cohort |>
    tidyr::expand_grid(time = obs_times) |>
    dplyr::transmute(id, time, evid = 0L, amt = NA_real_,
                     cmt = NA_character_, WT, dose_mg, treatment)
  dplyr::bind_rows(doses, obs) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

cohort <- dplyr::bind_rows(
  make_phase1_cohort(50, 3, id_offset =   0L) |> dplyr::mutate(treatment = "3 mg"),
  make_phase1_cohort(50, 6, id_offset =  50L) |> dplyr::mutate(treatment = "6 mg"),
  make_phase1_cohort(50, 9, id_offset = 100L) |> dplyr::mutate(treatment = "9 mg")
)
events_pk <- build_events(cohort)
stopifnot(!anyDuplicated(unique(events_pk[, c("id", "time", "evid")])))
```

# Simulation

```{r sim-pk}
sim_pk <- rxode2::rxSolve(mod_pk, events = events_pk,
                          keep = c("treatment", "dose_mg")) |>
  as.data.frame()
```

```{r sim-typical}
mod_pk_typ <- rxode2::zeroRe(mod_pk)
sim_pk_typ <- rxode2::rxSolve(mod_pk_typ, events = events_pk,
                              keep = c("treatment", "dose_mg")) |>
  as.data.frame()
```

# Replicate published figures

## Figure 2 - PK profile

```{r figure-2}
sim_pk |>
  dplyr::filter(time > 0, Cc > 0) |>
  dplyr::group_by(treatment, time) |>
  dplyr::summarise(Q05 = quantile(Cc, 0.05),
                   Q50 = quantile(Cc, 0.50),
                   Q95 = quantile(Cc, 0.95),
                   .groups = "drop") |>
  ggplot(aes(time, Q50, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.2, colour = NA) +
  geom_line() +
  scale_y_log10() +
  labs(x = "Time after first dose (h)", y = "PF-04878691 plasma concentration (ng/mL)",
       title = "Replicates Figure 2 of Jones 2011 (PK VPC)",
       caption = "Median and 5-95 percentile ribbon, 50 subjects per dose group.")
```

## Day 1 vs Day 11 Cmax

The paper reports that Cmax on day 11 is up to three times higher than
on day 1 because of the time-varying clearance.

```{r cmax-d1-vs-d11}
window_cmax <- function(df, start_h, end_h) {
  df |>
    dplyr::filter(time > start_h, time <= end_h) |>
    dplyr::group_by(id, treatment) |>
    dplyr::summarise(cmax = max(Cc, na.rm = TRUE), .groups = "drop")
}
d1  <- window_cmax(sim_pk, start_h =   0, end_h =  72)
d11 <- window_cmax(sim_pk, start_h = 240, end_h = 312)
ratio <- dplyr::inner_join(d1, d11, by = c("id", "treatment"),
                           suffix = c("_d1", "_d11")) |>
  dplyr::mutate(ratio = cmax_d11 / cmax_d1)
ratio |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(median_ratio = median(ratio),
                   q05 = quantile(ratio, 0.05),
                   q95 = quantile(ratio, 0.95),
                   .groups = "drop") |>
  knitr::kable(digits = 2,
               caption = "Simulated Cmax(day 11) / Cmax(day 1). Paper text reports 'up to three times higher'.")
```

## Figure 4 - OAS time course

```{r sim-oas}
sim_oas <- rxode2::rxSolve(mod_oas, events = events_pk,
                           keep = c("treatment", "dose_mg")) |>
  as.data.frame()
```

```{r figure-4}
sim_oas |>
  dplyr::filter(time >= 0) |>
  dplyr::group_by(treatment, time) |>
  dplyr::summarise(Q05 = quantile(oas, 0.05),
                   Q50 = quantile(oas, 0.50),
                   Q95 = quantile(oas, 0.95),
                   .groups = "drop") |>
  ggplot(aes(time, Q50, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.2, colour = NA) +
  geom_line() +
  labs(x = "Time after first dose (h)", y = "OAS fold change",
       title = "Replicates Figure 4 of Jones 2011 (OAS VPC)",
       caption = "Median and 5-95 percentile ribbon, 50 subjects per dose group.")
```

```{r oas-peak-fold}
sim_oas |>
  dplyr::group_by(id, treatment) |>
  dplyr::summarise(oas_peak = max(oas, na.rm = TRUE), .groups = "drop") |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(median_peak = median(oas_peak),
                   q05 = quantile(oas_peak, 0.05),
                   q95 = quantile(oas_peak, 0.95),
                   .groups = "drop") |>
  knitr::kable(digits = 2,
               caption = "Simulated peak OAS fold change per dose group. Paper observed OAS increases >= 8-fold from baseline in 3/6 individuals at 3 mg, 6/6 at 6 mg, and 5/6 at 9 mg.")
```

## Figure 6 - Lymphocyte time course

```{r sim-lymph}
sim_lymph <- rxode2::rxSolve(mod_lymph, events = events_pk,
                             keep = c("treatment", "dose_mg")) |>
  as.data.frame()
```

```{r figure-6}
sim_lymph |>
  dplyr::filter(time >= 0) |>
  dplyr::group_by(treatment, time) |>
  dplyr::summarise(Q05 = quantile(lymph, 0.05),
                   Q50 = quantile(lymph, 0.50),
                   Q95 = quantile(lymph, 0.95),
                   .groups = "drop") |>
  ggplot(aes(time, Q50, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.2, colour = NA) +
  geom_line() +
  geom_hline(yintercept = 500, linetype = "dashed", colour = "grey50") +
  geom_hline(yintercept = 800, linetype = "dotted", colour = "grey50") +
  labs(x = "Time after first dose (h)", y = "Absolute lymphocyte count (cells/uL)",
       title = "Replicates Figure 6 of Jones 2011 (lymphocyte VPC)",
       caption = "Dashed line = grade 3 lymphopenia (< 500 cells/uL), dotted line = grade 2 (< 800 cells/uL). 50 subjects per dose.")
```

## Figure 10 - Viral load chain (Jones 2011 4-week HCV simulation)

The Figure 10 simulation extends dosing to 4 weeks twice-weekly and
predicts a 1-2 log drop in viral load at the 6 mg dose. The chain
combines the PF-04878691 PK + OAS sub-models with the CPG-10101 / HCV
OAS-viral-load relationship.

```{r cohort-vl}
make_4wk_cohort <- function(n_per_dose, dose_mg, wt_median = 79,
                            wt_sd = 10, id_offset = 0L) {
  ids <- id_offset + seq_len(n_per_dose)
  tibble::tibble(id = ids, dose_mg = dose_mg,
                 WT = pmax(50, rnorm(n_per_dose, wt_median, wt_sd)))
}

dose_times_4wk <- seq(0, 24 * 25, by = 84)        # twice weekly for 4 weeks (8 doses)
obs_times_4wk  <- seq(0, 24 * 35, by = 6)         # 35 days = 4 weeks + tail

build_vl_events <- function(cohort, dose_times = dose_times_4wk,
                            obs_times = obs_times_4wk) {
  doses <- cohort |>
    tidyr::expand_grid(time = dose_times) |>
    dplyr::transmute(id, time, evid = 1L, amt = dose_mg,
                     cmt = "depot", WT, dose_mg, treatment)
  obs <- cohort |>
    tidyr::expand_grid(time = obs_times) |>
    dplyr::transmute(id, time, evid = 0L, amt = NA_real_,
                     cmt = NA_character_, WT, dose_mg, treatment)
  dplyr::bind_rows(doses, obs) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

cohort_vl <- dplyr::bind_rows(
  make_4wk_cohort(50, 3,   id_offset =   0L) |> dplyr::mutate(treatment = "3 mg BIW"),
  make_4wk_cohort(50, 6,   id_offset =  50L) |> dplyr::mutate(treatment = "6 mg BIW"),
  make_4wk_cohort(50, 4.5, id_offset = 100L) |> dplyr::mutate(treatment = "3 mg BIW then 4.5 mg BIW")
)
events_vl <- build_vl_events(cohort_vl)

sim_vl <- rxode2::rxSolve(mod_vl, events = events_vl,
                          keep = c("treatment", "dose_mg")) |>
  as.data.frame()
```

```{r figure-10}
sim_vl |>
  dplyr::filter(time >= 0) |>
  dplyr::group_by(treatment, time) |>
  dplyr::summarise(Q05 = quantile(vload, 0.05),
                   Q50 = quantile(vload, 0.50),
                   Q95 = quantile(vload, 0.95),
                   .groups = "drop") |>
  ggplot(aes(time / 24, Q50, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.2, colour = NA) +
  geom_line() +
  geom_hline(yintercept = 7.3 - 1, linetype = "dashed", colour = "grey50") +
  labs(x = "Time after first dose (day)", y = "HCV viral RNA (log10 copies/mL)",
       title = "Replicates Figure 10 of Jones 2011 (viral-load chain simulation)",
       caption = "Dashed line = 1 log10 drop from baseline (POC criterion threshold).")
```

```{r vl-trough}
sim_vl |>
  dplyr::filter(time >= 24 * 27) |>     # end of treatment window, day 27+
  dplyr::group_by(id, treatment) |>
  dplyr::summarise(vload_min = min(vload, na.rm = TRUE), .groups = "drop") |>
  dplyr::mutate(log_drop = 7.3 - vload_min) |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(median_drop = median(log_drop),
                   q05 = quantile(log_drop, 0.05),
                   q95 = quantile(log_drop, 0.95),
                   .groups = "drop") |>
  knitr::kable(digits = 2,
               caption = "Simulated maximum log10 viral-load drop from baseline (Day 27+). Paper Figure 10 predicts 1-2 log10 drop at doses > 6 mg.")
```

# PKNCA validation (PK)

The paper does not tabulate Cmax / AUC / half-life values, but it
reports a terminal half-life of 12-16 h. Run PKNCA on the simulated PK
profiles and inspect the terminal half-life.

```{r pknca-setup}
sim_nca <- sim_pk |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

# Guarantee a time=0 row per (id, treatment): for extravascular pre-dose
# Cc = 0 is the correct value.
sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |>
    dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

dose_df <- events_pk |>
  dplyr::filter(evid == 1L) |>
  dplyr::transmute(id, time, amt, treatment) |>
  dplyr::filter(time == 0)        # use single-dose interval for half-life

# Only keep concentrations from day 1 (single-dose interval) so the
# terminal phase is captured before re-dosing perturbs the profile.
conc_d1 <- sim_nca |> dplyr::filter(time <= 72)

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

intervals <- data.frame(
  start      = 0,
  end        = 72,
  cmax       = TRUE,
  tmax       = TRUE,
  auclast    = TRUE,
  half.life  = TRUE
)

nca_res <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj, dose_obj,
                                          intervals = intervals))
nca_summary <- as.data.frame(nca_res$result) |>
  dplyr::filter(PPTESTCD %in% c("cmax", "tmax", "auclast", "half.life")) |>
  dplyr::group_by(treatment, PPTESTCD) |>
  dplyr::summarise(median = median(PPORRES, na.rm = TRUE),
                   q05    = quantile(PPORRES, 0.05, na.rm = TRUE),
                   q95    = quantile(PPORRES, 0.95, na.rm = TRUE),
                   .groups = "drop")
```

```{r pknca-table}
published <- tibble::tribble(
  ~treatment, ~half.life,
  "3 mg",     14,
  "6 mg",     14,
  "9 mg",     14
)

simulated_wide <- nca_summary |>
  dplyr::transmute(treatment, parameter = PPTESTCD, value = median) |>
  tidyr::pivot_wider(names_from = parameter, values_from = value)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated     = nca_res,
  reference     = published,
  by            = "treatment",
  units         = c(cmax = "ng/mL", tmax = "h", auclast = "h*ng/mL",
                    half.life = "h"),
  tolerance_pct = 30
)

knitr::kable(cmp,
             caption = paste("Simulated NCA on day-1 interval (0-72 h)",
                             "vs. published terminal half-life range (12-16 h, midpoint 14 h).",
                             "* differs from reference by >30%."),
             align   = c("l", "l", "r", "r", "r"))
```

# Assumptions and deviations

- **Time-varying CL reparameterisation.** Jones 2011 Table 1 reports
  CLF, CL0 and DEG; the model files encode the equivalent canonical
  decomposition `CL(t) = CL_SS + CL_TIME * exp(-kdeg * t)` with
  `CL_SS = CLF = 1.7 L/h/kg`, `CL_TIME = CL0 - CLF = 1.8 L/h/kg`, and
  `kdeg = DEG = 0.24 1/h`. The two parameterisations are mathematically
  identical and produce identical PK time-courses; the canonical form
  matches the registered `lcl_ss` / `lcl_time` / `lkdeg` pattern used by
  Gibiansky 2014, Lu 2019, Wu 2024 and other time-varying-CL extractions.
- **Per-kg disposition parameters.** Jones 2011 normalised doses by body
  weight before fitting, so all clearance / volume parameters carry per-kg
  units. The model files require body weight `WT` as a covariate and
  recover per-subject disposition values inside `model()` as
  `cl_typ = exp(lcl + eta) * WT` (and similarly for Vc, Q, Vp). The paper
  itself reports "covariates were not included in the modelling"; WT is a
  scaling factor, not an estimated covariate effect.
- **Cc unit for PD effects.** The paper reports PK concentrations in
  ng/mL (LLOQ 0.1 ng/mL); the OAS / lymphocyte / viral-load Tables 2-4
  parameter values (slope, gamma, VO50) carry units consistent with Cc
  in ng/mL. The PK model files compute `Cc = 1000 * central / vc` so the
  scale is correct in every coupled model (`central / vc` alone has
  units mg/L = ug/mL = 1000x too low). Without this multiplier the
  simulated drug effects vanish because the typical Cmax 1-5 ng/mL would
  be 0.001-0.005 ug/mL.
- **Lymphocyte baseline-unit typo in paper Table 3.** Table 3 prints
  "pg/mL" for the BASE row, but lymphocytes were measured by
  immunophenotyping as an absolute count. The numerical value 1890 is
  the typical absolute lymphocyte count in cells/uL (within the
  reference range 1000-4500 cells/uL); the model files treat `rbase` as
  cells/uL.
- **IIV-on-Imax-and-VO50 vs IIV-on-Imax-and-gamma in the viral-load
  model.** Jones 2011 Results state "IIV was modelled as a variance
  covariance matrix for Imax and g [gamma]" and "the CV% for all
  parameters were reasonable (< 40%) except for IIV of g, 64.3%".
  Table 4 however labels its omega rows as OM2 = IIV Imax (0.29) and
  OM3 = IIV VO50 (0.25; %CV 64). The encoded model files follow the
  Table 4 labels (IIV on Imax and VO50) because the table is the more
  precise numerical source; the text-vs-table discrepancy is documented
  here without alteration. The numerical variances 0.29 and 0.25 are
  the reported point estimates regardless of which parameter they
  belong to, and no off-diagonal covariance is reported in Table 4 (so
  the variance-covariance matrix mentioned in the text is encoded as a
  diagonal matrix).
- **No off-diagonal IIV covariances.** Across all four tables only the
  viral-load model is described in the Results text as carrying a
  variance covariance matrix; no off-diagonal entries are reported in
  any of the four tables. All IIVs are encoded as independent omegas.
- **OAS fold-change interpretation in the viral-load layer.** The
  viral-load model's input is the OAS deviation from baseline expressed
  as a fold change above 1: `oas_fc_above = oas / rbase_oas - 1`. At
  the OAS baseline the deviation is zero and the viral load equals
  BASE = 7.3 log10 copies/mL; above baseline, viral load decreases by
  up to |Imax| = 2.7 log10 copies/mL. This interpretation reconciles
  paper Table 4 (BASE = 7.3 reproduces typical HCV-patient pre-dose
  viral load) with the OAS model's baseline of 0.96 fold change, and
  matches the paper's quantitative prediction of a 1-2 log10 reduction
  at 6 mg twice-weekly (verified above in the Figure 10 chunk).
- **Cross-cohort viral-load extrapolation.** The viral-load file
  combines PF-04878691 healthy-volunteer PK + OAS (Tables 1, 2) with the
  CPG-10101 HCV-patient OAS-viral-load relationship (Table 4) under the
  paper's two explicit assumptions: no PK or biomarker difference
  between healthy volunteers and HCV patients, and TLR7 / TLR9
  agonists produce comparable OAS-driven antiviral effects. These are
  paper assumptions, not nlmixr2lib assumptions; downstream users
  should weigh them before applying the model to other scenarios.
