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"

adm1 <- list(time=seq(0,66,by=12), amount=100)
adm2 <- list(time=seq(12,78,by=24), amount=200)
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)
g1 <- list(treatment=adm1, size=6)
g2 <- list(treatment=adm2, size=3)

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) +
        geom_point(data=res1$y, aes(x=time, y=y), colour="red", size=1)+ 
        facet_wrap( ~ id))

2 An ODE based model

The ODE based model Rmodel/modelMlxt_pk2.txt


[LONGITUDINAL]
input = {ka, V, Cl, k12, k21, k13, k31}

EQUATION:
k = Cl/V
ddt_Ad = -ka*Ad
ddt_Ac = ka*Ad - k*Ac + k21*Ap - k12*Ac + k31*Aq - k13*Ac
ddt_Ap = -k21*Ap + k12*Ac 
ddt_Aq = -k31*Aq + k13*Ac
C = Ac/V

can also be implemented in R:

library(deSolve)

modelR_pk2 <- function(parameter,dose,time)
{  
  MMmod <- function(Time, State, Pars) {
    with(as.list(c(State, Pars)), {
      k <- Cl/V
      ddt_Ad = -ka*Ad
      ddt_Ac = ka*Ad - k*Ac + k21*Ap - k12*Ac + k31*Aq - k13*Ac
      ddt_Ap = -k21*Ap + k12*Ac 
      ddt_Aq = -k31*Aq + k13*Ac
      return(list(c(ddt_Ad, ddt_Ac, ddt_Ap, ddt_Aq)))
    })
  }
  
  yini  <- c(Ad=0, Ac=0, Ap=0, Aq=0)
  
  t <- time[[1]]
  t <- sort(c(t,dose$data$time))
  out   <- lsode(yini, t, MMmod, parameter, events=dose)
  i0 <- which(diff(t)==0)
  out <- out[-i0,]
  
  C <- data.frame(time=out[,"time"], C=out[,"Ac"]/parameter["V"])
  r <- list(C=C)
  return(r)
}
pk.model <- "Rmodel/modelR_pk2.R"
# pk.model <- "Rmodel/modelMlxt_pk2.txt"

adm1 <- list(time=seq(1,144,by=24), amount=10, target="Ac")
adm2 <- list(time=seq(7,144,by=12), amount=4,  target="Ad")
adm <- list(adm1, adm2)
C <- list(name="C", time=seq(0,200, by=0.2))

p = c(ka=1, V=0.1, Cl=0.01, k12=0.2, k21=0.2, k13=0.3, k31=0.3)

res2 <- simulx(model     = pk.model, 
               output    = C,
               treatment = adm,
               parameter = p)

print(ggplot(data=res2$C)  + geom_line(aes(x=time, y=C), 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,
                treatment = adm,
                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)

```