The R codes and the Mlxtran model files can be found here

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

Introduction

simulx can be easily integrated in a workflow:

  1. modeling of the theophylline PK data is performed first with Monolix,

  2. simulation of PK data is then performed with simulx, using the parameters estimated by Monolix.


In this example we show how to simulate a model where the parameters have been estimated with Monolix. The Modelling and Simulation (M&S) workflow allows us to directly simulate the model just by providing the Monolix project as an input to simulx. Without any other inputs simulx will simulate the exact study design underlying the input data for parameter estimation. simulx gives us the flexibility to change any of the study design variables. We will demonstrate this by first simulating the original study design and then by changing observations, dosing regimens, number of patients and the patient covariates. Finally, we will show how to simulate a dose-concentration curve and how to do Monte Carlo simulations to assess how the number of patients affects the SE error of the concentration levels.


Examples

Define the project to be used for the simulations (relative paths)

project.file <- 'monolixRuns/theophylline_project.mlxtran'

1. Simulate a trial using the original design

Simulate model with individual parameters sampled from the estimated population distribution

  • Covariate WEIGHT is taken from data set
  • Dosage regimen and observation times are taken from data set
sim.res1  <- simulx(project = project.file)

plot the results ‘’as usual’’

print(ggplot(data=sim.res1$y1) + 
        geom_point(aes(x=time, y=y1, colour=id)) +
        geom_line(aes(x=time, y=y1, colour=id)) +
        scale_x_continuous("Time") + scale_y_continuous("Concentration"))


2. Remove the residual error

  • Covariate WEIGHT is taken from data set
  • Dosage regimen and observation times are taken from data set
  • and switch off the residual error by setting b = 0
sim.param <- c(b=0)
out  <- list(name = 'Cc', time = seq(0, 25, by=0.1))
sim.res2  <- simulx(project   = project.file,
                    output    = out,
                    parameter = sim.param)

print(ggplot() + 
        geom_point(data=sim.res2$y1, aes(x=time, y=y1, colour=id)) +
        geom_line(data=sim.res2$Cc, aes(x=time, y=Cc, colour=id)) +
        scale_x_continuous("Time") + scale_y_continuous("Concentration"))


3. Simulate a trial with N individuals

  • The designs and the weights of these N patients are sampled from the original dataset
N <- 50
sim.res3  <- simulx(project = project.file,
                    group   = list(size = N))

print(ggplot(data=sim.res3$y1) + 
        geom_point(aes(x=time, y=y1, colour=id)) +
        geom_line(aes(x=time, y=y1, colour=id)) +
        scale_x_continuous("Time") + scale_y_continuous("Concentration") +
        theme(legend.position="none"))


4. Define new observation times

  • Covariate WEIGHT and dosage regimen are taken from data set
out  <- list(name = 'y1', time = (0:12))

sim.res4  <- simulx(project = project.file,
                    output = out)

print(ggplot(data=sim.res4$y1) + 
        geom_point(aes(x=time, y=y1, colour=id)) +
        geom_line(aes(x=time, y=y1, colour=id)) +
        scale_x_continuous("Time") + scale_y_continuous("Concentration"))


5. Output the predicted concentrations, the individual parameters and the covariates

  • Covariate WEIGHT and dosage regimen are taken from data set
out1  <- list(name = 'y1', time = seq(0, 24, by=2))
out2  <- list(name = 'Cc', time = seq(0, 24, by=0.1))
out3  <- list(name = c('V', 'Cl', 'WEIGHT'))

sim.res5  <- simulx(project = project.file,
                    output  = list(out1, out2, out3))

print(sim.res5$parameter)
##    id        V       Cl WEIGHT
## 1   1 32.99802 4.898956   79.6
## 2   2 32.63807 3.875025   72.4
## 3   3 28.69185 2.621612   70.5
## 4   4 35.09865 3.101235   72.7
## 5   5 27.93620 1.629150   54.6
## 6   6 36.72302 2.674799   80.0
## 7   7 30.95468 5.271731   64.6
## 8   8 34.88231 3.673450   70.5
## 9   9 36.42052 3.723640   86.4
## 10 10 28.39625 2.811044   58.2
## 11 11 30.57689 2.425945   65.0
## 12 12 25.83735 2.642265   60.5
print(ggplot() + 
        geom_point(data=sim.res5$y1,aes(x=time, y=y1, colour=id)) +
        geom_line(data=sim.res5$Cc,aes(x=time, y=Cc, colour=id)) +
        scale_x_continuous("Time") + scale_y_continuous("Concentration"))