Jones_2011_PF04878691_oas <- function() {
  description <- paste(
    "Coupled PK + indirect-response pharmacodynamic model for the",
    "2',5'-oligoadenylate synthetase (OAS) gene-expression fold change",
    "during oral PF-04878691 (TLR7 agonist) administration in healthy adult",
    "volunteers (Jones 2011 BJCP). PK is the two-compartment time-varying",
    "clearance model inherited from the companion file",
    "Jones_2011_PF04878691.R (Table 1; all PK structural parameters and",
    "IIVs fixed at the published Table 1 values so the PK forcing function",
    "is the published popPK profile). Drug stimulates the production of OAS",
    "through a power function: dOAS/dt = kin * (1 + slope * Cc^gamma) -",
    "kout * OAS, with baseline OAS fold change rbase = kin / kout so that",
    "kin = rbase * kout (Jones 2011 Methods / Tables 2 OAS model). The",
    "typical Emax indirect-response model could not adequately identify the",
    "parameters given the limited number of dose levels studied, so the",
    "Emax * Cc^gamma / (EC50^gamma + Cc^gamma) drug effect was replaced",
    "with the power function slope * Cc^gamma (Methods 'Population PK-OAS",
    "and PK-lymphocyte models'). The OAS observation is unitless (fold",
    "change from baseline)."
  )
  reference <- paste(
    "Jones HM, Chan PLS, van der Graaf PH, Webster R. Use of modelling and",
    "simulation techniques to support decision making on the progression of",
    "PF-04878691, a TLR7 agonist being developed for hepatitis C.",
    "Br J Clin Pharmacol. 2012;73(1):77-92.",
    "doi:10.1111/j.1365-2125.2011.04047.x.",
    "Companion PK file Jones_2011_PF04878691.R supplies the upstream",
    "two-compartment time-varying clearance PK structure."
  )
  vignette <- "Jones_2011_PF04878691_HCV"
  paper_specific_compartments <- c("oas")
  paper_specific_etas <- c("etalrbase", "etalkout")
  paper_specific_residual_sds <- c("propSd_oas")
  units <- list(time = "h", dosing = "mg", concentration = "ng/mL")

  # 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 = "PF-04878691", units = "mg", specimen = "administration site", verified = FALSE),
    central     = list(analyte = "PF-04878691", units = "mg", specimen = "plasma", verified = FALSE),
    peripheral1 = list(analyte = "PF-04878691", units = "mg", specimen = "plasma", verified = FALSE),
    oas         = list(analyte = "OAS gene-expression fold change", units = "mg", specimen = "not applicable", verified = FALSE)
  )

  covariateData <- list(
    WT = list(
      description        = "Body weight",
      units              = "kg",
      type               = "continuous",
      reference_category = NULL,
      notes              = "Required scaling covariate for the inherited per-kg PK parameters; see Jones_2011_PF04878691.R. Population median 79 kg.",
      source_name        = "WT"
    )
  )

  population <- list(
    species        = "human",
    n_subjects     = 24L,
    n_studies      = 1L,
    age_range      = "21-55 years",
    age_median     = "34 years",
    weight_range   = "57-97 kg",
    weight_median  = "79 kg",
    sex_female_pct = 8,
    race_ethnicity = "Not tabulated in Jones 2011.",
    disease_state  = "Healthy adult volunteers (multiple-dose escalation Phase 1 study; ClinicalTrials.gov NCT00810758).",
    dose_range     = "PF-04878691 administered orally at 3, 6, or 9 mg twice weekly (days 1, 4, 8, 11) for 2 weeks; n = 6 active per dose cohort plus n = 2 placebo per cohort.",
    regions        = "Not specified.",
    notes          = "OAS RNA expression measured by reverse-transcription PCR (Jones 2011 reference [18]); OAS fold change was determined by the maximal observed median change from baseline at a specific time point using an internal control housekeeping gene as a control. Sequential population fit: individual EBE PK parameters from the Table 1 PK model were used as the exposure driver for the OAS fit."
  )

  ini({
    # ----------------------------------------------------------------------
    # PK structural parameters from the companion Jones_2011_PF04878691.R
    # PK file (Table 1). All four structural PK theta entries are wrapped
    # in fixed() because the OAS model was fit sequentially using the
    # individual EBE PK parameters from the Table 1 popPK fit -- the PK
    # parameters are inherited, not re-estimated. Reparameterised from the
    # paper's (CLF, CL0, DEG) into the canonical (cl_exp_inf, cl_time, kdeg)
    # form per the companion PK file.
    # ----------------------------------------------------------------------
    lcl_exp_inf       <- fixed(log(1.7));    label("Steady-state apparent clearance per kg body weight (CL_SS = paper CLF, L/h/kg)")              # Table 1 (CLF = 1.7 L/h/kg)
    lcl_exp_component  <- fixed(log(1.8));    label("Initial offset of the time-varying clearance component per kg (CL_TIME0 = CL0 - CLF, L/h/kg)") # Derived from Table 1 (CL0 = 3.5, CLF = 1.7)
    lcl_exp_kdes     <- fixed(log(0.24));   label("Exponential decay rate of the time-varying clearance component (paper DEG, 1/h)")              # Table 1 (DEG = 0.24 1/h)
    lvc       <- fixed(log(3.3));    label("Apparent central volume of distribution per kg body weight (Vc, L/kg)")                        # Table 1 (Vc = 3.3 L/kg)
    lq        <- fixed(log(0.74));   label("Apparent intercompartmental clearance per kg body weight (Q, L/h/kg)")                         # Table 1 (Q = 0.74 L/h/kg)
    lvp       <- fixed(log(21));     label("Apparent peripheral volume of distribution per kg body weight (Vp, L/kg)")                     # Table 1 (Vp = 21 L/kg)
    lka       <- fixed(log(0.078));  label("First-order absorption rate constant (ka, 1/h)")                                               # Table 1 (ka = 0.078 1/h)

    # PK IIVs from Table 1 -- also fixed because the OAS fit used EBE PK
    # parameters from the upstream PK fit.
    etalcl_exp_inf ~ fixed(0.067)                                                                                                                  # Table 1 (IIV CLF = 0.067)
    etalka ~ fixed(0.19)                                                                                                                   # Table 1 (IIV ka  = 0.19)

    # ----------------------------------------------------------------------
    # OAS indirect-response PD parameters (Jones 2011 Table 2). The
    # baseline OAS fold change rbase together with kout fix the production
    # rate via kin = rbase * kout. Drug stimulates production through the
    # power function slope * Cc^gamma (Methods 'Population PK-OAS' replaces
    # the Emax indirect response by a power function due to limited dose
    # range). gamma is the sigmoidicity exponent on Cc.
    # ----------------------------------------------------------------------
    lkout    <- log(0.034);  label("OAS first-order elimination rate constant kout (1/h)")                                                  # Table 2 (kout = 0.034 1/h, %CV 35)
    lslope   <- log(3.5);    label("Slope of the drug stimulation of OAS production per (ng/mL)^gamma (paper SLP, units 1 / (ng/mL)^gamma)") # Table 2 (SLP = 3.5, %CV 18)
    lrbase   <- log(0.96);   label("OAS baseline fold change (rbase = kin/kout, unitless fold change)")                                     # Table 2 (BASE = 0.96 fold change, %CV 6.7)
    lgamma   <- log(1.6);    label("Sigmoidicity exponent on Cc in the OAS stimulation power function (gamma, unitless)")                   # Table 2 (gamma = 1.6, %CV 10)

    # OAS IIVs from Table 2 (multiplicative exponential random effects on
    # log scale per Methods 'IIV ... multiplicative exponential random
    # effects'). Table 2 reports IIV on kout (omega = 1.7) and BASE
    # (omega = 0.18); no IIV on SLP or gamma in the final model.
    etalkout  ~ 1.7                                                                                                                         # Table 2 (IIV kout = 1.7,  %CV 97)
    etalrbase ~ 0.18                                                                                                                        # Table 2 (IIV BASE = 0.18, %CV 43)

    # OAS residual error (Methods 'proportional error model' for OAS).
    # Table 2 reports residual error = 0.19 (%CV 17).
    propSd_oas <- 0.19; label("OAS proportional residual SD (fraction of fold change)")                                                     # Table 2 (residual error = 0.19, %CV 17)
  })

  model({
    # ----------------------------------------------------------------------
    # Inherited PK (deterministic + per-subject log-normal IIV on cl_exp_inf
    # and ka, exactly as Table 1 reports). See Jones_2011_PF04878691.R for
    # the standalone PK model rationale.
    # ----------------------------------------------------------------------
    cl_exp_inf    <- exp(lcl_exp_inf + etalcl_exp_inf) * WT
    cl_exp_component <- exp(lcl_exp_component) * WT
    cl_exp_kdes  <- exp(lcl_exp_kdes)
    vc       <- exp(lvc) * WT
    q        <- exp(lq)  * WT
    vp       <- exp(lvp) * WT
    ka       <- exp(lka + etalka)

    cl  <- cl_exp_inf + cl_exp_component * exp(-cl_exp_kdes * time)
    kel <- cl / vc
    k12 <- q  / vc
    k21 <- q  / vp

    d/dt(depot)       <- -ka * depot
    d/dt(central)     <-  ka * depot - kel * central - k12 * central + k21 * peripheral1
    d/dt(peripheral1) <-                                k12 * central - k21 * peripheral1

    Cc <- 1000 * central / vc   # PF-04878691 plasma concentration in ng/mL

    # ----------------------------------------------------------------------
    # OAS indirect-response PD (Jones 2011 Methods 'Population PKPD
    # analysis', OAS model). Drug stimulates production through a power
    # function; baseline kin / kout balance sets the OAS state at rbase.
    #     dOAS/dt = kin * (1 + slope * Cc^gamma) - kout * OAS
    #     kin     = rbase * kout
    #     OAS(0)  = rbase
    # ----------------------------------------------------------------------
    rbase     <- exp(lrbase + etalrbase)
    kout      <- exp(lkout  + etalkout)
    slope_oas <- exp(lslope)
    gamma_oas <- exp(lgamma)

    kin_oas <- rbase * kout

    oas(0) <- rbase

    d/dt(oas) <- kin_oas * (1 + slope_oas * Cc^gamma_oas) - kout * oas

    oas ~ prop(propSd_oas)
  })
}
