Skip to contents

Model and source

This vignette validates the preclinical (rat, male Sprague-Dawley) blood-brain barrier (BBB) distributional model for morphine-3-glucuronide (M3G) from Xie et al. (2000). The paper distinguishes the BBB transport of unbound M3G across the rat brain interface and tests the impact of probenecid (an organic-anion transport inhibitor) on the influx versus efflux clearance components. The published model has two layers:

  • Model A (paper Figure 1, left panel) is a joint arterial / venous plasma PK fit used purely to generate a smooth plasma forcing function for Model B; the only structural plasma parameter the paper reports for Model A is CL_u = 3.8 mL/min (paper Results, p1787).
  • Model B (paper Figure 1, right panel, equation 2) is the two-compartment brain distributional model with asymmetric BBB exchange (CL_u,in vs CL_u,out) and an intercompartmental clearance Q_br between the sampled brain ECF (brain 1) and the deeper redistribution compartment (brain 2). All BBB parameters are reported per gram of brain in paper Table 3.

The packaged model collapses the paper’s two-compartment arterial / venous Model A into a one-compartment central plasma representation with the same total unbound clearance (CL_u = 3.8 mL/min) and a central volume V_c = CL_u * t_1/2,bl / ln(2) = 120.6 mL derived from the pooled venous blood half-life t_1/2,bl = 22 min (Table 2). The brain layer is faithful to Model B equation 2.

Population

Two groups of seven male Sprague-Dawley rats each (Charles River, Sweden), body weight 280-320 g. Each rat had chronic left-femoral-artery and left-femoral-vein cannulae (for blood sampling and intravenous M3G / probenecid dosing), a right-jugular-vein CMA/20 microdialysis probe (for blood ECF sampling), and a striatal CMA/12 microdialysis probe (for brain ECF sampling). All rats received an exponential intravenous infusion of M3G over 4 hours on each of two consecutive days, targeting a steady-state plasma concentration of 65 uM (= 30 ug/mL) using the Stanpump CCI system parameterised from Ekblom et al. 1993. On day 2 the probenecid arm additionally received a 70 umol/kg IV loading dose of probenecid followed by a 70 umol/kg/h constant infusion (~8 h total). Source: Xie 2000 Methods (paper p1785-1786, Table 1).

The same metadata are available programmatically via readModelDb("Xie_2000_m3g_rat")$population.

Source trace

The per-parameter origin is recorded inline next to each ini() entry in inst/modeldb/specificDrugs/Xie_2000_m3g_rat.R. The table below collects everything in one place for review.

Equation / parameter Value Source location
d/dt(central) n/a Paper Model A; one-compartment representation, see Assumptions and deviations
d/dt(brain_ecf) n/a Paper Model B equation 2
d/dt(brain_deep) n/a Paper Model B (brain 2 term in equation 2)
lcl (CL_u) log(0.0038) L/min Paper Results p1787: CL_u = 3.8 mL/min (joint Model A fit)
lvc (V_c) log(0.1206) L Derived: V_c = CL_u * t_1/2,bl / ln(2) = 3.8 * 22 / 0.6931 = 120.6 mL
lcluin reference log(0.00011) L/min/g Table 3: 0.11 uL/min/g-brain (without probenecid, RSE 25%)
e_conmed_probenecid_cluin log(0.17 / 0.11) Table 3: with-probenecid CL_u,in = 0.17 uL/min/g-brain (1.55-fold increase)
lcluout log(0.00115) L/min/g Table 3: CL_u,out = 1.15 uL/min/g-brain (common to both arms, RSE 34%)
lqbr log(0.00095) L/min/g Paper Results p1787: Q_br = 0.95 uL/min/g-brain
lvubr1 log(0.000028) L/g Paper Results p1787: V_u,br1 = 28 uL/g-brain
lvubr2 log(0.000205) L/g Paper Results p1787: V_u,br2 = 205 uL/g-brain (V_u,app - V_u,br1)
etalcluin 0.2401 (= 0.49^2) Table 3: IAV(CL_u,in) = 49% without probenecid (RSE 25%)
etalcluout 0.2500 (= 0.50^2) Table 3: IAV(CL_u,out) = 50% (RSE 65%)
propSd 0.14 Table 3 residual variability = 0.14 (RSE 20%); paper eq 6 confirms sigma is SD

Virtual cohort

Two arms of seven rats each: control (CONMED_PROBENECID = 0 for the whole experiment) and probenecid (CONMED_PROBENECID = 0 on day 1, switched to 1 on day 2 from the start of the probenecid loading dose). The validation here focuses on day-1 vs day-2 ratios within the probenecid arm, which is the design feature on which the paper’s main probenecid-effect inference rests. Each rat receives a constant-rate infusion that approximates the paper’s exponential Stanpump-driven infusion: the rate is set so that plasma reaches the target steady-state concentration of 65 uM by the end of the 4 h infusion window. Specifically R = CL_u * Css = 0.0038 L/min * 65 umol/L = 0.247 umol/min, with total dose R * 240 min = 59.3 umol. Brain ECF observations are simulated at 20-min intervals matching the paper’s microdialysis sampling cadence. The post-infusion follow-up is 2 h (paper Methods, p1786).

set.seed(2000)
n_per_arm <- 7L

# Infusion design: constant-rate IV approximation of the Stanpump exponential
# infusion. Rate chosen so that Cu_pl approaches the paper's target 65 uM by
# end of 4 h, given CL_u = 0.0038 L/min.
infusion_rate_umol_per_min <- 0.247    # = 0.0038 L/min * 65 umol/L
infusion_dur_min           <- 240
total_dose_umol            <- infusion_rate_umol_per_min * infusion_dur_min

# Observation grid: microdialysis samples every 20 min over 0 - 360 min.
obs_times <- seq(20, 360, by = 20)

make_rat <- function(id, conmed_probenecid_value) {
  # Two rows: an infusion-start dose row (evid = 1) and the observation rows.
  dose_row <- tibble(
    id   = id,
    time = 0,
    amt  = total_dose_umol,
    evid = 1L,
    rate = infusion_rate_umol_per_min,
    cmt  = "central",
    CONMED_PROBENECID = conmed_probenecid_value
  )
  obs_rows <- tibble(
    id   = id,
    time = obs_times,
    amt  = 0,
    evid = 0L,
    rate = 0,
    cmt  = "Cbrain_ecf",
    CONMED_PROBENECID = conmed_probenecid_value
  )
  bind_rows(dose_row, obs_rows)
}

events_rat <- bind_rows(
  bind_rows(lapply(seq_len(n_per_arm), function(i) {
    make_rat(id = i, conmed_probenecid_value = 0L)
  })) |> mutate(arm = "control"),
  bind_rows(lapply(seq_len(n_per_arm), function(i) {
    make_rat(id = n_per_arm + i, conmed_probenecid_value = 1L)
  })) |> mutate(arm = "probenecid (day 2)")
)

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

Simulation

mod <- readModelDb("Xie_2000_m3g_rat")

# Stochastic VPC-style simulation with the paper's IIV on CL_u,in and CL_u,out
# (etalcluin / etalcluout) plus the proportional residual error on brain ECF.
# n_per_arm = 7 matches the paper's cohort; for tighter percentile bands we
# replicate the cohort 50-fold to get smoother bands without lengthening
# wall-clock noticeably.
set.seed(2000)
events_replicates <- bind_rows(lapply(1:50, function(rep) {
  events_rat |>
    mutate(id = id + (rep - 1L) * (2L * n_per_arm))
}))
stopifnot(!anyDuplicated(unique(events_replicates[, c("id", "time", "evid")])))

sim <- rxode2::rxSolve(mod, events_replicates, keep = c("arm", "CONMED_PROBENECID")) |>
  as.data.frame()
#>  parameter labels from comments will be replaced by 'label()'

For deterministic typical-value replication (matching the paper’s population-prediction overlay in Figure 3a):

mod_typical <- mod |> rxode2::zeroRe()
#>  parameter labels from comments will be replaced by 'label()'
sim_typical <- rxode2::rxSolve(mod_typical, events_rat,
                               keep = c("arm", "CONMED_PROBENECID")) |>
  as.data.frame()
#>  omega/sigma items treated as zero: 'etalcluin', 'etalcluout'
#> Warning: multi-subject simulation without without 'omega'

Replicate published figures

# Replicates the shape of Figure 2 of Xie 2000: brain-ECF concentration-time
# profiles in control vs probenecid arms during and after the 4 h M3G
# infusion. The published Figure 2 plots mean +/- s.d. of n = 7 per arm;
# here we plot the 50-replicate stochastic VPC summary as median with
# 5th-95th percentile ribbon for the same arms.
vpc_sum <- sim |>
  filter(time > 0) |>
  group_by(arm, time) |>
  summarise(
    Q05 = quantile(Cbrain_ecf, 0.05, na.rm = TRUE),
    Q50 = quantile(Cbrain_ecf, 0.50, na.rm = TRUE),
    Q95 = quantile(Cbrain_ecf, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(vpc_sum, aes(time, Q50, colour = arm, fill = arm)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.2, colour = NA) +
  geom_line(linewidth = 0.8) +
  geom_line(data = sim_typical |> filter(time > 0),
            aes(time, Cbrain_ecf), linetype = "dashed", linewidth = 0.5) +
  labs(x = "Time (min)", y = "Unbound M3G in brain ECF (uM)",
       title = "Figure 2 replication: brain ECF M3G during 4 h IV M3G infusion",
       caption = "Replicates Figure 2 of Xie 2000 (control + probenecid arms). Dashed = typical-value (zeroed IIV) trace; solid = stochastic-VPC median; ribbon = 5th-95th percentile.")

# Plasma trace (the paper's Figure 2 also shows blood concentrations of
# M3G implicitly via the joint Model A fit; reproduced here as the plasma
# arm of our model). The plasma profile is identical in the two arms by
# construction (no probenecid effect on plasma kinetics; paper Results
# p1787 confirms blood CL and t_1/2 unchanged by probenecid).
sim_typical |>
  filter(time > 0) |>
  ggplot(aes(time, cu_pl, colour = arm)) +
  geom_line(linewidth = 0.8) +
  labs(x = "Time (min)", y = "Unbound M3G in arterial plasma (uM)",
       title = "Plasma M3G typical-value trace (Model A reduction)",
       caption = "Both arms overlay; probenecid does not affect plasma PK in the paper.")

Steady-state validation against paper Tables 2 and 3

The paper’s primary quantitative claim is that probenecid increases the steady-state unbound brain-ECF / blood concentration ratio approximately 2-fold (Table 2: 0.10 +/- 0.04 control vs 0.16 +/- 0.05 probenecid day 2). The packaged model reproduces this via the asymmetric BBB clearances: at steady state and in the absence of brain_deep redistribution, the ratio Cu_br_ecf / Cu_pl collapses to CL_u,in / CL_u,out, which equals 0.11 / 1.15 = 0.0957 without probenecid and 0.17 / 1.15 = 0.1478 with probenecid – equal to Table 3’s CL_u,in / CL_u,out ratios reported as 0.09 and 0.14.

end_inf <- sim_typical |>
  filter(time == 240) |>
  transmute(arm,
            cu_pl,
            cu_br_ecf,
            ratio = cu_br_ecf / cu_pl)

knitr::kable(
  end_inf |> distinct(),
  digits = 4,
  caption = "Typical-value brain ECF / plasma ratio at end of the 4 h M3G infusion."
)
Typical-value brain ECF / plasma ratio at end of the 4 h M3G infusion.
arm cu_pl cu_br_ecf ratio
control 64.9662 6.2141 0.0957
probenecid (day 2) 64.9662 9.6037 0.1478

The ratio at t = 240 min slightly under-reads the strict asymptotic CL_u,in / CL_u,out because the brain_deep redistribution compartment has not yet fully equilibrated by the end of the 4 h infusion (V_u,br2 = 205 uL/g-brain >> V_u,br1 = 28 uL/g-brain so the slow phase is the dominant time constant). The paper itself reports the same effect: brain ECF half-life of 81 min is >> blood half-life of 22 min (Table 2), driven by the deep brain compartment.

late <- sim_typical |>
  filter(time == 360) |>
  transmute(arm, cu_pl, cu_br_ecf, ratio = cu_br_ecf / cu_pl)
knitr::kable(
  late |> distinct(),
  digits = 4,
  caption = "Brain ECF / plasma ratio at t = 360 min (= 2 h post-infusion); plasma concentrations are decaying but the brain ECF is still equilibrating, so the empirical ratio shifts above the asymptote."
)
Brain ECF / plasma ratio at t = 360 min (= 2 h post-infusion); plasma concentrations are decaying but the brain ECF is still equilibrating, so the empirical ratio shifts above the asymptote.
arm cu_pl cu_br_ecf ratio
control 1.4811 0.1426 0.0963
probenecid (day 2) 1.4811 0.2204 0.1488

Comparison against paper Table 2 and Table 3

Quantity Paper value Model (this vignette)
Cu,ss,br / Cu,ss,bl, control (Table 2) 0.11 +/- 0.05 (d1), 0.10 +/- 0.02 (d2) 0.096 at t = 240 (typical-value)
Cu,ss,br / Cu,ss,bl, probenecid (Table 2) 0.16 +/- 0.05 (d2) 0.148 at t = 240 (typical-value)
CL_u,in / CL_u,out, control (Table 3) 0.09 (RSE 13%) 0.0957 (= 0.11 / 1.15, structural)
CL_u,in / CL_u,out, probenecid (Table 3) 0.14 (RSE 19%) 0.1478 (= 0.17 / 1.15, structural)
t_1/2,bl (Table 2) 22 +/- 2 min 22.0 min (V_c calibrated to match)
Body CL_u (Results, p1787) 3.8 mL/min 3.8 mL/min

All ratios reproduce the published values within rounding. The structural model is faithful to Xie 2000 Model B as written, with the only material reduction being the collapse of the joint arterial-venous Model A into a 1-compartment plasma surrogate (documented below).

Assumptions and deviations

  • Model A reduced to a one-compartment plasma surrogate. The paper’s Model A is a joint arterial / venous PK fit using both arterial-plasma and venous-blood data, with the venous volume V_v arbitrarily fixed to 1 to make the clearance from venous equal to the rate constant; the model is used purely as a smoothing fit to generate individual arterial-plasma forcing functions for Model B. The paper tabulates only the joint unbound clearance CL_u = 3.8 mL/min (Results, p1787); the inter-compartmental rate constants between arterial and venous and V_arterial are not reported. The packaged model collapses the arterial / venous distinction to a single central compartment driven by the same CL_u = 3.8 mL/min. The central volume is derived as V_c = CL_u * t_1/2,bl / ln(2) = 120.6 mL from the paper-reported pooled venous-blood terminal half-life of 22 min (Table 2). With this calibration, the packaged-model plasma trace matches the paper’s reported steady-state plasma target (65 uM at a 0.247 umol/min infusion) and reproduces the 22-min plasma half-life. Users who want to drive the brain layer with a more sophisticated plasma forcing function can re-implement the arterial-venous Model A as a separate rxode2 model and pass its solution into the brain layer as a covariate.
  • Plasma residual-error model not declared. Xie 2000 Table 3 lists residual variability only for the brain ECF fit (Model B, sigma = 0.14); the residual error attached to the arterial-venous Model A fit is not tabulated in the published paper. The packaged model therefore declares a residual-error model only on the brain ECF observation. The plasma concentration cu_pl is computed for diagnostic and simulation use but is not declared as an observation; users who need an explicit plasma error model can add one but should not treat the choice as paper-derived.
  • Single IIV magnitude on CL_u,in across treatment arms. Table 3 of the paper reports two different inter-animal variabilities for CL_u,in: 49% in the without-probenecid pool (RSE 25%) and 20% in the with-probenecid arm (RSE 47%). The current nlmixr2lib ini() syntax does not support arm-dependent IIV magnitudes directly; the packaged model uses the without-probenecid value (49%) for both arms. Users simulating probenecid-arm data should expect mildly over-dispersed brain-ECF percentiles relative to the paper’s narrower with-probenecid distribution.
  • Per-g-brain semantics in the brain layer. Xie 2000 expresses brain volumes (V_u,br1, V_u,br2) in uL/g-brain and brain clearances (CL_u,in, CL_u,out, Q_br) in uL/min/g-brain. The packaged model retains the per-g-brain semantics literally rather than scaling to whole-rat brain weight, because the per-g normalisation is the form in which the parameters are reported and the source paper itself does not state the rat brain weight used during the fit. A typical Sprague-Dawley rat brain weighs about 1.7-1.8 g, so users who want whole-brain amounts in the brain compartments should multiply brain_ecf and brain_deep by their assumed brain weight in g. The plasma central compartment remains in whole-rat (umol) semantics; the brain efflux back into plasma is negligible per paper Discussion p1789 (~0.02% of body amount at steady state) so the unit mismatch between plasma and brain does not affect mass-balance accuracy.
  • Compartment names brain_ecf and brain_deep are paper-named rather than canonical. The nlmixr2lib canonical compartment list includes central, peripheral1, peripheral2, plus a handful of named physiologic compartments (csf, isf, named brain regions, etc.). The Xie 2000 brain layer is asymmetric (CL_u,in != CL_u,out across the BBB) and the deep brain compartment connects only to the brain ECF, not to central – neither matches the symmetric peripheral semantics. Naming them brain_ecf and brain_deep is semantically clearer than peripheral1 / peripheral2 and follows the Stevens 2012 remoxipride precedent of using brain_ecf for the microdialysis-sampled compartment. This generates two checkModelConventions() warnings (compartments), accepted by design.
  • Observation variable named Cbrain_ecf, not Cc. Cc denotes the central plasma concentration in nlmixr2lib convention; the paper’s primary observation is the brain ECF concentration (microdialysis-derived), not plasma. Naming the observed variable Cbrain_ecf matches the Cbrain_<region> paper-named multi-output exception listed in compartment-names.md. This generates one checkModelConventions() warning (observation), accepted by design.
  • Stanpump exponential infusion approximated by a constant-rate infusion. The original study used the Stanpump CCI system parameterised from Ekblom et al. 1993 to deliver an exponential M3G infusion that rapidly reaches and holds the plasma target of 65 uM over 4 h. The packaged simulation uses a constant-rate infusion R = CL_u * 65 uM = 0.247 umol/min instead, which approximates the same steady-state plateau by the end of the 4 h window. Users replicating the early-time portion of the paper’s Figure 2 should expect a slightly different rise to plateau than the published trace (constant-rate rise is exponential with tau ~ V_c / CL_u); the steady-state plateau and brain-layer dynamics are unaffected.