#' 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.