Simpson_2013_artesunate <- function() {
  description <- "In vitro (P. falciparum). Sigmoid Emax inhibition model of artesunate effect on hypoxanthine uptake by clinical Plasmodium falciparum isolates from the Thai-Myanmar border (Shoklo Malaria Research Unit, 1993-2005), with pfmdr1 genotype covariate effects on EC50. The 'subject' in the NLME framework is a parasite isolate (n=474 isolates with artesunate data). STIM_ARTESUNATE_NM is the per-well drug concentration in the in vitro hypoxanthine-uptake-inhibition assay; the model has no PK and no time evolution. E0 and Emax are fixed per Simpson 2013 Table 3 footnote."
  reference <- paste(
    "Simpson JA, Jamsen KM, Anderson TJC, Zaloumis S, Nair S, Woodrow C, White NJ, Nosten F, Price RN. (2013).",
    "Nonlinear Mixed-Effects Modelling of In Vitro Drug Susceptibility and Molecular Correlates of",
    "Multidrug Resistant Plasmodium falciparum.",
    "PLoS ONE 8(7):e69505.",
    "doi:10.1371/journal.pone.0069505.",
    sep = " "
  )
  vignette <- "Simpson_2013_artesunate"
  units <- list(
    time          = "(unused; each record is one drug-well concentration in a 24+18 h hypoxanthine-uptake-inhibition assay)",
    dosing        = "(no PK dosing; STIM_ARTESUNATE_NM is the in vitro applied well concentration in nM)",
    concentration = "(observation is normalised hypoxanthine uptake fraction E, range 0-1)"
  )

  covariateData <- list(
    STIM_ARTESUNATE_NM = list(
      description        = "Applied artesunate concentration in the in vitro hypoxanthine-uptake-inhibition assay well (nM). Per-record covariate. Doubling-dilution series 0.044 to 87.0 nM, plus drug-free control well (CONC = 0). Drives the sigmoid Emax inhibition of normalised hypoxanthine uptake.",
      units              = "nM",
      type               = "continuous",
      reference_category = NULL,
      notes              = "Source paper used a doubling-dilution series from 87.0 nM down to 0.044 nM (Methods, In vitro Drug Assay), plus drug-free controls. New canonical entry registered alongside this model in inst/references/covariate-columns.md.",
      source_name        = "C"
    ),
    PFMDR1_86Y = list(
      description        = "Plasmodium falciparum pfmdr1 codon-86 tyrosine mutant indicator (1 = single-copy pfmdr1 with 86Y mutation, Simpson 2013 Genotype 2; 0 = otherwise). Time-fixed per isolate.",
      units              = "(binary)",
      type               = "binary",
      reference_category = "0 = Simpson 2013 Genotype 1 single-copy WT 86N/1042N (when all four PFMDR1 indicators are 0)",
      notes              = "Parasite-genome indicator (not host pharmacogenetics). Mutually exclusive with the other three PFMDR1 indicators in the Simpson 2013 cohort. Genotype 2 prevalence for artesunate: 24 of 474 isolates (5.1%) per Table 3.",
      source_name        = "X1"
    ),
    PFMDR1_1042D = list(
      description        = "Plasmodium falciparum pfmdr1 codon-1042 aspartate mutant indicator (1 = single-copy pfmdr1 with 1042D mutation, Simpson 2013 Genotype 3; 0 = otherwise). Time-fixed per isolate.",
      units              = "(binary)",
      type               = "binary",
      reference_category = "0 = Simpson 2013 Genotype 1 single-copy WT (when all four PFMDR1 indicators are 0)",
      notes              = "Mutually exclusive with the other three PFMDR1 indicators. Genotype 3 prevalence for artesunate: 24 of 474 isolates (5.1%) per Table 3.",
      source_name        = "X2"
    ),
    PFMDR1_CN2 = list(
      description        = "Plasmodium falciparum pfmdr1 double-copy amplification indicator (1 = two copies of pfmdr1, all WT 86N/1042N, Simpson 2013 Genotype 4; 0 = otherwise). Time-fixed per isolate.",
      units              = "(binary)",
      type               = "binary",
      reference_category = "0 = Simpson 2013 Genotype 1 single-copy WT (when all four PFMDR1 indicators are 0)",
      notes              = "Mutually exclusive with the other three PFMDR1 indicators. Genotype 4 prevalence for artesunate: 123 of 474 isolates (25.9%) per Table 3.",
      source_name        = "X3"
    ),
    PFMDR1_CN3PLUS = list(
      description        = "Plasmodium falciparum pfmdr1 triple-or-more-copy amplification indicator (1 = three or more copies of pfmdr1, all WT 86N/1042N, Simpson 2013 Genotype 5; 0 = otherwise). Time-fixed per isolate.",
      units              = "(binary)",
      type               = "binary",
      reference_category = "0 = Simpson 2013 Genotype 1 single-copy WT (when all four PFMDR1 indicators are 0)",
      notes              = "Mutually exclusive with the other three PFMDR1 indicators. Genotype 5 prevalence for artesunate: 69 of 474 isolates (14.6%) per Table 3.",
      source_name        = "X4"
    )
  )

  population <- list(
    species        = "in vitro (Plasmodium falciparum clinical isolates from the Shoklo Malaria Research Unit, western Thai-Myanmar border, 1993-2005)",
    n_subjects     = 474,
    n_studies      = 1,
    age_range      = "(not applicable; parasite isolates)",
    weight_range   = "(not applicable; parasite isolates)",
    sex_female_pct = NA_real_,
    disease_state  = "acute P. falciparum malaria",
    dose_range     = "in vitro doubling-dilution series 0.044 to 87.0 nM artesunate, plus drug-free control well (Methods, In vitro Drug Assay)",
    regions        = "western Thai-Myanmar border (Shoklo Malaria Research Unit clinics, Mae Sot, Thailand)",
    n_isolates_artesunate = 474,
    pfmdr1_distribution = "Genotype 1 (single-copy WT 86N/1042N): 234 of 474 (49.4%); Genotype 2 (single-copy 86Y): 24 of 474 (5.1%); Genotype 3 (single-copy 1042D): 24 of 474 (5.1%); Genotype 4 (double-copy WT): 123 of 474 (25.9%); Genotype 5 (triple+ copy WT): 69 of 474 (14.6%). Table 3.",
    observations_per_isolate = "median 22 observations per isolate across all drugs, range 7 to 44 (Results paragraph 1).",
    notes = "n=490 total isolates available across the four-drug study, with artesunate data on 474 isolates (Results paragraph 1; Table 3). Artesunate has the steepest concentration-effect curve in Simpson 2013 (slope gamma = 5.86 vs 2.7-4.1 for the other drugs) and the largest between-isolate variability in slope (SD 0.63 log_e units vs 0.41 for the other three; Table 1). The triple-or-more-copy WT amplification group has the largest relative shift on EC50 for artesunate (127% increase, Table 3). Drug-susceptibility assay = hypoxanthine uptake inhibition (Methods, In vitro Drug Assay; previously described in reference [13])."
  )

  ini({
    # ---------------------------------------------------------------------
    # Sigmoid Emax inhibition (Methods, Eq. 1):
    #   E = Emax - (Emax - E0) * C^gamma / (C^gamma + EC50^gamma)
    # Emax and E0 fixed per Table 3 footnote (#Emax fixed to 0.98;
    # #E0 fixed to 0.01).
    # ---------------------------------------------------------------------
    e0   <- fixed(0.01); label("Minimum normalised hypoxanthine uptake at maximal drug inhibition (fraction, unitless)") # Table 3 footnote (#E0 fixed to 0.01)
    emax <- fixed(0.98); label("Maximum normalised hypoxanthine uptake at zero drug (fraction, unitless)") # Table 3 footnote (#Emax fixed to 0.98)

    # ---------------------------------------------------------------------
    # Population EC50 and slope (Genotype 1 reference). EC50 from Table 3
    # row "Artesunate" Genotype 1 column: 2.3 nM (95% CI 2.1, 2.6). Note
    # this differs slightly from the Table 1 NLME row artesunate EC50 of
    # 2.58 nM (95% reference range 0.44 to 15.32), which is the
    # population-mean EC50 from the no-covariate base model; the Table 3
    # value 2.3 nM is the Genotype 1 (WT reference) population mean from
    # the genotype-covariate model and is the right starting point here.
    # Slope from Table 1 NLME row artesunate: 5.86 (95% reference range
    # 1.70 to 20.10). Slope covariates dropped per skill sidecar.
    # ---------------------------------------------------------------------
    lec50  <- log(2.3);  label("Population log-EC50 for the WT reference parasite (log of nM)") # Table 3 Artesunate Genotype 1 reference
    lgamma <- log(5.86); label("Population log-slope gamma of the sigmoid inhibition curve (log of unitless)") # Table 1 NLME row artesunate

    # ---------------------------------------------------------------------
    # pfmdr1 genotype covariate effects on EC50.
    # ---------------------------------------------------------------------
    e_pfmdr1_86y_ec50     <- -0.17;  label("Proportional effect of single-copy 86Y mutant on EC50 (fraction)")  # Table 3 Artesunate Genotype 2 percent change -17 (-39, 6)
    e_pfmdr1_1042d_ec50   <-  0.38;  label("Proportional effect of single-copy 1042D mutant on EC50 (fraction)") # Table 3 Artesunate Genotype 3 percent change 38 (-27, 102)
    e_pfmdr1_cn2_ec50     <-  0.63;  label("Proportional effect of double-copy WT on EC50 (fraction)")           # Table 3 Artesunate Genotype 4 percent change 63 (35, 92)
    e_pfmdr1_cn3plus_ec50 <-  1.27;  label("Proportional effect of triple-or-more-copy WT on EC50 (fraction)")   # Table 3 Artesunate Genotype 5 percent change 127 (74, 169)

    # ---------------------------------------------------------------------
    # Between-isolate variances. Artesunate is the only drug where the
    # slope variance differs from the other three (0.63 vs 0.41 SD on the
    # log_e scale, Table 1) -- it is the drug with the steepest curve and
    # the most isolate-to-isolate variability in slope.
    # ---------------------------------------------------------------------
    etalec50  ~ 0.67;    # Table 3 footnote: Between-isolate variance for EC50 artesunate = 0.67 (SE 0.048)
    etalgamma ~ 0.3969;  # Table 1 NLME artesunate slope SD = 0.63 log_e units; variance = 0.63^2 = 0.3969

    # ---------------------------------------------------------------------
    # Residual error.
    # ---------------------------------------------------------------------
    propSd <- sqrt(0.025);  label("Proportional residual SD on normalised uptake fraction E (unitless)") # Table 3 footnote: proportional variance 0.025 (SE 0.0037) artesunate; SD = sqrt(0.025) = 0.1581
    addSd  <- sqrt(0.0007); label("Additive residual SD on normalised uptake fraction E (unitless)")     # Table 3 footnote: additive variance 0.0007 (SE 0.0002) artesunate; SD = sqrt(0.0007) = 0.0265
  })

  model({
    # ---------------------------------------------------------------------
    # Per-record applied artesunate concentration (nM); 0 for the
    # drug-free control well.
    # ---------------------------------------------------------------------
    conc <- STIM_ARTESUNATE_NM

    # ---------------------------------------------------------------------
    # Genotype multiplier on EC50.
    # ---------------------------------------------------------------------
    ec50_geno <- 1 + e_pfmdr1_86y_ec50     * PFMDR1_86Y +
                     e_pfmdr1_1042d_ec50   * PFMDR1_1042D +
                     e_pfmdr1_cn2_ec50     * PFMDR1_CN2 +
                     e_pfmdr1_cn3plus_ec50 * PFMDR1_CN3PLUS

    # ---------------------------------------------------------------------
    # Individual structural parameters.
    # ---------------------------------------------------------------------
    ec50  <- exp(lec50 + etalec50) * ec50_geno
    gamma <- exp(lgamma + etalgamma)

    # ---------------------------------------------------------------------
    # Sigmoid Emax inhibition (Methods Eq. 1).
    # ---------------------------------------------------------------------
    effect <- emax - (emax - e0) * conc^gamma / (conc^gamma + ec50^gamma)

    # ---------------------------------------------------------------------
    # Observation model.
    # ---------------------------------------------------------------------
    effect ~ add(addSd) + prop(propSd)
  })
}
