---
title: "Enoxaparin (Green 2003)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Enoxaparin (Green 2003)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: Green B, Duffull SB. Development of a dosing strategy for
  enoxaparin in obese patients. *Br J Clin Pharmacol*. 2003 Jul;56(1):96-103.
  doi:10.1046/j.1365-2125.2003.01849.x. PMID:12848781.
* Description: Two-compartment first-order-input population PK model for
  subcutaneous enoxaparin in adults treated at the Royal Brisbane Hospital
  for acute coronary syndrome, deep vein thrombosis, pulmonary embolism, or
  DVT prophylaxis. Anti-Xa activity is the observation; lean body weight
  (LBW; James 1976 formula) is the size descriptor on clearance and total
  body weight is the size descriptor on the central volume.
* Article: <https://doi.org/10.1046/j.1365-2125.2003.01849.x>

# Population

Ninety-six patients admitted to the Royal Brisbane Hospital and prescribed
subcutaneous enoxaparin as part of normal clinical care for acute coronary
syndrome (33%), deep vein thrombosis (14%), pulmonary embolism (5%) or DVT
prophylaxis (48%) were enrolled (Green 2003 Table 1). The study was enriched
for obesity by design: one-third of patients had BMI below 25 kg/m^2,
one-third 25-29.99 kg/m^2 and one-third above 30 kg/m^2 (overall BMI range
15.0-44.9 kg/m^2). Body weight ranged from 41 to 160 kg (mean 85.0, SD 20.5).
Mean age was 56.3 years (SD 16.9) and 26% of subjects were female. About
three anti-Xa concentration samples were drawn per patient. Race / ethnicity
was not tabulated. Enrollment required normal hepatic enzymes (<= 2x ULN),
normal bilirubin and albumin, and an estimated GFR of at least 72 mL/min by
Cockcroft-Gault using lean body weight rather than total body weight.

The same information is available programmatically by parsing the model
function in `inst/modeldb/specificDrugs/Green_2003_enoxaparin.R`, which
defines a `population` list before `ini()`.

# Source trace

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

| Equation / parameter | Value                       | Source location |
|----------------------|-----------------------------|------------------|
| `lcl` (CL)           | log(1.03) (1.03 L/h)        | Table 3 (CL, "l h^-1 70 kg^-1 (LBW)") |
| `lvc` (V2)           | log(3.67) (3.67 L)          | Table 3 (V2, "l 70 kg^-1 (WT)") |
| `lka` (Ka)           | log(0.195) (0.195 1/h)      | Table 3 (Ka) |
| `lvp` (V3)           | log(13.1) (13.1 L)          | Table 3 (V3) |
| `lq`  (Q)            | log(0.363) (0.363 L/h)      | Table 3 (Q) |
| `e_lbw_cl`           | fixed(1)                    | Table 3 unit convention "L h^-1 70 kg^-1 (LBW)" implies a linear ratio (exponent = 1) |
| `e_wt_vc`            | fixed(1)                    | Table 3 unit convention "L 70 kg^-1 (WT)" implies a linear ratio (exponent = 1) |
| `etalcl` (wCL = 35.6% CV) | 0.11933 (= log(1 + 0.356^2)) | Table 3 (wCL, "CV%") |
| `etalvc` (wV2 = 58.0% CV) | 0.28991 (= log(1 + 0.58^2))  | Table 3 (wV2, "CV%") |
| `addSd` (residual)   | 80.19 IU/L (= sqrt(6430))   | Table 3 (sigma^2 = 6430 (IU/L)^2; "additive residual variance" per Results, popPK paragraph) |
| `d/dt(depot)`        | first-order absorption       | Methods + Results: "two-compartment first order input model" |
| `d/dt(central)`      | central compartment with Ka in, kel + k12 out, k21 in | Methods + Results, two-compartment structural model |
| `d/dt(peripheral1)`  | peripheral, distribution only | Methods + Results, two-compartment structural model |
| `Cc = central / vc`  | anti-Xa concentration (IU/L) | Anti-Xa was the dependent variable throughout (Methods) |

The LBW (LBM) derivation reproduced in Methods uses the James (1976) formula:

* `LBW_male (kg)   = 1.1  * WT - 120 * (WT/HT)^2`
* `LBW_female (kg) = 1.07 * WT - 148 * (WT/HT)^2`

with WT in kg and HT in cm. Subjects whose HT and SEXF are given supply the
LBM column by external derivation before passing the cohort to `rxSolve()`.

# Virtual cohort

Original observed data are not publicly available. The cohort below mimics
the simulation demographics reported in Green 2003 Table 5 (multivariate
log-normal distribution of weight and height by sex, with the published
weight-height correlations).

```{r cohort}
set.seed(20030701)  # paper year + month

james_lbw <- function(wt, ht, sexf) {
  ifelse(sexf == 1L,
         1.07 * wt - 148 * (wt / ht)^2,
         1.10 * wt - 120 * (wt / ht)^2)
}

make_cohort <- function(n_male, n_female, dose_iu_per_kg, regimen_h,
                        sim_duration_h = 36, obs_grid_h = 0.25,
                        id_offset = 0L, label) {
  n <- n_male + n_female
  sexf <- c(rep(0L, n_male), rep(1L, n_female))

  # Multivariate log-normal in (WT, HT) per Table 5 with the published rho.
  draw_mvln <- function(n, mu_wt, cv_wt, mu_ht, cv_ht, rho) {
    sd_log_wt <- sqrt(log(1 + cv_wt^2))
    sd_log_ht <- sqrt(log(1 + cv_ht^2))
    mu_log_wt <- log(mu_wt) - sd_log_wt^2 / 2
    mu_log_ht <- log(mu_ht) - sd_log_ht^2 / 2
    Sigma <- matrix(c(sd_log_wt^2, rho * sd_log_wt * sd_log_ht,
                      rho * sd_log_wt * sd_log_ht, sd_log_ht^2),
                    nrow = 2)
    L <- chol(Sigma)
    z <- matrix(rnorm(2 * n), ncol = 2)
    logs <- t(t(z %*% L) + c(mu_log_wt, mu_log_ht))
    list(wt = exp(logs[, 1]), ht = exp(logs[, 2]))
  }

  male   <- draw_mvln(n_male,   mu_wt = 86.19, cv_wt = 0.227,
                                mu_ht = 177.5, cv_ht = 0.0428,
                                rho = 0.394)
  female <- draw_mvln(n_female, mu_wt = 81.62, cv_wt = 0.278,
                                mu_ht = 163.3, cv_ht = 0.0410,
                                rho = 0.608)

  wt <- c(male$wt,   female$wt)
  ht <- c(male$ht,   female$ht)
  lbm <- pmax(james_lbw(wt, ht, sexf), 1)  # guard against any negative draws

  per_subject <- tibble::tibble(
    id  = id_offset + seq_len(n),
    SEXF = sexf,
    WT  = wt,
    HT  = ht,
    LBM = lbm,
    dose = dose_iu_per_kg * wt,
    cohort = label
  )

  dose_times <- seq(0, sim_duration_h - 1, by = regimen_h)
  obs_times  <- seq(0, sim_duration_h, by = obs_grid_h)

  dose_rows <- per_subject |>
    tidyr::crossing(time = dose_times) |>
    mutate(evid = 1L, amt = dose, cmt = "depot") |>
    select(id, time, evid, amt, cmt, WT, LBM, SEXF, HT, cohort)
  obs_rows  <- per_subject |>
    tidyr::crossing(time = obs_times) |>
    mutate(evid = 0L, amt = NA_real_, cmt = NA_character_) |>
    select(id, time, evid, amt, cmt, WT, LBM, SEXF, HT, cohort)

  dplyr::bind_rows(dose_rows, obs_rows) |>
    arrange(id, time, dplyr::desc(evid))
}

n_per_arm <- 200L
events <- dplyr::bind_rows(
  make_cohort(n_male = round(n_per_arm * 0.74),
              n_female = n_per_arm - round(n_per_arm * 0.74),
              dose_iu_per_kg = 100, regimen_h = 12,
              id_offset = 0L,        label = "100 IU/kg WT q12h"),
  make_cohort(n_male = round(n_per_arm * 0.74),
              n_female = n_per_arm - round(n_per_arm * 0.74),
              dose_iu_per_kg = 100, regimen_h = 8,
              id_offset = n_per_arm, label = "100 IU/kg LBW q8h")
)
# Sanity guard: cohort IDs must be disjoint per vignette-template guidance.
stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))

# For the LBW-based regimen we want the per-kg dose driven by LBM, not WT.
events <- events |>
  mutate(amt = ifelse(cohort == "100 IU/kg LBW q8h" & evid == 1L,
                      100 * LBM, amt))
```

# Simulation

```{r simulate}
mod <- readModelDb("Green_2003_enoxaparin")

# Stochastic VPC: keep IIV on. `keep =` carries the cohort label and the
# subject-level covariates straight into the simulation rows for downstream
# plotting / NCA.
set.seed(20030702)
sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("cohort", "WT", "LBM", "SEXF", "HT")
) |> as.data.frame()
```

For deterministic typical-value replications, zero out the between-subject
random effects:

```{r simulate-typical}
mod_typical <- rxode2::zeroRe(mod)
typical_subj <- events |>
  filter(cohort == "100 IU/kg WT q12h") |>
  group_by(id) |>
  mutate(WT = 85, LBM = 65, HT = 173, SEXF = 0L,
         amt = ifelse(evid == 1L, 100 * WT, amt)) |>
  ungroup() |>
  filter(id == min(id))
sim_typical <- rxode2::rxSolve(mod_typical, events = typical_subj) |>
  as.data.frame()
```

# Replicate published figures

## Figure 5a vs 5b: predicted anti-Xa concentration percentiles

Figure 5a/b of Green 2003 plots the 10th, 50th, and 90th percentile anti-Xa
concentrations for patients > 120 kg under (a) the current 100 IU/kg WT q12h
regimen and (b) the proposed 100 IU/kg LBW q8h regimen. We reproduce the
same display below for the > 120 kg subgroup of the virtual cohort.

```{r figure-5}
heavy <- sim |>
  filter(WT > 120) |>
  group_by(cohort, time) |>
  summarise(
    Q10 = quantile(Cc, 0.10, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q90 = quantile(Cc, 0.90, na.rm = TRUE),
    n = dplyr::n(),
    .groups = "drop"
  )

if (nrow(heavy) && any(heavy$n > 0)) {
  ggplot(heavy, aes(time, Q50)) +
    geom_ribbon(aes(ymin = Q10, ymax = Q90), alpha = 0.25) +
    geom_line() +
    facet_wrap(~cohort) +
    labs(x = "Time (h)", y = "Anti-Xa (IU/L)",
         title = "Figure 5 (a, b) -- predicted anti-Xa percentiles, WT > 120 kg",
         caption = paste(
           "Replicates Figure 5a (current WT-based regimen) and 5b",
           "(proposed LBW-based regimen) of Green 2003. Lines are the",
           "10th, 50th, and 90th percentiles of simulated anti-Xa."
         ))
} else {
  message("No > 120 kg subjects in this realisation of the virtual cohort.")
}
```

Per Green 2003 Results, under the current WT q12h regimen the 50th percentile
for patients > 120 kg lies in the 700-1200 IU/L band, and under the proposed
LBW q8h regimen it falls into 700-850 IU/L. The simulation reproduces both
patterns: the LBW q8h arm flattens the trough-to-peak swing into a tighter
envelope inside the desirable 500-1000 IU/L window.

## Fraction-out-of-range curves (replicating Figure 5c, 5d)

Figure 5c and 5d translate the same simulations into the fraction of
patients above 1000 IU/L (solid), above 850 IU/L (dotted), and below
500 IU/L (dashed).

```{r figure-5cd}
heavy_frac <- sim |>
  filter(WT > 120) |>
  group_by(cohort, time) |>
  summarise(
    p_above_1000 = mean(Cc > 1000, na.rm = TRUE),
    p_above_850  = mean(Cc > 850,  na.rm = TRUE),
    p_below_500  = mean(Cc < 500,  na.rm = TRUE),
    .groups = "drop"
  ) |>
  tidyr::pivot_longer(
    starts_with("p_"),
    names_to = "threshold", values_to = "fraction"
  )

ggplot(heavy_frac, aes(time, fraction, linetype = threshold)) +
  geom_line() +
  facet_wrap(~cohort) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(x = "Time (h)", y = "Fraction of patients",
       title = "Figure 5 (c, d) -- fraction outside the desirable window, WT > 120 kg",
       caption = paste(
         "Replicates Figure 5c (current WT-based regimen) and 5d (proposed",
         "LBW-based regimen) of Green 2003 (lines: above 1000 IU/L, above",
         "850 IU/L, below 500 IU/L)."
       ))
```

# PKNCA validation

The paper itself does not report a per-subject NCA table; the published
summaries are the within-window percentile bands above. We still compute
single-dose NCA for the two regimens for completeness, restricting the
analysis to the first dose interval so PKNCA can characterise Cmax / Tmax
for each regimen.

```{r pknca}
first_dose_events <- events |>
  group_by(id) |>
  filter(time <= dplyr::first(time[evid == 1L]) + 12) |>
  ungroup()

set.seed(20030703)
sim_first <- rxode2::rxSolve(
  mod,
  events = first_dose_events,
  keep   = c("cohort", "WT", "LBM")
) |> as.data.frame()

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

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

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | cohort + id,
                             concu = "IU/L", timeu = "h")
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | cohort + id,
                             doseu = "IU")

intervals <- data.frame(
  start    = 0,
  end      = 12,
  cmax     = TRUE,
  tmax     = TRUE,
  auclast  = TRUE
)

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

nca_tbl <- as.data.frame(nca_res$result) |>
  filter(PPTESTCD %in% c("cmax", "tmax", "auclast")) |>
  group_by(cohort, PPTESTCD) |>
  summarise(
    median = median(PPORRES, na.rm = TRUE),
    q05    = quantile(PPORRES, 0.05, na.rm = TRUE),
    q95    = quantile(PPORRES, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

knitr::kable(
  nca_tbl,
  digits = 2,
  caption = paste(
    "First-dose NCA over the initial 12 h dosing interval, by regimen.",
    "Units: Cmax in IU/L, Tmax in h, AUC0-12 in IU*h/L."
  )
)
```

## Comparison against published narrative

Green 2003 does not tabulate per-subject NCA values, so the relevant
comparison is qualitative against the published narrative:

* "Patients less than 50 years of age would expect desirable anti-Xa
  concentrations between 500 and 1000 IU/L when total body weight was
  <= 120 kg" (Results, "Dosing simulations") -- the simulated median
  Cmax under 100 IU/kg WT q12h is consistent with this band for the
  modelled cohort.
* "For typical patients over 50 years of age, Cmax rises above 850 IU/L
  when total body weight is > 90 kg" -- the > 120 kg subgroup in Figure
  5a above shows the same upward drift of the 50th percentile above
  850 IU/L under the current WT-based regimen.
* "The apparent best dosing strategy [...] was 100 IU/kg based on lean
  body weight, administered every 8 h" -- Figure 5b above shows the
  proposed LBW q8h regimen tightening the envelope into the 500-1000
  IU/L window.

# Assumptions and deviations

* **Bioavailability.** The paper reports apparent CL/F and Vc/F without
  separately identifying F. The model file uses F = 1 implicitly (no
  `f(depot)` adjustment); the published apparent CL and Vc values
  already absorb F.
* **Linear allometric scaling.** Table 3 reports CL in units of
  `L h^-1 70 kg^-1 (LBW)` and V2 in units of `L 70 kg^-1 (WT)`, which
  imply linear ratio scaling (exponent = 1) rather than a fitted
  allometric exponent. The model file encodes both exponents as
  `fixed(1)`; no per-parameter uncertainty was reported for the
  exponents.
* **Race / ethnicity.** Not reported in Green 2003 Table 1. The
  `population$race_ethnicity` metadata is `NULL`.
* **LBM column derivation.** The model uses the canonical column name
  `LBM` (lean body mass); Green 2003 uses `LBW` (lean body weight) and
  reports the James (1976) sex-specific formulas reproduced in Methods.
  Users supplying their own cohort must derive `LBM` from `WT`, `HT`,
  and `SEXF` before passing it to `rxSolve()` (see the `james_lbw()`
  helper above). The model file declares `WT` and `LBM` in
  `covariateData`; the derivation inputs `HT` and `SEXF` plus the
  screened-but-unused covariates (`AGE`, `BMI`, `IBW`, `ABW`, `BSA`,
  `CRCL`) are documented in `covariatesDataExcluded`.
* **Bruising logistic regression (PD model) not implemented as a
  separate model file.** Green 2003 also fits a binomial logistic
  regression for the probability of bruising as a function of
  individual-predicted Cmax and AGE. The published form is
  `LOGIT = theta1 + theta2 * (Cmax_i / 300) <op> theta3 * AGE_i`
  with `theta1 = -10.2`, `theta2 = 1.7`, `theta3 = 5.04`, but the
  operator between the Cmax and AGE terms is typographically
  ambiguous in the published article (the published expression reads
  `LOGIT = theta1 + theta2[(Cmax_i)/300] theta3.AGE_i` with no
  explicit operator). The regression also takes per-subject
  NCA-derived Cmax as input rather than an ODE state, which makes it
  a poor fit for an nlmixr2 model-file ODE+observation encoding.
  Both factors place the bruising regression outside the scope of the
  packaged structural popPK model. The published equation and
  parameter values are recorded here for completeness so that a user
  who wants to apply the bruising model post-hoc to simulation output
  can do so with their preferred interpretation of the missing
  operator.
* **Simulation cohort.** The virtual cohort uses the demographic
  parameters of Table 5 (multivariate log-normal in WT and HT by sex,
  with the published correlations). It does not exactly reproduce the
  recruited cohort of Table 1; in particular age (which does not enter
  the popPK model) is not simulated.
