Zierhut_2008_osteoprotegerin <- function() {
  description <- "Population PK/PD model for Fc-osteoprotegerin (Fc-OPG, AMG 162 / AMGN-0007 precursor) in healthy postmenopausal women (Zierhut 2008): two-peripheral-compartment IV/SC PK with parallel linear and Michaelis-Menten elimination from the central compartment, first-order absorption from the SC depot with a logistic-style bioavailability F = FSC / (1 + FSC), and an indirect-response biomarker turnover model for urinary N-telopeptide (uNTX) where Fc-OPG inhibits bone-resorption-driven NTX synthesis via an Imax = 1, IC50 sigmoidal Hill term."
  reference <- paste(
    "Zierhut ML, Gastonguay MR, Martin SW, Vicini P, Bekker PJ,",
    "Holloway D, Leese PT, Peterson MC. (2008).",
    "Population PK-PD model for Fc-osteoprotegerin in healthy",
    "postmenopausal women. J Pharmacokinet Pharmacodyn 35(4):379-99.",
    "doi:10.1007/s10928-008-9093-5 (PMID 18633695).",
    "DDMORE Foundation Model Repository: DDMODEL00000233 (scenario 4).",
    sep = " "
  )
  vignette <- "Zierhut_2008_osteoprotegerin"
  units <- list(time = "h", dosing = "mg", concentration = "ng/mL", biomarker = "nM BCE/mM creatinine per hr")

  ddmore_id    <- "DDMODEL00000233"
  replicate_of <- NULL

  # Issue #482: what each ODE state holds, in what amount units, in what
  # biological matrix. analyte/specimen proposed by a local model from the
  # model description; units derived from the units block. verified = FALSE
  # means NOT checked against the source paper.
  compartmentData <- list(
    depot       = list(analyte = "Fc-OPG", units = "mg", specimen = "administration site", verified = FALSE),
    central     = list(analyte = "Fc-OPG", units = "mg", specimen = "plasma", verified = FALSE),
    peripheral1 = list(analyte = "Fc-OPG", units = "mg", specimen = "plasma", verified = FALSE),
    peripheral2 = list(analyte = "Fc-OPG", units = "mg", specimen = "plasma", verified = FALSE),
    effect      = list(analyte = "uNTX synthesis inhibition", units = "mg", specimen = "not applicable", verified = FALSE)
  )

  covariateData <- list(
    ROUTE_IV = list(
      description        = "Indicator for intravenous (IV) dosing route",
      units              = "(binary)",
      type               = "binary",
      reference_category = "0 (SC)",
      notes              = "Per-subject dosing-route indicator (1 = IV cohort, 0 = SC cohort) used to switch the PK observation residual SD between the IV-arm value (`CcpropSdIV`) and the SC-arm value (`CcpropSdSC`). Source column `IV` in the DDMODEL00000233 MDL `dataObj` (where the .mdl keys it as `1 in CMT as SC, 2 in CMT as CENT` for the dosing-target encoding and uses `IV == 2` to select `ADDIV` in the residual error piecewise; the canonical 0/1 register form is used here, with `ROUTE_IV = 1` selecting `CcpropSdIV`). Note that the dosing route is also encoded by the rxode2 `cmt` event column (cmt = central for IV, cmt = depot for SC); `ROUTE_IV` is the additional per-subject covariate the residual-error switch needs.",
      source_name        = "IV"
    )
  )

  population <- list(
    n_subjects     = NA_integer_,
    n_studies      = NA_integer_,
    age_range      = NA_character_,
    weight_range   = NA_character_,
    sex_female_pct = 100,
    disease_state  = "Healthy postmenopausal women in a multiple-ascending-dose Phase 1 program of single IV and SC Fc-osteoprotegerin (Fc-OPG, an Fc-fused recombinant osteoprotegerin construct; the precursor program to denosumab). Detailed cohort sizes, weight, and age distributions are not reproduced in the DDMORE bundle for DDMODEL00000233 and the original Zierhut 2008 publication is not on disk in this worktree.",
    dose_range     = "Single IV (3 mg/kg in the bundle's simulation) or SC (3 mg/kg = 210 mg for a 70 kg subject in the bundle's vignette) Fc-OPG, with the .mdl piecewise residual error capturing distinct IV vs SC assay variability. The bundle's `Simulated_opg.txt` is generated by `Vignette_opg.R` for a single 3 mg/kg SC dose (amt = 210 mg, cmt = 1, ROUTE_IV = 0) in 50 typical subjects under `mrgsolve::zero_re(sigma)` (typical-value simulation, no IIV / no residual error).",
    regions        = NA_character_,
    notes          = "Population demographic detail (n_subjects, weight range, age range, race) is not reproduced in the DDMORE Foundation Model Repository bundle for DDMODEL00000233. The `Model_Accommodations.txt` file ships only the citation and asserts 'There are no model differences' relative to Zierhut 2008. The original publication (Zierhut et al. 2008, J Pharmacokinet Pharmacodyn) is not on disk in this worktree and a paper-side cross-check of the parameter values is therefore not in scope of this extraction; see the validation vignette's Errata section for the bundle-only provenance trail."
  )

  ini({
    # Structural PK / PD parameter values are taken directly from the DDMODEL00000233
    # MDL `parObj` STRUCTURAL block (Fc_opg_uNTx_PKPD.mdl lines 32-46), mirrored without
    # change in the mrgsolve [PARAM] block of opg.cpp (lines 17-30). The bundle does not
    # ship a NONMEM Output_real_*.lst (the executable is mrgsolve, not NONMEM), so
    # there is no FINAL PARAMETER ESTIMATE block; the Model_Accommodations.txt asserts
    # "There are no model differences" relative to Zierhut 2008, so these are the
    # publication-derived point values that drive the bundle's simulation rather than
    # refitted final estimates.

    # Structural - PK
    lka  <- log(0.0131); label("First-order SC absorption rate ka (1/h)")                # TVKA   = 0.0131  in .mdl STRUCTURAL
    lcl  <- log(168);    label("Linear (first-order) clearance CL (mL/h)")               # TVCL   = 168     in .mdl STRUCTURAL
    lvc  <- log(2800);   label("Central volume of distribution Vc (mL)")                 # TVVC   = 2800    in .mdl STRUCTURAL
    lvp  <- log(443);    label("First peripheral volume Vp (mL)")                        # TVVP1  = 443     in .mdl STRUCTURAL
    lvp2 <- log(269);    label("Second peripheral volume Vp2 (mL)")                      # TVVP2  = 269     in .mdl STRUCTURAL
    lq   <- log(15.5);   label("First inter-compartmental clearance Q (mL/h)")           # TVQ1   = 15.5    in .mdl STRUCTURAL
    lq2  <- log(3.02);   label("Second inter-compartmental clearance Q2 (mL/h)")         # TVQ2   = 3.02    in .mdl STRUCTURAL
    lvmax <- log(13300); label("Michaelis-Menten Vmax for nonlinear elimination (ng/h)") # TVVMAX = 13300   in .mdl STRUCTURAL
    lkm  <- log(6.74);   label("Michaelis-Menten Km for nonlinear elimination (ng/mL)")  # TVKM   = 6.74    in .mdl STRUCTURAL
    lfdepot <- log(0.0719); label("SC bioavailability THETA (FSC, used as F = FSC / (1 + FSC); raw FSC is unitless)")  # TVFSC  = 0.0719  in .mdl STRUCTURAL

    # Structural - PD biomarker (uNTX) turnover
    lksyn <- log(0.864);  label("Biomarker (uNTX) zero-order synthesis rate ksyn (nM BCE/mM Cr per h)")  # TVKSYN = 0.864   in .mdl STRUCTURAL
    lkdeg <- log(0.0204); label("Biomarker (uNTX) first-order degradation rate kdeg (1/h)")              # TVKDEG = 0.0204  in .mdl STRUCTURAL
    lic50 <- log(5.38);   label("Half-maximal Fc-OPG concentration suppressing uNTX synthesis IC50 (ng/mL)") # TVIC50 = 5.38    in .mdl STRUCTURAL

    # Inter-individual variability. The DDMORE bundle ships two parallel encodings
    # of the IIV: an MDL `parObj` VARIABILITY block (.mdl lines 47-63) where each
    # PPV is annotated as either a variance (`var = ...`) or a standard deviation
    # (`sd = ...`) via the `Normal(mean = 0, var | sd = ...)` form in the
    # RANDOM_VARIABLE_DEFINITION block (.mdl lines 107-117); and an mrgsolve
    # `[OMEGA]` block in opg.cpp (lines 54-66) that copies the same numeric
    # values but interprets them all as variances (mrgsolve's `[OMEGA]`
    # convention). The two encodings are numerically inconsistent for the
    # `sd = ...` PPVs (PPV_KA, PPV_FSC), but the bundle's `Simulated_opg.txt`
    # is generated by `Vignette_opg.R` running mrgsolve on opg.cpp, so the
    # as-simulated values are variances. Matching the mrgsolve interpretation
    # below preserves self-consistency with the shipped simulated dataset
    # (verified by reproducing the bundle's per-time-point median Cc and NTX
    # within ~10% across 200 IIV-included subjects); see the validation
    # vignette's Errata section for a fuller discussion.
    #
    #   PPV_CL  = 0.0391   variance (mrgsolve [OMEGA])
    #   PPV_VC  = 0.0102   variance
    #   PPV_VP1 = 0.0144   variance
    #   PPV_VP2 = 0.0333   variance
    #   PPV_Q1  = 0.0379   variance
    #   PPV_KA  = 0.0457   variance (.mdl annotates as `sd` but bundle simulates
    #                                with mrgsolve [OMEGA] = variance)
    #   PPV_FSC = 0.263    variance (same caveat as PPV_KA)
    etalcl     ~ 0.0391
    etalvc     ~ 0.0102
    etalvp     ~ 0.0144
    etalvp2    ~ 0.0333
    etalq      ~ 0.0379
    etalka     ~ 0.0457
    etalfdepot ~ 0.263

    # Block IIV on (KSYN, KDEG, IC50). The .mdl declares this as a 3x3 covariance
    # matrix via `PPV_KSYN_KDEG_IC50 : {matrixValue = triangle([...], 3, true)}`,
    # using the lower-triangular row-by-row layout
    #   row 1: 0.281                       (var KSYN)
    #   row 2: 0.0867 0.0325               (cov KSYN-KDEG, var KDEG)
    #   row 3: 0.000  0.000   1.18         (cov KSYN-IC50, cov KDEG-IC50, var IC50)
    # which gives a high KSYN-KDEG correlation (rho = 0.0867 / sqrt(0.281 * 0.0325)
    # = 0.908) and zero correlation of either with IC50. nlmixr2's correlated-IIV
    # `eta1 + eta2 + eta3 ~ c(...)` form takes the same lower-triangular ordering.
    etalksyn + etalkdeg + etalic50 ~ c(0.281,
                                       0.0867, 0.0325,
                                       0,      0,      1.18)

    # Residual error. The bundle ships two parallel encodings: the .mdl
    # OBSERVATION block (.mdl lines 156-161) declares the additive /
    # proportional parameters as SD coefficients of a unit-variance Gaussian
    # (`eps ~ Normal(mean = 0, var = 1)`, .mdl lines 151-154), while the
    # mrgsolve `[SIGMA]` block in opg.cpp (lines 68-74) copies the same
    # numeric values but interprets them as variances (mrgsolve `[SIGMA]`
    # convention). The bundle's `Simulated_opg.txt` is generated under
    # `mrgsolve::zero_re(sigma)`, so the two encodings produce the same
    # simulated dataset (sigma is zeroed); however the as-encoded simulator
    # interpretation is mrgsolve = variance, matching the IIV interpretation
    # above. SDs are therefore taken as the square root of the listed values.
    #
    # PK observation:
    #   .mdl: type is additiveError, trans is ln,
    #         additive = piecewise{{ ADDIV when IV==2; otherwise ADDSC }}
    #   nlmixr2: log-additive (= proportional in linear space) with a route-switched
    #            SD; ROUTE_IV = 1 selects sqrt(ADDIV), ROUTE_IV = 0 selects sqrt(ADDSC).
    CcpropSdIV <- 0.13892;  label("Log-scale residual SD on Cc, IV cohort (= sqrt(0.0193); proportional CV ~ 14%)")    # sqrt(ADDIV  = 0.0193) in .mdl VARIABILITY (selected when IV==2)
    CcpropSdSC <- 0.85616;  label("Log-scale residual SD on Cc, SC cohort (= sqrt(0.7330); proportional CV ~ 86%)")    # sqrt(ADDSC  = 0.7330) in .mdl VARIABILITY (selected otherwise)

    # PD observation:
    #   .mdl: type is combinedError1, additive = PDADD, proportional = PDPROP
    #   The combined1 form (linear sum of additive + proportional SDs) is the
    #   `add() + prop() + combined1()` triple in nlmixr2.
    propSd_NTX <- 0.20174;  label("Proportional residual SD on NTX (= sqrt(0.0407); fraction)")                       # sqrt(PDPROP = 0.0407) in .mdl VARIABILITY
    addSd_NTX  <- 4.5497;   label("Additive residual SD on NTX (= sqrt(20.7); nM BCE/mM creatinine per h)")          # sqrt(PDADD  = 20.7  ) in .mdl VARIABILITY
  })

  model({
    # ------------------------------------------------------------------
    # 1. Individual PK / PD parameters (.mdl `mdlObj` INDIVIDUAL_VARIABLES,
    #    lines 120-134; equivalent [MAIN] block in opg.cpp lines 35-52).
    #    Q2, Vmax, and Km have no IIV in the .mdl source.
    # ------------------------------------------------------------------
    ka     <- exp(lka     + etalka)
    cl     <- exp(lcl     + etalcl)
    vc     <- exp(lvc     + etalvc)
    vp     <- exp(lvp     + etalvp)
    vp2    <- exp(lvp2    + etalvp2)
    q      <- exp(lq      + etalq)
    q2     <- exp(lq2)
    fdepot <- exp(lfdepot + etalfdepot)
    ksyn   <- exp(lksyn   + etalksyn)
    kdeg   <- exp(lkdeg   + etalkdeg)
    ic50   <- exp(lic50   + etalic50)
    vmax   <- exp(lvmax)
    km     <- exp(lkm)

    # ------------------------------------------------------------------
    # 2. Concentration in central compartment.
    #    State amounts are in mg (matches the .mdl dosing target: AMT in mg
    #    on CMT 1 (SC) or CMT 2 (CENT)); volumes are in mL. Multiplying
    #    central(mg) / vc(mL) by 1e6 gives concentration in ng/mL, which
    #    is the unit used by Km, IC50, and the published Cc reporting.
    #    Equivalent to opg.cpp's `#define CP (CENT/(VC/1000000.0))`.
    # ------------------------------------------------------------------
    Cc <- (central / vc) * 1e6

    # ------------------------------------------------------------------
    # 3. Nonlinear elimination clearance term CLNL = Vmax / (Cc + Km).
    #    Vmax is in ng/h, Cc and Km in ng/mL, so CLNL is in mL/h, summing
    #    correctly with the linear CL and inter-compartmental Q / Q2 terms
    #    in the central-compartment ODE below.
    # ------------------------------------------------------------------
    clnl <- vmax / (Cc + km)

    # ------------------------------------------------------------------
    # 4. ODE system. Reproduces the .mdl MODEL_PREDICTION DEQ block
    #    (lines 140-147) and the equivalent [ODE] block in opg.cpp
    #    (lines 76-82). Compartment ordering follows the .mdl/.cpp
    #    [CMT] block: depot = 1 (SC), central = 2 (CENT), peripheral1 =
    #    3 (P1), peripheral2 = 4 (P2), effect = 5 (NTX).
    # ------------------------------------------------------------------
    d/dt(depot)       <- -ka * depot
    d/dt(central)     <-  ka * depot -
                          (cl + q + q2 + clnl) * central / vc +
                          q  * peripheral1 / vp +
                          q2 * peripheral2 / vp2
    d/dt(peripheral1) <-  central * q  / vc - peripheral1 * q  / vp
    d/dt(peripheral2) <-  central * q2 / vc - peripheral2 * q2 / vp2

    # Indirect-response uNTX turnover. Imax is implicitly 1 (.mdl line
    # 146: `KSYN*(1.0 - CP/(IC50+CP))` -- the inhibition fraction reaches
    # unity as Cc -> infinity), and the baseline uNTX is the steady-state
    # ratio ksyn / kdeg (.mdl line 137: `NTX_0 = KSYN / KDEG`).
    ntx0 <- ksyn / kdeg
    effect(0)    <- ntx0
    d/dt(effect) <-  ksyn * (1 - Cc / (ic50 + Cc)) - kdeg * effect

    # ------------------------------------------------------------------
    # 5. SC-depot bioavailability. The .mdl uses a logistic-style
    #    transformation F = FSC / (1 + FSC) (line 138), where FSC is
    #    sampled log-normally with IIV. With the population FSC =
    #    0.0719 this maps to a typical absolute SC F of
    #    0.0719 / 1.0719 = 0.0671 (~6.7%). f(depot) only affects
    #    SC dosing rows (cmt = depot); IV doses go directly to
    #    central with F = 1 by rxode2 default.
    # ------------------------------------------------------------------
    f(depot) <- fdepot / (1 + fdepot)

    # ------------------------------------------------------------------
    # 6. Observations and route-conditional residual error.
    #    PK: log-additive (= proportional) error with SD switched between
    #        the IV and SC cohorts by the per-subject ROUTE_IV indicator.
    #        Same shape as the Cirincione_2017_exenatide multi-study
    #        residual-error switch.
    #    PD: combined1 (linear sum of additive + proportional SDs).
    # ------------------------------------------------------------------
    NTX <- effect

    CcpropSd <- CcpropSdIV * ROUTE_IV + CcpropSdSC * (1 - ROUTE_IV)
    Cc  ~ lnorm(CcpropSd)
    NTX ~ add(addSd_NTX) + prop(propSd_NTX) + combined1()
  })
}
