Bioequivalence and Bioavailability Forum • Simulation framework

Simulation framework [NCA / SHAM]

posted by mittyri – Russia, 2019-04-17 01:21 – Posting: # 20179
Views: 2,462

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

—
Kind regards,
Mittyri

Complete thread:

Metrics for absorption Helmut 2019-04-15 12:08
- Metrics for absorption nobody 2019-04-15 12:15
  - Metrics for absorption ElMaestro 2019-04-15 12:23
    - Metrics for absorption nobody 2019-04-15 12:39
  - Metrics for absorption Helmut 2019-04-15 13:02
    - Metrics for absorption nobody 2019-04-15 13:48
      - Metrics for absorption Helmut 2019-04-15 19:31
        
        Metrics for absorption nobody 2019-04-16 09:34
        
        Metrics for absorption nobody 2019-04-16 15:15
        
        Metrics for absorption Helmut 2019-04-16 17:01
        
        Metrics for absorption nobody 2019-04-17 11:15
        
        Simulation frameworkmittyri 2019-04-17 01:21
        
        Simulation framework Helmut 2019-04-17 13:31
        
        Conclusion mittyri 2019-04-17 15:42
        
        Conclusion & beyond Helmut 2019-04-17 17:05