Multiple peaks: fallback to linear trapezoidal [NCA / SHAM]

posted by Helmut Homepage – Vienna, Austria, 2013-05-22 22:25 (3963 d 13:03 ago) – Posting: # 10620
Views: 17,657

Hi Yung-jin!

❝ However is there any special consideration when using lin-up/log-down trapezoidal rule to calculate AUCs if there is more than one absorption peak?


Not that I am aware of. If values increase again the algo falls back to the linear trapezoidal. Try this clumsy code (a bimodal release formulation with two lag-times, two absorption rate constants, and a common elimination):

#set.seed(29081957) # keep this line to compare results only
require(truncnorm)
n       <- 48        # no of profiles
t       <- c(0,0.25,0.5,0.75,1,1.25,1.5,2,2.5,3,3.5,
             4,4.25,4.5,4.75,5,5.25,5.5,6,6.5,7,7.75,9,12,16,24)
D       <- 40
D.f     <- c(0.5, 0.5) # dose fractions of components
# Index ".d" denotes default values, ".c" CV
F.d     <- 0.5      # fraction absorbed (BA)
F.c     <- 0.1      # ±0.1, uniform
V.d     <- 1        # volume of distribution
V.c     <- 0.1      # CV 10%, lognormal
k01.1.d <- 1.4      # absorption rate constant (1st component)
k01.2.d <- 0.7      # absorption rate constant (2nd component)
k01.c   <- 0.2      # CV 20%, lognormal
k10.d   <- log(2)/3 # elimination rate constant (common to both components)
k10.c   <- 0.2      # CV 20%, lognormal
tlag1.d <- 0.5      # 1st lag-time
tlag1.l <- 0        # lower truncation 1
tlag1.u <- 1.5      # upper truncation 1
tlag2.d <- 4        # 2nd lag-time
tlag2.l <- 3        # lower truncation 2
tlag2.u <- 6        # upper truncation 2
tlag.c  <- 0.5      # CV 50%, truncated normal
tmax.th <- max((log(k01.1.d/k10.d)/(k01.1.d-k10.d)+tlag1.d),
               (log(k01.2.d/k10.d)/(k01.2.d-k10.d)+tlag2.d-tlag1.d))
Cmax.th <- F.d*D*D.f[1]/V.d*k01.1.d/(k01.1.d-k10.d)*
             (exp(-k10.d*(tmax.th-tlag1.d))-exp(-k01.1.d*(tmax.th-tlag1.d)))+
           F.d*D*D.f[2]/V.d*k01.2.d/(k01.2.d-k10.d)*
             (exp(-k10.d*(tmax.th-tlag2.d))-exp(-k01.2.d*(tmax.th-tlag2.d)))
x       <- seq(min(t), max(t), length.out=((max(t)-min(t))*4+1))
C1      <- F.d*D*D.f[1]/V.d*k01.1.d/(k01.1.d-k10.d)*
             (exp(-k10.d*(x-tlag1.d))-exp(-k01.1.d*(x-tlag1.d)))
C1[C1<0]<- 0        # set negative values of 1st component to zero
C2      <- F.d*D*D.f[2]/V.d*k01.2.d/(k01.2.d-k10.d)*
            (exp(-k10.d*(x-tlag2.d))-exp(-k01.2.d*(x-tlag2.d)))
C2[C2<0]<- 0        # set negative values of 2nd component to zero
C       <- C1+C2    # add components
AErr    <- 0.10     # analytical error CV 10%, lognormal
LLOQ    <- 0.01*Cmax.th # 1% of theoretical Cmax
split.screen(c(2,1))
for (s in 1:2) {
  screen(s)
  if (s == 1){
    plot(x, C, type="l", lwd=3, col="red", ylim=c(0, Cmax.th*1.75),
      main="linear trapezoidal", cex.main=1.1,
      xlab="time", ylab="concentration")
  } else {
    plot(x, C, type="l", lwd=3, col="red", ylim=c(0, Cmax.th*1.75),
      main="lin-up/log-down trapezoidal", cex.main=1.1,
      xlab="time", ylab="concentration")
  }
  rug(t, side=3)
  abline(v=0, lty=3)
  abline(h=0, lty=3)
}
AUCt.th <- NULL; AUCt.lin <- NULL; AUCt.linlog <- NULL
for (i in 1:n) {
  F      <- runif(n=1, min=F.d-F.c, max=F.d+F.c)
  V      <- rlnorm(n=1, meanlog=log(V.d)-0.5*log(V.c^2+1),
              sdlog=sqrt(log(V.c^2+1)))
  k01.1  <- rlnorm(n=1, meanlog=log(k01.1.d)-0.5*log(k01.c^2+1),
              sdlog=sqrt(log(k01.c^2+1)))
  k01.2  <- rlnorm(n=1, meanlog=log(k01.2.d)-0.5*log(k01.c^2+1),
              sdlog=sqrt(log(k01.c^2+1)))
  k10    <- rlnorm(n=1, meanlog=log(k10.d)-0.5*log(k10.c^2+1),
              sdlog=sqrt(log(k10.c^2+1)))
  tlag1  <- rtruncnorm(n=1, a=tlag1.l, b=tlag1.u, mean=tlag1.d, sd=tlag.c)
  tlag2  <- rtruncnorm(n=1, a=tlag2.l, b=tlag2.u, mean=tlag2.d, sd=tlag.c)
  C0.1   <- F*D*D.f[1]/V  # fictive iv concentration at t=0
  B1.1   <- C0.1*k01.1/(k01.1-k10)*exp(k10*tlag1)
  B2.1   <- C0.1*k01.1/(k01.1-k10)*exp(k01.1*tlag1)
  C1     <- B1.1*exp(-k10*t)-B2.1*exp(-k01.1*t)
  C1[C1<0] <- 0
  C0.2   <- F*D*D.f[2]/V  # fictive iv concentration at t=0
  B1.2   <- C0.2*k01.2/(k01.2-k10)*exp(k10*tlag2)
  B2.2   <- C0.2*k01.2/(k01.2-k10)*exp(k01.2*tlag2)
  C2     <- B1.2*exp(-k10*t)-B2.2*exp(-k01.2*t)
  C2[C2<0] <- 0
  C <- C1+C2
  for (j in 1:length(t)) { # analytical error (multiplicative 10% C)
    if (j==1)
      AErr1 <- rnorm(n=1, mean=0, sd=abs(C[j]*AErr))
    else
      AErr1 <- c(AErr1, rnorm(n=1, mean=0, sd=abs(C[j]*AErr)))
  }
  C <- C+AErr1 # add analytical error
  C[C<LLOQ] <- NA  # assign NAs to Cs below LLOQ
  tfirst <- (t[!is.na(C)])[1]              # 1st C >= LLOQ
  t0     <- t[length(head((t[is.na(C)])))] # last LLOQ before tfirst
  tlast  <- rev(t[!is.na(C)])[1]           # last C >= LLOQ
  screen(1, FALSE)
  plot(t, C, type="l", axes=FALSE,
    ylim=c(0, Cmax.th*1.75), xlab="", ylab="") # linear spaghetti plots
  screen(2, FALSE)
  for (j in 1:(length(t)-1)) {
    if (is.na(C[j]) == FALSE & is.na(C[j+1]) == FALSE) {
      if (C[j+1] >= C[j]) {
        lines(x=c(t[j],t[j+1]), y=c(C[j],C[j+1])) # increasing segments
      } else {
        tk <- NULL; Ck <- NULL # interpolated segments
        for (k in seq(t[j], t[j+1], length.out=10)) {
          tk <- c(tk, k)
          Ck <- c(Ck, exp(log(C[j])+abs((k-t[j])/(t[j+1]-t[j]))*(log(C[j+1])-log(C[j]))))
        }
        lines(tk, Ck) # decreasing segments
      }
    }
  }
  Cwo.err1 <- function(x) {B1.1*exp(-k10*x)-B2.1*exp(-k01.1*x)}
  Cwo.err2 <- function(x) {B1.2*exp(-k10*x)-B2.2*exp(-k01.2*x)}
  AUCt.1 <- integrate(Cwo.err1, lower=tlag1, upper=max(t))
  AUCt.2 <- integrate(Cwo.err2, lower=tlag2, upper=max(t))
  AUCt <- AUCt.1$value+AUCt.2$value # AUCt without analytical error
  AUCt.th <- c(AUCt.th, AUCt)
  pAUC.lin <- 0 # linear trapezidal rule
  for (j in seq(which(t == t0), which(t == tlast)-1)){
    if (t[j] == t0) { # 1st triangle
      pAUC.lin <- pAUC.lin + 0.5*(t[j+1]-t[j])*C[j+1]
    } else {
      pAUC.lin <- pAUC.lin + 0.5*(t[j+1]-t[j])*(C[j+1]+C[j])
    }
  }
  AUCt.lin <- c(AUCt.lin, pAUC.lin)
  pAUC.linlog <- 0 # lin-up/log-down trapezidal rule
  for (j in seq(which(t == t0), which(t == tlast)-1)){
    if (t[j] == t0) { # 1st triangle
      pAUC.linlog <- pAUC.linlog + 0.5*(t[j+1]-t[j])*C[j+1]
      } else {
      if (C[j+1] >= C[j]) { # linear trap.
        pAUC.linlog <- pAUC.linlog + 0.5*(t[j+1]-t[j])*(C[j+1]+C[j])
      } else { # loglinear trap.
        pAUC.linlog <- pAUC.linlog + (t[j+1]-t[j])*(C[j+1]-C[j])/(log(C[j+1]/C[j]))
      }
    }
  }
  AUCt.linlog <- c(AUCt.linlog, pAUC.linlog)
}
close.screen(all = TRUE)
par(ask=T)
RE.AUCt.lin <- 100*(AUCt.lin-AUCt.th)/AUCt.th
RE.AUCt.linlog <- 100*(AUCt.linlog-AUCt.th)/AUCt.th
boxplot(RE.AUCt.lin, RE.AUCt.linlog, ylab="% RE",
  ylim=c(-max(RE.AUCt.lin, RE.AUCt.linlog)*1.1, max(RE.AUCt.lin, RE.AUCt.linlog)*1.1),
  main="Bias of algos", cex.main=1.1,
  names=c("linear", "lin-up/log-down"))
abline(h=0, lty=3)
cat("\nCVtotal; linear trap.:", round(sd(AUCt.lin)/mean(AUCt.lin),4),
  "– lin-up/log-down trap.:", round(sd(AUCt.linlog)/mean(AUCt.linlog),4), "\n\n")
wilcox.test(RE.AUCt.lin, mu=0, conf.int=T)
wilcox.test(RE.AUCt.linlog, mu=0, conf.int=T)
wilcox.test(RE.AUCt.lin, RE.AUCt.linlog, conf.int=T)
par(ask=F)


No big deal, but the lin-up/log-down should perform better if there are missings.

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes

Complete thread:

UA Flag
Activity
 Admin contact
22,957 posts in 4,819 threads, 1,636 registered users;
96 visitors (0 registered, 96 guests [including 6 identified bots]).
Forum time: 10:28 CET (Europe/Vienna)

With four parameters I can fit an elephant,
and with five I can make him wiggle his trunk.    John von Neumann

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5