R script: analytical.R
Mlxtran code: model/analytical.txt
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
Mlxtran
or PharmML script),Results are returned as a list of data frames.
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)))