inst/examples/ex.9.PK.2.comp.oral.md.ode.compiled.R

In PopED: Population (and Individual) Optimal Experimental Design

#' An implementation of a two compartment model with oral absorption using ODEs.
library(PopED)
library(deSolve)

#' Define the ODE system
PK.2.comp.oral.ode <- function(Time, State, Pars){
  with(as.list(c(State, Pars)), {    
    dA1 <- -KA*A1 
    dA2 <- KA*A1 + A3* Q/V2 -A2*(CL/V1+Q/V1)
    dA3 <- A2* Q/V1-A3* Q/V2
    return(list(c(dA1, dA2, dA3)))
  })
}

#' define the initial conditions and the dosing
ff.PK.2.comp.oral.md.ode <- function(model_switch, xt, parameters, poped.db){
  with(as.list(parameters),{
    A_ini <- c(A1=0, A2=0, A3=0)
    times_xt <- drop(xt)
    dose_times = seq(from=0,to=max(times_xt),by=TAU)
    eventdat <- data.frame(var = c("A1"), 
                           time = dose_times,
                           value = c(DOSE), method = c("add"))
    times <- sort(c(times_xt,dose_times))
    out <- ode(A_ini, times, PK.2.comp.oral.ode, parameters, events = list(data = eventdat))#atol=1e-13,rtol=1e-13)
    y = out[, "A2"]/(V1/Favail)
    y=y[match(times_xt,out[,"time"])]
    y=cbind(y)
    return(list(y=y,poped.db=poped.db))
  })
}

#' parameter definition function 
#' names match parameters in function ff
fg <- function(x,a,bpop,b,bocc){
  parameters=c( CL=bpop[1]*exp(b[1]),
                V1=bpop[2],
                KA=bpop[3]*exp(b[2]),
                Q=bpop[4],
                V2=bpop[5],
                Favail=bpop[6],
                DOSE=a[1],
                TAU=a[2])
  return( parameters ) 
}


#' create poped database
poped.db <- create.poped.database(ff_fun="ff.PK.2.comp.oral.md.ode",
                                  fError_fun="feps.add.prop",
                                  fg_fun="fg",
                                  groupsize=20,
                                  m=1,      #number of groups
                                  sigma=c(prop=0.1^2,add=0.05^2),
                                  bpop=c(CL=10,V1=100,KA=1,Q= 3.0, V2= 40.0, Favail=1),
                                  d=c(CL=0.15^2,KA=0.25^2),
                                  notfixed_bpop=c(1,1,1,1,1,0),
                                  xt=c( 48,50,55,65,70,85,90,120), 
                                  minxt=0,
                                  maxxt=144,
                                  discrete_xt = list(0:144),
                                  a=c(DOSE=100,TAU=24),
                                  maxa=c(DOSE=1000,TAU=24),
                                  mina=c(DOSE=0,TAU=8),
                                  discrete_a = list(DOSE=seq(0,1000,by=100),TAU=8:24))


#' plot intial design just PRED
plot_model_prediction(poped.db,model_num_points = 500)

#' plot intial design with BSV and RUV in model
plot_model_prediction(poped.db,IPRED=T,DV=T)

#' how long does one evaluation of the FIM take? 
tic(); (eval <- evaluate_design(poped.db)); toc()


#' To make optimization more reasonable you can use compiled code
#' 
#' compile and load the two_comp_oral_CL.c code.
#' To set this up see the 
#' "R Package deSolve, Writing Code in Compiled Languages" 
#' vingette in the deSolve documentation
#' 
#' make sure your working directory contains the c code
system("R CMD SHLIB two_comp_oral_CL.c")
dyn.load(paste("two_comp_oral_CL", .Platform$dynlib.ext, sep = ""))

#' define the initial conditions and the dosing
ff.PK.2.comp.oral.md.ode.compiled <- function(model_switch, xt, parameters, poped.db){
  with(as.list(parameters),{
    A_ini <- c(A1=0, A2=0, A3=0)
    times_xt <- drop(xt)
    dose_times = seq(from=0,to=max(times_xt),by=TAU)
    eventdat <- data.frame(var = c("A1"), 
                           time = dose_times,
                           value = c(DOSE), method = c("add"))
    times <- sort(c(times_xt,dose_times))
    out <- ode(A_ini, times, func = "derivs", parms = parameters,
               #jacfunc = "jac", # not really needed, speed up is minimal if this is defined or not.
               dllname = "two_comp_oral_CL",
               initfunc = "initmod", nout = 1, outnames = "Sum",
               events = list(data = eventdat))#atol=1e-13,rtol=1e-13))  
    y = out[, "A2"]/(V1/Favail)
    y=y[match(times_xt,out[,"time"])]
    y=cbind(y)
    return(list(y=y,poped.db=poped.db))
  })
}

#' create poped database
poped.db.compiled <- create.poped.database(poped.db, ff_fun="ff.PK.2.comp.oral.md.ode.compiled")

##  create plot of model without variability 
plot_model_prediction(poped.db.compiled,model_num_points = 500)

##  create plot of model with variability 
plot_model_prediction(poped.db.compiled,IPRED=T,DV=T,model_num_points = 500)

#' how long does one evaluation of the FIM take? 
tic(); (eval_compiled <- evaluate_design(poped.db.compiled)); toc()

#' very small differences in computation value
#' but a large difference in computation time (8 times faster with the compiled code)
(eval_compiled$ofv-eval$ofv)/eval$ofv

#' making optimization times more resonable
#output <- poped_optim(poped.db.compiled,opt_xt=T, opt_a=T, parallel=T)