Simulation framework [NCA / SHAM]

posted by mittyri – Russia, 2019-04-16 23:21 (680 d 05:42 ago) – Posting: # 20179
Views: 4,726

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,355 posts in 4,458 threads, 1,493 registered users;
online 20 (0 registered, 20 guests [including 4 identified bots]).
Forum time: Thursday 05:03 UTC (Europe/Vienna)

The rise of biometry in this 20th century,
like that of geometry in the 3rd century before Christ,
seems to mark out one of the great ages or critical periods
in the advance of the human understanding.    R.A. Fisher

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