R script: analytical.R

Mlxtran code: model/analytical.txt

$$\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 Introduction

Functions of time are implemented in a block EQUATION of the section [LONGITUDINAL] of a script Mlxtran.

We consider parametric functions of time of the form $$f(t;\phi)$$, where $$\phi$$ is a vector of structural parameters.

The [LONGITUDINAL] section starts with the definition of the vector of input parameters $$\phi$$ in a list input.

Then, the simulx script should provide

• the model which is used (i.e.Â the name of the Mlxtran or PharmML script),
• the values of the components of $$\phi$$,
• the list of desired outputs and the times at which they should be evaluated.

Results are returned as a list of data frames.

# 2 Example

We consider in this example a model where two functions $$f_1$$ and $$f_2$$ depending on a vector of parameters $$\phi=(ka,V,k)$$ are defined:

\begin{aligned} f_1(t) &= \frac{D}{V}e^{-k\, t} \\ f_2(t) &= \frac{D\, k_a}{V(k_a - k)}(e^{-k\, t} - e^{-k_a \, t}) \end{aligned} where $$D=100$$ is given.

This model is implemented in the model file model/analytical.txt:

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

EQUATION:
D  = 100
f1 = D/V*exp(-k*t)
f2 = D*ka/(V*(ka-k))*(exp(-k*t) - exp(-ka*t))


We will start evaluating only $$f_1$$ every 0.1h between time 0 and time 25, with $$k_a=0.5$$, $$V=10$$ and $$k=0.2$$:

f <- list(name = 'f1',  time = seq(0, 25, by=0.1))

p <- c(ka = 0.5, V = 10, k = 0.2)

res <- simulx(model     = 'model/analytical.txt',
parameter = p,
output    = f)

res is a list with element f1 which is a data frame with 251 rows and 2 columns:

res$f1[1:5,] ## time f1 ## 1 0.0 10.000000 ## 2 0.1 9.801987 ## 3 0.2 9.607894 ## 4 0.3 9.417645 ## 5 0.4 9.231163 res$f1[248:251,]
##     time         f1
## 248 24.7 0.07154598
## 249 24.8 0.07012928
## 250 24.9 0.06874063
## 251 25.0 0.06737947

we can plot this data using ggplot

print(ggplot(data=res$f1) + geom_line(aes(x=time, y=f1)))  We way now want evaluate both $$f_1$$ and $$f_2$$ at the same time points. We can still define the output in a single list, but the field name has now two elements f1 and f2: f <- list(name = c('f1','f2'), time = seq(0, 25, by=0.1)) res <- simulx(model = 'model/analytical.txt', parameter = p, output = f) print(ggplot() + geom_line(data=res$f1, aes(x=time, y=f1), color="black") +
geom_line(data=res$f2, aes(x=time, y=f2), color="red") + ylab('concentration')) It can be useful to reshape the data before plotting it with ggplot. We merge the two outputs into a unique data frame r with three columns: time, f (which takes the values f1 or f2) and value. library(reshape2) r <- merge(res$f1,res$f2) r <- melt(r , id = 'time', variable.name = 'f', value.name="concentration") r[c(1:4),] ## time f concentration ## 1 0.0 f1 10.000000 ## 2 0.1 f1 9.801987 ## 3 0.2 f1 9.607894 ## 4 0.3 f1 9.417645 r[c(250:253),] ## time f concentration ## 250 24.9 f1 0.06874063 ## 251 25.0 f1 0.06737947 ## 252 0.0 f2 0.00000000 ## 253 0.1 f2 0.48282081 This format is now suitable for ggplot: print(ggplot(r, aes(time,concentration)) + geom_line(aes(colour = f),size=1) + guides(colour=guide_legend(title=NULL)) + theme(legend.position=c(.9, .8))) It is also possible to evaluate $$f_1$$ and $$f_2$$ at different time points. The two outputs should be then defined in two separate lists: f1 <- list(name = 'f1', time = seq(0, 15, by=0.1)) f2 <- list(name = 'f2', time = seq(10, 25, by=1)) res <- simulx(model = 'model/analytical.txt', parameter = p, output = list(f1,f2)) plot(ggplot() + geom_line(data=res$f1, aes(x=time, y=f1), color="black") +
geom_line(data=res$f2, aes(x=time, y=f2), color="red") + ylab('concentration')) We may again prefer to reshape the data before plotting it with ggplot r <- merge(res$f1,res$f2,all=TRUE) r <- melt(r , id = 'time', variable.name = 'f', value.name="concentration") r=r[(!is.na(r$concentration)),]
print(ggplot(r, aes(time,concentration)) + geom_line(aes(colour = f),size=1) +
guides(colour=guide_legend(title=NULL)) + theme(legend.position=c(.9, .8)))