We show how to simulate PK data using the theophylline project, where the parameters have been estimated with Monolix.

Remark: in this project, the amount is given in mg per kilo.

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

Without any other inputs simulx will simulate the exact study design underlying the input data for parameter estimation

res1  <- simulx(project = project.file)
names(res1)
## [1] "y1"         "treatment"  "originalId" "population"
print(res1$treatment)
##    id time amount
## 1   1    0    320
## 2   2    0    320
## 3   3    0    320
## 4   4    0    320
## 5   5    0    320
## 6   6    0    320
## 7   7    0    320
## 8   8    0    320
## 9   9    0    320
## 10 10    0    320
## 11 11    0    320
## 12 12    0    320
print(res1$population)
##       ka_pop        V_pop  beta_V_lw70       Cl_pop beta_Cl_lw70 
##      1.51895     33.11292      0.74870      2.80800      0.51409 
##     omega_ka      omega_V     omega_Cl            b 
##      0.64446      0.09442      0.23628      0.15909

plot the results ‘’as usual’’

print(ggplot(data=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"))


Some parameters values can be modified, e.g. switch off the residual error by setting b = 0, and some additional outputs can be defined, e.g. the predicted concentration \(Cc\), the individual PK parameters and the covariates, .

sim.param <- c(b=0.0)
out1  <- list(name = 'Cc', time = seq(0, 25, by=0.1))
outp <- c("ka","V", "Cl","WEIGHT")
res2  <- simulx(project   = project.file,
                output    = list(out1,outp),
                parameter = sim.param)
names(res2)
## [1] "y1"         "Cc"         "parameter"  "treatment"  "originalId"
## [6] "population"
print(res2$population)
##       ka_pop        V_pop  beta_V_lw70       Cl_pop beta_Cl_lw70 
##      1.51895     33.11292      0.74870      2.80800      0.51409 
##     omega_ka      omega_V     omega_Cl            b 
##      0.64446      0.09442      0.23628      0.00000
head(res2$parameter)
##   id        ka        V       Cl WEIGHT
## 1  1 2.5026398 37.32696 2.002180   79.6
## 2  2 1.8167616 38.17348 1.899981   72.4
## 3  3 1.0505773 31.00704 2.204756   70.5
## 4  4 1.8912941 30.95903 2.097024   72.7
## 5  5 0.4238495 28.06020 2.357542   54.6
## 6  6 0.5895110 38.97677 3.386650   80.0
print(ggplot() + 
        geom_point(data=res2$y1, aes(x=time, y=y1, colour=id)) +
        geom_line(data=res2$Cc, aes(x=time, y=Cc, colour=id)) +
        scale_x_continuous("Time") + scale_y_continuous("Concentration"))

Observation times can also be redefined:

out2  <- list(name = 'y1', time = seq(1, 25, by=2))
res3  <- simulx(project   = project.file,
                output    = list(out1, out2),
                parameter = sim.param)

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

Remove the simulated concentrations from the outputs,

out0 <- list(name='y1', time='none')
res4  <- simulx(project   = project.file,
                output    = list(out1,out0))
names(res4)
## [1] "Cc"         "treatment"  "originalId" "population"

Use the estimated individual parameters (mode = EBEs),

sim.param <- "mode"
res5  <- simulx(project   = project.file,
                output    = out1,
                parameter = sim.param)

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

Use the estimated individual parameters and set b=0

sim.param <- list("mode",c(b=0))
res6  <- simulx(project   = project.file,
                output    = out1,
                parameter = sim.param)

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

Remove treatement from the outputs,

res7  <- simulx(project   = project.file,
                settings  = list(out.trt=F))
names(res7)
## [1] "y1"         "originalId" "population"


A new administration schedule can be specified

adm   <- list(time = c(0,6), amount = c(6, 3))

res8 <- simulx(project = project.file, treatment = adm, output = list(out1, out2))

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