R script: analytical.R

Mlxtran code: model/analytical.txt


# 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)))