Simulation framework [NCA / SHAM]

posted by mittyri – Russia, 2019-04-17 01:21 (588 d 17:34 ago) – Posting: # 20179
Views: 4,367

Hi Helmut,

» If t½ of active metabolite > t½ of parent, assess only the metabolite.

Could you please explain a little bit? When did I miss that good old times?

ready for simulation:

library(ggplot2)
# input paraemeters
Nsub <- 1000 # number of subjects to simulate
D        <- 400
ka       <- 1.39  # 1/h
ka.omega <- 0.1
Vd       <- 1     # L
Vd.omega <- 0.2
CL       <- 0.347 # L/h
CL.omega <- 0.15
t<- c(seq(0, 1, 0.25), seq(2,6,1), 8,10,12,16,24) # some realistic sequence
ratio <- 2^(seq(-3,3,0.2)) # ratios of ka(T)/ka/R)

# helper functions
C.sd <- function(F=1, D, Vd, ka, ke, t) {
  if (!identical(ka, ke)) { # common case ka != ke
    C <- F*D/Vd*(ka/(ka - ke))*(exp(-ke*t) - exp(-ka*t))
  } else {                  # equal input & output
    C <- F*D/Vd*ke*t*exp(-ke*t)
  }
  return(C)
}
AUCcalc <- function(t,C){
  linlogflag <- C[-length(C)] <= C[-1]
  AUCsegments <- ifelse(linlogflag,
                   diff(t)*(C[-1]+C[-length(C)])/2,
                   (C[-length(C)] - C[-1])*diff(t)/(log(C[-length(C)]) - log(C[-1])))
  return(sum(AUCsegments))
}

AbsorptionDF <- function(D, ka, Vd, CL,t,ratio){
  # Reference
  ke       <- CL/Vd
  C        <- C.sd(D=D, Vd=Vd, ka=ka, ke=ke, t=t)
  tmax     <- t[C == max(C)][1]
  Cmax     <- C.sd(D=D, Vd=Vd, ka=ka, ke=ke, t=tmax)
  AUC.t    <- AUCcalc(t, C)
  t.1      <- t[which(t <= tmax)]
  t.cut    <- max(t.1)
  C.1      <- C[which(t <= t.cut)]
  pAUC     <- AUCcalc(t.1, C.1)
  Cmax.AUC <- Cmax/AUC.t
 
  # Tests
  ka.t  <- ka*ratio                           # Tests' ka
  res   <- data.frame(kaR=ka, kaT_kaR=ratio, kaT=signif(ka.t, 5),
                      Cmax=NA, Cmax.r=NA, pAUC=NA, pAUC.r=NA,
                      Cmax_AUC=NA, Cmax_AUC.r=NA)
 
  for (j in seq_along(ratio)) {
    # full internal precision, 4 significant digits for output
    C.tmp    <- C.sd(D=D, Vd=Vd, ka=ka.t[j], ke=ke, t=t)
    if (!identical(ka.t[j], ke)) { # ka != ke
      tmax.tmp <- log(ka.t[j]/ke)/(ka.t[j] - ke)
    } else {                       # ka = ke
      tmax.tmp <- 1/ke
    }
    Cmax.tmp <- C.sd(D=D, Vd=Vd, ka=ka.t[j], ke=ke, t=tmax.tmp)
    res[j, "Cmax"]   <- signif(Cmax.tmp, 4)
    res[j, "Cmax.r"] <- signif(Cmax.tmp/Cmax, 4)
    AUC.t.tmp <- AUCcalc(t,C.tmp)
    t.1.tmp   <- t[which(t <= t.cut)]
    C.1.tmp   <- C.tmp[which(t <= t.cut)] # cut at tmax of R!
    pAUC.tmp  <- AUCcalc(t.1.tmp, C.1.tmp)
    res[j, "pAUC"]       <- signif(pAUC.tmp, 4)
    res[j, "pAUC.r"]     <- signif(pAUC.tmp/pAUC, 4)
    res[j, "Cmax_AUC"]   <- signif(Cmax.tmp/AUC.t.tmp, 4)
    res[j, "Cmax_AUC.r"] <- signif((Cmax.tmp/AUC.t.tmp)/Cmax.AUC, 4)
  }
  return(res)
}

SubjectsDF <- data.frame()
for(isub in 1:Nsub){
  # sampling parameters
  ka.sub       <- ka * exp(rnorm(1, sd = sqrt(ka.omega)))
  Vd.sub       <- Vd * exp(rnorm(1,sd = sqrt(Vd.omega)))
  CL.sub       <- CL * exp(rnorm(1,sd = sqrt(CL.omega)))
  DF.sub <- cbind(Subject = isub, V = Vd.sub, CL = CL.sub, AbsorptionDF(D, ka.sub, Vd.sub, CL.sub, t, ratio))
  SubjectsDF <- rbind(SubjectsDF, DF.sub)
}

SubjectsDFstack <-
  reshape(SubjectsDF[, -c(2,3,4,6,7,9,11)],
        direction = 'long', varying = 3:5, v.names = "ratio", timevar = "metric", times = names(SubjectsDF)[3:5]) # hate this one!

ggplot(SubjectsDFstack, aes(x=kaT_kaR, y=ratio, color=factor(metric)) ) +
  theme_bw() +
  geom_point(size=.3) +
  geom_smooth(method = 'loess', se = FALSE) +
  stat_density_2d(data = subset(SubjectsDFstack, metric == unique(SubjectsDFstack$metric)[1]), geom = "raster", aes(alpha = ..density..), fill = "#F8766D" , contour = FALSE) +
  stat_density_2d(data = subset(SubjectsDFstack, metric == unique(SubjectsDFstack$metric)[2]), geom = "raster", aes(alpha = ..density..), fill = "#6daaf8" , contour = FALSE) +
  stat_density_2d(data = subset(SubjectsDFstack, metric == unique(SubjectsDFstack$metric)[3]), geom = "raster", aes(alpha = ..density..), fill = "#6df876" , contour = FALSE) +
  scale_alpha(range = c(0, 0.7)) +
  scale_x_continuous(trans='log2') +
  scale_y_continuous(trans='log')


[image]

Kind regards,
Mittyri

Complete thread:

Activity
 Admin contact
21,210 posts in 4,426 threads, 1,481 registered users;
online 3 (0 registered, 3 guests [including 3 identified bots]).
Forum time: Wednesday 17:55 CET (Europe/Vienna)

You can’t fix by analysis
what you bungled by design.    Richard J. Light, Judith D. Singer, John B. Willett

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