Simulation framework [NCA / SHAM]

posted by mittyri – Russia, 2019-04-16 23:21 (592 d 06:00 ago) – Posting: # 20179
Views: 4,402

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,216 posts in 4,427 threads, 1,481 registered users;
online 7 (0 registered, 7 guests [including 7 identified bots]).
Forum time: Sunday 05:21 UTC (Europe/Vienna)

Half the harm that is done in this world
Is due to people who want to feel important.    T. S. Eliot

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