---
title: "Liposomal amphotericin B (Hong 2006)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Liposomal amphotericin B (Hong 2006)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Hong_2006_amphotericinB_liposomal"))$reference`
* Description: `r rxode2::rxode(readModelDb("Hong_2006_amphotericinB_liposomal"))$description`
* Article: [Antimicrob Agents Chemother 2006;50(3):935-942](https://doi.org/10.1128/AAC.50.3.935-942.2006) (open access)

# Population

The model was developed from a prospective, single-centre, open-label
observational pharmacokinetic study of 39 paediatric oncology patients
at the Children's Hospital at Westmead (Sydney, NSW, Australia) (Hong
2006 Table 1). Patients were all neutropenic and febrile, with a mix of
diagnoses (15 acute lymphoblastic leukaemia, 4 acute myeloid leukaemia,
20 other), and 16 of the 39 had received a bone marrow transplant (12
allogeneic, 4 autologous). Liposomal amphotericin B (AmBisome) was
administered once daily as a 1-h intravenous infusion through a central
venous catheter at doses ranging from 0.8 to 5.9 mg/kg/day. A total of
637 plasma concentration observations were collected over 48 dosing
courses (mean 13.3 samples per course, range 3-30). Baseline
demographics (Hong 2006 Table 1, mean (SD), median, range): age 7.1
years (5.1), 6.5, 0.17-17; weight 28.8 kg (19.8), 21.1, 6.1-84.1; height
119.2 cm (33.7), 118.5, 61.5-190.0; body surface area 0.96 m^2 (0.46),
0.82, 0.32-2.11; serum creatinine 74 umol/L (108), 47, 19-659; gamma
glutamyl transferase 105 U/L (211), 39, 11-1398. The cohort was 67%
male. The model was fit in NONMEM 5.1.1 with first-order conditional
estimation, with lognormal IIV on CL, V1, Q, and V2, additional
between-occasion variability (BOV) on CL and V1, and a combined
additive plus proportional residual-error model. Robustness and accuracy
were verified by a bootstrap of 1000 successful replicates in Wings for
NONMEM.

The same information is available programmatically via
`readModelDb("Hong_2006_amphotericinB_liposomal")$population`.

# Source trace

Every numeric value in `ini()` carries an in-file comment pointing to
the Hong 2006 source location. The table below collects them in one
place for review.

| Equation / parameter           | Value         | Source location                                                |
|--------------------------------|---------------|----------------------------------------------------------------|
| `lcl` (CL)                     | 0.44 L/h      | Table 3 row "CL (l h-1)" (RSE 27%); Table 4 row "theta CL"      |
| `lvc` (V1)                     | 3.12 L        | Table 3 row "V1 (l)" (RSE 40%); Table 4 row "theta V1"          |
| `lq`  (Q)                      | 0.73 L/h      | Table 3 row "Q (l h-1)" (RSE 18%); Table 4 row "theta Q"        |
| `lvp` (V2)                     | 18.0 L        | Table 3 row "V2 (l)" (RSE 40%); Table 4 row "theta V2"          |
| `e_wt_cl` (theta_5)            | 0.0152 /kg    | Table 3 row "theta_5"; Table 4 row "theta CL ~ WT" (RSE 12%)    |
| `e_wt_vc` (theta_6)            | 0.0241 /kg    | Table 3 row "theta_6"; Table 4 row "theta V1 ~ WT" (RSE 29%)    |
| `etalcl`  (omega^2 CL = 0.0108)| 0.0108        | Table 3 CL IIV 10%; Table 4 omega^2 CL                          |
| `etalvc`  (omega^2 V1 = 0.0450)| 0.0450        | Table 4 bootstrap omega^2 V1 (original failed; see notes below) |
| `etalq`   (omega^2 Q = 0.593)  | 0.593         | Table 3 Q IIV 77%; Table 4 omega^2 Q                            |
| `etalvp`  (omega^2 V2 = 0.546) | 0.546         | Table 3 V2 IIV 74%; Table 4 omega^2 V2                          |
| `etalcl_bov_mult` (pi^2 = 0.209)| 0.209        | Table 3 CL BOV 46%; Table 4 pi^2 CL                             |
| `etalvc_bov_mult` (pi^2 = 0.319)| 0.319        | Table 3 V1 BOV 56%; Table 4 pi^2 V1                             |
| `propSd` (27% proportional)    | 0.27          | Table 3 row "Proportional (%)" (RSE 7%); Table 4 sigma^2 prop   |
| `addSd`  (0.034 mg/L additive) | 0.034         | Table 3 row "Additive (mg l-1)" (RSE 24%); Table 4 sigma^2 add  |
| WT centering (21 kg)           | 21            | Methods text after Table 2 footnote: WTmedian = 21 kg           |
| 2-cmt IV zero-order infusion   | n/a           | Abstract; Results "Population pharmacokinetic modeling" para 1  |
| Combined add + prop residual   | n/a           | Methods Eq. (combined error model); Results para 2              |
| Lognormal IIV on CL/V1/Q/V2    | n/a           | Methods "Development of base population pharmacokinetic model"  |
| BOV on CL and V1 (BLOCK(2) SAME)| n/a          | Methods Eq. (BOV definition); Results para 2                    |

IIV / BOV variance derivation. Hong 2006 reports IIV in Table 3 as a
percentage; the corresponding numbers in Table 4 are `omega^2` (the
variance on the log scale). The "%" values in Table 3 therefore
represent the log-scale SD (`omega`) expressed as a percent, so
`omega^2 ~ ("%"/100)^2`:

* CL: 10% IIV -> `omega^2 = 0.0108` (Table 4)
* V1: failed-to-estimate in Table 3; bootstrap `omega^2 = 0.0450` (Table 4)
* Q:  77% IIV -> `omega^2 = 0.593`  (Table 4)
* V2: 74% IIV -> `omega^2 = 0.546`  (Table 4)
* CL BOV: 46% -> `pi^2 = 0.209`     (Table 4)
* V1 BOV: 56% -> `pi^2 = 0.319`     (Table 4)

# Virtual cohort

Original observed data are not publicly available. The cohort below
covers four scenarios bracketing the paper's body-weight range: a
typical patient at the cohort median weight (WT = 21 kg), a small
infant (WT = 10 kg), a school-aged child (WT = 35 kg), and an
adolescent (WT = 60 kg). All scenarios receive a 1-h IV infusion of
liposomal amphotericin B once daily for 7 doses at 3 mg/kg/day (the
mid-range prophylactic regimen described in Hong 2006 Methods and
Results).

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

n_sub <- 200L

build_arm <- function(label, wt_kg, dose_mg_per_kg_per_day, id_offset) {
  ids <- id_offset + seq_len(n_sub)

  dose_amt_mg <- wt_kg * dose_mg_per_kg_per_day
  infusion_h  <- 1
  dose_times  <- seq(0, 144, by = 24)        # seven doses, once daily

  dose_rows <- tidyr::expand_grid(id = ids, time = dose_times) |>
    mutate(
      evid     = 1L,
      amt      = dose_amt_mg,
      cmt      = "central",
      rate     = dose_amt_mg / infusion_h,    # 1-h IV infusion
      cohort   = label,
      WT       = wt_kg,
      dose_mgkg = dose_mg_per_kg_per_day
    )

  obs_times <- c(seq(0, 24, by = 0.25),
                 seq(25, 168, by = 2),
                 seq(170, 240, by = 4))
  obs_rows <- tidyr::expand_grid(id = ids, time = obs_times) |>
    mutate(
      evid     = 0L,
      amt      = 0,
      cmt      = NA_character_,
      rate     = 0,
      cohort   = label,
      WT       = wt_kg,
      dose_mgkg = dose_mg_per_kg_per_day
    )

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

events <- bind_rows(
  build_arm("infant_10kg",      10, 3,    0L),
  build_arm("typical_21kg",     21, 3,  200L),
  build_arm("schoolage_35kg",   35, 3,  400L),
  build_arm("adolescent_60kg",  60, 3,  600L)
)

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

# Simulation

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

sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("cohort", "WT", "dose_mgkg")
) |> as.data.frame()
```

For typical-value comparisons against the Hong 2006 Table 3 point
estimates, also simulate with the random effects zeroed:

```{r simulate-typical}
mod_typical <- mod |> rxode2::zeroRe()

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

# Concentration-time profile (typical-value, 1.8 mg/kg)

Hong 2006 Figure 2 shows the representative concentration-time profile
of a patient treated with 1-h IV infusions of L-AmB at 1.8 mg/kg per
day over a single 24-h dosing interval. To reproduce that figure, we
simulate a typical 21-kg patient (cohort median) at the 1.8 mg/kg
regimen, then plot the first 24 hours after dose 1.

```{r figure-2, fig.width = 7, fig.height = 4.5}
events_fig2 <- build_arm("fig2_1.8mgkg", 21, 1.8, 0L)
sim_fig2    <- rxode2::rxSolve(mod_typical, events = events_fig2,
                               keep = c("cohort", "WT", "dose_mgkg")) |>
  as.data.frame()

sim_fig2 |>
  filter(time <= 24, id == min(id)) |>
  ggplot(aes(time, Cc)) +
  geom_line(linewidth = 0.9) +
  labs(
    x = "Time after first dose (h)",
    y = "L-AmB plasma concentration (mg/L)",
    title = "Replicates Hong 2006 Figure 2: typical-value Cc over one 24-h interval",
    subtitle = "1-h IV infusion of L-AmB at 1.8 mg/kg in a typical 21 kg patient"
  )
```

# Concentration-time profile by weight cohort

```{r profile-cohort, fig.width = 7, fig.height = 4.5}
sim |>
  group_by(cohort, time) |>
  summarise(
    Q05 = quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(time, Q50, colour = cohort, fill = cohort)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.2, colour = NA) +
  geom_line() +
  scale_y_log10() +
  labs(
    x = "Time (h)",
    y = "Simulated L-AmB Cc (mg/L, log scale)",
    title = "Simulated L-AmB by weight cohort, 3 mg/kg/day x 7 doses",
    caption = "Lines: median; ribbons: 5th-95th percentile across 200 simulated subjects."
  )
```

# PKNCA validation (steady-state, last dosing interval)

The block below characterises steady-state Cmax, Cmin (= Ctau), AUC0-tau,
average concentration, and the inferred terminal half-life from the
typical-value time course over the seventh dosing interval (`tau = 24
h`). The treatment grouping is `cohort`, matching the four
weight-stratified scenarios.

```{r pknca}
last_dose_time <- 144  # seventh dose at t = 144 h; tau = 24 h

sim_nca <- sim_typical |>
  filter(!is.na(Cc), time >= last_dose_time, time <= last_dose_time + 24) |>
  mutate(time_in_tau = time - last_dose_time) |>
  select(id, time = time_in_tau, Cc, cohort)

dose_df <- events |>
  filter(evid == 1, time == last_dose_time) |>
  select(id, time, amt, cohort) |>
  mutate(time = 0)

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

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

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

nca_summary <- summary(nca_res)
knitr::kable(
  nca_summary,
  caption = "Simulated steady-state NCA parameters (typical-value, seventh dose interval, tau = 24 h) by weight cohort."
)
```

## Terminal half-life (typical-value)

Hong 2006 reports a secondary terminal half-life `t1/2,beta = 59.4 SD
36.5 h` derived from individual Bayesian posterior estimates. For a
typical 2-compartment IV model with `CL = 0.44 L/h, V1 = 3.12 L, Q =
0.73 L/h, V2 = 18.0 L`, the two macro-rate constants `alpha` and `beta`
satisfy `alpha + beta = kel + k12 + k21` and `alpha * beta = kel *
k21`, giving `beta ~= 0.0142 1/h` and `t1/2,beta ~= 48.8 h` for the
typical (21 kg) subject. The figure below evaluates the post-infusion
log-linear decline numerically from the stopped-dose simulation.

```{r tbeta-check}
event_singledose <- build_arm("typical_21kg_single", 21, 3, 0L) |>
  filter(time <= 200, !(evid == 1 & time > 0))      # one dose only
sim_single <- rxode2::rxSolve(mod_typical, events = event_singledose,
                              keep = c("WT")) |>
  as.data.frame()

# Estimate terminal slope on natural log scale from the 72-200 h window
tail_fit <- sim_single |>
  filter(id == min(id), time >= 72, time <= 200, Cc > 0) |>
  mutate(lnCc = log(Cc))
fit <- lm(lnCc ~ time, data = tail_fit)
lambda_z <- -coef(fit)[["time"]]
tbeta_h <- log(2) / lambda_z

knitr::kable(
  tibble(
    quantity = c(
      "Simulated terminal slope lambda_z (1/h)",
      "Simulated terminal half-life t1/2,beta (h)",
      "Published t1/2,beta (mean SD 1/h)",
      "Analytical beta from typical THETAs (1/h)",
      "Analytical t1/2,beta from typical THETAs (h)"
    ),
    value = c(
      sprintf("%.4f", lambda_z),
      sprintf("%.1f", tbeta_h),
      "59.4 +/- 36.5",
      "0.0142",
      "48.8"
    )
  ),
  caption = "Simulated vs published terminal half-life for the typical 21 kg subject."
)
```

## Steady-state volume of distribution (typical-value)

Hong 2006 reports a secondary `Vss = 1.14 +/- 0.68 L/kg` from individual
Bayesian estimates. For the typical 21 kg subject the model gives `Vss
= (V1 + V2) / WT = (3.12 + 18.0) / 21 = 1.006 L/kg`, comparable to the
published mean.

```{r vss-check}
vc_typ <- 3.12
vp_typ <- 18.0
wt_typ <- 21
vss_lkg <- (vc_typ + vp_typ) / wt_typ
knitr::kable(
  tibble(
    quantity = c(
      "Simulated Vss for typical 21 kg subject (L/kg)",
      "Published Vss (mean SD 1/kg)"
    ),
    value = c(
      sprintf("%.3f", vss_lkg),
      "1.14 +/- 0.68"
    )
  ),
  caption = "Simulated vs published Vss for the typical 21 kg subject."
)
```

## Cmax,ss / MIC target attainment

Hong 2006 demonstrated that `Cmax,ss / MIC` discriminates between
complete and partial clinical responses to L-AmB in patients with
proven invasive fungal infection (Hong 2006 Table 5 and Figure 3). The
authors then used the population PK model to explore whether
weight-stratified dosing achieves the proposed target of `Cmax,ss /
MIC > 50` for canonical MICs of 0.5 mg/L (Candida) and 1.0 mg/L
(Aspergillus). The table below reproduces that calculation using the
typical-value `Cmax,ss` at the seventh dose interval for the four
weight cohorts at 3 mg/kg/day.

```{r cmax-mic-target}
cmax_ss <- sim_typical |>
  filter(time >= last_dose_time, time <= last_dose_time + 24) |>
  group_by(cohort, WT, dose_mgkg) |>
  summarise(Cmax_ss_mgL = max(Cc, na.rm = TRUE), .groups = "drop") |>
  mutate(
    Cmaxss_over_MIC_Candida    = round(Cmax_ss_mgL / 0.5, 1),
    Cmaxss_over_MIC_Aspergillus = round(Cmax_ss_mgL / 1.0, 1)
  )

knitr::kable(
  cmax_ss,
  digits = 2,
  caption = "Typical-value Cmax,ss and Cmax,ss/MIC ratios by weight cohort at 3 mg/kg/day, day 7. Hong 2006 proposed Cmax,ss/MIC >= 50 as the antifungal efficacy target."
)
```

# Assumptions and deviations

- **Between-occasion variability (BOV) encoded as IIV on multiplicative
  anchors.** Hong 2006 estimates BOV on CL (46%) and V1 (56%) via
  NONMEM's `$OMEGA BLOCK(2) SAME` (a single magnitude shared across all
  dosing occasions for a given subject). nlmixr2lib has no canonical
  pattern for occasion-keyed IOV without an `OCC` covariate column, so
  the BOV variance is carried here as ordinary IIV on fixed-at-1
  multiplicative anchors `lcl_bov_mult` and `lvc_bov_mult` (following
  the Bellanti 2015 deferoxamine precedent). The total log-scale
  variance for the typical subject is therefore `etalcl + etalcl_bov_mult
  = 0.0108 + 0.209 = 0.220` (CL) and `etalvc + etalvc_bov_mult = 0.045
  + 0.319 = 0.364` (V1). The model preserves the population-level
  spread but loses the within-subject between-occasion drift that the
  paper estimates (a single simulated subject does not switch
  parameter values between courses).

- **IIV on V1 from the bootstrap.** Hong 2006 Table 3 leaves the IIV
  V1 column blank because the original NONMEM fit did not estimate a
  stable `omega^2` for V1 (the Discussion notes that BOV exceeded IIV
  on V1, and NONMEM had difficulty disentangling them). The bootstrap
  mean from 1000 successful replicates (`omega^2 V1 = 0.0450`, Table
  4) is used here as the best available point estimate. Treat this
  number as a lower bound on the V1 IIV; the paper does not commit to
  a published value.

- **Weight is time-varying across dosing courses.** Hong 2006 Methods
  states that covariate values were updated per course; the model
  encodes WT as time-varying. For the simulations above we hold WT
  fixed within each cohort (matched to the median weight for that
  band). Users simulating real per-patient time-varying WT should
  supply the column at every record.

- **No erratum identified.** A Crossref search for "Hong 2006 liposomal
  amphotericin erratum / corrigendum" returned no relevant results in
  2026-05-31; the model uses the values as published in the original
  paper without correction.

- **Race / ethnicity not reported.** Hong 2006 does not report race /
  ethnicity for the cohort (the study was conducted at a single centre
  in Sydney). No race-based stratification or covariate is included.

- **Renal-function markers not retained.** The paper screened
  serum creatinine and other renal markers as covariates but did not
  retain them because renal elimination is a minor pathway for L-AmB
  (Hong 2006 Discussion). The model has no renal-function covariate.

- **Sex not retained.** Sex did not significantly affect L-AmB PK in
  the Mann-Whitney U test (Hong 2006 Results). The model has no sex
  covariate.
