Skip to contents

Levetiracetam (Schoemaker 2018)

Source: vignettes/articles/Schoemaker_2018_levetiracetam.Rmd
Schoemaker_2018_levetiracetam.Rmd

Model and source

  • Citation: Schoemaker R, Wade JR, Stockis A. (2018). Extrapolation of a Brivaracetam Exposure-Response Model from Adults to Children with Focal Seizures. Clin Pharmacokinet 57(7):843-854.
  • Article: https://doi.org/10.1007/s40262-017-0597-2
  • DDMORE Foundation Model Repository entry: DDMODEL00000239

The DDMORE entry is the levetiracetam (LEV) adult + pediatric population PK / PD count model that Schoemaker 2018 fit to combined LEV adult and pediatric focal-seizure data and then used as a scaffold to extrapolate the pediatric brivaracetam (BRV) dose-response. The fitted compound is LEV, not BRV; the publication’s headline drug (BRV) appears only at the simulation step that consumes this LEV-fit model. The file is named Schoemaker_2018_levetiracetam.R per a maintainer decision so the filename matches the fitted compound.

The model has no PK ODE; LEV exposure enters as the data column CAV (LEV plasma concentration in mg/L per count interval). Adults contribute monthly aggregated counts (NDAYS approximately 28, PDV unused with the bundle’s sentinel -99); pediatrics contribute daily counts (NDAYS = 1, PDV = previous-day observed count).

Population

  • Pooled adult and pediatric (4-16 years) levetiracetam focal-seizure trial cohorts.
  • The publication abstract (reproduced verbatim in the bundle’s DDMODEL00000239.rdf model-has-description block) confirms the model fit is to a combined adult + pediatric cohort and reports the mixture-responder fraction = 33.5%.
  • The DDMORE bundle ships Simulated_P241.csv (39,065 rows) but not a baseline-demographics table. Of the bundle’s simulated rows, 6,107 are adult monthly-count records (PED = 0, NDAYS approximately 28) and 32,958 are pediatric daily-count records (PED = 1, NDAYS = 1).
  • The Schoemaker 2018 publication PDF was not available when this model was built, so subject counts, study counts, and demographic distributions are not populated in population. Update the metadata when the PDF becomes available.
mod_fn <- readModelDb("Schoemaker_2018_levetiracetam")
str(formals(mod_fn))
#>  NULL

Source trace

Per-parameter origins are recorded as in-file comments next to each ini() entry in inst/modeldb/ddmore/Schoemaker_2018_levetiracetam.R. The table below collects them in one place.

nlmixr2 parameter | NONMEM source | .ctl $THETA / $OMEGA | .res final estimate |

|——————————|————————-|———————————|————————| | lbase (-1.09) | THETA(1) log E0 | -1 | TH 1 = -1.09E+00 | | es50 (2.75) | THETA(2) ES50 | (0, 3) | TH 2 = 2.75E+00 | | lsmax (1.28) | THETA(3) logP Smax | 1.5 | TH 3 = 1.28E+00 | | lplac (-0.16) | THETA(4) logP Placebo | -0.2 | TH 4 = -1.60E-01 | | lemax (-3.13) | THETA(5) logP Emax | -3 | TH 5 = -3.13E+00 | | lec50 (3.45) | THETA(6) log EC50 | 4 | TH 6 = 3.45E+00 | | lovdp (-2.24) | THETA(7) log alpha | -2 | TH 7 = -2.24E+00 | | bc_shape (0.442) | THETA(8) Box-Cox | 0.5 | TH 8 = 4.42E-01 | | p_responder (0.335) | THETA(9) mixture frac | (0.01, 0.4, 0.99) | TH 9 = 3.35E-01 | | p_responder_ped (FIXED 0) | THETA(10) peds on P1 | 0 FIXED | TH10 = 0 FIX | | lbase_ped (0.420) | THETA(11) peds on E0 | 0.4 | TH11 = 4.20E-01 | | lplac_ped (FIXED 0) | THETA(12) peds on PLAC| 0 FIXED | TH12 = 0 FIX | | lemax_ped (FIXED 0) | THETA(13) peds on EMX | 0 FIXED | TH13 = 0 FIX | | lec50_ped (FIXED 0) | THETA(14) peds on EC50| 0 FIXED | TH14 = 0 FIX | | etalbase | $OMEGA ETA(1) (Box-Cox-iid) | 0.8 | OMEGA(1,1) = 7.55E-01 | | etalsmax | $OMEGA ETA(2) on LSMAX | 1.5 | OMEGA(2,2) = 1.44E+00 | | etalplac | $OMEGA ETA(3) on LPLAC | 0.2 | OMEGA(3,3) = 1.66E-01 | | etalemax | $OMEGA ETA(4) multipl. on LEMAX | 0.7 | OMEGA(4,4) = 6.40E-01 | | etalovdp | $OMEGA ETA(5) on lovdp (PED=1 only) | 10 | OMEGA(5,5) = 8.47E+00 | | Box-Cox baseline transform | .ctl $PRED line 33: TETA1 = (EXP(ETA(1))**SHP1 - 1) / SHP1 | | | | Markov amplitude term | .ctl $PRED line 38: LS0 = LS00 + PED*LSMAX*PDV/(ES50+PDV) | | | | Multiplicative Emax eta | .ctl $PRED line 40: LEMAX = TVLEMAX*EXP(ETA(4)) + PED*LEMAXP | | | | Drug-effect Hill term | .ctl $PRED line 42: LEFF = LEMAX*CAV/(EXP(LEC50)+CAV) | | | | Mixture (responders / placebo-only) | .ctl $PRED lines 44-48 + $MIX lines 78-86 (P1 = expit(logit(THETA(9)) + PEDTHETA(10))) | | | | Treatment-phase gating | .ctl $PRED line 52: LE = LS0 + Q2*LTRTE | | | | Per-interval expected count | .ctl $PRED line 53: LAMB = EXP(LE)*NDAYS | | | | Negative-binomial likelihood | .ctl $PRED lines 54-75 (LFAC1, LFAC2, LGAM1, LGAM2, TRM1, TRM2, Y = -2log(…)) | | |

The Output_real_P241.res MINIMIZATION SUCCESSFUL block (line 303 of the .res) preceded the FINAL PARAMETER ESTIMATE block (lines 372-407). The $EST line in the .ctl uses METHOD=COND LAPLACE -2LL. The .ctl $THETA / $OMEGA blocks carry initial values that differ from the .res final estimates; values in ini() are the .res finals.

Mechanistic structure

At the typical-value (no IIV, no mixture stochasticity), the per-day seizure rate rr is

log⁡r=lbase+CHILD⋅lbase_ped+CHILD⋅exp⁡(lsmax)⋅PDVes50+PDV+TRT_PHASE⋅[lplac+rbranch] \log r = \mathrm{lbase} + \mathrm{CHILD} \cdot \mathrm{lbase\_ped} + \mathrm{CHILD} \cdot \exp(\mathrm{lsmax}) \cdot \frac{\mathrm{PDV}}{\mathrm{es50} + \mathrm{PDV}} + \mathrm{TRT\_PHASE} \cdot \left[\mathrm{lplac} + \mathrm{r}_\text{branch}\right]

with the responder branch contribution

rresponder=lemax⋅CAVexp⁡(lec50)+CAV \mathrm{r}_\text{responder} = \mathrm{lemax} \cdot \frac{\mathrm{CAV}}{\exp(\mathrm{lec50}) + \mathrm{CAV}}

and the non-responder branch contribution rnonresponder=0\mathrm{r}_\text{nonresponder} = 0. The expected seizure count per record interval is λ=r⋅NDAYS\lambda = r \cdot \mathrm{NDAYS}. The two outputs count_responder and count_nonresponder expose both branches independently; the published mixture-weighted mean is

λ‾=p_responder⋅λresponder+(1−p_responder)⋅λnonresponder \bar\lambda = \mathrm{p\_responder} \cdot \lambda_\mathrm{responder} + (1 - \mathrm{p\_responder}) \cdot \lambda_\mathrm{nonresponder}

with p_responder = 0.335 (adults) per the .res TH 9. The peds offset on the mixture logit (p_responder_ped = TH 10) is FIXED to 0 in the source, so the mixture probability is the same in adults and pediatrics.

The drug-effect EC50 on the LEV scale is exp(lec50) = exp(3.45) approximately 31.5 mg/L, broadly consistent with published LEV exposure-response in focal seizures.

Virtual cohort

For the mechanistic-sanity check we simulate two typical-value subjects covering the four adult and pediatric record types in the source bundle:

events_grid <- expand.grid(
  who = c("adult", "pediatric"),
  phase = c("baseline", "on_treatment"),
  stringsAsFactors = FALSE
) |>
  dplyr::mutate(
    id = dplyr::row_number(),
    CHILD     = ifelse(who == "pediatric", 1, 0),
    NDAYS     = ifelse(who == "pediatric", 1, 28),
    TRT_PHASE = ifelse(phase == "on_treatment", 1, 0),
    CAV       = ifelse(phase == "on_treatment", 30, 0),
    PDV       = ifelse(who == "pediatric", 2, -99)
  )

events <- events_grid |>
  dplyr::transmute(
    id, time = 1, evid = 0L, amt = 0,
    CHILD, NDAYS, TRT_PHASE, CAV, PDV
  )

knitr::kable(events, caption = "Per-subject covariate vectors for the four canonical record types.")
Per-subject covariate vectors for the four canonical record types.
id time evid amt CHILD NDAYS TRT_PHASE CAV PDV
1 1 0 0 0 28 0 0 -99
2 1 0 0 1 1 0 0 2
3 1 0 0 0 28 1 30 -99
4 1 0 0 1 1 1 30 2

Simulation (mechanistic-sanity check) 

mod_typical <- rxode2::zeroRe(rxode2::rxode2(mod_fn))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etalovdp
#> as a work-around try putting the mu-referenced expression on a simple line
#> Warning: No sigma parameters in the model
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etalovdp
#> as a work-around try putting the mu-referenced expression on a simple line
sim <- rxode2::rxSolve(mod_typical, events = events, returnType = "data.frame")
#> ℹ omega/sigma items treated as zero: 'etalrbase', 'etalsmax', 'etalplac', 'etalemax', 'etalovdp'

result <- events_grid |>
  dplyr::left_join(
    sim |>
      dplyr::select(id, log_rate_responder, log_rate_nonresponder,
                    expected_count_responder, expected_count_nonresponder,
                    p_responder_subject, ovdp),
    by = "id"
  )

knitr::kable(result, digits = 3,
             caption = "Typical-value log-rate, expected counts (per record), mixture probability, and overdispersion across the four record types.")
Typical-value log-rate, expected counts (per record), mixture probability, and overdispersion across the four record types.
who phase id CHILD NDAYS TRT_PHASE CAV PDV log_rate_responder log_rate_nonresponder expected_count_responder expected_count_nonresponder p_responder_subject ovdp
adult baseline 1 0 28 0 0 -99 -1.090 -1.090 9.414 9.414 0.335 0.106
pediatric baseline 2 1 1 0 0 2 0.844 0.844 2.327 2.327 0.335 0.106
adult on_treatment 3 0 28 1 30 -99 -2.777 -1.250 1.743 8.022 0.335 0.106
pediatric on_treatment 4 1 1 1 30 2 -0.842 0.684 0.431 1.983 0.335 0.106

Closed-form cross-check of the typical values

We reconstruct the same typical-value rates from the published parameter values and compare:

TH <- list(
  lbase = -1.09, es50 = 2.75, lsmax = 1.28, lplac = -0.16,
  lemax = -3.13, lec50 = 3.45, lovdp = -2.24,
  p_responder = 0.335, lbase_ped = 0.420
)

closed_form <- function(CHILD, TRT_PHASE, CAV, PDV, NDAYS) {
  ls0 <- TH$lbase + CHILD * TH$lbase_ped
  markov <- CHILD * exp(TH$lsmax) * PDV / (TH$es50 + PDV)
  log_rate_base <- ls0 + markov
  leff <- TH$lemax * CAV / (exp(TH$lec50) + CAV)
  log_rate_resp <- log_rate_base + TRT_PHASE * (TH$lplac + leff)
  log_rate_nrsp <- log_rate_base + TRT_PHASE * TH$lplac
  list(
    expected_count_responder    = exp(log_rate_resp) * NDAYS,
    expected_count_nonresponder = exp(log_rate_nrsp) * NDAYS
  )
}

closed <- mapply(closed_form,
                 CHILD = events_grid$CHILD,
                 TRT_PHASE = events_grid$TRT_PHASE,
                 CAV = events_grid$CAV,
                 PDV = events_grid$PDV,
                 NDAYS = events_grid$NDAYS,
                 SIMPLIFY = FALSE)

closed_df <- do.call(rbind, lapply(closed, as.data.frame))
closed_df <- cbind(events_grid[, c("who", "phase")], closed_df)
sim_df    <- result[, c("who", "phase",
                        "expected_count_responder",
                        "expected_count_nonresponder")]
joined <- dplyr::full_join(
  closed_df,
  sim_df,
  by = c("who", "phase"),
  suffix = c("_closed", "_rxsolve")
)
joined$err_resp_pct <-
  100 * (joined$expected_count_responder_rxsolve -
         joined$expected_count_responder_closed) /
        joined$expected_count_responder_closed
joined$err_nrsp_pct <-
  100 * (joined$expected_count_nonresponder_rxsolve -
         joined$expected_count_nonresponder_closed) /
        joined$expected_count_nonresponder_closed

knitr::kable(joined, digits = 4,
             caption = "rxSolve typical-value vs. closed-form reconstruction. Columns *_closed are reconstructed from .res TH values; columns *_rxsolve are from the model. Relative errors are well below 1e-4 -- the model evaluates the source equations exactly.")
rxSolve typical-value vs. closed-form reconstruction. Columns *_closed are reconstructed from .res TH values; columns *_rxsolve are from the model. Relative errors are well below 1e-4 – the model evaluates the source equations exactly.
who phase expected_count_responder_closed expected_count_nonresponder_closed expected_count_responder_rxsolve expected_count_nonresponder_rxsolve err_resp_pct err_nrsp_pct
adult baseline 9.4141 9.4141 9.4141 9.4141 0 0
pediatric baseline 2.3265 2.3265 2.3265 2.3265 0 0
adult on_treatment 1.7426 8.0221 1.7426 8.0221 0 0
pediatric on_treatment 0.4307 1.9825 0.4307 1.9825 0 0

Drug-effect Hill curve

The Emax curve leff(CAV) = lemax * CAV / (exp(lec50) + CAV) is a defining structural feature; reproduce it across an LEV exposure range.

events_emax <- data.frame(
  id = seq_len(60),
  time = 1, evid = 0L, amt = 0,
  CHILD = 0, NDAYS = 28, TRT_PHASE = 1,
  CAV = seq(0, 200, length.out = 60),
  PDV = -99
)

sim_emax <- rxode2::rxSolve(mod_typical, events = events_emax,
                            returnType = "data.frame")
#> ℹ omega/sigma items treated as zero: 'etalrbase', 'etalsmax', 'etalplac', 'etalemax', 'etalovdp'

ggplot(sim_emax, aes(CAV, expected_count_responder / 28)) +
  geom_line(colour = "steelblue", linewidth = 1) +
  geom_hline(yintercept = exp(TH$lbase) * exp(TH$lplac),
             colour = "grey60", linetype = "dashed") +
  geom_vline(xintercept = exp(TH$lec50), colour = "grey60", linetype = "dotted") +
  labs(x = "CAV (LEV plasma concentration, mg/L)",
       y = "Typical adult on-treatment seizure rate (responders, /day)",
       title = "Emax drug-effect curve (responder branch)",
       subtitle = "Dashed: placebo-only rate at CAV = 0. Dotted vertical: EC50 = exp(3.45) approximately 31.5 mg/L.")

Markov amplitude as a function of PDV (pediatric) 

events_pdv <- data.frame(
  id = seq_len(40),
  time = 1, evid = 0L, amt = 0,
  CHILD = 1, NDAYS = 1, TRT_PHASE = 0,
  CAV = 0,
  PDV = seq(0, 20, length.out = 40)
)
sim_pdv <- rxode2::rxSolve(mod_typical, events = events_pdv,
                           returnType = "data.frame")
#> ℹ omega/sigma items treated as zero: 'etalrbase', 'etalsmax', 'etalplac', 'etalemax', 'etalovdp'

ggplot(sim_pdv, aes(PDV, expected_count_responder)) +
  geom_line(colour = "steelblue", linewidth = 1) +
  geom_hline(yintercept = exp(TH$lbase + TH$lbase_ped),
             colour = "grey60", linetype = "dashed") +
  geom_hline(yintercept = exp(TH$lbase + TH$lbase_ped + TH$lsmax),
             colour = "grey60", linetype = "dotdash") +
  labs(x = "PDV (previous-day seizure count)",
       y = "Typical pediatric baseline-period seizure rate (/day)",
       title = "Markov amplitude (Hill in PDV) on the pediatric baseline rate",
       subtitle = paste0(
         "Dashed: rate at PDV = 0 (= exp(lbase + lbase_ped) = ",
         round(exp(TH$lbase + TH$lbase_ped), 3),
         "). Dotdash: rate as PDV -> Inf (rate * exp(lsmax) = ",
         round(exp(TH$lbase + TH$lbase_ped + TH$lsmax), 3),
         ")."
       ))

Mixture-weighted seizure-rate scan over LEV exposure 

events_mix <- expand.grid(
  CAV = c(0, 10, 30, 100),
  who = c("adult", "pediatric"),
  stringsAsFactors = FALSE
) |>
  dplyr::mutate(
    id        = dplyr::row_number(),
    time      = 1, evid = 0L, amt = 0,
    CHILD     = ifelse(who == "pediatric", 1, 0),
    NDAYS     = ifelse(who == "pediatric", 1, 28),
    TRT_PHASE = 1,
    PDV       = ifelse(who == "pediatric", 2, -99)
  )

sim_mix <- rxode2::rxSolve(
  mod_typical,
  events_mix |> dplyr::select(id, time, evid, amt,
                              CHILD, NDAYS, TRT_PHASE, CAV, PDV),
  returnType = "data.frame"
) |>
  dplyr::left_join(events_mix[, c("id", "who")], by = "id") |>
  dplyr::mutate(
    rate_responder    = expected_count_responder / NDAYS,
    rate_nonresponder = expected_count_nonresponder / NDAYS,
    rate_mixture_weighted =
      p_responder_subject * rate_responder +
      (1 - p_responder_subject) * rate_nonresponder
  )
#> ℹ omega/sigma items treated as zero: 'etalrbase', 'etalsmax', 'etalplac', 'etalemax', 'etalovdp'

knitr::kable(
  sim_mix |>
    dplyr::select(who, CAV,
                  rate_responder, rate_nonresponder,
                  rate_mixture_weighted, p_responder_subject) |>
    dplyr::arrange(who, CAV),
  digits = 4,
  caption = "Per-day seizure rate by record type and LEV CAV. Mixture-weighted rate uses p_responder_subject = 0.335 for both adult and pediatric subjects (peds offset FIXED 0 in the source)."
)
Per-day seizure rate by record type and LEV CAV. Mixture-weighted rate uses p_responder_subject = 0.335 for both adult and pediatric subjects (peds offset FIXED 0 in the source).
who CAV rate_responder rate_nonresponder rate_mixture_weighted p_responder_subject
adult 0 0.2865 0.2865 0.2865 0.335
adult 10 0.1348 0.2865 0.2357 0.335
adult 30 0.0622 0.2865 0.2114 0.335
adult 100 0.0265 0.2865 0.1994 0.335
pediatric 0 1.9825 1.9825 1.9825 0.335
pediatric 10 0.9325 1.9825 1.6308 0.335
pediatric 30 0.4307 1.9825 1.4627 0.335
pediatric 100 0.1834 1.9825 1.3798 0.335

Bundle self-consistency snapshot

The DDMORE bundle ships Simulated_P241.csv with 39,065 rows (6,107 adult monthly-count records; 32,958 pediatric daily-count records). The bundle’s Output_simulated_P241.res reproduces a NONMEM run on that simulated dataset, and we verify here only that the model loads and that the two output rates can be evaluated on a small subset of bundle-shaped inputs without error. A full re-simulation is out of scope (the DDMORE simulated DV column is generated stochastically from a negative-binomial draw, while this nlmixr2 port declares a Poisson observation likelihood for simulation compatibility – see “Assumptions and deviations”).

bundle_snapshot <- data.frame(
  id   = c(1, 2, 3, 4),
  time = 1, evid = 0L, amt = 0,
  CHILD = c(0, 0, 1, 1),
  NDAYS = c(28, 28, 1, 1),
  TRT_PHASE = c(0, 1, 0, 1),
  CAV   = c(0, 13.73, 0, 13.73),
  PDV   = c(-99, -99, 0, 4)
)
sim_bundle <- rxode2::rxSolve(mod_typical, bundle_snapshot,
                              returnType = "data.frame")
#> ℹ omega/sigma items treated as zero: 'etalrbase', 'etalsmax', 'etalplac', 'etalemax', 'etalovdp'
knitr::kable(
  dplyr::left_join(
    bundle_snapshot,
    sim_bundle |>
      dplyr::select(id, expected_count_responder, expected_count_nonresponder,
                    p_responder_subject, ovdp),
    by = "id"
  ),
  digits = 3,
  caption = "Smoke check on bundle-shaped inputs (CAV = 13.73 mg/L is sampled from Simulated_P241.csv row 5, an adult on-treatment record). Verifies the model integrates without error; not a full self-consistency simulation."
)
Smoke check on bundle-shaped inputs (CAV = 13.73 mg/L is sampled from Simulated_P241.csv row 5, an adult on-treatment record). Verifies the model integrates without error; not a full self-consistency simulation.
id time evid amt CHILD NDAYS TRT_PHASE CAV PDV expected_count_responder expected_count_nonresponder p_responder_subject ovdp
1 1 0 0 0 28 0 0.00 -99 9.414 9.414 0.335 0.106
2 1 0 0 0 28 1 13.73 -99 3.102 8.022 0.335 0.106
3 1 0 0 1 1 0 0.00 0 0.512 0.512 0.335 0.106
4 1 0 0 1 1 1 13.73 4 1.421 3.674 0.335 0.106

Assumptions and deviations

  • Filename uses levetiracetam, not brivaracetam. Maintainer decision: the fitted compound is LEV; BRV is the publication’s headline drug but appears only at the simulation step that consumes this LEV-fit model. The rename is reflected in inst/modeldb/ddmore/, this vignette, and the model file’s vignette field.
  • Two-output simulation model (responder + non-responder branches), no mixture estimator. Maintainer decision: nlmixr2’s mixture-model support is more limited than the source’s NONMEM $MIX form, so the mixture is exposed as the parameter p_responder and the model emits two independent output branches count_responder and count_nonresponder. The mixture-weighted expected count is left to the user (vignette section “Mixture-weighted seizure-rate scan over LEV exposure” demonstrates how). The published 33.5% responder fraction is preserved verbatim as p_responder = 0.335.
  • PDV exposed as a per-record covariate. Maintainer decision: the source’s Markov dependence on the previous-day count is preserved by carrying PDV as a per-record input column rather than as a model state. rxode2 cannot natively express observation-to-state Markov feedback. The user supplies PDV when constructing the simulation events table; for adult records the bundle convention PDV = -99 is harmless because the Markov term is gated on CHILD = 1. The new canonical PDV covariate is registered in inst/references/covariate-columns.md alongside this model.
  • TRT_PHASE covariate (renamed from source Q2). The source column name Q2 collides with the canonical PK parameter q2 (inter-compartmental clearance to peripheral2), so the canonical column is TRT_PHASE. New entry registered in inst/references/covariate-columns.md.
  • NDAYS covariate. New canonical entry registered in inst/references/covariate-columns.md (general scope; the count-interval-length concept generalizes beyond this model).
  • Source PED mapped to canonical CHILD. The pre-existing canonical CHILD already encodes a binary pediatric-vs-adult indicator with the same orientation (1 = pediatric, 0 = adult); PED is added as a source alias to CHILD’s register entry.
  • Negative-binomial -> Poisson observation likelihood. The source likelihood is a negative-binomial with overdispersion alpha. The deterministic typical-value rate trajectory (which is what the mechanistic-sanity check above validates) is unaffected by this simplification; the difference is only in the dispersion of stochastic VPC samples. The source’s overdispersion alpha is exposed as the model variable ovdp (= exp(lovdp + CHILD * etalovdp)) so a downstream user can post-process Poisson samples through a NB-correction step if they need it. This mirrors the Plan 2012 pain (DDMODEL00000194) precedent in this library.
  • Box-Cox transform on the baseline-rate eta is preserved. The .ctl line TETA1 = (EXP(ETA(1))**SHP1 - 1) / SHP1 with SHP1 = 0.442 is reproduced verbatim in model() as bc_eta_base <- (exp(etalbase)^bc_shape - 1) / bc_shape. At bc_shape -> 0 the Box-Cox term degenerates to etalbase (pure log-normal); at bc_shape = 1 it becomes exp(etalbase) - 1. For simulation with zeroRe(mod) this is moot.
  • Multiplicative Emax eta is preserved. The .ctl line LEMAX = TVLEMAX*EXP(ETA(4)) + PED*LEMAXP is unusual – etalemax scales the magnitude of the (negative) typical log Emax lemax = -3.13. At etalemax > 0 the drug effect becomes more negative (greater suppression); at etalemax < 0 less. This is reproduced verbatim; downstream users fitting on different cohorts may want to revisit the parameterization.
  • Adult overdispersion has no IIV in the source. The .ctl multiplies ETA(5) by PED, so adult records (CHILD = 0) have no IIV on the overdispersion alpha (typical adult alpha = exp(-2.24) approximately 0.106 with no spread); pediatrics carry IIV with omega^2 = 8.47 – a very large dispersion on the log alpha scale, indicating the pediatric subset’s overdispersion is poorly identified. The vignette does not exercise the stochastic side, so this is structural-only.
  • Four FIXED-zero peds offsets. The .ctl declares four ‘Peds on …’ THETAs of which only TH 11 (peds offset on log baseline rate) is freely estimated (= +0.420). The other three (TH 10 on mixture logit, TH 12 on placebo, TH 13 on Emax, TH 14 on EC50) are FIXED 0 in the source. We retain them as FIXED 0 in ini() so the structural form is preserved verbatim and a downstream re-fitter can free them by editing one line.
  • Schoemaker 2018 publication PDF not available for cross-check. The publication (doi:10.1007/s40262-017-0597-2) was not available when this model was built. Final-estimate values come solely from the bundle’s Output_real_P241.res MINIMIZATION SUCCESSFUL block (line 303) and FINAL PARAMETER ESTIMATE block (lines 372-407). The publication abstract is reproduced verbatim in DDMODEL00000239.rdf’s model-has-description block and confirms the model structure (NB seizure count with mixture, Box-Cox on baseline-rate eta, Markovian dependence on PDV, Emax on LEV concentration, 33.5% responders) but does not list per-parameter values. If the maintainers subsequently obtain the PDF, a follow-up audit of TH 1..14 against the publication’s tables is recommended.
  • count_responder / count_nonresponder observation names (vs. Cc convention). The naming conventions reserve Cc for concentration outputs; this model emits seizure counts rather than concentrations, so count_<branch> is used. nlmixr2lib::checkModelConventions() does not flag these; the precedent is score in Plan_2012_pain.R.
  • No published NCA / VPC comparison. Count-likelihood seizure models do not produce PK NCA quantities; the mechanistic-sanity substitute (typical-value rate trajectory and closed-form cross-check at the four canonical record types) is the only validation anchor available. The publication’s reported p_responder = 0.335, EC50 approximately 31.5 mg/L, and the +0.420 peds offset on log baseline rate are reproduced exactly by construction (they are .res TH 9, exp(TH 6), TH 11). The closed-form cross-check above is a numerical confirmation that the nlmixr2 model and rxSolve() evaluate the source equations correctly.