Bioequivalence and Bioavailability Forum

Hi Martin,

❝ I just come accross a topic relating to calculation of AUC from 0 to tau at steady state

Where?

❝ Should the concentration for time point zero (for calculation of AUC from 0 to tau) be somehow related to the pre-dose concentrations …

Yes.

❝ … or set to zero where I think setting to zero would be in-line with the principle of superposition (e.g. have a look here).

Why on earth would you ignore the measured pre-dose concentration and impute zero instead?
Since I’m guilty for that one …

… I had the script to generate it. See the one on steroids at the end.

Model    C0   Cmax     Cτ AUC0.25     AUCτ     AUC∞
      0.000 60.959 20.839   1.556  962.047 1322.815
NCA (AUC by lin-up/log-down trapezoidal)
Dose     C0   Cmax     Cτ AUC0.25     AUCτ     AUC∞
   1  0.000 60.959 20.839   1.508  960.939 1321.727
   2 20.839 77.648 26.048   6.681 1231.522
   3 26.048 81.953 27.351   7.974 1299.168
   4 27.351 83.030 27.676   8.297 1316.079
 Inf 27.785 83.389 27.785   8.405 1321.716
C0,ss forced to zero        4.932 1318.243

The lin-up/log-down trapezoidal performs nicely, since 1321.716 (AUC_ss,τ) ≈ 1321.727 (AUC_SD,∞). The latter matches the 1322.815 obtained by the model pretty well. The plot shows also why the trapezoidal rule underestimates increasing concentrations (the thin black line is the model).

If you force the pre-dose concentration to zero, you fall into a trap (see there for a similar case but the other way ’round).
Get the concentrations after the 4^th dose with superposition(tau, t, 4)
Correct (trapezoid)$$\small{AUC_{0-0.25}=0.5\times(0.25-0)\times(27.351+39.026)=8.297}$$Wrong (triangle)$$\small{AUC_{0-0.25}=0.5\times(0.25-0)\times39.026=4.878}$$Check

trap.linlog(t, superposition(tau, t, 4), 0, tau)$pAUC[2]
# [1] 8.297082
trap.linlog(t, c(0, superposition(tau, t, 4)[2:length(t)]), 0, tau)$pAUC[2]
# [1] 4.878223

---

one.comp <- function(f, D, V, k01, k10, t, p.t, t.last) { # no lagtime
  if (!isTRUE(all.equal(k01, k10))) { # common: k01 ≠ k10
    C    <- f*D*k01/(V*(k01-k10))*(exp(-k10*t)-exp(-k01*t))
    tmax <- log(k01/k10)/(k01-k10)
    Cmax <- f*D*k01/(V*(k01-k10))*(exp(-k10*tmax)-exp(-k01*tmax))
    AUC  <- f*D/V/k10
    x    <- micro.macro(f, D, V, k01, k10)
    pAUC <- (x$C[["A"]]-x$C[["A"]]*exp(-x$E[["alpha"]]*p.t))/x$E[["alpha"]]+
            (x$C[["B"]]-x$C[["B"]]*exp(-x$E[["beta"]]*p.t))/x$E[["beta"]]
    AUCt <- (x$C[["A"]]-x$C[["A"]]*exp(-x$E[["alpha"]]*t.last))/x$E[["alpha"]]+
            (x$C[["B"]]-x$C[["B"]]*exp(-x$E[["beta"]]*t.last))/x$E[["beta"]]
  } else {                            # flip-flop
    k    <- k10
    C    <- f*D/V*k*t*exp(-k*t)
    tmax <- 1/k
    Cmax <- f*D/V*k*tmax*exp(-k*tmax)
    AUC  <- f*D/V/k
    pAUC <- AUCt <- NA # no idea
  }
  res <- list(C = C, tmax = tmax, Cmax = Cmax,
              AUC = AUC, pAUC = pAUC, AUCt = AUCt)
  return(res)
}
micro.macro <- function(f, D, V, k01, k10) {
  # Coefficients (C) and exponents (E)
  C <- f*D*k01/(V*(k01-k10))
  C <- setNames(c(-C, +C), c("A", "B"))
  E <- setNames(c(k01, k10), c("alpha", "beta"))
  macro <- list(C = C, E = E)
  return(macro)
}
trap.linlog <- function(t, C, t.first, t.last) {
  pAUC.linlog <- 0
  for (j in seq(which(t == t.first), which(t == t.last))) {
    if (t[j] == t.first) {  # first triangle if 0
      pAUC.linlog <- c(pAUC.linlog, 0.5*(t[j]-t[j-1])*C[j])
      } else {
      if (C[j] >= C[j-1]) { # linear up
        pAUC.linlog <- c(pAUC.linlog, 0.5*(t[j]-t[j-1])*(C[j]+C[j-1]))
      } else {              # log down
        pAUC.linlog <- c(pAUC.linlog,
                         (t[j]-t[j-1])*(C[j]-C[j-1])/
                         (log(C[j]/C[j-1])))
      }
    }
  }
  AUC.linlog <- cumsum(pAUC.linlog)
  res <- list(pAUC = pAUC.linlog, AUC = tail(AUC.linlog, 1))
  return(res)
}
lz <- function(t, C) {
  m <- lm(log(C) ~ t)
  return(as.numeric(-coef(m)[2]))
}
superposition <- function(tau, t, N) {
  x <- micro.macro(f, D, V, k01, k10)
  if (!N == Inf) { # Nth dose: Gibaldi/Perrier (3.18)
    f.alpha <- (1-exp(-N*x$E[["alpha"]]*tau))/
                 (1-exp(-x$E[["alpha"]]*tau))
    f.beta  <- (1-exp(-N*x$E[["beta"]]*tau))/
                 (1-exp(-x$E[["beta"]]*tau))
    C <- f.alpha*x$C[["A"]]*exp(-x$E[["alpha"]]*t) +
         f.beta*x$C[["B"]]*exp(-x$E[["beta"]]*t)
  } else {         # Steady state: Gibaldi/Perrier (3.19)
    C <- x$C[["A"]]*exp(-x$E[["alpha"]]*t)/
           (1-exp(-x$E[["alpha"]]*tau)) +
         x$C[["B"]]*exp(-x$E[["beta"]]*t)/
           (1-exp(-x$E[["beta"]]*tau))
  }
  return(C)
}
# input #
f        <- 2/3
D        <- 400
V        <- 3.49
t12.01   <- 1
t12.10   <- 12
k01      <- log(2)/t12.01
k10      <- log(2)/t12.10
tau      <- 24
doses    <- 4
t        <- c(0, 0.25, 0.6, 1, 1.5, 1.83, 2.23, 2.71, 3.3,
              3.910868183, 4.89, 5.95, 7.24, 8.77798079, 10.9,
              13.53, 16.8, 20.84355645, 24)           # crazy but good looking
res      <- one.comp(f, D, V, k01, k10, t, t[2], tau) # model
t.first  <- t[!is.na(res$C)][1]
t.last   <- rev(t[!is.na(res$C)])[1]
lambda.z <- lz(t = tail(t, 3), tail(res$C, 3))
C0.sup   <- Ctau.sup <- setNames(rep(NA, doses-1), paste0("dose.", 2:doses))
Cmax.sup <- AUC.sup  <- setNames(rep(NA, doses-1), paste0("dose.", 2:doses))
pAUC.sup <- setNames(rep(NA, doses-1), paste0("dose.", 2:doses))
for (dose in 2:doses) {
  C                <- superposition(tau, dose)
  C0.sup[dose-1]   <- C[1]
  Ctau.sup[dose-1] <- tail(C, 1)
  Cmax.sup[dose-1] <- max(C)
  x                <- trap.linlog(t, C, 0, tau)
  AUC.sup[dose-1]  <- x$AUC
  pAUC.sup[dose-1] <- x$pAUC[2]
}
C.ss     <- superposition(tau, Inf)
AUC.ss   <- trap.linlog(t, C.ss, 0, tau)$AUC
C.ss.0   <- c(0, C.ss[2:length(C.ss)])
AUC.ss.0 <- trap.linlog(t, C.ss.0, 0, tau)$AUC
heading  <- paste("   C0   Cmax     C\u03C4 AUC0.25     AUC\u03C4     AUC\u221E")
txt <- paste("\nModel", heading, "\n",
             sprintf("%10.3f", head(res$C, 1)), sprintf("%6.3f", res$Cmax),
             sprintf("%6.3f", tail(res$C, 1)), sprintf("%7.3f", res$pAUC),
             sprintf("%8.3f", res$AUCt), sprintf("%8.3f", res$AUC))
txt <- paste(txt, "\nNCA (AUC by lin-up/log-down trapezoidal)\nDose ", heading,
             "\n  ", 1, sprintf("%6.3f", head(res$C, 1)),
             sprintf("%6.3f", max(res$C)), sprintf("%6.3f", tail(res$C, 1)),
             sprintf("%7.3f", trap.linlog(t, res$C, 0, tail(t, 1))$pAUC[2]),
             sprintf("%8.3f", trap.linlog(t, res$C, 0, tail(t, 1))$AUC),
             sprintf("%8.3f",
             trap.linlog(t, res$C, 0, tau)$AUC+tail(res$C, 1)/lambda.z))
for (j in seq_along(C0.sup)) {
  txt <- paste(txt, "\n  ", j+1,
               sprintf("%6.3f", C0.sup[j]), sprintf("%6.3f", Cmax.sup[j]),
               sprintf("%6.3f", Ctau.sup[j]), sprintf("%7.3f", pAUC.sup[j]),
               sprintf("%8.3f", AUC.sup[j]))
}
txt <- paste(txt, "\n Inf",
             sprintf("%6.3f", C.ss[1]), sprintf("%6.3f", max(C.ss)),
             sprintf("%6.3f", tail(C.ss, 1)),
             sprintf("%7.3f", trap.linlog(t, C.ss, 0, tail(t, 1))$pAUC[2]),
             sprintf("%8.3f", AUC.ss))
txt <- paste(txt, "\nC0,ss forced to zero",
             sprintf("%12.3f", trap.linlog(t, C.ss.0, 0, tail(t, 1))$pAUC[2]),
             sprintf("%8.3f", AUC.ss.0), "\n")
cat(txt, "\n")
# Steady state plot #
t1 <- t[1:2]
C1 <- C.ss[1:2]
C0 <- C.ss.0[1:2]
tm <- seq(0, t[2], length.out = 501)
Cm <- superposition(tau, tm, Inf)
dev.new(width = 4.5, height = 4.5)
op <- par(no.readonly = TRUE)
par(mar = c(3.4, 3.1, 0, 0.1), cex.axis = 0.9)
plot(t1, C1, type = "n", ylim = c(0, max(C1)),
     xlab = "", ylab = "", axes = FALSE)
mtext(expression(italic(t)*" (h)"), 1, line = 2.5)
mtext(expression(italic(C)[ss]*" (µg/mL)"), 2, line = 2)
grid(nx = NA, ny = NULL)
polygon(c(t1, rev(t1)), c(rep(0, length(C1)), rev(C1)),
        border = "#87CEFA80", col = "#87CEFA80")
polygon(c(t1, rev(t1)), c(rep(0, length(C0)), rev(C0)),
        border = "#FF808080", col = "#FF808080")
lines(tm, Cm) # model
lines(t1, C0, col = "red", lwd = 2)
lines(t1, C1, col = "blue", lwd = 2)
points(t1[1], C0[1], cex = 1.25, pch = 19, col = "red")
points(t1, C1, cex = 1.25, pch = 19, col = "blue")
axis(1, t1)
axis(2, las = 1)
box()
par(op)

PS: I derived the definite integral \(\small{\int_{0}^{t}C(t)\,\textrm{d}t}\) based on the micro-constants \(\small{A,\alpha,B,\beta}\) back in 1983 because I couldn’t find any reference. Ugly. If you can simplify it, let me know.