Shen_2013_belatacept <- function() {
  description <- "PK/PD model for belatacept (CTLA-4/IgG1 fusion protein, selective T-cell co-stimulation blocker) in adult kidney transplant recipients (Shen 2013). The PK side is a one-compartment IV-infusion model derived from the paper's noncompartmental analysis (Table 1, 10 mg/kg substudy, n = 10): typical clearance and volume for a 70 kg adult are set so the model reproduces the reported geometric-mean CL, Vss, AUC over a 4-week dosing interval, and ~8-9 day terminal half-life. The PD side is the inhibitory Emax model of Eq. 2 (Section 3.2, n = 62 in the phase II corticosteroid-avoidance substudy IM103034) describing free CD86 receptor expression on peripheral-blood monocytes (MESF) as E0 - Emax * Cc / (EC50 + Cc); CD86 receptor occupancy is derived as 100 * (E0 - freeCD86) / E0. Belatacept exhibited linear PK across 5-10 mg/kg with relatively low between-subject variability; the full population PK with body-weight covariates was published separately by Zhou et al. (2012) and is not refit here."
  reference <- "Shen J, Townsend R, You X, Shen Y, Zhan P, Zhou Z, Geng D, Wu D, McGirr N, Soucek K, Proszynski E, Pursley J, Masson E. Pharmacokinetics, pharmacodynamics, and immunogenicity of belatacept in adult kidney transplant recipients. Clin Drug Investig. 2014;34(2):117-126 (published online 12 November 2013). doi:10.1007/s40261-013-0153-2"
  vignette <- "Shen_2013_belatacept"
  units <- list(
    time          = "day",
    dosing        = "mg",
    concentration = "ug/mL (belatacept Cc); MESF (free CD86 receptor); percent (CD86 receptor occupancy)"
  )

  # Issue #482: what each ODE state holds, in what amount units, in what
  # biological matrix. Derived mechanically; verified = FALSE means it has
  # NOT been checked against the source paper.
  compartmentData <- list(
    central = list(analyte = "belatacept", units = "mg", specimen = "plasma", verified = FALSE)
  )

  covariateData <- list()

  population <- list(
    species        = "human",
    n_subjects     = 24L,
    n_pd_subjects  = 62L,
    n_studies      = 5L,
    age_range      = NA,
    weight_range   = NA,
    sex_female_pct = NA,
    race_ethnicity = NA,
    disease_state  = "De novo adult kidney transplant recipients (KTRs) receiving kidneys from living, standard-criteria, or extended-criteria donors. All subjects received basiliximab induction plus mycophenolate mofetil and a corticosteroid taper.",
    dose_range     = "Belatacept 5 or 10 mg/kg administered as a 30-min IV infusion. Initial phase: 10 mg/kg loading doses; maintenance phase: 5 mg/kg every 4 weeks. Two phase III regimens (LI and MI) differ only in loading-phase frequency (see Fig. 1 of Shen 2013).",
    regions        = NA,
    notes          = "Pharmacokinetic analyses pool the open-label 10 mg/kg PK study (n = 10 completing week 16; intensive sampling around week 12 dose) with the 5 mg/kg PK substudy of the phase II long-term-extension study IM103100 (n = 14; intensive sampling after first substudy dose). Pharmacodynamic analyses use the phase II corticosteroid-avoidance study IM103034 (n = 62 receiving the MI regimen). Trough-concentration tables draw on the phase III BENEFIT (n = 197-208 LI; n = 162-202 MI) and BENEFIT-EXT (n = 95-150 LI; n = 89-155 MI) studies. Body weight was the only significant covariate identified by Zhou et al. (2012) on CL and V; that covariate effect is not encoded in this model because Shen 2013 itself reports only NCA point estimates and the Zhou paper is the appropriate source for a covariate-aware popPK."
  )

  ini({
    # ------------------------------------------------------------------
    # Structural PK parameters. Values derived from Shen 2013 Table 1
    # (10 mg/kg column, n = 10) for a typical 70 kg adult:
    #   CL  = 0.47 mL/h/kg * 70 kg * 24 h/day / 1000 = 0.7896 L/day
    #   Vss = 0.11 L/kg * 70 kg                     = 7.7 L
    # The 1-compartment approximation (Vc == Vss) is consistent with
    # the paper's observation that Vss is approximately equal to the
    # vascular volume and that exposure is approximately dose-proportional
    # between 5 and 10 mg/kg (Section 3.1, Results).
    # ------------------------------------------------------------------
    lcl <- log(0.79); label("Belatacept clearance CL (L/day) for typical 70 kg adult")   # Shen 2013 Table 1, 10 mg/kg column (0.47 mL/h/kg)
    lvc <- log(7.7);  label("Belatacept central volume Vc (L) for typical 70 kg adult")  # Shen 2013 Table 1, 10 mg/kg column (Vss = 0.11 L/kg)

    # ------------------------------------------------------------------
    # PD parameters from the inhibitory Emax model of Shen 2013 Eq. 2
    # (Section 3.2, IM103034 corticosteroid-avoidance n = 62):
    #   R(C) = E0 - Emax * C / (EC50 + C)
    # Point estimates (with 95% CI) reported in Section 3.2.
    # ------------------------------------------------------------------
    le0   <- log(9820); label("Baseline free CD86 receptor level E0 (MESF)")                     # Shen 2013 Section 3.2 (95% CI 8319-11320)
    lemax <- log(9144); label("Maximal decrease in free CD86 receptors Emax (MESF)")             # Shen 2013 Section 3.2 (95% CI 7662-10619)
    lec50 <- log(2.4);  label("Belatacept concentration for 50% of Emax on CD86 EC50 (ug/mL)")   # Shen 2013 Section 3.2 (95% CI 1.2-3.5)

    # ------------------------------------------------------------------
    # IIV - diagonal. omega^2 = log(1 + CV^2) from the NCA-derived CV%
    # in Shen 2013 Table 1 (10 mg/kg column). The reported CV% lumps
    # between-subject variability with residual variability from the
    # assay; popPK refits in Zhou et al. (2012) attribute the bulk to
    # BSV. No IIV is set on the PD parameters because Shen 2013 reports
    # 95% CIs around the point estimates (parameter uncertainty) but
    # does not tabulate between-subject variances; the paper states only
    # that a non-linear mixed-effect approach with random effect,
    # exponential variance, and exponential spatial correlation was used
    # (Section 2.6).
    # ------------------------------------------------------------------
    etalcl ~ 0.07037  # log(1 + 0.27^2); Shen 2013 Table 1 CL CV% = 27% (10 mg/kg)
    etalvc ~ 0.08618  # log(1 + 0.30^2); Shen 2013 Table 1 Vss CV% = 30% (10 mg/kg)

    # ------------------------------------------------------------------
    # Residual error. Shen 2013 does not separately report PK residual
    # error magnitude (NCA-based reporting only). The PD analysis used
    # exponential heteroscedastic variance with exponential within-
    # subject correlation but the variance estimate is not tabulated.
    # Values below are placeholders that yield plausible scatter for
    # simulation; document in the vignette's Assumptions and deviations.
    # ------------------------------------------------------------------
    propSd          <- 0.10; label("Proportional residual error on belatacept Cc (fraction)")     # placeholder; Shen 2013 reports NCA CV% only
    propSd_freeCD86 <- 0.15; label("Proportional residual error on free CD86 (fraction)")         # placeholder; Shen 2013 Section 2.6 (exponential variance, magnitude not reported)
  })

  model({
    # ------------------------------------------------------------------
    # 1. Individual PK parameters (no covariate effects in this model;
    #    Zhou et al. 2012 is the covariate-aware popPK reference).
    # ------------------------------------------------------------------
    cl <- exp(lcl + etalcl)
    vc <- exp(lvc + etalvc)

    # ------------------------------------------------------------------
    # 2. PD typical-value parameters. No eta on the PD parameters
    #    (see ini() note).
    # ------------------------------------------------------------------
    e0   <- exp(le0)
    emax <- exp(lemax)
    ec50 <- exp(lec50)

    # ------------------------------------------------------------------
    # 3. ODE system. IV infusion or bolus enters `central` directly
    #    via the user's event table (cmt = central; rate or dur set on
    #    the dose row for the 30-min infusion).
    # ------------------------------------------------------------------
    kel <- cl / vc
    d/dt(central) <- -kel * central

    # ------------------------------------------------------------------
    # 4. Observations.
    #    Cc       - belatacept serum concentration (ug/mL).
    #    freeCD86 - free CD86 receptor expression on monocytes (MESF),
    #               inhibitory Emax of Shen 2013 Eq. 2.
    #    ro86     - CD86 receptor occupancy by belatacept (%), derived
    #               from the Emax fit: 100 * (E0 - freeCD86) / E0.
    # ------------------------------------------------------------------
    Cc       <- central / vc
    freeCD86 <- e0 - emax * Cc / (ec50 + Cc)
    ro86     <- 100 * (e0 - freeCD86) / e0

    Cc       ~ prop(propSd)
    freeCD86 ~ prop(propSd_freeCD86)
  })
}
