Simulation framework [NCA / SHAM]

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

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,162 posts in 4,410 threads, 1,476 registered users;
online 2 (0 registered, 2 guests [including 2 identified bots]).
Forum time: Tuesday 04:55 CEST (Europe/Vienna)

In theory, there is no difference between theory and practice.
But, in practice, there is.    Jan L.A. van de Snepscheut

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