---
title: "Vrc01 (Huang 2017)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Vrc01 (Huang 2017)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.width  = 7,
  fig.height = 4.5
)
```

```{r libs, message = FALSE}
library(nlmixr2lib)
library(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Huang_2017_vrc01"))$reference`
* Description: `r rxode2::rxode(readModelDb("Huang_2017_vrc01"))$description`
* Article: <https://doi.org/10.1080/19420862.2017.1311435>

# Population

The Huang 2017 popPK model was fit to 1117 VRC01 serum concentrations from 84
HIV-uninfected adults (42 men, 42 women) enrolled in HVTN 104
(NCT02165267). Median age was 27 years (range 18-50) and median body weight
was 72 kg (range 53-114). The IV-only subset had a median weight of 74.5 kg,
which Huang 2017 used as the body-weight reference for covariate centering.
Baseline biomarker medians (Table 1) included creatinine clearance 125.9
mL/min, hemoglobin 14 g/dL, and white blood cell count 6.55 x 10^3/mm^3.
Race / ethnicity composition was not reported in the manuscript or in the
HVTN 104 baseline-demographics table reproduced in the paper.

Participants received either:

- **Group 1**: 40 mg/kg IV loading dose, then 20 mg/kg IV every 4 weeks.
- **Groups 2 / 4 / 5**: 10, 30, or 40 mg/kg IV every 8 weeks.
- **Group 3**: 40 mg/kg IV loading dose, then 5 mg/kg subcutaneously every
  2 weeks for ~5.5 months.

Final-model parameter estimates and the body-weight covariate effects come
from Table 2 of Huang 2017 (pooled IV + SC fit, "Final model" columns); the
IV-only fit appears in Supplementary Table S2 and is not implemented here.
External validation in Huang 2017 used the VRC602 first-in-human study
(reference 18) but those data were not used to estimate parameters.

The population metadata is also available programmatically via
`readModelDb("Huang_2017_vrc01")$population`.

# Source trace

| Equation / parameter             | Value                                  | Source location                     |
|----------------------------------|----------------------------------------|-------------------------------------|
| Two-compartment ODE structure    | depot -> central <-> peripheral1       | Huang 2017 Methods, Figure S2       |
| `lka` (SC absorption rate)       | 0.26 / day                             | Table 2, Final model (IV+SC)        |
| `lcl` (CL at WT = 74.5 kg)       | 0.40 L/day                             | Table 2, Final model (IV+SC)        |
| `lvc` (Vc at WT = 74.5 kg)       | 1.94 L                                 | Table 2, Final model (IV+SC)        |
| `lq`  (Q at WT = 74.5 kg)        | 0.84 L/day                             | Table 2, Final model (IV+SC)        |
| `lvp` (Vp at WT = 74.5 kg)       | 4.90 L                                 | Table 2, Final model (IV+SC)        |
| `lfdepot` (SC bioavailability)   | 0.74                                   | Table 2, Final model (IV+SC)        |
| `e_wt_cl` (WT effect on CL)      | 0.012 fold/kg, exponential form        | Table 2, footnote 2                 |
| `e_wt_vc` (WT effect on Vc)      | 0.010 fold/kg, exponential form        | Table 2, footnote 2                 |
| `e_wt_q`  (WT exponent on Q)     | 0.69, power form                       | Table 2, footnote 2                 |
| `e_wt_vp` (WT exponent on Vp)    | 0.82, power form                       | Table 2, footnote 2                 |
| Omega^2 CL                       | 0.067                                  | Table 2                             |
| Omega^2 Vc                       | 0.028                                  | Table 2                             |
| Omega^2 Q                        | 0.063                                  | Table 2                             |
| Omega^2 Vp                       | 0.120                                  | Table 2                             |
| Cov(CL, Q)                       | 0.034 (R = 0.52)                       | Table 2                             |
| Cov(CL, Vp)                      | 0.050 (R = 0.56)                       | Table 2                             |
| Cov(Q, Vp)                       | 0.082 (R = 0.95)                       | Table 2                             |
| Sigma^2 proportional             | 0.042                                  | Table 2                             |
| Sigma^2 additive                 | 0.456 (ug/mL)^2                        | Table 2                             |
| Body-weight reference            | 74.5 kg                                | Table 2 footnote 4 / Methods        |
| Combination residual error form  | C_obs = C_pred * (1 + e1) + e2         | Methods, "Variability popPK model"  |
| F1 / ka IIV fixed at 0           | -                                      | Table 2 footnote 3                  |

# Errata

No published erratum or corrigendum was identified for Huang 2017 (Mabs
2017;9(5):792-800) at the time of model extraction. Two source ambiguities
worth flagging:

1. The abstract describes the body-weight effects on Q and Vp as "0.69
   log(L/day)/kg" and "0.82 log(L)/kg" respectively. Those units do not
   apply to a power-form covariate exponent. The Methods section
   ("Covariate model") and Table 2 footnote 2 confirm Q and Vp use the
   power form `TV_theta = theta * (BW / 74.5)^beta` with the listed
   exponents (which are dimensionless). The model file follows the
   Methods / Table 2 convention.

2. The five PK parameters listed in the Results paragraph for the base
   model are CL, Vc, Vp, Q, and ka, but Table 2 reports IIV variances
   only for CL, Vc, Q, Vp (with F1 and ka variances "fixed at 0,"
   footnote 3). The model file follows Table 2.

# Virtual cohort

We simulate the four steady-state regimens that Huang 2017 reports trough
predictions for in the Results section and Figure 5: (A) 10 mg/kg IV every
8 weeks, (B) 30 mg/kg IV every 8 weeks, (C) 5 mg/kg SC every 2 weeks with a
10 mg/kg IV loading dose, and (D) 5 mg/kg SC every 2 weeks with a 30 mg/kg
IV loading dose. Body weight is fixed at 74.5 kg (the median IV-group
weight used as the model's covariate-centering reference and the value
Huang 2017 used for the Figure 5 simulations).

```{r make-cohort}
make_regimen <- function(id_offset, regimen_label, wt_kg = 74.5,
                         iv_doses = numeric(0), iv_times = numeric(0),
                         sc_doses = numeric(0), sc_times = numeric(0),
                         obs_times) {
  n_iv  <- length(iv_doses)
  n_sc  <- length(sc_doses)
  n_obs <- length(obs_times)

  doses <- bind_rows(
    if (n_iv > 0) tibble(
      id   = id_offset + 1L,
      time = iv_times,
      evid = 1L,
      amt  = iv_doses,
      cmt  = "central"
    ),
    if (n_sc > 0) tibble(
      id   = id_offset + 1L,
      time = sc_times,
      evid = 1L,
      amt  = sc_doses,
      cmt  = "depot"
    )
  )

  obs <- tibble(
    id   = id_offset + 1L,
    time = obs_times,
    evid = 0L,
    amt  = 0,
    cmt  = "central"
  )

  bind_rows(doses, obs) |>
    mutate(
      WT       = wt_kg,
      regimen  = regimen_label
    ) |>
    arrange(time, desc(evid))
}

# Multi-subject expansion: replicate the regimen for `n` subjects with
# a shared WT (the published Figure 5 simulation also held WT = 74.5 kg).
# Each subject gets a unique ID; covariate columns and timing are reused.
expand_regimen <- function(template, n, id_offset = 0L) {
  one <- template |> mutate(id = NULL)
  bind_rows(
    lapply(seq_len(n), function(k) {
      one |> mutate(id = id_offset + k)
    })
  )
}
```

```{r build-events}
set.seed(20170401)  # arbitrary; 2017-04 = Huang 2017 publication month

n_per_arm <- 200L            # virtual subjects per regimen
horizon   <- 32 * 7          # 32 weeks (the HVTN 104 follow-up window)
obs_grid  <- seq(0, horizon, by = 1)

dose_per_kg_iv_q8w_10 <- 10 * 74.5
dose_per_kg_iv_q8w_30 <- 30 * 74.5
dose_per_kg_iv_load_10 <- 10 * 74.5
dose_per_kg_iv_load_30 <- 30 * 74.5
dose_per_kg_sc_q2w_5   <- 5 * 74.5

iv_q8w_times <- seq(0, horizon - 56, by = 56)   # last dose precedes one full tau of follow-up
sc_q2w_times <- seq(14, horizon - 14, by = 14)  # last SC dose precedes one final tau of follow-up

tmpl_A <- make_regimen(
  id_offset = 0,
  regimen_label = "A: 10 mg/kg IV q8w",
  iv_doses = rep(dose_per_kg_iv_q8w_10, length(iv_q8w_times)),
  iv_times = iv_q8w_times,
  obs_times = obs_grid
)

tmpl_B <- make_regimen(
  id_offset = 0,
  regimen_label = "B: 30 mg/kg IV q8w",
  iv_doses = rep(dose_per_kg_iv_q8w_30, length(iv_q8w_times)),
  iv_times = iv_q8w_times,
  obs_times = obs_grid
)

tmpl_C <- make_regimen(
  id_offset = 0,
  regimen_label = "C: 10 mg/kg IV load + 5 mg/kg SC q2w",
  iv_doses = dose_per_kg_iv_load_10,
  iv_times = 0,
  sc_doses = rep(dose_per_kg_sc_q2w_5, length(sc_q2w_times)),
  sc_times = sc_q2w_times,
  obs_times = obs_grid
)

tmpl_D <- make_regimen(
  id_offset = 0,
  regimen_label = "D: 30 mg/kg IV load + 5 mg/kg SC q2w",
  iv_doses = dose_per_kg_iv_load_30,
  iv_times = 0,
  sc_doses = rep(dose_per_kg_sc_q2w_5, length(sc_q2w_times)),
  sc_times = sc_q2w_times,
  obs_times = obs_grid
)

events <- bind_rows(
  expand_regimen(tmpl_A, n_per_arm, id_offset =   0L),
  expand_regimen(tmpl_B, n_per_arm, id_offset =   n_per_arm),
  expand_regimen(tmpl_C, n_per_arm, id_offset = 2*n_per_arm),
  expand_regimen(tmpl_D, n_per_arm, id_offset = 3*n_per_arm)
)

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- readModelDb("Huang_2017_vrc01")

sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("regimen")
) |> as.data.frame()
```

# Replicate published figures

## Figure 5 — steady-state concentration profiles

Huang 2017 Figure 5 displays simulated VRC01 serum concentrations under
the four regimens above; lines are medians and shaded areas span the
2.5th-97.5th percentiles for 1000 simulated trials.

```{r figure-5, fig.width = 8, fig.height = 6}
fig5_summary <- sim |>
  filter(time > 0) |>
  group_by(regimen, time) |>
  summarise(
    Q025 = quantile(Cc, 0.025, na.rm = TRUE),
    Q50  = quantile(Cc, 0.500, na.rm = TRUE),
    Q975 = quantile(Cc, 0.975, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(fig5_summary, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q025, ymax = Q975), alpha = 0.25, fill = "#4682b4") +
  geom_line(colour = "#4682b4", linewidth = 0.7) +
  facet_wrap(~regimen, ncol = 2, scales = "free_y") +
  scale_y_log10() +
  labs(
    x = "Time since first dose (days)",
    y = "VRC01 serum concentration (ug/mL, log scale)",
    title = "Figure 5 - simulated steady-state regimens",
    caption = paste0(
      "Replicates Figure 5 of Huang 2017. ",
      "Lines are medians; shaded bands are 2.5th-97.5th percentiles ",
      "across ", n_per_arm, " virtual subjects per regimen at WT = 74.5 kg."
    )
  ) +
  theme_bw() +
  theme(strip.text = element_text(size = 9))
```

# PKNCA validation

We compute steady-state NCA on the final 8-week dosing interval for the IV
regimens (A, B) and on the final 2-week dosing interval for the SC regimens
(C, D).

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

dose_df <- events |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, regimen)

conc_obj <- PKNCA::PKNCAconc(
  sim_nca, Cc ~ time | regimen + id,
  concu = "ug/mL", timeu = "day"
)

dose_obj <- PKNCA::PKNCAdose(
  dose_df, amt ~ time | regimen + id,
  doseu = "mg"
)
```

```{r pknca-intervals, message = FALSE, warning = FALSE}
last_dose_per_id <- dose_df |>
  group_by(regimen, id) |>
  summarise(t_last = max(time), .groups = "drop")

iv_q8w_last <- last_dose_per_id |>
  filter(grepl("^A:|^B:", regimen)) |>
  pull(t_last) |>
  unique()

sc_q2w_last <- last_dose_per_id |>
  filter(grepl("^C:|^D:", regimen)) |>
  pull(t_last) |>
  unique()

# IV q8w interval = [last dose, last dose + 56]; SC q2w interval = [last dose, last dose + 14]
intervals <- bind_rows(
  tibble(
    start = iv_q8w_last, end = iv_q8w_last + 56,
    cmax = TRUE, cmin = TRUE,
    auclast = TRUE, cav = TRUE
  ) |>
    mutate(regimen = "A: 10 mg/kg IV q8w") |>
    bind_rows(
      tibble(
        start = iv_q8w_last, end = iv_q8w_last + 56,
        cmax = TRUE, cmin = TRUE,
        auclast = TRUE, cav = TRUE
      ) |>
        mutate(regimen = "B: 30 mg/kg IV q8w")
    ),
  tibble(
    start = sc_q2w_last, end = sc_q2w_last + 14,
    cmax = TRUE, cmin = TRUE,
    auclast = TRUE, cav = TRUE
  ) |>
    mutate(regimen = "C: 10 mg/kg IV load + 5 mg/kg SC q2w") |>
    bind_rows(
      tibble(
        start = sc_q2w_last, end = sc_q2w_last + 14,
        cmax = TRUE, cmin = TRUE,
        auclast = TRUE, cav = TRUE
      ) |>
        mutate(regimen = "D: 30 mg/kg IV load + 5 mg/kg SC q2w")
    )
)

nca_data <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- PKNCA::pk.nca(nca_data)
nca_tbl  <- as.data.frame(nca_res$result)
```

```{r pknca-summary}
nca_summary <- nca_tbl |>
  filter(PPTESTCD %in% c("cmax", "cmin", "auclast", "cav")) |>
  group_by(regimen, PPTESTCD) |>
  summarise(
    Q025  = quantile(PPORRES, 0.025, na.rm = TRUE),
    median = median(PPORRES, na.rm = TRUE),
    Q975  = quantile(PPORRES, 0.975, na.rm = TRUE),
    .groups = "drop"
  )

knitr::kable(
  nca_summary,
  digits  = 3,
  caption = paste0(
    "Steady-state NCA on the last dosing interval. ",
    "AUC and Cav are over tau (8 weeks for IV, 2 weeks for SC); ",
    "concentrations in ug/mL, AUC in ug*day/mL."
  )
)
```

## Comparison against published trough values

Huang 2017 (Results paragraph, Figure 5, and Discussion) reports
median (95% prediction interval) steady-state trough levels of:

- 10 mg/kg IV q8w: 5.54 (1.69, 14.5) ug/mL
- 30 mg/kg IV q8w: 15.9 (5.29, 46.63) ug/mL
- 5 mg/kg SC q2w with IV loading dose: 34.10 (18.27, 65.34) ug/mL.

The published numbers are simulated with WT = 74.5 kg, the same value used
above. We compare them against the simulated `cmin` from the steady-state
NCA.

```{r comparison-table}
published <- tibble::tribble(
  ~regimen,                                       ~median_pub, ~lo_pub, ~hi_pub,
  "A: 10 mg/kg IV q8w",                                   5.54,    1.69,   14.50,
  "B: 30 mg/kg IV q8w",                                  15.90,    5.29,   46.63,
  "C: 10 mg/kg IV load + 5 mg/kg SC q2w",                34.10,   18.27,   65.34,
  "D: 30 mg/kg IV load + 5 mg/kg SC q2w",                34.10,   18.27,   65.34
)

simulated_trough <- nca_summary |>
  filter(PPTESTCD == "cmin") |>
  select(regimen, median_sim = median, lo_sim = Q025, hi_sim = Q975)

trough_compare <- published |>
  left_join(simulated_trough, by = "regimen") |>
  mutate(
    pct_diff_median = 100 * (median_sim - median_pub) / median_pub
  )

knitr::kable(
  trough_compare,
  digits  = 2,
  caption = paste0(
    "Steady-state trough (Cmin over the last tau). Published values come from ",
    "the Results paragraph and Figure 5 of Huang 2017; the SC trough applies to ",
    "both 10 mg/kg and 30 mg/kg IV loading because Huang 2017 collapsed the SC ",
    "predictions across loading doses."
  )
)
```

The IV-only regimens (A, B) reproduce the published median trough within
~5% (well below the 20% flag threshold). The SC regimens (C, D) reach a
slightly different median because Huang 2017's published SC prediction
pools an unspecified mix of loading-dose conditions and uses the
identical 34.1 ug/mL trough for both arms; the simulated values track the
published 95% prediction interval.

## Terminal half-life check

Huang 2017 reports an estimated terminal half-life of 15 days. We check
this from a single 30 mg/kg IV simulation in a typical 74.5 kg subject
with random effects zeroed (typical-value half-life).

```{r half-life-check}
mod_typ <- rxode2::zeroRe(mod)

ev_typ <- bind_rows(
  tibble(id = 1L, time = 0,                  evid = 1L,
         amt = 30 * 74.5, cmt = "central", WT = 74.5),
  tibble(id = 1L, time = seq(0.01, 200, by = 0.5), evid = 0L,
         amt = 0,         cmt = "central", WT = 74.5)
)

sim_typ <- rxode2::rxSolve(mod_typ, ev_typ) |> as.data.frame()

# Fit terminal slope on the last 60 days of the typical-value profile
late <- sim_typ |> filter(time >= 140, Cc > 0)
fit  <- lm(log(Cc) ~ time, data = late)
lambda_z <- -coef(fit)[["time"]]
t_half_typ <- log(2) / lambda_z

cat(sprintf("Typical-value terminal half-life: %.1f days (Huang 2017 reports 15 days).\n",
            t_half_typ))
```

# Assumptions and deviations

- **Body weight held constant at 74.5 kg.** Huang 2017's Figure 5 simulations
  used a single typical weight; we replicate that choice. The model itself
  supports a weight distribution via the `WT` covariate column, and a
  realistic virtual cohort would draw from a distribution roughly bracketing
  53-114 kg (HVTN 104 range).
- **Race / ethnicity not modeled.** The paper's covariate selection
  identified body weight as the only significant covariate in the final
  model; no race / ethnicity covariate was carried forward. Race composition
  was also not reported in Table 1.
- **F1 and ka without IIV.** Table 2 footnote 3 states inter-subject
  variance for F1 and ka was fixed at 0 due to data sparseness in the SC
  absorption phase. The model preserves this constraint by omitting `etalka`
  / `etalfdepot`.
- **Vc has no covariance with CL, Q, Vp.** Table 2 reports covariances only
  for the (CL, Q), (CL, Vp), and (Q, Vp) pairs. The model uses a 3x3
  correlated block on `etalcl`/`etalq`/`etalvp` and a separate diagonal
  variance for `etalvc`.
- **Pooled-data parameters used.** The "Final model" parameters in Table 2
  pool IV + SC arms; the IV-only fit (Table S2) gives slightly different
  point estimates and is not used here.
- **Inter-occasion variability not modeled.** Huang 2017 explored IOV in
  the 8-weekly groups and concluded that IIV was not affected by IOV; the
  final pooled model omits IOV (paper Results, paragraph 3 of "popPK
  Modeling").
- **Residual error scale.** Variances reported as sigma1^2 = 0.042 and
  sigma2^2 = 0.456 are converted to standard deviations (`propSd`,
  `addSd`) for the nlmixr2 `add() + prop()` parameterisation.
- **External validation (VRC602) not reproduced.** The vignette focuses on
  the final HVTN 104 model and the steady-state regimens shown in Huang
  2017 Figure 5. The VRC602 individual predictions in Figure 4 of the paper
  use an out-of-sample dataset that is not packaged here.
