---
title: "Recombinant Factor VIII / von Willebrand factor QSS-TMDD (Vargas Christensen 2019) -- rat"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Recombinant Factor VIII / von Willebrand factor QSS-TMDD (Vargas Christensen 2019) -- rat}
  %\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(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
library(PKNCA)
```

# Model and source

* Citation: Vargas Christensen I, Loftager M, Rode F, Morck Nielsen H,
  Kreilgaard M, Larsen MS. Impact of capacity-limited binding on recombinant
  factor VIII and von Willebrand factor pharmacokinetics in hemophilia A
  rats. *J Thromb Haemost* 2019;17(6):964-974.
* Article: [doi:10.1111/jth.14441](https://doi.org/10.1111/jth.14441)

# Population

The model was developed from 30 hemophilia A rats (mean age 11 weeks, range
8-14; median body weight 300 g, range 200-550) bred at Novo Nordisk A/S
(Maaloev, Denmark). Rats received a single IV bolus dose of recombinant
Factor VIII (rFVIII) into the tail vein at one of four dose levels: 17.5 IU/kg
(n = 8, 5 F / 3 M), 100 IU/kg (n = 6, 4 F / 2 M), 1000 IU/kg (n = 8, 5 F / 3 M),
or 5000 IU/kg (n = 8, 5 F / 3 M). Plasma samples were collected over 48 hours;
the 2-min post-dose samples were excluded from the analysis. Total rFVIII was
measured by a chromogenic FVIII activity assay (Chromogenix, calibrated against
a Turoctocog alfa SRM standard; LLOQ 120 mU/mL). Total VWF was measured by a
novel luminescent oxygen channeling immunoassay (LOCI) calibrated against a
literature-derived 14.94 nmol/L human VWF standard (LLOQ 0.3 nmol/L); no rat
VWF calibrator was available. The rFVIII:VWF complex was measured by a second
LOCI with no calibrator, so its dependent variable is the raw signal in counts
per second (cps).

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

# Source trace

Per-parameter origin is recorded as an in-file comment next to each `ini()`
entry in
`inst/modeldb/specificDrugs/VargasChristensen_2019_rfviii_rat.R`. The table
below collects them in one place for review. All parameter estimates come from
Table 1 of the paper; the NONMEM control stream lives in Data S1 (not on disk).

| Parameter / equation                         | Value (RSE)                    | Source location                    |
|----------------------------------------------|--------------------------------|------------------------------------|
| CL (`lcl`)                                   | 28.7 mL/h (23%) -> 0.0287 L/h  | Table 1                            |
| V  (`lvc`)                                   | 9.9 mL (4%) -> 0.0099 L        | Table 1                            |
| k12 (`lk12`)                                 | 0.032 1/h (21%)                | Table 1                            |
| k21 (`lk21`)                                 | 0.0488 1/h (24%)               | Table 1                            |
| Kss (`lkss`)                                 | 0.142 nmol/L (30%)             | Table 1                            |
| R0 (`lrbase`)                                | 11.5 nmol/L (6%)               | Table 1                            |
| kdeg (`lkdeg`)                               | 0.031 1/h (21%)                | Table 1                            |
| KM,comp (`lkm_comp`)                         | 0.871 nmol/L (39%)             | Table 1                            |
| Vmax,comp (`lvmax_comp`)                     | 0.653 1/h (16%)                | Table 1                            |
| base_cps (`lbase_cps`)                       | 298 cps (3%)                   | Table 1                            |
| gamma (`lhill_cps`)                          | 2.24 (3%)                      | Table 1                            |
| cps50 (`lcps50`)                             | 0.228 nmol/L (15%)             | Table 1                            |
| max_cps (`lmax_cps`)                         | 127,000 cps (13%)              | Table 1                            |
| Allometric exponent on CL (`e_wt_cl`)        | 0.75 (fixed)                   | Table 1 + Results                  |
| IIV CL (%CV)                                 | 57.1% (25); [shrinkage 13]     | Table 1                            |
| IIV R0 (%CV)                                 | 25.6% (13); [shrinkage 7]      | Table 1                            |
| IIV Vmax,comp (%CV)                          | 36.3% (23); [shrinkage 36]     | Table 1                            |
| Proportional error rFVIII (%CV)              | 27.6% (10)                     | Table 1                            |
| Proportional error VWF (%CV)                 | 32.5% (9)                      | Table 1                            |
| Proportional error rFVIII:VWF signal (%CV)   | 35.5% (8)                      | Table 1                            |
| QSS TMDD algebra                             | Complex, Cfree, Rfree closed-form | Methods (Mager & Krzyzanski 2005) |
| Complex Michaelis-Menten elimination         | k_e,comp = Vmax,comp*Cfree/(KM,comp+Cfree) | Methods Eq. (unnumbered)  |
| LOCI 4-PL signal function                    | base_cps + (max_cps-base_cps)*complex^gamma/(cps50^gamma+complex^gamma) | Methods Eq. (unnumbered) |
| VWF steady-state derivation                  | ksyn = kdeg * R0               | Methods                            |
| Reference weight for allometry               | 0.3 kg (median cohort weight; assumed) | Methods (cohort weights)   |

# Loading the model

```{r load}
mod <- readModelDb("VargasChristensen_2019_rfviii_rat")
mod_typical <- rxode2::zeroRe(mod)
```

# Virtual cohort

To replicate the dose-ranging visual predictive checks of Figure 3, we simulate
40 typical-value rats at each dose level (dose in IU/kg is converted to nmol
per rat using the paper's median body weight of 300 g and the FVIII activity
specific activity of 5000 IU per mg of rFVIII protein, with a molecular weight
of 170 kDa for B-domain-deleted rFVIII; the exact conversion is only needed to
produce concentration-time profiles on the paper's nmol/L scale and does not
affect any of the fitted parameters).

```{r cohort}
median_wt_kg <- 0.3
mg_per_iu    <- 1 / 5000                     # rFVIII specific activity (approx)
mw_rfviii    <- 170e3                        # BDD-rFVIII molecular weight (g/mol)
iu_per_kg_to_nmol_per_rat <- function(dose_iu_per_kg, wt_kg = median_wt_kg) {
  dose_iu   <- dose_iu_per_kg * wt_kg
  dose_mg   <- dose_iu * mg_per_iu
  dose_nmol <- dose_mg / mw_rfviii * 1e6     # mg / (g/mol) * 1e6 = nmol
  dose_nmol
}

dose_levels_iu_per_kg <- c(17.5, 100, 1000, 5000)
dose_levels_nmol       <- vapply(dose_levels_iu_per_kg, iu_per_kg_to_nmol_per_rat, numeric(1L))
data.frame(dose_iu_per_kg = dose_levels_iu_per_kg, dose_nmol_per_rat = dose_levels_nmol)
```

# Simulation helper

```{r helper}
make_events <- function(dose_nmol, wt_kg = median_wt_kg, tmax = 48, dt = 0.25) {
  times   <- seq(0, tmax, by = dt)
  # Multi-output model: dose goes to the central ODE state; observations
  # sit on any of the three algebraic observables (Cc / totalVWF /
  # lociComplex are auto-injected as observation slots by rxode2 after
  # the ODE states), and rxSolve returns all three as columns at each
  # observation time.
  obs     <- data.frame(id = 1L, time = times, evid = 0L, amt = 0,
                        cmt = "Cc", WT = wt_kg)
  dose_row <- data.frame(id = 1L, time = 0, evid = 1L, amt = dose_nmol,
                         cmt = "central", WT = wt_kg)
  dplyr::bind_rows(dose_row, obs)
}
```

# Replicating Figure 3 -- typical-value dose-ranging profiles

Figure 3 of the paper shows the visual predictive checks of total rFVIII,
total VWF, and the rFVIII:VWF LOCI signal over 0-24 h for each dose level. We
replicate the typical-value trajectories below (VPC bands would require sampling
IIV; see `Stochastic simulation with IIV` further down).

```{r figure3-simulate}
sim_typ <- lapply(seq_along(dose_levels_iu_per_kg), function(k) {
  ev  <- make_events(dose_levels_nmol[k])
  sol <- as.data.frame(rxode2::rxSolve(mod_typical, ev))
  sol$dose_iu_per_kg <- dose_levels_iu_per_kg[k]
  sol
}) |> dplyr::bind_rows() |>
  dplyr::mutate(dose_label = factor(paste0(dose_iu_per_kg, " IU/kg"),
                                    levels = paste0(dose_levels_iu_per_kg, " IU/kg")))
```

```{r figure3a-rfviii}
ggplot(sim_typ, aes(time, Cc, colour = dose_label)) +
  geom_line(linewidth = 0.7) +
  scale_y_log10() +
  labs(x = "Time (h)", y = "Total rFVIII (nmol/L, log scale)",
       colour = "Dose",
       title = "Total rFVIII plasma concentration vs. time -- typical value",
       caption = "Replicates Figure 3A of Vargas Christensen 2019 (typical-value trajectories, 0-24 h subset).") +
  coord_cartesian(xlim = c(0, 24))
```

```{r figure3b-vwf}
ggplot(sim_typ, aes(time, totalVWF, colour = dose_label)) +
  geom_line(linewidth = 0.7) +
  geom_hline(yintercept = 11.5, linetype = "dotted", colour = "grey40") +
  labs(x = "Time (h)", y = "Total VWF (nmol/L)",
       colour = "Dose",
       title = "Total VWF plasma concentration vs. time -- typical value",
       caption = paste0(
         "Replicates Figure 3B of Vargas Christensen 2019 (typical-value trajectories, 0-24 h subset).\n",
         "Dotted line at 11.5 nmol/L = estimated baseline R0. High-dose rFVIII transiently ",
         "reduces total VWF then rebounds."
       )) +
  coord_cartesian(xlim = c(0, 24))
```

```{r figure3c-loci}
ggplot(sim_typ, aes(time, lociComplex, colour = dose_label)) +
  geom_line(linewidth = 0.7) +
  scale_y_log10() +
  labs(x = "Time (h)", y = "rFVIII:VWF LOCI signal (cps, log scale)",
       colour = "Dose",
       title = "rFVIII:VWF complex LOCI signal vs. time -- typical value",
       caption = "Replicates Figure 3C of Vargas Christensen 2019 (typical-value trajectories, 0-24 h subset).") +
  coord_cartesian(xlim = c(0, 24))
```

# Replicating Figure 2 -- rFVIII:VWF complex half-life vs. free rFVIII

Figure 2 of the paper shows how the estimated elimination half-life of the
rFVIII:VWF complex depends on the unbound plasma concentration of rFVIII. At
Cfree = 0 the complex has no direct elimination pathway (k_e,comp = 0) and its
effective loss reduces to the free-VWF turnover (`ln(2) / kdeg` = ~22 h). At
saturating Cfree the half-life approaches `ln(2) / Vmax,comp` = ~1 h.

```{r figure2}
cfree_grid <- 10 ^ seq(-3, 1.5, length.out = 200)
vmax_comp  <- 0.653
km_comp    <- 0.871
kdeg       <- 0.031
k_e_comp   <- vmax_comp * cfree_grid / (km_comp + cfree_grid)
# Effective complex loss combines direct k_e,comp with the free-VWF turnover
# that drains the complex through QSS re-equilibration.
effective_k <- k_e_comp + kdeg
thalf       <- log(2) / effective_k

fig2 <- data.frame(cfree = cfree_grid, thalf = thalf)

ggplot(fig2, aes(cfree, thalf)) +
  geom_line(linewidth = 0.7) +
  scale_x_log10() +
  labs(x = "Unbound rFVIII (nmol/L, log scale)",
       y = "Estimated rFVIII:VWF complex half-life (h)",
       title = "Complex half-life vs. unbound rFVIII -- typical value",
       caption = "Replicates Figure 2 of Vargas Christensen 2019 (analytical, typical-value).")
```

# PKNCA validation of total rFVIII

Population NCA of the simulated total rFVIII concentration profile provides a
model-vs.-model sanity check: single-dose Cmax scales with dose, AUC scales
with dose (though less than proportionally at high doses because saturation of
VWF binding increases the unbound fraction and CL of unbound rFVIII), and
apparent half-life is dose-dependent because of the same nonlinearity. The
paper does not report NCA-derived parameters, so no side-by-side comparison
table is rendered.

```{r pknca}
pk_input <- sim_typ |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::mutate(id = as.integer(dose_iu_per_kg),  # one virtual subject per arm
                treatment = as.character(dose_iu_per_kg))

dose_input <- data.frame(id        = as.integer(dose_levels_iu_per_kg),
                         treatment = as.character(dose_levels_iu_per_kg),
                         time      = 0,
                         amt       = dose_levels_nmol)

conc_obj <- PKNCA::PKNCAconc(pk_input, Cc ~ time | treatment + id,
                             concu = "nmol/L", timeu = "h")
dose_obj <- PKNCA::PKNCAdose(dose_input, amt ~ time | treatment + id,
                             doseu = "nmol")
intervals <- data.frame(start = 0, end = Inf,
                         cmax = TRUE, tmax = TRUE, half.life = TRUE,
                         auclast = TRUE, aucinf.obs = TRUE)
data_obj <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- PKNCA::pk.nca(data_obj)
nca_summary <- as.data.frame(nca_res$result) |>
  dplyr::mutate(PPTESTCD = as.character(PPTESTCD)) |>
  dplyr::filter(PPTESTCD %in% c("cmax", "tmax", "half.life", "auclast", "aucinf.obs")) |>
  dplyr::select(treatment, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  dplyr::rename(
    "Dose (IU/kg)"         = treatment,
    "Cmax (nmol/L)"        = cmax,
    "Tmax (h)"             = tmax,
    "AUC0-last (h*nmol/L)" = auclast,
    "AUC0-inf (h*nmol/L)"  = aucinf.obs,
    "t1/2 (h)"             = half.life
  )
knitr::kable(nca_summary,
             caption = "PKNCA-derived NCA parameters of typical-value simulated total rFVIII.",
             digits  = 3)
```

# Stochastic simulation with IIV

The IIV magnitudes reported in Table 1 (CL 57.1% CV, R0 25.6% CV, Vmax,comp
36.3% CV; the last with 36% shrinkage) can be sampled to produce a stochastic
cohort. The example below simulates 40 rats per dose level. Given the cap on
per-arm cohort sizes documented in the skill (200 subjects per arm), 40 rats
per arm is well within budget and reproduces the between-subject variability
shown in the paper's VPC bands qualitatively.

```{r vpc-simulate}
set.seed(42)
per_arm     <- 40L
sim_stoch <- lapply(seq_along(dose_levels_iu_per_kg), function(k) {
  # Reuse a single set of observation times for all IDs in this arm.
  ev1 <- make_events(dose_levels_nmol[k])
  ev  <- lapply(seq_len(per_arm), function(i) {
    e <- ev1
    e$id <- i
    e
  }) |> dplyr::bind_rows()
  sol <- as.data.frame(rxode2::rxSolve(mod, ev))
  sol$dose_iu_per_kg <- dose_levels_iu_per_kg[k]
  sol
}) |> dplyr::bind_rows() |>
  dplyr::mutate(dose_label = factor(paste0(dose_iu_per_kg, " IU/kg"),
                                    levels = paste0(dose_levels_iu_per_kg, " IU/kg")))

sim_stoch_summary <- sim_stoch |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::group_by(dose_label, time) |>
  dplyr::summarise(
    median = stats::median(Cc),
    lo     = stats::quantile(Cc, 0.05),
    hi     = stats::quantile(Cc, 0.95),
    .groups = "drop"
  )

ggplot(sim_stoch_summary, aes(time, median, colour = dose_label, fill = dose_label)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.2, colour = NA) +
  geom_line(linewidth = 0.7) +
  scale_y_log10() +
  labs(x = "Time (h)", y = "Total rFVIII (nmol/L, log scale)",
       colour = "Dose", fill = "Dose",
       title = "Total rFVIII -- median (5%-95%) of a 40-rat/arm stochastic cohort",
       caption = "Ribbon shows the 5-95% range across simulated rats; complements Figure 3A of the paper.") +
  coord_cartesian(xlim = c(0, 24))
```

# Assumptions and deviations

* **Reference weight for allometry.** The paper reports that CL scales with body
  weight at fixed allometric exponent 0.75 but does not print the reference
  weight explicitly. We use 0.3 kg (the reported cohort median weight); if the
  NONMEM code in Data S1 uses a different reference (e.g. 0.25 kg standard
  laboratory rat), simulated typical CL scales by
  `(0.3 / actual_ref)^0.75`, but the overall shape of every profile
  (including all nonlinearities) is unchanged.
* **IU-to-nmol dose conversion.** Doses were reported in IU/kg. The Model
  disposition and observations are computed in nmol / nmol/L. We use a
  molecular weight of 170 kDa (B-domain-deleted rFVIII) and a specific
  activity of 5000 IU per mg of rFVIII protein to convert IU/kg to nmol per
  rat. Neither the exact molecular weight nor the specific activity is
  specified in the paper; using different assumed values changes only the
  x-axis scaling of dose (and the corresponding concentration axis linearly)
  without altering the ratio of the four dose-arm profiles.
* **Sex distribution.** The paper reports 19/30 female rats (63.3%); no sex
  effect was retained. Neither sex nor age was carried into the model.
* **VWF calibrator.** The paper's own Discussion flags that the estimated R0
  (11.5 nmol/L) is likely a distortion caused by using a human VWF preparation
  as calibrator in a rat assay -- endogenous rat VWF may be higher. This is
  a discussion-level qualitative caution; the R0 estimate is left at the
  reported value.
* **rFVIII:VWF complex calibrator.** No calibrator was available for the
  rFVIII:VWF complex LOCI assay, so the LOCI signal is the observation
  dependent variable rather than a concentration. Absolute complex
  concentrations predicted by the model (used inside `lociComplex`) are model-
  interpreted only and should not be treated as directly measured.
* **VPC vs. typical-value figures.** The paper's Figure 3 shows VPCs (data +
  simulation-based 95% CIs). The typical-value replications above show only
  the median trajectory; the stochastic-cohort section shows the 5-95%
  variability under sampled IIV.
* **Complex effective-half-life plot (Figure 2).** The paper's Figure 2
  reports "estimated elimination half-life of rFVIII:VWF complex" as a
  function of the "estimated unbound concentration of rFVIII". We reproduce
  it by combining the direct Michaelis-Menten k_e,comp with the free-VWF
  kdeg (the two mechanisms that drain the complex under QSS re-equilibration).
  The exact analytic Figure 2 that the paper used may differ in detail;
  either form gives ~22 h at Cfree = 0 and ~1 h at saturating Cfree, matching
  the paper's stated range.
