---
title: "Tacrine (Holford 1992)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Tacrine (Holford 1992)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5)
library(nlmixr2lib)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Holford_1992_tacrine"))$reference`
* Description: `r rxode2::rxode(readModelDb("Holford_1992_tacrine"))$description`
* Article: <https://doi.org/10.1073/pnas.89.23.11471>
* Companion methodology paper: <https://doi.org/10.1073/pnas.89.23.11466>

# Population

The packaged parameters come from a pooled analysis of two clinical trials of identical
design enrolling patients with probable Alzheimer's disease and assessing cognitive
function with the ADAS-cog total subscale (0-70 score). The pooled cohort had 909
patients with 5253 ADAS-cog observations split across the United States protocol 970-01
(n = 632) and the French protocol 970-04 (n = 277). Observation periods ran up to
five months per subject. Active treatment cycled through 0 (placebo), 40 mg/day, and
80 mg/day tacrine in three randomized titration orderings.

Demographic detail (age range, sex split, baseline ADAS-cog distribution) is described
in the companion methodology paper (Holford and Peace 1992, PNAS 89:11466-11470), which
was not on disk at extraction time; the on-disk `Results and validation` paper reports
only that the population mean ideal body weight is 60 kg (used as the size-normalisation
reference inside `model()`). The same information is available programmatically via
the model's `population` metadata:

```{r population}
rxode2::rxode(readModelDb("Holford_1992_tacrine"))$population
```

# Source trace

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lS0` (log baseline ADAS-cog) | log(28.7) | Table 2, Disease class: S_o units = 28.7 (SE 0.44, CV 37.7%) |
| `lalpha` (log progression rate) | log(6.17) | Table 2, Disease class: a units/year = 6.17 (SE 1.27, CV 208%) |
| `lbeta_a` (log absolute tacrine potency) | log(2.99) | Table 2, Pharmacodynamic class: beta_a units/80 mg/day = -2.99 (SE 0.67, CV 126%) -- magnitude stored, sign applied in model() |
| `lbeta_p` (log absolute placebo potency) | log(1.42) | Table 2, Pharmacodynamic class: beta_p units = -1.42 (SE 0.20, CV 128%) -- magnitude stored, sign applied in model() |
| `lt12eqa` | log(20.9) | Table 2, Pharmacokinetic class: t1/2,eq,a days = 20.9 (SE 6.0) -- tacrine effect-compartment equilibration |
| `lt12eqp` | log(1.58) | Table 2, Pharmacokinetic class: t1/2,eq,p days = 1.58 (SE 0.56) -- placebo on-treatment build-up |
| `lt12elp` | log(61.0) | Table 2, Pharmacokinetic class: t1/2,el,p days = 61.0 (SE 28.6) -- placebo off-treatment elimination, US reference |
| `lt12tolp` | log(13.5) | Table 2, Pharmacokinetic class: t1/2,tol,p day = 13.5 (SE 3.4) -- placebo tolerance development |
| `e_region_france_s0` | 0.08 | Table 2, Scale class: FS0_4 = 1.08 (SE 0.03), encoded as fractional shift = 0.08 |
| `e_region_france_betap` | 0.76 | Table 2, Scale class: F_beta-p_4 = 1.76 (SE 0.25), encoded as 0.76 |
| `e_region_france_t12elp` | 1.78 | Table 2, Scale class: F_t1/2,el,p_4 = 2.78 (SE 1.09), encoded as 1.78 |
| `etalS0 + etalalpha + etalbeta_a` block | c(0.1331, 0.1746, 1.6730, 0.1102, 0.6562, 0.9508) | Table 2 Population CV column + Results paragraph: CV 37.7 / 208 / 126 %; correlations rho(S0,beta_a) = 0.31, rho(S0,alpha) = 0.37, rho(beta_a,alpha) = 0.52 |
| `etalbeta_p` (diagonal) | 0.9702 | Table 2 Population CV column: CV 128% |
| `propSd` | 0.105 | Table 2, Error class: SD ADASC = 3.14 (SE 0.08); interpreted as proportional residual at typical ADAS-cog ~30, propSd = 3.14 / 30 ~= 0.105 (see Assumptions below) |
| Disease progression `ADAS_cog = S0 + alpha*t + ...` | n/a | Results paragraph + Discussion ("linear disease progression model", "an effect shifting the disease progression curve (offset model)") |
| Tacrine effect compartment `d/dt(effect1) = keqa * (Dr_norm - effect1)` | n/a | Discussion ("response may be proportional to the average steady-state concentration") + Data paragraph ("Tacrine clearance was calculated from the patient's size covariate divided by the mean value ... 60 kg for IBW") + Results ("3-week equilibration half-time") |
| Placebo effect compartment `d/dt(effect2) = TRT_PHASE * keqp * (1 - effect2) - (1 - TRT_PHASE) * kelp * effect2` | n/a | Field-standard reconstruction (companion paper not on disk); see Assumptions and deviations |
| Tolerance compartment `d/dt(effect3) = ktolp * (TRT_PHASE - effect3)` | n/a | Field-standard reconstruction (companion paper not on disk); see Assumptions and deviations |

# Virtual cohort

The original observed ADAS-cog records are not publicly available. The figures below
use a virtual cohort sized to the source paper's pooled tacalll4 dataset (909 patients)
with the US / France split, a representative titration sequence, and an IBW distribution
plausible for an elderly Alzheimer's cohort.

```{r cohort}
set.seed(1992)

make_subject_rows <- function(id, region_france, ibw) {
  # Representative titration sequence approximating tacseq114 (placebo -> 40 mg -> 80 mg)
  # over a five-month observation window.
  # Phase boundaries (days):
  #   0 - 28   : placebo  (DOSE = 0,  TRT_PHASE = 1)
  #   28 - 56  : 40 mg/day (DOSE = 40, TRT_PHASE = 1)
  #   56 - 140 : 80 mg/day (DOSE = 80, TRT_PHASE = 1)
  #   140-150  : washout (DOSE = 0,   TRT_PHASE = 0)
  #
  # Observations are placed every 14 days and at phase transitions so the
  # dose / phase indicators update correctly.
  obs_times <- sort(unique(c(seq(0, 150, by = 14), c(28, 56, 140))))
  dose_for <- function(t) ifelse(t < 28, 0, ifelse(t < 56, 40, ifelse(t < 140, 80, 0)))
  trt_for  <- function(t) ifelse(t < 140, 1L, 0L)
  data.frame(
    id            = id,
    time          = obs_times,
    evid          = 0L,
    amt           = 0,
    DOSE          = dose_for(obs_times),
    TRT_PHASE     = trt_for(obs_times),
    IBW           = ibw,
    REGION_FRANCE = as.integer(region_france)
  )
}

# Build cohort: ~70% US, ~30% France, IBW ~ Normal(60, 10) truncated to [40, 90].
n_total <- 200L
n_us    <- 140L
n_fr    <- n_total - n_us
ibw_draws <- pmin(pmax(rnorm(n_total, mean = 60, sd = 10), 40), 90)

events <- bind_rows(
  lapply(seq_len(n_us), function(i) {
    make_subject_rows(id = i, region_france = 0L, ibw = ibw_draws[i])
  }),
  lapply(seq_len(n_fr), function(i) {
    make_subject_rows(id = n_us + i, region_france = 1L, ibw = ibw_draws[n_us + i])
  })
)
events <- bind_rows(events)
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))

events |> dplyr::count(REGION_FRANCE)
range(events$IBW)
```

# Simulation

```{r simulate}
mod <- rxode2::rxode(readModelDb("Holford_1992_tacrine"))
sim <- rxode2::rxSolve(mod, events = events,
                       keep = c("DOSE", "TRT_PHASE", "REGION_FRANCE", "IBW"))
sim <- as.data.frame(sim)
head(sim)
```

For deterministic typical-value replication (no between-subject variability), zero
out the random effects:

```{r simulate-typical}
mod_typical <- mod |> rxode2::zeroRe()
sim_typical <- rxode2::rxSolve(mod_typical, events = events,
                               keep = c("DOSE", "TRT_PHASE", "REGION_FRANCE", "IBW"))
sim_typical <- as.data.frame(sim_typical)
```

# Replicate published quantities

## Disease-progression rate (untreated arm)

The source paper reports a linear disease-progression rate of 6.17 ADAS-cog units/year
in the absence of any treatment effect. Confirm this on an untreated typical-value
simulation (DOSE = 0, TRT_PHASE = 0 throughout) by extracting the slope of ADAS_cog
vs. time:

```{r progression-rate}
untreated_grid <- data.frame(
  id            = 1L,
  time          = seq(0, 365, by = 7),
  evid          = 0L,
  amt           = 0,
  DOSE          = 0,
  TRT_PHASE     = 0L,
  IBW           = 60,
  REGION_FRANCE = 0L
)
sim_untreated <- rxode2::rxSolve(mod_typical, events = untreated_grid) |>
  as.data.frame()

slope_per_year <- with(sim_untreated, coef(lm(ADAS_cog ~ time))[["time"]]) * 365.25
slope_per_year   # should be ~6.17 units/year
```

## Tacrine delay-of-progression metric

The source paper reports a "delay" of 177 days at 80 mg/day, defined as the
postponement of disease progression delivered by full-effect tacrine at the highest
dose. Algebraically `delay = |beta_a| / alpha_day = 2.99 / (6.17/365.25) ~= 177 days`.
Confirm from the packaged model:

```{r delay}
inst <- rxode2::rxode(readModelDb("Holford_1992_tacrine"))$theta
beta_a_abs <- exp(inst[["lbeta_a"]])
alpha_yr   <- exp(inst[["lalpha"]])
delay_days <- beta_a_abs / (alpha_yr / 365.25)
delay_days    # should be ~177
```

## Steady-state tacrine effect on 80 mg/day at IBW = 60

At long-time steady state on DOSE = 80 mg/day with IBW = 60 kg the tacrine effect
compartment reaches effect1 = 1 (since `Dr_norm = (DOSE / 80) * (60 / IBW) = 1`), so the
tacrine contribution to ADAS_cog is exactly `beta_a = -2.99` units. Verify with a
long-time simulation in a single typical-value subject held on 80 mg/day:

```{r ss-effect}
ss_grid <- data.frame(
  id            = 1L,
  time          = seq(0, 365, by = 7),
  evid          = 0L,
  amt           = 0,
  DOSE          = 80,
  TRT_PHASE     = 1L,
  IBW           = 60,
  REGION_FRANCE = 0L
)
sim_ss <- rxode2::rxSolve(mod_typical, events = ss_grid) |> as.data.frame()
# Drug effect = ADAS_cog - (S0 + alpha_day*t + placebo*effect2*(1-effect3))
# Simpler: extract effect1 at the longest time point.
tail(sim_ss[, c("time", "effect1", "effect2", "effect3", "ADAS_cog")], 3)
```

The reported `effect1` at long time should approach 1.0 (tacrine effect compartment at
steady state with `DOSE = 80, IBW = 60`).

## Placebo build-up and tolerance time-course

Reproduce the qualitative placebo-response shape: build-up over the first ~1.58 days
of treatment, tolerance development over ~13.5 days flattening the placebo benefit,
and slow ~61-day decline if treatment is removed.

```{r placebo-shape}
placebo_grid <- data.frame(
  id            = 1L,
  time          = seq(0, 365, by = 1),
  evid          = 0L,
  amt           = 0,
  DOSE          = 0,
  # On placebo for the first 180 days, then off treatment (washout) for the rest.
  TRT_PHASE     = ifelse(seq(0, 365, by = 1) < 180, 1L, 0L),
  IBW           = 60,
  REGION_FRANCE = 0L
)
sim_pla <- rxode2::rxSolve(mod_typical, events = placebo_grid) |> as.data.frame()

placebo_long <- sim_pla |>
  dplyr::select(time, effect2, effect3) |>
  tidyr::pivot_longer(c(effect2, effect3), names_to = "state", values_to = "value")

ggplot(placebo_long, aes(time, value, colour = state)) +
  geom_line() +
  geom_vline(xintercept = 180, linetype = "dashed", colour = "grey50") +
  labs(x = "Time (days)", y = "Compartment value (unitless)",
       colour = "State",
       title = "Placebo dynamics: build-up, tolerance, and post-washout decline",
       caption = "Dashed line at day 180 marks the end of placebo treatment.")
```

## France vs US placebo-response contrast

Reproduce the source paper's observation that the French cohort has a 76 percent
larger placebo response and a 2.78x longer placebo washout half-time than the US
cohort. Compare typical-value trajectories under the same placebo schedule:

```{r france-vs-us}
make_grid <- function(rf) {
  data.frame(
    id            = 1L,
    time          = seq(0, 365, by = 1),
    evid          = 0L,
    amt           = 0,
    DOSE          = 0,
    TRT_PHASE     = ifelse(seq(0, 365, by = 1) < 180, 1L, 0L),
    IBW           = 60,
    REGION_FRANCE = rf
  )
}

sim_us  <- rxode2::rxSolve(mod_typical, events = make_grid(0L)) |>
  as.data.frame() |> mutate(region = "US")
sim_fr  <- rxode2::rxSolve(mod_typical, events = make_grid(1L)) |>
  as.data.frame() |> mutate(region = "France")
sim_both <- bind_rows(sim_us, sim_fr)

ggplot(sim_both, aes(time, ADAS_cog, colour = region)) +
  geom_line() +
  geom_vline(xintercept = 180, linetype = "dashed", colour = "grey50") +
  labs(x = "Time (days)", y = "ADAS-cog (units)",
       colour = "Cohort",
       title = "Typical-value placebo-arm trajectories: US vs France",
       caption = "Dashed line at day 180 marks the end of placebo treatment.")
```

The peak placebo benefit and the post-washout decline rate differ between cohorts in
the direction the source paper describes (larger placebo benefit in France, slower
washout).

## Population-variability VPC during titration

Show the population-variability ADAS-cog trajectory under the titration sequence
encoded in the virtual cohort (placebo -> 40 mg -> 80 mg -> washout over 150 days),
split by US vs France:

```{r vpc}
vpc <- sim |>
  dplyr::filter(!is.na(ADAS_cog)) |>
  dplyr::mutate(cohort = ifelse(REGION_FRANCE == 1L, "France", "US")) |>
  dplyr::group_by(cohort, time) |>
  dplyr::summarise(
    Q05 = quantile(ADAS_cog, 0.05, na.rm = TRUE),
    Q50 = quantile(ADAS_cog, 0.50, na.rm = TRUE),
    Q95 = quantile(ADAS_cog, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(vpc, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25) +
  geom_line() +
  facet_wrap(~cohort) +
  labs(x = "Time (days)", y = "ADAS-cog (units)",
       title = "Population variability VPC: titration placebo -> 40 -> 80 mg/day",
       caption = "Ribbon = 5th-95th percentile; line = median.")
```

# Validation against published quantities

| Quantity | Source value | Simulated value | Match? |
|----------|--------------|-----------------|--------|
| Disease-progression rate (untreated) | 6.17 ADAS-cog units/year | `r round(slope_per_year, 2)` units/year | algebraic identity (typical-value `alpha`) |
| Delay at 80 mg/day | 177 days | `r round(delay_days, 1)` days | algebraic identity (`|beta_a| / alpha_day`) |
| Steady-state tacrine effect at 80 mg/day, IBW = 60 | -2.99 ADAS-cog units | `r round(-tail(sim_ss$effect1, 1) * beta_a_abs, 2)` units | numerical (effect1 -> 1 at long t) |
| Placebo on-treatment half-time | 1.58 days | governed by `lt12eqp` (typical-value, no IIV) | structural identity |
| Placebo off-treatment half-time (US) | 61.0 days | governed by `lt12elp` (typical-value, no IIV) | structural identity |
| Placebo tolerance half-time | 13.5 days | governed by `lt12tolp` (typical-value, no IIV) | structural identity |
| Placebo potency, France vs US | beta_p France ~ 76% larger | `r round(100 * 0.76, 1)`% larger by construction | structural identity |
| Placebo washout half-time, France vs US | 2.78x longer | `r round(2.78, 2)`x longer by construction | structural identity |

The agreement is exact for the algebraic-identity rows because the packaged model is a
straight transcription of the source parameter table; verification mostly confirms that
the time-conversion and dose-rate-scaling steps in `model()` were done correctly.

# Assumptions and deviations

* **Companion paper not on disk; placebo ODE form is a field-standard reconstruction.**
  The lead Holford 1992 PNAS 89:11471-11475 `Results and validation` paper supplies all
  parameter values (Table 2) but defers the structural model definition to the companion
  methodology paper, Holford and Peace 1992 PNAS 89:11466-11470, which was not on disk at
  extraction time and could not be retrieved through the acquisition script (PNAS direct
  download is Cloudflare-blocked; PMC PDF download is bot-protected). The placebo
  ODE form encoded here -- asymmetric on/off effect compartment with rate `keqp` during
  treatment and rate `kelp` after treatment, plus a multiplicative tolerance compartment
  driven by `ktolp` -- is the field-standard reconstruction consistent with the lead
  paper's prose ("equilibration", "elimination", "tolerance" half-times described as
  three distinct phenomena) and with how the parameters are reported (separate half-time
  estimates, separate SEs, separate units in Table 2 "Pharmacokinetic class" rows). When
  the companion paper becomes available a second pass should verify the exact ODE form
  and the precise coupling between the tolerance and placebo-effect compartments.

* **Tacrine clearance unobserved; effect-compartment input uses IBW-normalised dose
  rate.** The source paper explicitly says "Clearance could not be estimated directly,
  but potential relationships between body size and tacrine clearance were explored."
  The IBW-corrected clearance was the best size descriptor over total body weight or
  height (Table 1, AObj 2.69 worse without size correction, AObj 0.2 worse with height).
  The unobservable absolute clearance is absorbed into beta_a (defined at IBW = 60 kg and
  DOSE = 80 mg/day), leaving a relative dose-rate input `Dr_norm = (DOSE / 80) *
  (60 / IBW)` that drives the effect compartment toward 1 at steady state on
  80 mg/day, IBW = 60.

* **Residual error interpretation.** The source paper reports "Error SD ADASC 3.14"
  (SE 0.08) in Table 2 and describes the residual as proportional ("Proportional error
  better" in Table 1's reduced-model comparison). The reported SD is interpreted as the
  implied residual SD on the ADAS-cog scale at the typical observation (ADAS-cog ~30),
  giving `propSd ~ 3.14 / 30 = 0.105` in NONMEM proportional-residual convention
  Y = F * (1 + EPS). If the source actually meant an additive residual with SD 3.14 in
  ADAS-cog units, the encoding here would over-report variability at high ADAS-cog and
  under-report at low ADAS-cog. The vignette's VPC plots use the proportional form;
  re-encoding as additive (`ADAS_cog ~ add(addSd)` with `addSd = 3.14`) is a single-line
  change in `model()` if a future user wants to compare.

* **Beta_a and beta_p sign convention.** Both potency parameters are negative in the
  source paper (tacrine and placebo reduce ADAS-cog, i.e., improve cognitive function).
  Log-normal IIV is well-defined for positive quantities; the model stores the absolute
  magnitudes (`lbeta_a = log(2.99)`, `lbeta_p = log(1.42)`) and applies the negative sign
  multiplicatively in `model()`. This preserves the lognormal IIV shape (typical |beta|
  with CV 126 / 128 percent) and keeps the parameter values strictly negative across the
  population. The alternative -- additive IIV directly on the negative value, which
  would let sign flip when an extreme draw crosses zero -- is not what the source paper
  reports (it explicitly says the model failed when beta_a or beta_p was fixed to zero,
  and the inter-individual variability is described as proportional).

* **Protocol scale factors entered as fractional shifts.** The source reports the France
  scale factors multiplicatively (FS0_4 = 1.08, F_beta-p_4 = 1.76, F_t1/2el-p_4 = 2.78).
  The covariate-effect-name convention in nlmixr2lib is `e_<cov>_<param>`, with
  fractional shifts entering as `(1 + e * REGION_FRANCE)`, so the source `F_x4` values
  enter the model file as `e_region_france_x = F_x4 - 1`. The numerical model output is
  identical; only the bookkeeping differs.

* **Renal-function and age/sex covariates were tested and rejected by the source
  paper.** Holford 1992 Table 1 shows AObj 0.55 for female-sex on beta_a, AObj 0.063 for
  per-year-of-age on beta_a, and no improvement with creatinine-clearance scaling of
  tacrine clearance. None of these covariates are carried in the packaged model.

* **Observation name `ADAS_cog` rather than `Cc`.** The nlmixr2lib `checkModelConventions()`
  lint warns that single-output observations should be named `Cc` (the convention for PK
  concentration outputs). Holford 1992 is a disease-progression PD model with no
  concentration output; the cognitive endpoint is the ADAS-cog total score (0-70,
  unitless). Following the same paper-named-PD pattern used by `Conrado_2014_alzheimer`
  (`ADAS_NORM`) and `Delor_2013_alzheimer` (`SCORE_CDR_SOB`), the observation is kept as
  `ADAS_cog`. The `units$concentration` field documents this as a unitless ADAS-cog score
  rather than a mass/volume concentration; the lint warning on the units field is the
  same downstream of that choice.

* **Population demographic detail is partial.** The on-disk source paper reports only the
  pooled n (909), the US n (632), the France n (277), and the population mean IBW
  (60 kg); it does not tabulate age range / median, sex split, race / ethnicity, or
  weight range. Those details live in the companion methodology paper (Holford and
  Peace 1992 PNAS 89:11466-11470) which was not on disk. The virtual cohort here uses
  a plausible IBW distribution centered at 60 kg.

# Errata

* No errata for the lead paper are listed on the PNAS landing page or in PubMed
  (`Holford 1992 erratum`, `Holford Peace 1992 corrigendum` searched at extraction
  time). No erratum is listed for the companion methodology paper either.
