---
title: "MAb PBPK (Shah 2012)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{MAb PBPK (Shah 2012)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: Shah DK, Betts AM. Towards a platform PBPK model to characterize the
  plasma and tissue disposition of monoclonal antibodies in preclinical species
  and human. J Pharmacokinet Pharmacodyn. 2012;39(1):67-86.
  doi:10.1007/s10928-011-9232-2
* Article: <https://doi.org/10.1007/s10928-011-9232-2>

This vignette validates the Shah & Betts 2012 platform PBPK model
parameterized for the human (71 kg male) reference subject. The packaged
file in `inst/modeldb/pharmacokinetics/Shah_2012_mAb_PBPK.R` carries the
full 93-state ODE system (15 tissues x {vascular plasma, vascular blood
cells, endosomal unbound mAb, endosomal FcRn-bound mAb, endosomal free
FcRn, interstitial} plus central plasma, central blood cells, and lymph
node).

# Population

The platform model was fit simultaneously to 52 plasma and tissue
concentration profiles spanning four species (mouse, rat, monkey,
human) and several mAbs (control IgG, the murine anti-platelet mAb 7E3,
MOPC21, the chimeric anti-CEA cT84.66, the human anti-HIV OST577, and
adalimumab). Four system parameters (FcRn endosomal concentration,
pinocytosis rate per endosomal volume, endosomal degradation rate, and
the lymph-node-to-plasma transfer scaling C_LNLF) were estimated from
the pooled cross-species data; tissue volumes, plasma and blood-cell
flows, vascular reflection coefficients, FR (the fraction of FcRn-bound
mAb recycled to vascular space), and the FcRn-IgG association /
dissociation rate constants were fixed from prior literature (Shah &
Betts 2012, Tables 1-4 and text p.73).

The human-specific validation cohort is the 5 mg/kg single-dose IV arm
of Weisman et al. 2003 (n=6 adults with rheumatoid arthritis on
concomitant methotrexate, plasma adalimumab concentration vs time, used
in Shah & Betts 2012 Fig 8c).

```{r read-model}
mod <- readModelDb("Shah_2012_mAb_PBPK")
```

# Source trace

The packaged model is a literal transcription of the equations and
parameter values from Shah & Betts 2012. Per-parameter source comments
appear inline in the model file; the table below summarizes the
provenance.

| Equation / parameter | Value (human) | Source |
|----------------------|---------------|--------|
| Eq 1: central plasma mass balance | derivation | Shah & Betts 2012 Eq 1 |
| Eq 2: central blood cells | derivation | Eq 2 |
| Eq 3: lymph node | derivation | Eq 3 |
| Eq 4: typical-tissue vascular plasma | derivation | Eq 4 |
| Eq 5: tissue vascular blood cells | derivation | Eq 5 |
| Eq 6: endosomal unbound mAb | derivation | Eq 6 |
| Eq 7: endosomal FcRn-bound mAb | derivation | Eq 7 |
| Eq 8: endosomal free FcRn | derivation | Eq 8 |
| Eq 9: interstitial mAb (Ag terms zeroed) | derivation | Eq 9 |
| Eq 11: liver vascular plasma (portal inflows) | derivation | Eq 11 |
| Eq 12: liver vascular blood cells (portal inflows) | derivation | Eq 12 |
| FcRn endosomal concentration | 4.98e-5 mol/L | Table 6 |
| CLup per endosomal volume | 0.0366 L/h/L | Table 6 |
| Kdeg (endosomal mAb degradation) | 42.9 1/h | Table 6 |
| C_LNLF | 9.1 (unitless) | Table 6 |
| Kon (FcRn-IgG association, human) | 5.59e8 1/M/h | text p.73 |
| Koff (FcRn-IgG dissociation, human) | 23.9 1/h | text p.73 |
| FR (recycle fraction to vascular space) | 0.715 | text p.73 (Garg & Balthasar 2007) |
| sigma_IS (lymph reflection, all tissues) | 0.2 | text p.73 |
| Tissue volumes (V_total, V_plasma, V_BC, V_E, V_IS, V_C) | per Table 4 | Table 4 |
| Tissue plasma flows (PLQ_i) | per Table 4 | Table 4 |
| Tissue blood-cell flows (BCQ_i) | per Table 4 | Table 4 |
| Tissue vascular reflection coefficients (sigma_V) | 0.85 to 0.99 by tissue | text p.73 |
| Lymph flows (L_i = PLQ_i / 500) | per tissue | text p.73 (Swartz 2001) |
| Lymph node return flow (L_LymphNode) | 3.670 L/h | Table 4 row Ly. Node |

## Units in the ODE system

| Quantity | Units |
|----------|-------|
| State amounts (mAb in vp, bc, eu, eb, is) | nmol |
| State amounts (free FcRn in fr_*) | nmol (1:1 stoichiometry with mAb) |
| Concentrations derived as state / volume | nmol/L (= nM) |
| Time | h |
| Volumes | L |
| Plasma / blood-cell / lymph flows | L/h |
| Kon | 1/(M h) (converted internally to 1/(nM h) via x 1e-9) |
| Koff, Kdeg | 1/h |
| CLup_per_v | L/h per L of endosomal volume |
| FcRn baseline | mol/L (converted internally to nM via x 1e9) |

Doses are supplied to the `plasma` compartment in nmol. To translate a
mass dose: `dose_nmol = dose_mg / MW_g_per_mol * 1e6` (e.g., adalimumab
MW 144,190 g/mol gives 5 mg/kg x 71 kg = 2,462 nmol).

# Mass-balance / flux check

With endosomal degradation switched off (Kdeg -> 0), the platform
model is a closed system: total mAb mass should be conserved. This
catches structural ODE errors that would not be visible in a
numerical-only check.

```{r mass-balance}
mod_no_deg <- mod |> ini(lkdeg = log(1e-12))

dose_nmol <- 5 * 71 / 144190 * 1e6
ev <- et(amt = dose_nmol, cmt = "plasma", time = 0) |>
  et(seq(0, 200, by = 4))
sim_no_deg <- as.data.frame(rxSolve(mod_no_deg, events = ev))

organ_suffix <- c("heart", "lung", "muscle", "skin", "adipose",
                  "bone", "brain", "kidney", "liver",
                  "small_intestine", "large_intestine",
                  "pancreas", "thymus", "spleen", "other")
igg_cols <- c(
  "plasma", "bcc", "lnode",
  paste0("vp_", organ_suffix),
  paste0("bc_", organ_suffix),
  paste0("eu_", organ_suffix),
  paste0("eb_", organ_suffix),
  paste0("is_", organ_suffix)
)
sim_no_deg$total_igg <- rowSums(sim_no_deg[, igg_cols])

balance_summary <- sim_no_deg |>
  filter(time %in% c(0, 24, 100, 200)) |>
  transmute(
    time_h          = time,
    total_igg_nmol  = round(total_igg, 3),
    pct_of_dose     = round(100 * total_igg / dose_nmol, 4)
  )
knitr::kable(balance_summary,
             caption = paste0("Total mAb mass with Kdeg = 0. ",
                              "Conservation should hold to machine precision."))
```

The total mAb mass is held at the dose value across the simulation
window when Kdeg is zero, confirming that the convective flows, the
FcRn binding kinetics, the recycling pathway, and the lymph-node
return are all balanced in the implementation.

## Free + bound FcRn conservation per tissue

In each tissue endosomal space, free FcRn + FcRn-bound mAb totals to
the initial free FcRn (equal to `fcrn_M * V_endo_tissue`). Verify in
the heart endosome:

```{r fcrn-conservation}
fcrn_init_heart <- 4.98e-5 * 1e9 * 0.00171  # nmol; v_he_e = 0.00171 L
fcrn_check <- sim_no_deg |>
  filter(time %in% c(0, 24, 100, 200)) |>
  transmute(
    time_h            = time,
    fr_heart          = round(fr_heart, 6),
    eb_heart          = round(eb_heart, 6),
    sum_fr_eb_heart   = round(fr_heart + eb_heart, 6),
    expected          = round(fcrn_init_heart, 6)
  )
knitr::kable(fcrn_check,
             caption = "FcRn mass balance in the heart endosome.")
```

# Adalimumab 5 mg/kg IV - replicates Shah & Betts 2012 Figure 8c

The human plasma profile in Shah & Betts Fig 8c is from Weisman et al.
2003. We simulate the typical-value plasma trajectory and compare its
shape and approximate magnitudes to the published figure.

```{r adalimumab-sim}
ev_long <- et(amt = dose_nmol, cmt = "plasma", time = 0) |>
  et(seq(0, 1008, by = 4))  # 6 weeks of follow-up

sim <- as.data.frame(rxSolve(mod, events = ev_long))

# Convert plasma nM -> mg/L using adalimumab MW for mass-unit comparisons
mw_ada <- 144190  # g/mol
sim$Cc_mg_per_L <- sim$Cc * mw_ada * 1e-6
```

```{r fig-8c, fig.cap = "Replicates Shah & Betts 2012 Figure 8c: human plasma adalimumab after 5 mg/kg IV. Solid line is the platform PBPK model simulation."}
sim |>
  filter(time > 0) |>
  ggplot(aes(time / 24, Cc)) +
  geom_line() +
  scale_y_log10() +
  labs(
    x = "Time (days)",
    y = "Plasma adalimumab (nM)",
    title = "Human plasma adalimumab after 5 mg/kg IV",
    caption = "Replicates Shah & Betts 2012 Figure 8c (panel C)."
  )
```

```{r profile-table}
profile <- sim |>
  filter(time %in% c(0, 8, 24, 72, 168, 336, 504, 672, 1008)) |>
  transmute(
    time_h  = time,
    time_d  = round(time / 24, 1),
    Cc_nM   = round(Cc, 2),
    Cc_mg_per_L = round(Cc_mg_per_L, 3)
  )
knitr::kable(profile,
             caption = "Simulated human plasma adalimumab over 6 weeks.")
```

# PKNCA validation

PKNCA-derived NCA parameters from the simulated typical profile.

```{r pknca}
nca_input <- sim |>
  filter(!is.na(Cc)) |>
  mutate(id = 1L, treatment = "5 mg/kg IV") |>
  select(id, time, Cc, treatment)

dose_df <- data.frame(
  id        = 1L,
  time      = 0,
  amt       = dose_nmol,
  treatment = "5 mg/kg IV"
)

conc_obj <- PKNCA::PKNCAconc(nca_input, Cc ~ time | treatment + id)
dose_obj <- PKNCA::PKNCAdose(dose_df,   amt ~ time | treatment + id)

intervals <- data.frame(
  start      = 0,
  end        = 1008,
  cmax       = TRUE,
  tmax       = TRUE,
  auclast    = TRUE,
  aucinf.obs = TRUE,
  half.life  = 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)
knitr::kable(
  nca_summary[, c("PPTESTCD", "PPORRES")],
  caption = "Simulated NCA parameters for human plasma adalimumab (5 mg/kg IV)."
)
```

## Comparison against published values

The Shah & Betts 2012 paper reports a median percent-prediction-error of
16.0% for the human plasma data set (Table 5), the best of the four
species considered. Specific NCA parameters are not tabulated in the
paper.

Published terminal half-life of adalimumab in healthy and rheumatoid
arthritis subjects after IV dosing is roughly 10-20 days (Weisman 2003,
Nestorov 2014). The simulated terminal half-life from the typical-value
PBPK trajectory is in the lower part of this range, consistent with the
fitted Kdeg = 42.9 1/h being a population-average across multiple mAbs
and species (and the human cohort dominating the half-life signal less
than the longer-lived monkey and rat data).

# Assumptions and deviations

* **No antigen-binding terms.** The packaged model implements the
  nonspecific-mAb / control-IgG version of the Shah & Betts 2012
  framework. Antigen-specific terms (cell-membrane binding,
  antibody-antigen complex internalization, etc., paper Eqs 9-10
  Ag-related terms) are dropped. For tumor or target-antigen
  applications the additional `Kon_Ag`, `Koff_Ag`, `Kint`, and tumor /
  cellular-space ODEs would have to be added. The four estimated
  system parameters (FcRn, CLup, Kdeg, C_LNLF) and the binding-rate
  constants come from the multi-species fit and would not need to be
  re-estimated.
* **Human (71 kg) parameter set only.** The paper additionally reports
  rodent and primate parameter sets (Tables 1-3 plus species-specific
  Kon / Koff). Those are not in this file. To extend, copy the
  function, swap the `v_X_*`, `q_X`, `bcq_X`, `sv_X`, `l_lnode`, `lkon`,
  and `lkoff` numeric values, and rename the function and vignette per
  the file-naming convention.
* **Lung plasma flow derived for mass balance.** Shah & Betts 2012
  Table 4 reports `Q_lu = 181.913 L/h` and `BCQ_lu = 148.838 L/h`, but
  these values do not equal `sum_X Q_X + L_lu` and `sum_X BCQ_X` from
  the same table (about 1.8 to 2.0% high). Using the table values
  literally creates a slow but non-trivial mass leak at the
  lung-arterial junction (lung outputs more mass per unit time than
  the tissues can absorb). To preserve total mAb mass conservation
  the model file derives `Q_lu = sum_X Q_X / (1 - 1/500)`, `L_lu =
  Q_lu / 500`, and `BCQ_lu = sum_X BCQ_X` from the tissue-row values.
  Numerically this is a ~1.8% reduction relative to the Table 4 lung
  flows; the dynamics are essentially unchanged but the closed-system
  mass balance verified above is satisfied exactly.
* **C_LNLF carried as a metadata parameter only.** The paper estimates
  `C_LNLF = 9.1` as the proportionality constant between the lymph-
  node-to-plasma flow `L_LymphNode` and "the plasma flow of the given
  species" (Shah & Betts 2012 p.73). The species-specific
  `L_LymphNode` values are tabulated directly in the rows of Tables
  1-4 (3.670 L/h for human). The model file uses the table value
  `L_LymphNode` as the lymph node return flow and exposes
  `lclnlf` as a separately-estimated parameter for traceability, but
  `lclnlf` does not enter the ODE system.
* **No IIV, no residual error.** Shah & Betts 2012 fits 52 mean
  digitized plasma and tissue profiles with a single point estimate
  per parameter; the variance model (Eq 14) operates on the residual
  between observed and model-predicted concentrations and the paper
  does not report `sigma_intercept` / `sigma_slope`. The packaged
  model is a typical-value mechanistic simulator. For
  estimation use, IIV blocks could be added to `lclup`, `lkdeg`, etc.
  guided by the Table 6 CV% values (CV%=3.48 on `clup`, 15.7 on
  `kdeg`, 11.1 on `fcrn`, > 50 on `clnlf`).
* **Compartment naming deviation from `naming-conventions.md`.** The
  PBPK structure has 93 compartments (15 anatomical tissues x
  {vascular plasma, vascular blood cells, endosomal unbound mAb,
  endosomal FcRn-bound mAb, endosomal free FcRn, interstitial} plus
  central plasma, central blood cells, lymph node). These do not map
  onto the standard `central` / `peripheral1` / `depot` / `effect`
  vocabulary; `checkModelConventions("Shah_2012_mAb_PBPK")` flags
  every PBPK compartment as a non-canonical name. The naming used in
  this file (`vp_<tissue>`, `bc_<tissue>`, `eu_<tissue>`,
  `eb_<tissue>`, `fr_<tissue>`, `is_<tissue>`, plus `plasma`, `bcc`,
  `lnode`) follows the paper's symbolic conventions
  (C^V, C^BC, C^E_unbound, etc.) and is the most readable mapping
  available. No convention change in the rest of nlmixr2lib is
  implied.
