## Simulation framework [NCA / SHAM]

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?

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')

Kind regards,
Mittyri