R script: Rmodel.R

$$\newcommand{\esp}[1]{\mathbb{E}\left(#1\right)} \newcommand{\var}[1]{\mbox{Var}\left(#1\right)} \newcommand{\deriv}[1]{\dot{#1}(t)} \newcommand{\prob}[1]{ \mathbb{P}\!(#1)} \newcommand{\eqdef}{\mathop{=}\limits^{\mathrm{def}}} \newcommand{\by}{\boldsymbol{y}} \newcommand{\bc}{\boldsymbol{c}} \newcommand{\bpsi}{\boldsymbol{\psi}} \def\pmacro{\texttt{p}} \def\like{{\cal L}} \def\llike{{\cal LL}} \def\logit{{\rm logit}} \def\probit{{\rm probit}} \def\one{{\rm 1\!I}} \def\iid{\mathop{\sim}_{\rm i.i.d.}} \def\simh0{\mathop{\sim}_{H_0}} \def\df{\texttt{df}} \def\res{e} \def\xomega{x} \newcommand{\argmin}[1]{{\rm arg}\min_{#1}} \newcommand{\argmax}[1]{{\rm arg}\max_{#1}} \newcommand{\Rset}{\mbox{\mathbb{R}}} \def\param{\theta} \def\setparam{\Theta} \def\xnew{x_{\rm new}} \def\fnew{f_{\rm new}} \def\ynew{y_{\rm new}} \def\nnew{n_{\rm new}} \def\enew{e_{\rm new}} \def\Xnew{X_{\rm new}} \def\hfnew{\widehat{\fnew}} \def\degree{m} \def\nbeta{d} \newcommand{\limite}[1]{\mathop{\longrightarrow}\limits_{#1}} \def\ka{k{\scriptstyle a}} \def\ska{k{\scriptscriptstyle a}} \def\kel{k{\scriptstyle e}} \def\skel{k{\scriptscriptstyle e}} \def\cl{C{\small l}} \def\Tlag{T\hspace{-0.1em}{\scriptstyle lag}} \def\sTlag{T\hspace{-0.07em}{\scriptscriptstyle lag}} \def\Tk{T\hspace{-0.1em}{\scriptstyle k0}} \def\sTk{T\hspace{-0.07em}{\scriptscriptstyle k0}} \def\thalf{t{\scriptstyle 1/2}} \newcommand{\Dphi}[1]{\partial_\pphi #1} \def\asigma{a} \def\pphi{\psi} \newcommand{\stheta}{{\theta^\star}} \newcommand{\htheta}{{\widehat{\theta}}}$$

# 1 A basic PK model

A PK model is implemented in Rmodel/modelMlxt_pk1.txt:

[LONGITUDINAL]
input = {V, k, a}

EQUATION:
C = pkmodel(V, k)

DEFINITION:
y = {distribution=normal, prediction=C, sd=a}

[INDIVIDUAL]
input = {V_pop, omega_V, k_pop, omega_k, w}

EQUATION:
V_pred = V_pop*(w/70)

DEFINITION:
V = {distribution = lognormal, prediction = V_pred, sd = omega_V}
k = {distribution = lognormal, prediction = k_pop,  sd = omega_k}


This model is implemented as a R function in Rmodel/modelR_pk1.R

Show the R code

Then, any of these two models can be used with Simulx

pk.model <- "Rmodel/modelR_pk1.R"
# pk.model <- "Rmodel/modelMlxt_pk1.txt"

y <- list(name="y", time=seq(18, 80, by=6))
C <- list(name="C", time=seq(0,100, by=0.5))
ind <- list(name=c("V","k"))
p <- c(V_pop=10, omega_V=0.3, w=50, k_pop=0.1, omega_k=0.2, a=2)

res1 <- simulx(model     = pk.model,
output    = list(C,y,ind),
parameter = p,
group     = list(g1,g2))

print(res1$parameter) ## group id V k ## 1 1 1 9.451768 0.11216190 ## 2 1 2 10.871362 0.07699689 ## 3 1 3 9.346003 0.09052976 ## 4 1 4 6.052929 0.12801316 ## 5 1 5 4.980617 0.12138887 ## 6 1 6 7.492797 0.07191773 ## 7 2 7 6.925656 0.09012701 ## 8 2 8 12.061768 0.08847391 ## 9 2 9 13.229024 0.10122134 print(ggplot() + geom_line( data=res1$C, aes(x=time, y=C), colour="black", size=0.5) +

# 3 Comparing computational times between R and Mlxtran

Using R instead of Mlxtran can be extremly slow, even with a model which is not based on ODEs

adm <- list(time=seq(0,66,by=12), amount=100)
C <- list(name="C", time=seq(0,100, by=0.5))
ind <- list(name=c("w","V","k"))
p <- c(w_pop=70, omega_w=10, V_pop=10, omega_V=0.3, k_pop=0.1, omega_k=0.2)
g <- list(size=1000)

ptm <- proc.time()
res3a <- simulx(model     = "Rmodel/modelR_pk3.R",
output    = list(C,ind),
parameter = p,
group     = g)
print(proc.time() - ptm)
##    user  system elapsed
##    8.72    0.28    9.08
print(prctilemlx(res3a$C)+ylim(c(0,25))) ptm <- proc.time() res3b <- simulx(model = "Rmodel/modelMlxt_pk3.txt", output = list(C,ind), parameter = p, treatment = adm, group = g) print(proc.time() - ptm) ## user system elapsed ## 2.18 0.15 3.97 print(prctilemlx(res3b$C)+ylim(c(0,25)))

# 4 A model with regression variables

modelR_reg1 <- function(parameter,regressor,time)
{
with(as.list(parameter),{
t <- time[[1]]
regt <- regressor[[1]]$time regC <- regressor[[1]]$value
C <- approx(regt,regC,t)$y E = Emax*C/(C+EC50) r <- list(C=data.frame(time=t, C=C), E=data.frame(time=t, E=E)) return(r) }) } reg <- data.frame(time=seq(0,50,by=5), Cin=exp(-0.1*seq(0,50,by=5))) out <- list(name=c('C','E'), time=seq(0,50, by=0.2)) res4 <- simulx( model = "modelR_reg1", parameter = c(Emax=100, EC50=0.3), regressor = reg, output = out) plot1 <- ggplot(data=res4$C) + geom_line(aes(x=time, y=C))
plot2 <- ggplot(data=res4\$E) + geom_line(aes(x=time, y=E))
gridExtra::grid.arrange(plot1, plot2, ncol=2)

`