---
title: "[177Lu]Lu-PSMA-I&T PSA dynamics (Siebinga 2024)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{[177Lu]Lu-PSMA-I&T PSA dynamics (Siebinga 2024)}
  %\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(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Siebinga_2024_lu177psmaIT"))$reference`
* Article: <https://doi.org/10.1186/s40658-024-00642-2>
* PubMed Central: <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC11043318/>

This is the sequel to `Siebinga_2023_lu177psma617` and the first published
PK/PD model that ties radioligand tumor exposure to a biochemical response
measure. The PK layer is a five-compartment first-order system fitted to 409
quantitative SPECT/CT scans; the PD layer is a single prostate-specific antigen
(PSA) state fitted sequentially to 566 PSA observations.

The structure (Figure 1 of the paper):

| Paper index | Compartment | nlmixr2lib state | Observed as |
|---|---|---|---|
| 1 | Central | `central` | `Cc` (MBq/L) |
| 2 | Salivary glands (lumped) | `salivary_gland` | `Asalivary_gland` (MBq) |
| 3 | Kidneys (lumped) | `kidney` | `Akidney` (MBq) |
| 4 | Tumor lesions (lumped) | `tumor` | `Atumor` (MBq) |
| 5 | Remaining tissues (lumped) | `other` | not observed |
| -- | Effect compartment | `effect` | not observed |
| -- | PSA | `PSA` | `PSA` (ug/L) |

Only the central compartment has a volume (`V1` = 10.3 L, fixed), so it is
observed as a concentration and every tissue as a decay-corrected radioactivity
amount. Salivary-gland uptake is the only saturable process. The PD layer is
Equation 12,

```
d(PSA)/dt = kG * PSA - (Edrug + Edelayed) * PSA
```

with a direct effect linear in the tumor activity concentration
(`Edrug = kD,direct * Ctumor`, Equation 9) and a delayed effect linear in an
effect-compartment concentration equilibrating at `ke0`.

```{r model}
mod <- readModelDb("Siebinga_2024_lu177psmaIT")
ui <- rxode2::rxode(mod)
ui
```

# Population

Seventy-six men with metastatic castration-resistant prostate cancer (mCRPC)
treated at the Netherlands Cancer Institute between September 2019 and January
2023 (IRBd21288). All had PSMA-positive lesions on diagnostic
[68Ga]Ga-PSMA PET/CT and no PSMA-negative lesions on contrast-enhanced CT, per
the VISION criteria. Median age 73 years (48-91), weight 79 kg (61-116),
Cockcroft-Gault GFR 80.9 mL/min (25.9-181), baseline PSA 260 ug/L (0.12-4896),
and segmented target-tumor volume 0.0443 L (0.000122-0.546) (Table 1).

Two schedules were used clinically: four cycles six weeks apart ("4 x 6",
n = 10, patients 1-10) and two cycles two weeks apart repeated after twelve
weeks ("2 x 2 - repeated after twelve weeks", n = 66). Median injected activity
was 7378 MBq (5605-7763) per cycle and patients received up to eight cycles.

```{r population}
str(ui$meta$population)
```

# Source trace

Every value in `ini()` and every equation in `model()`, with its location in
Siebinga 2024.

| Quantity | Model name | Value | Source |
|---|---|---|---|
| Renal excretion rate | `lkel` | 0.253 1/h | Table 2, k10 |
| Central volume (fixed) | `lvc` | 10.3 L | Table 2, V1, footnote c |
| Central -> salivary gland | `lkin_salivary_gland` | 0.0105 1/h | Table 2, k12 |
| Salivary gland -> central | `lkout_salivary_gland` | 0.0629 1/h | Table 2, k21 |
| Central -> kidney | `lkin_kidney` | 0.0321 1/h | Table 2, k13 |
| Kidney -> central | `lkout_kidney` | 0.0625 1/h | Table 2, k31 |
| Central -> tumor | `lkin_tumor` | 0.00967 1/h | Table 2, k14 |
| Tumor -> central | `lkout_tumor` | 0.0150 1/h | Table 2, k41 |
| Central -> remaining tissue | `lkin_other` | 0.275 1/h | Table 2, k15 |
| Remaining tissue -> central | `lkout_other` | 0.0247 1/h | Table 2, k51 |
| Salivary-gland Bmax | `lbmax` | 134 MBq | Table 2 |
| Tumor volume on tumor uptake | `e_tum_vol_kin_tumor` | 1.08 | Table 2, footnote a (Eq. 5) |
| Tumor volume on salivary uptake | `e_tum_vol_kin_salivary_gland` | 0.0910 | Table 2, footnote a (Eq. 5) |
| Cycle-2 uptake fraction | `e_cycle2_kin_tumor` | 0.731 | Table 2, footnote b (Eq. 8) |
| Cycle-3 uptake fraction | `e_cycle3_kin_tumor` | 0.498 | Table 2, footnote b (Eq. 8) |
| Cycle-4-7 uptake fraction | `e_cycle47_kin_tumor` | 0.436 | Table 2, footnote b (Eq. 8) |
| Allometric exponent, volume | `e_wt_vc` | 1 (fixed) | Methods, "Pharmacokinetic model development" |
| Allometric exponent, rates | `e_wt_k` | -0.25 (fixed) | Methods, as above |
| Baseline PSA (fixed) | `rbase` | 140 ug/L | Table 3, footnote a |
| Tumor volume on baseline PSA | `e_tum_vol_rbase` | 57.5 ug/L | Table 3, footnote b (Eq. 6) |
| PSA growth rate | `lkgrow` | 0.000408 1/h | Table 3, kG |
| Direct drug effect | `lkd_direct` | 0.00335 L/day/GBq | Table 3, kD,direct (see Errata) |
| Effect-compartment rate | `lke0` | 0.00128 1/h | Table 3, ke0 |
| Delayed drug effect | `lkd_delay` | 0.0000328 L/day/MBq | Table 3, kD,delay |
| Box-Cox shape on the kG eta | `boxcox_lkgrow` | -0.822 | Table 3 |
| IIV k10 / k13 / k14 / k15 / Bmax | `etalkel` etc. | 33.1 / 33.4 / 63.4 / 29.1 / 67.0 % CV | Table 2 |
| IOV on tumor uptake | `etaiov_kin_tumor_*` | 37.8 % CV | Table 2 |
| IIV baseline PSA / kG / kD,direct / kD,delay | `etarbase` etc. | 179 / 90.9 / 140 / 87.5 % CV | Table 3 |
| Residual error, central | `propSd`, `addSd` | 55.5 % CV, 9.57 MBq/L | Table 2 (Eq. 3) |
| Residual error, tissues | `propSd_A*` | 39.7 / 31.9 / 32.7 % CV | Table 2 (Eq. 4) |
| Residual error, PSA | `propSd_PSA` | 29.3 % CV | Table 3 (Eq. 4) |
| Saturable uptake | `fluxSalivaryGland` | -- | Equation 1 |
| Exponential IIV / IOV | -- | -- | Equation 2 |
| Power covariate model | -- | -- | Equation 5 |
| Fractional cycle effect | `cycleFactor` | -- | Equation 8 |
| Linear direct drug effect | `edrug` | -- | Equation 9 |
| PSA turnover | `d/dt(PSA)` | -- | Equation 12 |

IIV percentages are converted to log-normal variances as
`omega^2 = log(CV^2 + 1)`, which is the variance implied by the exponential
model of Equation 2.

# Virtual cohort

Original patient data are not publicly available. Two cohorts are built.

The **PK cohort** is a single typical patient (median weight 79 kg, median
segmented tumor volume 0.0443 L = 44300 mm^3) receiving one 7378 MBq cycle,
with a dense observation grid so that trapezoidal integration is accurate. It
is used for the deterministic structural checks.

The **PD cohort** is 200 subjects per dosing regimen (the paper simulated 2000;
200 per arm is the nlmixr2lib cap and is ample, because the reported endpoints
are cohort proportions). Weight is drawn normally and truncated to the Table 1
range; segmented tumor volume is drawn log-normally with the Table 1 median and
truncated to the Table 1 range, because the paper reports only median and range.
The **same** covariate draws and the **same** random seed are used for both
regimens, so the two arms share common random numbers and the between-regimen
difference is not swamped by Monte-Carlo noise.

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

n_per_arm <- 200L
week <- 7 * 24

# --- PK cohort: one typical patient, one cycle, dense grid -----------------
tgrid_pk <- sort(unique(c(seq(0, 4, by = 0.05), seq(4, 48, by = 0.5),
                          seq(48, 1008, by = 4))))

pk_endpoints <- tibble::tribble(
  ~cmt,             ~dvid,
  "central",        1L,
  "salivary_gland", 2L,
  "kidney",         3L,
  "tumor",          4L
)

events_pk <- pk_endpoints |>
  tidyr::crossing(time = tgrid_pk) |>
  mutate(id = 1L, amt = NA_real_, evid = 0L) |>
  bind_rows(tibble(id = 1L, time = 0, amt = 7378, evid = 1L,
                   cmt = "central", dvid = NA_integer_)) |>
  mutate(WT = 79, TUM_VOL = 44300, CYCLE = 1L) |>
  arrange(id, time, desc(evid))

# --- PD cohort: 200 subjects per regimen, shared covariates ----------------
subj <- tibble(
  id = seq_len(n_per_arm),
  WT = pmin(pmax(rnorm(n_per_arm, 79, 12), 61), 116),
  TUM_VOL = pmin(pmax(exp(rnorm(n_per_arm, log(44300), 1.6)), 122), 546000)
)

regimens <- list(
  "4 x 6"  = seq(0, by = 6 * week, length.out = 4),
  "2 x 2 - repeated" = c(0, 2 * week, 12 * week, 12 * week + 2 * week)
)

# Observe PSA weekly out to 32 weeks, and force a record exactly at baseline
# and at the 24-week response endpoint.
tgrid_pd <- sort(unique(c(seq(0, 32 * week, by = week), 24 * week)))

make_pd_events <- function(dose_times) {
  doses <- subj |>
    select(id, WT, TUM_VOL) |>
    tidyr::crossing(time = dose_times) |>
    mutate(amt = 7378, evid = 1L, cmt = "central", dvid = NA_integer_)
  obs <- subj |>
    select(id, WT, TUM_VOL) |>
    tidyr::crossing(time = tgrid_pd) |>
    mutate(amt = NA_real_, evid = 0L, cmt = "PSA", dvid = 5L)
  bind_rows(doses, obs) |>
    # CYCLE is time-varying: it counts how many administrations have occurred.
    mutate(CYCLE = pmax(1L, findInterval(time, dose_times - 1e-6))) |>
    arrange(id, time, desc(evid))
}

events_pd <- lapply(regimens, make_pd_events)

stopifnot(
  nrow(events_pk) > 0,
  vapply(events_pd, function(d) length(unique(d$id)) == n_per_arm, logical(1))
)
```

# Simulation

`useLinCmt = FALSE` is required: rxode2's automatic ODE-to-`linCmt` conversion
corrupts the `dvid` -> `cmt` mapping for multi-endpoint models like this one.

```{r simulate}
# Typical-value PK profile (all etas zeroed) for the structural checks.
mod_typical <- rxode2::zeroRe(mod)
sim_pk <- rxode2::rxSolve(
  mod_typical, events = events_pk,
  keep = c("WT", "TUM_VOL", "CYCLE"), useLinCmt = FALSE
) |>
  as.data.frame() |>
  distinct(time, .keep_all = TRUE) |>
  # A single-subject solve returns no `id` column; PKNCA needs one.
  mutate(id = 1L)

# Population PD simulation. Reseeding immediately before each solve gives the
# two regimens common random numbers over identical covariates.
solve_pd <- function(ev, label, seed = 20240426) {
  set.seed(seed)
  out <- rxode2::rxSolve(
    mod, events = ev, keep = c("WT", "TUM_VOL", "CYCLE"),
    useLinCmt = FALSE
  ) |>
    as.data.frame()
  stopifnot(length(unique(out$id)) == n_per_arm)   # rxSolve can silently drop subjects
  out$regimen <- label
  out
}

sim_pd <- bind_rows(Map(solve_pd, events_pd, names(events_pd)))
```

# Replicate published figures

## Figure 1 -- radioligand disposition in the four observed compartments

```{r fig-disposition}
sim_pk |>
  select(time, Cc, Asalivary_gland, Akidney, Atumor) |>
  tidyr::pivot_longer(-time, names_to = "channel", values_to = "value") |>
  filter(time <= 336) |>
  mutate(channel = recode(channel,
    Cc = "Central (MBq/L)", Asalivary_gland = "Salivary glands (MBq)",
    Akidney = "Kidneys (MBq)", Atumor = "Tumor (MBq)")) |>
  ggplot(aes(time, value)) +
  geom_line(linewidth = 0.7) +
  facet_wrap(~channel, scales = "free_y") +
  scale_x_continuous(breaks = seq(0, 336, by = 48)) +
  labs(x = "Time after injection (h)", y = NULL,
       title = "Typical patient, one 7378 MBq cycle",
       subtitle = "Reproduces the compartments of Figure 1 / Figure 2 of Siebinga 2024") +
  theme_bw()
```

## The cycle effect on tumor uptake

The paper's headline PK finding is that tumor uptake declines over cycles, to
73%, 50% and 44% of the cycle-1 rate in cycles 2, 3 and 4-7 (Table 2,
Equation 8). Reading the tumor uptake rate constant back out of the packaged
model confirms the encoding.

```{r cycle-effect}
cycle_check <- lapply(1:4, function(cy) {
  ev <- events_pk |> mutate(CYCLE = cy)
  rxode2::rxSolve(mod_typical, events = ev, useLinCmt = FALSE) |>
    as.data.frame() |>
    summarise(cycle = cy, kin_tumor = first(kin_tumor))
}) |>
  bind_rows() |>
  mutate(
    fraction_of_cycle1 = kin_tumor / kin_tumor[cycle == 1],
    published = c(1, 0.731, 0.498, 0.436)
  )

cycle_check |>
  dplyr::rename("Cycle" = cycle, "kin_tumor (1/h)" = kin_tumor,
                "Fraction of cycle 1" = fraction_of_cycle1,
                "Published fraction (Table 2)" = published) |>
  knitr::kable(digits = 5,
               caption = "Structural cycle effect on the tumor uptake rate constant.")

stopifnot(all(abs(cycle_check$fraction_of_cycle1 - cycle_check$published) < 1e-6))
```

## Structural effect of tumor volume on tumor uptake

The paper states that "a two-times higher volume resulted in a 2.11-fold
increased tumor uptake rate" (Results), which is the power exponent 1.08.

```{r tumvol-effect}
two_fold <- vapply(c(44300, 88600), function(tv) {
  ev <- events_pk |> mutate(TUM_VOL = tv)
  rxode2::rxSolve(mod_typical, events = ev, useLinCmt = FALSE) |>
    as.data.frame() |>
    (\(d) d$kin_tumor[1])()
}, numeric(1))

cat(sprintf("Doubling tumor volume multiplies kin_tumor by %.3f (paper: 2.11)\n",
            two_fold[2] / two_fold[1]))
stopifnot(abs(two_fold[2] / two_fold[1] - 2^1.08) < 1e-6)
```

## Figure 5 -- simulated PSA dynamics for the two dosing regimens

The paper plots the median simulated PSA with a 50% prediction interval for
both schedules and observes "only slight differences", with the
"2 x 2 - repeated" arm showing a more prominent decrease just after the second
cycle.

```{r fig5, fig.height = 4.5}
psa_summary <- sim_pd |>
  group_by(regimen, time) |>
  summarise(
    median = stats::median(PSA),
    lo = stats::quantile(PSA, 0.25),
    hi = stats::quantile(PSA, 0.75),
    .groups = "drop"
  )

ggplot(psa_summary, aes(time / week, median, fill = regimen, colour = regimen)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.2, colour = NA) +
  geom_line(linewidth = 0.8) +
  geom_vline(xintercept = 24, linetype = "dashed", linewidth = 0.3) +
  scale_x_continuous(breaks = seq(0, 32, by = 4)) +
  labs(x = "Time after start of treatment (weeks)", y = "PSA (ug/L)",
       colour = "Regimen", fill = "Regimen",
       title = "Simulated PSA dynamics (median and 50% prediction interval)",
       subtitle = "Replicates Figure 5 of Siebinga 2024; dashed line is the 24-week response endpoint") +
  theme_bw() + theme(legend.position = "top")
```

## Simulated response rates

The paper's quantitative simulation endpoints (Results, "Simulations"):
response is a >= 50% decrease in PSA at 24 weeks and stable disease is no
increase in PSA at 24 weeks.

```{r response-rates}
response <- sim_pd |>
  filter(time %in% c(0, 24 * week)) |>
  select(regimen, id, time, PSA) |>
  tidyr::pivot_wider(names_from = time, values_from = PSA,
                     names_prefix = "t") |>
  mutate(change = (.data[[paste0("t", 24 * week)]] - t0) / t0) |>
  group_by(regimen) |>
  summarise(
    response = 100 * mean(change <= -0.5),
    stable = 100 * mean(change <= 0),
    .groups = "drop"
  ) |>
  mutate(
    published_response = c(42.4, 44.7)[match(regimen, c("4 x 6", "2 x 2 - repeated"))],
    published_stable = c(54.8, 56.4)[match(regimen, c("4 x 6", "2 x 2 - repeated"))]
  )

response |>
  dplyr::rename(
    "Regimen" = regimen,
    "Response, simulated (%)" = response,
    "Stable disease, simulated (%)" = stable,
    "Response, published (%)" = published_response,
    "Stable disease, published (%)" = published_stable
  ) |>
  knitr::kable(digits = 1,
               caption = "Treatment response at 24 weeks vs. the published simulation (n = 2000 in the paper, 200 per arm here).")
```

Both proportions land within a few percentage points of the published values,
and the ordering (the "2 x 2 - repeated" arm slightly ahead) is reproduced. The
residual gap is consistent with the smaller cohort and with the covariate
distributions being reconstructed from medians and ranges rather than from the
original patient-level data.

# PKNCA validation

Two NCA passes. The first is a conventional NCA of the central-compartment
activity concentration, whose clearance is directly comparable to the renal
clearance the paper quotes. The second integrates the tumor activity
concentration to the cumulative tumor exposure that the paper's threshold
analysis uses.

```{r pknca-central}
conc_central <- sim_pk |>
  filter(!is.na(Cc)) |>
  select(id, time, Cc) |>
  mutate(treatment = "7378 MBq single cycle")

conc_central <- bind_rows(
  conc_central,
  conc_central |> distinct(id, treatment) |> mutate(time = 0, Cc = 0)
) |>
  distinct(id, treatment, time, .keep_all = TRUE) |>
  arrange(id, treatment, time)

dose_obj <- PKNCA::PKNCAdose(
  events_pk |> filter(evid == 1L) |>
    transmute(id, time, amt, treatment = "7378 MBq single cycle"),
  amt ~ time | treatment + id
)

nca_central <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  PKNCA::PKNCAconc(conc_central, Cc ~ time | treatment + id),
  dose_obj,
  intervals = data.frame(start = 0, end = Inf, cmax = TRUE, tmax = TRUE,
                         aucinf.obs = TRUE, cl.obs = TRUE, half.life = TRUE)
))

nca_central_df <- as.data.frame(nca_central) |>
  filter(PPTESTCD %in% c("cmax", "tmax", "aucinf.obs", "cl.obs", "half.life"))

nca_central_df |>
  select(PPTESTCD, PPORRES) |>
  dplyr::rename("NCA parameter" = PPTESTCD, "Simulated" = PPORRES) |>
  knitr::kable(digits = 3,
               caption = "NCA of the simulated central-compartment activity concentration (typical patient).")
```

```{r pknca-tumor}
conc_tumor <- sim_pk |>
  filter(!is.na(Ctumor)) |>
  select(id, time, Ctumor) |>
  mutate(treatment = "7378 MBq single cycle")

conc_tumor <- bind_rows(
  conc_tumor,
  conc_tumor |> distinct(id, treatment) |> mutate(time = 0, Ctumor = 0)
) |>
  distinct(id, treatment, time, .keep_all = TRUE) |>
  arrange(id, treatment, time)

nca_tumor <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  PKNCA::PKNCAconc(conc_tumor, Ctumor ~ time | treatment + id),
  dose_obj,
  intervals = data.frame(start = 0, end = 1008, auclast = TRUE, cmax = TRUE)
))

tumor_auc_cycle1 <- as.data.frame(nca_tumor) |>
  filter(PPTESTCD == "auclast") |>
  pull(PPORRES)

# The paper's threshold is the cumulative tumor AUC of the FIRST TWO cycles;
# cycle 2 carries the 0.731 uptake fraction of Table 2.
two_cycle_auc_mbq_h_per_ml <- tumor_auc_cycle1 * (1 + 0.731) / 1000

cat(sprintf(
  "Typical-patient tumor AUC, cycles 1+2: %.0f MBq*h/mL\nPublished response thresholds: 709.5 (nadir PSA) and 1188 (last PSA)\n",
  two_cycle_auc_mbq_h_per_ml))
```

The typical (median-tumor-volume, median-weight, no-random-effects) patient
lands **between the paper's two published thresholds**: just above the
709.5 MBq*h/mL required for a >= 50% decrease at the PSA *nadir*, and well below
the 1188 MBq*h/mL required for a >= 50% decrease at the *last* measurement after
two cycles. That is the expected placement. The published response rate at the
later, stricter endpoint is 42-45%, so the median patient should fall short of
it, while a majority transiently reach a 50% decrease at nadir. Whether an
individual clears either threshold is driven by the tumor-uptake random effects
(IIV 63.4% CV plus IOV 37.8% CV on `kin_tumor`).

```{r auc-threshold-check}
stopifnot(
  two_cycle_auc_mbq_h_per_ml > 709.5,
  two_cycle_auc_mbq_h_per_ml < 1188
)
```

# Comparison against published values

The paper reports no plasma NCA table, but the Discussion gives three derived
quantities that the NCA above reproduces directly, plus the tumor half-life.

The clearance of the central compartment is the one NCA parameter the paper
states outright: "Renal clearance was estimated 2.61 L/h" (Discussion). Because
renal excretion via `k10` is the model's only route of elimination, the
NCA-derived clearance must equal `k10 * V1` exactly, which makes this a sharp
check on the encoded `lkel` and `lvc`.

```{r nca-compare}
cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = data.frame(
    PPTESTCD = "cl.obs",
    PPORRES = nca_central_df$PPORRES[nca_central_df$PPTESTCD == "cl.obs"]
  ),
  reference = data.frame(cl.obs = 2.61),
  units = c(cl.obs = "L/h"),
  tolerance_pct = 20
)

knitr::kable(cmp, digits = 3,
             caption = paste("Simulated central-compartment clearance vs. the",
                             "renal clearance quoted in the Discussion of",
                             "Siebinga 2024."))
```

The terminal half-life returned by the NCA above is *not* `log(2) / k10`: the
slowest eigenvalue of the five-compartment system is the tumor efflux
(`k41` = 0.0150 1/h), so the central-compartment profile decays with the tumor
compartment in its terminal phase. The compartment-specific half-lives the
paper quotes are the individual efflux rate constants, tabulated next.

```{r structural-halflives}
tibble::tibble(
  Quantity = c("Organs at risk (salivary gland, k21)",
               "Organs at risk (kidney, k31)",
               "Tumor (k41)",
               "PSA growth rate kG (1/h)"),
  Model = c(log(2) / 0.0629, log(2) / 0.0625, log(2) / 0.0150, 0.000408),
  Published = c(11, 11, 46, 0.000408),
  `Published comparator` = c("~11 h (Discussion)", "~11 h (Discussion)",
                             "~46 h (Discussion)",
                             "0.000366 1/h in castration-resistant PCa (van Hasselt 2015)")
) |>
  knitr::kable(digits = 4,
               caption = "Structural half-lives and the PSA growth rate against the values quoted in the Discussion.")
```

# Assumptions and deviations

## Errata and value conflicts in the source

* **`kD,direct` is printed two ways.** The Results text states
  "K D, direct was estimated 0.000335 L*day-1*GBq-1 (40.1% RSE)" while Table 3
  gives 0.00335 (40.1% RSE, 95% CI 0.000961-0.006147) -- a factor of ten apart.
  **Table 3 is used.** Its three numbers are mutually consistent: an RSE of
  40.1% on 0.00335 implies a 95% CI of roughly 0.00072-0.0060, which is the
  interval printed; the same RSE on 0.000335 would imply an interval ten times
  narrower than the one printed. The text value is inconsistent with the
  interval printed beside it. The choice is also low-impact: because the
  delayed effect integrates to roughly ten times the direct effect, the
  simulated 24-week response rate moves by only about five percentage points
  between the two readings, and both remain compatible with the published
  42-45%.
* **The sign of the tumor-volume effect on salivary-gland uptake.** The Results
  narrative says "an increased tumor volume resulted in a lower salivary gland
  uptake rate (k12)", describing a tumor sink effect, but Table 2 prints the
  power exponent as +0.0910 with a 95% CI of 0.0350-0.544 that is entirely
  positive. **Table 2 is used** (a positive exponent), because the estimate,
  its RSE and the lower confidence bound all agree on the sign while the
  narrative is a single qualitative clause. The effect is weak either way: a
  ten-fold increase in tumor volume changes k12 by only 23%. Note also that the
  printed upper bound of 0.544 is not consistent with an RSE of 30.0%, which
  suggests a typographical error somewhere in that table row; the point
  estimate is unaffected.
* **Radioactive decay.** Methods say the PD driver is "the estimated
  concentration in the tumor compartment (taking into account radioactive
  decay)", while all radioactivity data were "corrected for decay to the time
  of injection". The model here takes the states to be decay-corrected (as in
  the predecessor `Siebinga_2023_lu177psma617`) and drives the PD directly from
  them, reading the parenthetical as a statement that decay was handled in the
  data correction. Encoding an *additional* explicit Lu-177 physical-decay term
  on every state was tested and rejected: it drops the simulated 24-week
  response rate to about 29% against the published 42-45%, far outside
  Monte-Carlo error, whereas the decay-corrected reading reproduces the
  published rates. This also keeps the model free of any constant not printed
  in the paper.
* **Supplement.** Additional file 1 (individual model fit results, Figures
  S1-S4) is not on disk. It contains individual prediction-vs-observation plots
  only; every parameter used here comes from Table 2, Table 3 or the Methods
  text, so nothing is missing from the extraction.

## Modelling assumptions

* **Allometric scaling.** Methods state that body weight was scaled "in
  relation to the median body weight, where exponents for volume of
  distribution (V) and rate constants (k) were set to 1 and -0.25". Both
  exponents are therefore encoded as `fixed()`, and the reference is the
  Table 1 median of 79 kg (the paper does not restate the centering value).
  `Bmax` is neither a volume nor a rate constant and the sentence names only V
  and k, so `Bmax` is left unscaled.
* **Covariate reference values.** The paper's covariate equations are all
  normalised to "the median value of the tested covariate" without restating
  it, so the Table 1 medians are used throughout: 79 kg for weight and
  0.0443 L (44300 mm^3) for segmented tumor volume.
* **Interoccasion variability.** nlmixr2 has no equivalent of NONMEM's
  `$OMEGA BLOCK(1) SAME`, so occasions 2-8 each carry their own eta with the
  variance fixed equal to the estimated occasion-1 variance. Eight occasions
  cover the maximum number of cycles in Table 1.
* **Box-Cox transformed eta.** The IIV on `kG` was Box-Cox transformed by the
  authors (Table 3 shape parameter -0.822). It is implemented in `model()` in
  the standard Petersson form, `phi = (exp(eta)^lambda - 1) / lambda` followed
  by `kG = kG_pop * exp(phi)`, because rxode2's `boxCox()` applies to residual
  error rather than to etas. The reported 90.9% CV is taken to be the variance
  of the untransformed eta, which is what NONMEM reports in `$OMEGA`.
* **Cycle counting.** `CYCLE` is a time-varying covariate that counts
  administrations. It carries both the structural Equation 8 cycle effect and
  the interoccasion-variability occasion index, and is decomposed into binary
  indicators inside `model()` rather than entering any equation as an integer.
* **Covariate distributions for the simulation.** The paper simulated 2000
  patients from "patient characteristic distributions of all included
  covariates" but does not publish those distributions. Weight is drawn
  normally and tumor volume log-normally, both matched to the Table 1 median
  and truncated to the Table 1 range. The simulated response rates are
  therefore expected to differ from the published ones by a few percentage
  points.
* **Screened but unused covariates.** Creatinine clearance (on `k10`) and age
  (on baseline PSA) were tested and not retained, and no coefficients are
  published for them; hematocrit is reported in Table 1 but never tested. All
  three are recorded in `covariatesDataExcluded` rather than `covariateData`.
* **Cycle effects on the salivary gland.** The Discussion mentions a cycle
  effect on salivary-gland uptake (84.6% in cycle 2, 98.1% in cycle 3), but
  Table 2 lists cycle effects only on `k14`, so the final model carries none on
  `k12` and neither does this encoding.
