---
title: "Ursodeoxycholic acid (Zuo 2016)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Ursodeoxycholic acid (Zuo 2016)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: Zuo P, Dobbins RL, O'Connor-Semmes RL, Young MA. A Systems
  Model for Ursodeoxycholic Acid Metabolism in Healthy and Patients With
  Primary Biliary Cirrhosis. CPT Pharmacometrics Syst Pharmacol.
  2016;5(8):418-426. doi:10.1002/psp4.12100
* Article: <https://doi.org/10.1002/psp4.12100>

The packaged model in
`inst/modeldb/specificDrugs/Zuo_2016_UDCA.R` is the 19-ODE
enterohepatic-recirculation systems model of Zuo et al. (Figure 1 +
Table 1) for the parent drug ursodeoxycholic acid (UDCA) and its two
hepatic conjugates, glycoursodeoxycholic acid (GUDCA) and
tauroursodeoxycholic acid (TUDCA). Compartments per analyte are
stomach (UDCA only), intestine, portal vein, blood (plasma), liver,
biliary system, and feces; all fluxes are first-order linear rate
constants in 1/hour. Two physiological modulators are encoded via
time-varying covariates: meal-induced and snack-induced scaling of the
biliary-to-intestine fluxes (gallbladder contraction surrogate), and a
disease-state indicator that switches the model from the healthy
parameter set to a primary biliary cirrhosis (PBC) adaptation by
multiplying the liver-to-biliary rate constants by paper-defined
disease-state factors.

The model is deterministic: there is no IIV and no residual error,
matching the source paper's status as a typical-value systems model fit
to mean published profiles from Xiang et al. 2011, Dilger et al. 2012,
and Hess et al. 2004.

# Population

UDCA is a minor bile acid in healthy human bile (< 1-5% of total
biliary bile acids) and is the FDA-approved oral treatment of PBC, an
autoimmune destruction of intrahepatic bile ducts. The model was
calibrated against mean concentration-time profiles from two healthy
adult cohorts and validated against a third:

* Calibration -- Xiang et al. (2011): n = 27 healthy adult volunteers,
  single 150 mg oral UDCA tablet, plasma UDCA + GUDCA + TUDCA sampled
  to 24 h with a post-dose warm meal at 4 h and two snacks at 7 h and
  10 h.
* Calibration -- Dilger et al. (2012): n = 11 healthy adults and n =
  11 PBC patients receiving 15 mg/kg/day oral UDCA for 3 weeks
  (modelled as 1000 mg/day for the cohort mean), with plasma and
  duodenal-bile sampling on day 22 (lunch / light snack / dinner at
  4 h / 7 h / 10 h after the morning dose).
* Validation -- Hess et al. (2004): n = 21 healthy adults receiving
  900 mg oral UDCA twice daily for 21 days, with plasma exposure
  (Cmax, Cavg, AUC0-4h) and integrated daily fecal output reported on
  day 21.

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

```{r read-model}
mod <- readModelDb("Zuo_2016_UDCA")
```

# Source trace

Per-parameter source comments appear inline in the model file. The
table below collects them in one place.

| Equation / parameter | Value | Source |
|----------------------|-------|--------|
| Stomach -> intestine K_SI | 16.61 1/h | Table 1 |
| Intestine -> portal K_IP,0 (UDCA) | 2.70 1/h | Table 1 |
| Intestine -> portal K_IP,1 (GUDCA, also TUDCA) | 0.32 1/h | Table 1 |
| Portal -> blood K_PB,0 (UDCA) | 0.61 1/h | Table 1 |
| Portal -> blood K_PB,1 (GUDCA, also TUDCA) | 0.36 1/h | Table 1 |
| Blood -> portal K_BP,0 (UDCA) | 6.94 1/h | Table 1 |
| Blood -> portal K_BP,1 (GUDCA, also TUDCA) | 0.66 1/h | Table 1 |
| Portal -> liver K_PL,0 (UDCA) | 0.82 1/h | Table 1 |
| Portal -> liver K_PL,1 (GUDCA) | 1.68 1/h | Table 1 |
| Portal -> liver K_PL,2 (TUDCA) | 3.98 1/h | Table 1 |
| Liver -> portal K_LP,1 (GUDCA) | 0.10 1/h | Table 1 |
| Liver -> portal K_LP,2 (TUDCA) | 0.03 1/h | Table 1 |
| Liver -> biliary K_LB,0 (UDCA) | 0.33 1/h | Table 1 |
| Liver -> biliary K_LB,1 (GUDCA, also TUDCA) | 0.32 1/h | Table 1 |
| Biliary -> intestine K_BI,0 (UDCA) | 0.64 1/h | Table 1 |
| Biliary -> intestine K_BI,1 (GUDCA, also TUDCA) | 0.13 1/h | Table 1 |
| Intestine -> feces K_IF,0 (UDCA) | 0.79 1/h | Table 1 |
| Intestine -> feces K_IF,1 (GUDCA, also TUDCA) | 0.21 1/h | Table 1 |
| UDCA -> GUDCA conjugation K_CONJ,1 | 54.61 1/h | Table 1 |
| GUDCA -> UDCA deconjugation K_DECONJ,1 | 0.76 1/h | Table 1 |
| UDCA -> TUDCA conjugation K_CONJ,2 | 3.83 1/h | Table 1 |
| TUDCA -> UDCA deconjugation K_DECONJ,2 | 0.17 1/h | Table 1 |
| Meal scaling E_meal on K_BI (1 h window) | 35.33 (unitless) | Table 1 |
| Snack scaling E_snack on K_BI (0.5 h window) | 9.53 (unitless) | Table 1 |
| Plasma volume V_plasma | 2.5 L | Methods, citing Hofmann 1983 |
| Bile volume V_bile | 0.075 L (75 mL) | Methods, citing Hofmann 1983 |
| UDCA molecular weight | 392.572 g/mol | Standard physical constant |
| GUDCA molecular weight | 449.626 g/mol | Standard physical constant |
| TUDCA molecular weight | 499.703 g/mol | Standard physical constant |
| Dose-dependent F (UDCA at 150 mg) | 0.66 | Methods (regression on log-dose; Walker 1992 / Crosignani 1991) |
| Dose-dependent F (UDCA at 1000 mg) | 0.31 | Methods (Walker 1992) |
| PBC K_LB,0 scaling | 0.10 (90% reduction) | Methods 'Model for PBC' / Figure 3 legend |
| PBC K_LB,1 scaling | 0.30 (70% reduction) | Methods 'Model for PBC' / Figure 3 legend |
| PBC K_LB,2 scaling | 0.10 (90% reduction) | Methods 'Model for PBC' / Figure 3 legend |
| Square-wave oral absorption duration | 0.5 h | Methods (dose enters stomach over 0.5 h) |
| Meal duration | 1 h | Methods (manually calibrated) |
| Snack duration | 0.5 h | Methods (manually calibrated) |

# Helpers: building event tables with meal / snack timing

Two helpers below build event tables with appropriate dosing rows and
time-varying meal / snack flags. They are vignette-local utilities (not
exported from the package).

```{r helpers}
# Build a vector of meal and snack interval edges (start and end times) for a
# given dosing-day schedule. Each meal / snack is represented as a (start, end)
# pair with the indicator switching to 1 at start and back to 0 at end.
make_food_intervals <- function(meal_times, snack_times,
                                meal_dur = 1, snack_dur = 0.5) {
  list(
    meal  = data.frame(start = meal_times,  end = meal_times  + meal_dur),
    snack = data.frame(start = snack_times, end = snack_times + snack_dur)
  )
}

# Evaluate MEAL_FLAG and SNACK_FLAG at an arbitrary time grid, given an
# intervals object produced by make_food_intervals().
food_flags_at <- function(times, intervals) {
  meal  <- vapply(times, function(t)
    as.integer(any(t >= intervals$meal$start  & t < intervals$meal$end)),
    integer(1))
  snack <- vapply(times, function(t)
    as.integer(any(t >= intervals$snack$start & t < intervals$snack$end)),
    integer(1))
  data.frame(MEAL_FLAG = meal, SNACK_FLAG = snack)
}

# Build an event table for one dose record per dosing day, with intra-day
# meal / snack timings replicated across N dosing days. Each dosing day is
# 24 h; t = 0 is the first dose, doses repeat at multiples of `dose_interval`.
make_chronic_events <- function(n_days,
                                dose_mg,
                                fracabs,
                                dis_pbc = 0L,
                                meal_times_relative,
                                snack_times_relative,
                                dose_interval = 24,
                                obs_per_day = 49,        # 30-min spacing
                                obs_dt_around_food = 0.25) {
  # Per-day meal / snack absolute times across the simulation horizon.
  day_starts <- (seq_len(n_days) - 1) * dose_interval
  meal_times  <- as.vector(outer(meal_times_relative,  day_starts, "+"))
  snack_times <- as.vector(outer(snack_times_relative, day_starts, "+"))
  intervals <- make_food_intervals(meal_times, snack_times)

  total_h <- n_days * dose_interval

  # Observation grid: uniform spacing + dense around food edges + meal/snack
  # boundaries to land MEAL_FLAG / SNACK_FLAG transitions on actual rows.
  uniform_grid <- seq(0, total_h, length.out = obs_per_day * n_days + 1)
  edge_times <- sort(unique(c(
    intervals$meal$start, intervals$meal$end,
    intervals$snack$start, intervals$snack$end,
    # add fine grid near food edges for cleaner plots
    intervals$meal$start  - obs_dt_around_food, intervals$meal$end  + obs_dt_around_food,
    intervals$snack$start - obs_dt_around_food, intervals$snack$end + obs_dt_around_food
  )))
  edge_times <- edge_times[edge_times >= 0 & edge_times <= total_h]
  obs_times  <- sort(unique(c(uniform_grid, edge_times)))

  # Dose rows: one per dosing day at the day start, dur = 0.5 h.
  dose_rows <- data.frame(
    id    = 1L,
    time  = day_starts,
    evid  = 1L,
    amt   = dose_mg,
    dur   = 0.5,
    cmt   = "stomach_udca",
    stringsAsFactors = FALSE
  )
  flags_dose <- food_flags_at(dose_rows$time, intervals)
  dose_rows$MEAL_FLAG  <- flags_dose$MEAL_FLAG
  dose_rows$SNACK_FLAG <- flags_dose$SNACK_FLAG
  dose_rows$FRACABS    <- fracabs
  dose_rows$DIS_PBC    <- dis_pbc

  # Observation rows.
  obs_rows <- data.frame(
    id    = 1L,
    time  = obs_times,
    evid  = 0L,
    amt   = 0,
    dur   = NA_real_,
    cmt   = NA_character_,
    stringsAsFactors = FALSE
  )
  flags_obs <- food_flags_at(obs_rows$time, intervals)
  obs_rows$MEAL_FLAG  <- flags_obs$MEAL_FLAG
  obs_rows$SNACK_FLAG <- flags_obs$SNACK_FLAG
  obs_rows$FRACABS    <- fracabs
  obs_rows$DIS_PBC    <- dis_pbc

  ev <- rbind(dose_rows, obs_rows)
  ev <- ev[order(ev$time, -ev$evid), ]
  ev$id <- 1L
  ev
}
```

# Replicate Figure 2 (top row) -- 150 mg UDCA single dose, healthy adults

The Xiang et al. 2011 single-dose study gave one 150 mg oral UDCA
tablet at t = 0 with a postdose warm meal at 4 h and snacks at 7 h and
10 h. Plasma UDCA, GUDCA, and TUDCA were measured to 24 h. The
fractional absorption at 150 mg is 0.66 per the paper's regression on
log dose.

```{r fig2-top, fig.cap = "Figure 2 (top row) replication: UDCA, GUDCA, and TUDCA plasma concentrations following a single 150 mg oral UDCA tablet (Xiang 2011 study design). Solid lines are the deterministic model prediction."}
ev_sd <- make_chronic_events(
  n_days               = 1,
  dose_mg              = 150,
  fracabs              = 0.66,
  meal_times_relative  = 4,
  snack_times_relative = c(7, 10)
)

sim_sd <- rxSolve(mod, ev_sd, keep = c("MEAL_FLAG", "SNACK_FLAG"))
sim_sd_long <- sim_sd |>
  as.data.frame() |>
  select(time, Cc_udca, Cc_gudca, Cc_tudca) |>
  pivot_longer(starts_with("Cc_"), names_to = "analyte", values_to = "C",
               names_prefix = "Cc_") |>
  mutate(analyte = toupper(analyte))

ggplot(sim_sd_long, aes(time, C, color = analyte)) +
  geom_line(linewidth = 1) +
  labs(x = "Time after morning dose (h)",
       y = "Plasma concentration (umol/L)",
       title = "150 mg UDCA single oral dose -- Xiang 2011 design",
       caption = "Replicates Figure 2 top row of Zuo 2016 (model only; observed digitised data not redistributed).") +
  scale_color_manual(values = c(UDCA = "#1f77b4", GUDCA = "#d62728", TUDCA = "#2ca02c"),
                     name = "Analyte") +
  theme_minimal()

# Summary metrics for cross-check against Xiang 2011 / Zuo Figure 2.
sim_sd_long |>
  group_by(analyte) |>
  summarise(Cmax      = round(max(C), 3),
            tmax_h    = round(time[which.max(C)], 2),
            C_24h     = round(C[which.min(abs(time - 24))], 3),
            .groups = "drop") |>
  knitr::kable(caption = "Simulated single-dose plasma summary metrics (deterministic typical-value).")
```

# Replicate Figure 2 (middle and bottom rows) -- chronic 1000 mg/day, healthy adults

The Dilger et al. 2012 study gave 1000 mg/day oral UDCA (modelling
the mean of the 15 mg/kg arm) for 21 days with lunch at +4 h, a light
snack at +7 h, and dinner at +10 h after the morning dose each day.
Plasma and duodenal bile were sampled on day 22. Fractional absorption
at 1000 mg is 0.31. To land near steady state without paying for a
full 21-day integration in the vignette time budget, we simulate
14 days and inspect the final 24-hour window.

```{r fig2-chronic, fig.cap = "Figure 2 (middle and bottom rows) replication: plasma and biliary UDCA / GUDCA / TUDCA concentration-time profiles on the 14th day of chronic 1000 mg/day oral UDCA (Dilger 2012 design)."}
ev_chronic <- make_chronic_events(
  n_days               = 14,
  dose_mg              = 1000,
  fracabs              = 0.31,
  meal_times_relative  = c(4, 10),   # lunch and dinner
  snack_times_relative = 7
)
sim_chronic <- rxSolve(mod, ev_chronic)
sim_chronic_df <- as.data.frame(sim_chronic)
# Keep the final 24-hour window (day 14 -> day 15) and shift to a 0-24 axis.
day14_start <- (14 - 1) * 24
sim_day14 <- sim_chronic_df |>
  filter(time >= day14_start, time <= day14_start + 24) |>
  mutate(time_in_day = time - day14_start)

# Plasma panel
sim_day14_plasma <- sim_day14 |>
  select(time_in_day, Cc_udca, Cc_gudca, Cc_tudca) |>
  pivot_longer(starts_with("Cc_"), names_to = "analyte", values_to = "C",
               names_prefix = "Cc_") |>
  mutate(matrix = "Plasma", analyte = toupper(analyte))

# Biliary panel
sim_day14_bile <- sim_day14 |>
  select(time_in_day, Cb_udca, Cb_gudca, Cb_tudca) |>
  pivot_longer(starts_with("Cb_"), names_to = "analyte", values_to = "C",
               names_prefix = "Cb_") |>
  mutate(matrix = "Biliary", analyte = toupper(analyte))

sim_day14_long <- bind_rows(sim_day14_plasma, sim_day14_bile)

ggplot(sim_day14_long, aes(time_in_day, C, color = analyte)) +
  geom_line(linewidth = 1) +
  facet_wrap(~ matrix, scales = "free_y") +
  labs(x = "Hours after morning dose on day 14 (h)",
       y = "Concentration (umol/L)",
       title = "Day-14 plasma and biliary profiles -- 1000 mg UDCA daily",
       caption = "Replicates Figure 2 middle (plasma) and bottom (biliary) rows of Zuo 2016 (chronic 1000 mg/day arm).") +
  scale_color_manual(values = c(UDCA = "#1f77b4", GUDCA = "#d62728", TUDCA = "#2ca02c"),
                     name = "Analyte") +
  theme_minimal()
```

# Comparison against Hess 2004 chronic dosing (validation)

Zuo et al. validated the calibrated model against Hess et al. (2004),
who gave 900 mg oral UDCA twice daily for 21 days and reported plasma
postprandial Cmax, AUC0-4h-derived Cavg, AUC0-4h, and average daily
fecal UDCA output. The paper's Table 2 best-model predictions are
reproduced here against the simulation.

```{r hess-validation}
ev_hess <- make_chronic_events(
  n_days               = 14,
  dose_mg              = 900,
  fracabs              = 0.31,           # 900 mg ~ 1000 mg per the paper's regression
  meal_times_relative  = c(4, 10),
  snack_times_relative = 7,
  dose_interval        = 12              # BID
)
sim_hess <- rxSolve(mod, ev_hess)
sim_hess_df <- as.data.frame(sim_hess)

day14_start <- (14 - 1) * 24
day14 <- sim_hess_df |>
  filter(time >= day14_start, time <= day14_start + 24) |>
  mutate(t_in_day = time - day14_start)

cmax_udca_uM    <- max(day14$Cc_udca)
cmax_udca_ugmL  <- cmax_udca_uM * 392.572 / 1000   # umol/L -> mg/L = ug/mL

# AUC0-4h after the morning dose, trapezoidal on a fine grid.
auc_grid <- day14 |> filter(t_in_day <= 4) |> arrange(t_in_day)
auc_0_4_uMmin <- sum(diff(auc_grid$t_in_day) *
                     (head(auc_grid$Cc_udca, -1) + tail(auc_grid$Cc_udca, -1)) / 2) * 60
auc_0_4_ugmLmin <- auc_0_4_uMmin * 392.572 / 1000

cavg_udca_uM    <- auc_0_4_uMmin / (4 * 60)
cavg_udca_ugmL  <- cavg_udca_uM * 392.572 / 1000

# Average daily fecal output (mg/day): difference in cumulative fecal mass
# across the final 24 h, summed across UDCA + GUDCA + TUDCA states (mmol)
# and converted to mg via each MW.
last_t <- max(day14$time)
prev_t <- min(day14$time)
feces_start <- day14 |> filter(time == prev_t) |> as.data.frame()
feces_end   <- day14 |> filter(time == last_t) |> as.data.frame()

# UDCA-only fecal output (the paper's Table 2 sums solid and aqueous phase
# UDCA; conjugates are not reported separately in Hess 2004).
fec_udca_mg <- (feces_end$feces_udca - feces_start$feces_udca) * 392.572

comparison <- tibble::tribble(
  ~Parameter,          ~Hess2004,            ~Zuo2016_bestmodel,   ~Simulation,
  "Cmax (ug/mL)",      "5.61 +/- 6.55",      "5.4 +/- 2.6",        sprintf("%.2f", cmax_udca_ugmL),
  "Cavg 0-4h (ug/mL)", "4.63 +/- 2.05",      "3.0 +/- 2.0",        sprintf("%.2f", cavg_udca_ugmL),
  "AUC0-4h (ug.min/mL)","1105 +/- 287",      "731 +/- 469",        sprintf("%.0f", auc_0_4_ugmLmin),
  "Fecal UDCA (mg/day)","13.5 + 136 = 149.5", "147 +/- 26",        sprintf("%.0f", fec_udca_mg)
)
knitr::kable(comparison,
             caption = paste("Hess 2004 reported values vs. Zuo 2016 best-model predictions",
                             "(paper Table 2) vs. this packaged-model deterministic simulation",
                             "(no IIV; +/- not applicable)."))
```

The deterministic simulation lies within the same order of magnitude
as the paper's best-model predictions for all four metrics. Differences
versus the paper's mean Table 2 values arise primarily because the
packaged model uses the single best-fit parameter set rather than the
bootstrap distribution, and because the vignette stops at day 14 rather
than running the full 21-day calibration horizon.

# Replicate Figure 3 -- adaptation to PBC patients

Zuo et al. adapted the healthy model to PBC by reducing the
liver-to-biliary rate constants: K_LB,0 by 90%, K_LB,1 by 70%, and
K_LB,2 by 90%. The packaged model encodes this via the `DIS_PBC`
covariate (0 = healthy, 1 = PBC). Below we re-run the day-14
1000 mg/day scenario in both healthy and PBC modes.

```{r fig3, fig.cap = "Figure 3 replication: day-14 plasma concentration-time profiles for the 1000 mg/day UDCA chronic regimen in healthy adults vs. PBC patients. Conjugate plasma levels rise 2- to 4-fold under PBC while UDCA itself changes less, matching the paper."}
ev_pbc <- make_chronic_events(
  n_days               = 14,
  dose_mg              = 1000,
  fracabs              = 0.31,
  dis_pbc              = 1L,
  meal_times_relative  = c(4, 10),
  snack_times_relative = 7
)
sim_pbc <- rxSolve(mod, ev_pbc)
sim_pbc_df <- as.data.frame(sim_pbc)

# Combine healthy day-14 from above and PBC day-14 from this run.
day14_pbc <- sim_pbc_df |>
  filter(time >= day14_start, time <= day14_start + 24) |>
  mutate(t_in_day = time - day14_start) |>
  select(t_in_day, Cc_udca, Cc_gudca, Cc_tudca) |>
  mutate(state = "PBC")

day14_healthy <- sim_chronic_df |>
  filter(time >= day14_start, time <= day14_start + 24) |>
  mutate(t_in_day = time - day14_start) |>
  select(t_in_day, Cc_udca, Cc_gudca, Cc_tudca) |>
  mutate(state = "Healthy")

day14_both <- bind_rows(day14_healthy, day14_pbc) |>
  pivot_longer(starts_with("Cc_"), names_to = "analyte", values_to = "C",
               names_prefix = "Cc_") |>
  mutate(analyte = toupper(analyte),
         state   = factor(state, levels = c("Healthy", "PBC")))

ggplot(day14_both, aes(t_in_day, C, color = state)) +
  geom_line(linewidth = 1) +
  facet_wrap(~ analyte, scales = "free_y") +
  labs(x = "Hours after morning dose on day 14 (h)",
       y = "Plasma concentration (umol/L)",
       title = "Day-14 plasma profiles: Healthy vs. PBC (1000 mg/day)",
       caption = "Replicates Figure 3 of Zuo 2016. K_LB,0 reduced 90%, K_LB,1 reduced 70%, K_LB,2 reduced 90% in PBC.") +
  scale_color_manual(values = c(Healthy = "#1f77b4", PBC = "#d62728")) +
  theme_minimal()

# Quantitative ratios of PBC vs Healthy Cmax for each analyte.
day14_both |>
  group_by(state, analyte) |>
  summarise(Cmax = max(C), .groups = "drop") |>
  pivot_wider(names_from = state, values_from = Cmax) |>
  mutate(ratio_PBC_over_Healthy = round(PBC / Healthy, 2)) |>
  knitr::kable(caption = "Plasma Cmax on day 14 -- PBC vs Healthy. Zuo 2016 reports 2-4x increase in conjugate levels under PBC with minimal change in UDCA itself.")
```

# Mass-balance / sanity check

For a single 150 mg dose, the total mmol entering the system across
the bioavailability hook (`f(stomach_udca) <- FRACABS / mw_udca`) is
`150 * 0.66 / 392.572 ~= 0.252 mmol`. At 24 hours after the dose, the
sum of all in-body amounts plus fecal output should match this value
(no other elimination pathways exist in the model).

```{r mass-balance}
sim_24 <- sim_sd |> as.data.frame() |> filter(time == 24) |> tail(1)

state_cols <- c("stomach_udca", "intestine_udca", "portal_udca", "blood_udca",
                "liver_udca", "biliary_udca",
                "intestine_gudca", "portal_gudca", "blood_gudca",
                "liver_gudca", "biliary_gudca",
                "intestine_tudca", "portal_tudca", "blood_tudca",
                "liver_tudca", "biliary_tudca")
feces_cols <- c("feces_udca", "feces_gudca", "feces_tudca")

in_body <- sum(unlist(sim_24[, state_cols]))
in_feces <- sum(unlist(sim_24[, feces_cols]))
dose_mmol <- 150 * 0.66 / 392.572

tibble::tibble(
  metric          = c("Dose absorbed (mmol)",
                      "In body at 24 h (mmol)",
                      "In feces at 24 h (mmol)",
                      "Sum (mmol)",
                      "Relative error (%)"),
  value           = c(sprintf("%.5f", dose_mmol),
                      sprintf("%.5f", in_body),
                      sprintf("%.5f", in_feces),
                      sprintf("%.5f", in_body + in_feces),
                      sprintf("%.3e", (in_body + in_feces - dose_mmol) / dose_mmol * 100))
) |>
  knitr::kable(caption = "Mass-balance check at t = 24 h for the 150 mg single-dose simulation.")
```

A relative error well below 0.01% confirms that the ODE topology
conserves mass across the seven compartments per analyte and the
hepatic conjugation / deconjugation cycle.

# Assumptions and deviations

This is a deterministic typical-value systems model with no IIV and no
residual error; the packaged model is intended for mechanistic
simulation, not for individual-level Bayesian fitting. The following
deviations from the standard nlmixr2lib conventions are intentional and
flagged here because `checkModelConventions()` warns on each.

* **Compartment naming.** The 19 compartments use mechanism-specific
  names (`stomach_udca`, `intestine_udca`, ..., `feces_tudca`) rather
  than the canonical `depot` / `central` / `peripheral1` / `peripheral2`
  / `effect` set. The canonical names cannot represent a system with
  seven anatomical compartments per analyte coupled through hepatic
  conjugation / deconjugation. The same precedent applies to the
  packaged `Shah_2012_mAb_PBPK` model, which uses
  tissue-and-sub-compartment names like `vp_li`, `eu_mu`, `bc_lu`.
* **Dose vs. concentration units.** `units$dosing = "mg"` and
  `units$concentration = "umol/L"` are dimensionally distinct (mass vs.
  amount-of-substance); they are bridged inside `model()` by the
  bioavailability hook `f(stomach_udca) <- FRACABS / mw_udca`, which
  applies both the dose-dependent fractional absorption and the
  mg-to-mmol unit conversion. `checkModelConventions()` warns on the
  unit mismatch because it cannot inspect the bioavailability hook.
* **No IIV, no residual error.** The source paper fit the model to mean
  published profiles using least-squares, not maximum-likelihood
  population fitting; only the typical-value parameter set + bootstrap
  CIs are reported. The vignette uses `rxSolve()` directly; downstream
  users who want to fit subject-level data would need to add IIV (and
  residual error) before passing the model to `nlmixr2()`.
* **PBC adaptation is a binary structural switch, not a continuous
  covariate.** The `DIS_PBC` indicator selects between two fixed
  parameter sets (healthy and PBC); intermediate disease severity is
  not represented. This matches the paper's hypothesis-testing scope.
* **Meal / snack windows are user-controlled.** `MEAL_FLAG` and
  `SNACK_FLAG` are time-varying covariates that the user supplies via
  the event table. The helpers in this vignette implement the Xiang
  2011 / Dilger 2012 schedules; alternative meal / snack timings can
  be tested by re-using `make_food_intervals()`.
* **Initial conditions = 0.** The paper used reported baseline
  concentrations in plasma and duodenal bile as initial conditions and
  set zero in other tissues. The packaged model leaves all initial
  conditions at zero; for chronic-dosing simulations this is reached
  through accumulation across dosing days. Single-dose simulations may
  therefore slightly underestimate the absolute concentrations at very
  early time points.
* **Fractional absorption supplied as a covariate.** The paper reports
  F as a function of dose via a Walker 1992 / Crosignani 1991
  regression; this is implemented in the packaged model as a `FRACABS`
  covariate that the user supplies per dose record, rather than as an
  estimable bioavailability parameter. Values 0.66 (150 mg), 0.31
  (1000 mg, 900 mg) are read off the paper.
* **No erratum applies.** A targeted search for corrections to Zuo et
  al. 2016 returned no errata as of the extraction date; the original
  publication is the canonical source.
