---
title: "Digoxin (Zhou 2010)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Digoxin (Zhou 2010)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode2(readModelDb("Zhou_2010_digoxin"))$reference`
* Description: `r rxode2::rxode2(readModelDb("Zhou_2010_digoxin"))$description`
* Article: [Acta Pharmacologica Sinica. 2010;31(6):753-758](https://doi.org/10.1038/aps.2010.51)
* Free full-text: [www.nature.com/aps](https://www.nature.com/articles/aps201051)

# Population

Zhou et al. (2010) developed a population PK model of digoxin in 119 older
Chinese inpatients (>= 60 years, 95% with congestive heart failure)
recruited retrospectively at the General Hospital of the Air Force, PLA
(Beijing). Subjects had been on oral digoxin (0.25 mg tablets, Sine
Pharmaceutical, Shanghai) for at least 7 days prior to a single PK draw,
and 173 sparse serum concentrations were collected at a median 22.9 h
(range 6-192 h) after the last dose. Patients with severe hepatic or
renal dysfunction and concentrations outside the assay range
(0.2 ug/L LLOQ, fluorescence polarization immunoassay; Abbott TDx-FLx)
were excluded.

Baseline demographics (Zhou 2010 Table 1):

| Variable                  | Mean  | Range            |
|---------------------------|-------|------------------|
| N                         | 119   |                  |
| Male : Female             | 69:50 |                  |
| Age (years)               | 71.0  | 60-88            |
| Weight (kg)               | 62.9  | 34-91            |
| ALT (U/L)                 | 29.3  | 3-383            |
| AST (U/L)                 | 33.9  | 8-402            |
| BUN (mmol/L)              | 9.0   | 2.5-28.9         |
| Serum creatinine (umol/L) | 126.8 | 36-686           |
| Albumin (g/L)             | 63.3  | 31-89            |
| Digoxin serum (ug/L)      | 1.11  | 0.07-4.45        |

Concomitant medications recorded: spironolactone (SPI, 32 patients);
calcium-channel blockers (nifedipine 53, diltiazem 1); nitrates (86);
propafenone (27). Of these, only SPI was retained in the final
covariate model (on Cl/F).

The same metadata is available programmatically via
`readModelDb("Zhou_2010_digoxin")$population`.

# Source trace

Per-parameter origin is recorded as an in-file comment next to each `ini()`
entry in `inst/modeldb/specificDrugs/Zhou_2010_digoxin.R`. The table below
collects them for review.

| Equation / parameter        | Value                  | Source location                                              |
|-----------------------------|------------------------|--------------------------------------------------------------|
| Cl/F structural             | 5.90 L/h               | Table 7: Cl/F = 5.90 (RSE 6.97%, 95% CI 5.09-6.71)           |
| Vd/F structural             | 550 L                  | Table 7: V/F = 550 (RSE 19.6%, 95% CI 338-762)               |
| Ka                          | 1.63 1/h (FIXED)       | Table 7: Ka = 1.63 (fixed); Methods (Results): "Ka and IIV were fixed to 1.63 h-1 and 0" |
| theta_SPI-Cl                | 0.412                  | Table 7: theta_SPI-Cl = 0.412 (RSE 26.5%, 95% CI 0.198-0.626) |
| theta_WT-Cl                 | 0.0101 1/kg            | Table 7: theta_WT-Cl = 0.0101 (RSE 120%, 95% CI -0.0136-0.0338) |
| theta_Cr-Cl                 | 0.00120 L/umol         | Table 7: theta_Cr-Cl = 0.00120 (RSE 45.8%, 95% CI 0.000122-0.00228) |
| IIV Cl/F                    | 49.0% CV               | Table 7: Inter-Indv 49.0% on Cl/F                            |
| IIV Vd/F                    | 94.3% CV               | Table 7: Inter-Indv 94.3% on Vd/F                            |
| IIV Ka                      | 0 (FIXED)              | Results section: "Ka and inter-individual variation were fixed to 1.63 h-1 and 0, respectively" |
| Residual additive SD        | 0.365 ug/L             | Table 7 / Abstract: "residual variability (SD) ... was 0.365 ug/L" |
| Proportional residual       | dropped in final model | Table 7: "Proportional, CV (%)" reported as "-"              |
| WT reference                | 62.9 kg                | Table 1 population mean                                      |
| CREAT (Cr) reference        | 126.8 umol/L           | Table 1 population mean                                      |
| Structure (1-cmt 1st-order) | n/a                    | Methods (Structural pharmacokinetic model): "a one-compartment open model" with NONMEM ADVAN2 |
| Cl/F covariate equation     | n/a                    | Discussion: "Cl/F = 5.9 x [1 - 0.412*SPI] x [1 - 0.0101*(WT - 62.9)] x [1 - 0.00120*(Cr - 126.8)] (L/h)" |

# Virtual cohort

Zhou 2010 Table 8 reports 8 external-validation patients receiving
0.125 mg digoxin orally; the table lists each patient's WT, dose, Cr,
SPI flag, individual-predicted (IPRE) and observed (Obs) concentrations.
The cohort below reproduces those 8 patients exactly so the typical-value
trajectories can be compared against the IPRE column from the source.

```{r cohort}
set.seed(20100501L)

# Zhou 2010 Table 8: 8 elderly heart-failure patients on long-term oral digoxin.
# The paper reports steady-state-leaning individual-predicted concentrations from
# the final model. Dose 0.125 mg/day for everyone; concentration column is
# steady-state-leaning at a single phlebotomy time. The actual phlebotomy time
# is not stated per patient; we simulate to a 24-h dosing window at steady state
# (day 14) and report the 22-h post-dose point as a comparable trough.
table8 <- tibble::tribble(
  ~id, ~WT,   ~dose_mg, ~CREAT, ~CONMED_SPIRON, ~IPRE_paper, ~Obs_paper, ~Percent_diff,
  1L,  58,    0.125,    250,    0,              1.03,        1.37,       24.8,
  2L,  70,    0.125,    125 * 1, 1,             0.99,        0.93,       -6.5,    # CREAT 87 per Table 8
  3L,  63,    0.125,    80,     0,              1.11,        1.36,       18.4,
  4L,  58,    0.125,    171,    1,              1.93,        1.99,       3.02,
  5L,  67.5,  0.125,    74,     0,              0.59,        0.51,       -15.7,
  6L,  86.5,  0.125,    65,     1,              2.28,        2.56,       10.9,
  7L,  55,    0.125,    57,     1,              1.20,        1.15,       -4.3,
  8L,  65,    0.125,    115,    0,              1.28,        1.58,       18.99
)
# Patch row 2 CREAT (the inline expression above was a placeholder so the value 87
# is documented in source); set the canonical value:
table8$CREAT[2] <- 87

# Build the event table: 14 days of 0.125 mg QD oral dosing for each subject,
# observation at every hour during day 14 (the steady-state window).
build_subject_events <- function(row) {
  dose_times <- seq(0, 13 * 24, by = 24)
  obs_times <- seq(13 * 24, 14 * 24, by = 1)
  ev <- rxode2::et(amt = row$dose_mg, time = dose_times[1], cmt = "depot")
  for (tt in dose_times[-1]) {
    ev <- rxode2::et(ev, amt = row$dose_mg, time = tt, cmt = "depot")
  }
  ev <- rxode2::et(ev, obs_times)
  out <- as.data.frame(ev)
  out$id <- row$id
  out$WT <- row$WT
  out$CREAT <- row$CREAT
  out$CONMED_SPIRON <- row$CONMED_SPIRON
  out
}

events <- bind_rows(lapply(seq_len(nrow(table8)),
                           function(i) build_subject_events(table8[i, ])))

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

The paper's individual-predicted concentrations (IPRE) are typical-value
predictions for the per-subject covariates with the inter-individual etas
zeroed out (NONMEM IPRED). The chunks below reproduce that behaviour with
`rxode2::zeroRe()` so the simulated `Cc` is comparable to Table 8's IPRE
column.

```{r simulate-typical}
mod <- readModelDb("Zhou_2010_digoxin")
mod_typical <- rxode2::zeroRe(mod)

sim_typical <- rxode2::rxSolve(
  mod_typical,
  events = events,
  keep   = c("WT", "CREAT", "CONMED_SPIRON")
) |>
  as.data.frame()
```

# Replicate Table 8 (external-validation cohort)

Table 8 of Zhou 2010 reports individual-predicted vs observed digoxin
serum concentrations for 8 elderly heart-failure patients. The exact
phlebotomy time per patient is not stated (the source describes
trough-leaning samples at median 22.9 h post-dose across the parent
cohort), so the simulated comparison takes the day-14 22-h post-dose
concentration as a representative trough.

```{r table8-comparison}
trough_window <- function(df) {
  # day-14 22 h post-dose corresponds to t = 13*24 + 22 = 334 h in the event grid
  df |>
    filter(time == 13 * 24 + 22) |>
    select(id, time, Cc, WT, CREAT, CONMED_SPIRON)
}

sim_trough <- trough_window(sim_typical)

t8_compare <- table8 |>
  select(id, WT, CREAT, CONMED_SPIRON, IPRE_paper, Obs_paper) |>
  left_join(sim_trough |> select(id, Cc_sim = Cc), by = "id") |>
  mutate(percent_diff_sim_vs_paper_IPRE = 100 * (Cc_sim - IPRE_paper) / IPRE_paper)

knitr::kable(
  t8_compare,
  digits  = 2,
  caption = "Zhou 2010 Table 8 patients: simulated typical-value day-14 22-h post-dose concentration vs the paper-reported IPRE."
)
```

```{r table8-plot}
t8_compare |>
  ggplot(aes(IPRE_paper, Cc_sim, label = id)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  geom_point(size = 3) +
  geom_text(nudge_x = 0.05, nudge_y = 0.05) +
  expand_limits(x = 0, y = 0) +
  labs(x = "Paper IPRE (ug/L)",
       y = "Simulated Cc at day 14 + 22 h (ug/L)",
       title = "Zhou 2010 Table 8: simulated vs reported individual predictions",
       caption = "Points labelled by Table 8 patient ID. Dashed line is unity.")
```

# Replicate Cl/F covariate surfaces

The Discussion of Zhou 2010 reproduces the final-model formula
`Cl/F = 5.9 x [1 - 0.412*SPI] x [1 - 0.0101*(WT - 62.9)] x [1 - 0.0012*(Cr - 126.8)]`.
The chunk below evaluates the formula analytically across the observed WT
and CREAT ranges so the multiplicative effects can be sanity-checked
against the in-file `ini()` coefficients.

```{r covariate-surfaces}
typical_cl <- function(WT, CREAT, CONMED_SPIRON,
                       cl0 = 5.9, theta_spi = 0.412,
                       theta_wt = 0.0101, theta_cr = 0.00120,
                       WT_ref = 62.9, CREAT_ref = 126.8) {
  cl0 *
    (1 - theta_spi * CONMED_SPIRON) *
    (1 - theta_wt  * (WT    - WT_ref)) *
    (1 - theta_cr  * (CREAT - CREAT_ref))
}

cov_grid <- expand.grid(
  WT = seq(40, 90, by = 5),
  CREAT = seq(40, 400, by = 20),
  CONMED_SPIRON = c(0, 1)
) |>
  as_tibble() |>
  mutate(cl_typical = typical_cl(WT, CREAT, CONMED_SPIRON),
         SPI_label = ifelse(CONMED_SPIRON == 1, "On spironolactone", "Off spironolactone"))

ggplot(cov_grid, aes(CREAT, cl_typical, colour = factor(WT), group = WT)) +
  geom_line() +
  facet_wrap(~ SPI_label) +
  scale_colour_viridis_d(name = "Weight (kg)") +
  labs(x = "Serum creatinine (umol/L)",
       y = "Typical Cl/F (L/h)",
       title = "Zhou 2010 final-model Cl/F surface",
       caption = "Reproduces the Discussion's final-model Cl/F formula across the observed WT (34-91 kg) and CREAT (36-686 umol/L) ranges.")
```

# PKNCA single-dose validation

Zhou 2010 does not report any formal NCA-derived parameters
(no Cmax / Tmax / AUC / half-life summaries appear in the article).
The chunk below runs PKNCA on a simulated single 0.125 mg oral dose
for the reference subject (WT = 62.9 kg, CREAT = 126.8 umol/L,
no spironolactone) as a sanity-check that the model produces
non-negative, dimensionally consistent serum concentrations with a
half-life close to the canonical digoxin half-life of 36-48 h
(`t1/2 = ln(2) * V / CL = ln(2) * 550 / 5.9 ~= 64.6 h` at the typical
subject; consistent with the upper end of the published range for
elderly patients).

```{r pknca}
single_dose_evt <- rxode2::et(amt = 0.125, cmt = "depot") |>
  rxode2::et(seq(0, 192, by = 0.25))
single_dose_evt <- as.data.frame(single_dose_evt)
single_dose_evt$id <- 1L
single_dose_evt$WT <- 62.9
single_dose_evt$CREAT <- 126.8
single_dose_evt$CONMED_SPIRON <- 0L
single_dose_evt$treatment <- "single 0.125 mg oral, typical subject"

sim_single <- rxode2::rxSolve(
  rxode2::zeroRe(mod),
  events = single_dose_evt,
  keep   = c("WT", "CREAT", "CONMED_SPIRON", "treatment")
) |>
  as.data.frame()

if (!"id" %in% names(sim_single)) sim_single$id <- 1L

sim_nca <- sim_single |>
  filter(!is.na(Cc)) |>
  select(id, time, Cc, treatment)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id)

dose_df <- single_dose_evt |>
  filter(evid == 1L) |>
  select(id, time, amt, treatment)

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id)

intervals <- data.frame(
  start      = 0,
  end        = Inf,
  cmax       = TRUE,
  tmax       = TRUE,
  aucinf.obs = TRUE,
  half.life  = TRUE
)

nca_data <- PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_res  <- PKNCA::pk.nca(nca_data)

knitr::kable(
  as.data.frame(nca_res$result),
  digits  = 4,
  caption = "Simulated single 0.125 mg oral PKNCA output for the Zhou 2010 reference subject (WT = 62.9 kg, CREAT = 126.8 umol/L, no spironolactone)."
)
```

# Assumptions and deviations

- **Residual error is additive-only in the final model.** Zhou 2010 Table 2
  (basic model) reports both proportional (CV) and additive (SD) residual
  terms, consistent with the basic-model equation
  `Y = YPRED + YPRED * ERR(1) + ERR(2)`. Tables 5 and 7 (full-regression and
  final models) both record the proportional CV as "-" with only an additive
  SD value (0.365 ug/L in the final model). This is read as the proportional
  term having been dropped during covariate-model refinement; the model file
  therefore implements `Cc ~ add(addSd)` with `addSd = 0.365`.
- **Ka and Ka IIV are fixed.** The paper's Results section (Basic model)
  states: "Because the absorption phase in the data had been completed, Ka
  and inter-individual variation were fixed to 1.63 h-1 and 0, respectively,
  in the next calculation to avoid the effect of Ka fluctuation on the
  stability of model." `lka` is therefore wrapped in `fixed()` and `etalka`
  is omitted from the IIV block.
- **Covariate-effect sign convention.** The final-model formula reported in
  the Discussion uses positive coefficients with `(1 - theta * cov)`
  multiplicative linear-deviation factors. The model file preserves the
  paper-reported positive signs of `theta_SPI-Cl = 0.412`,
  `theta_WT-Cl = 0.0101`, and `theta_Cr-Cl = 0.00120`; the minus signs
  appear explicitly in `model()` so the in-file `ini()` value matches the
  source table verbatim.
- **WT effect coefficient is statistically marginal but retained.** Table 7
  reports `theta_WT-Cl = 0.0101` with RSE 120% and 95% CI -0.0136 to 0.0338
  -- the confidence interval crosses zero. Zhou 2010 retained the WT-Cl
  term in the final model on the strength of the backward-elimination dOFV
  criterion (7.515 > 6.63 at the 0.01 level; Table 6); the model file
  therefore retains it as the published final model rather than dropping
  it on the 95% CI argument.
- **Reference covariate values are population means, not medians.** Table 1
  reports means (62.9 kg WT, 126.8 umol/L CREAT) and not medians; the
  Discussion formula and the model file use those means as the reference
  values, matching the source.
- **No NCA values published.** Zhou 2010 reports population Cl/F and Vd/F
  estimates but not Cmax / Tmax / AUC / half-life. The PKNCA section above
  is a sanity check on the implementation; it is not a comparison to a
  published NCA table.
- **Table 8 IPRE vs simulated typical-value.** Table 8 reports
  individual-predicted (IPRE) concentrations from the NONMEM fit. The
  vignette compares the simulated day-14 22-h post-dose concentration
  (zero-eta, per-subject covariates) against IPRE. The exact phlebotomy time
  is not reported per patient in Table 8 (only the cohort-level median
  22.9 h post-dose from Table 1), so the simulated time is an approximation;
  some scatter around the unity line is expected and is documented in the
  Table 8 comparison plot above.
- **Population field `assay`.** The paper reports fluorescence polarization
  immunoassay on Abbott TDx-FLx with a 0.2 ug/L LLOQ; this is recorded in
  `population$assay` rather than the canonical schema so the assay context
  is searchable via `readModelDb("Zhou_2010_digoxin")$population`.

# Errata

No errata or corrigenda were located in a search of the Acta Pharmacologica
Sinica journal site or PubMed against Zhou XD 2010 digoxin as of the model
extraction date.
