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

Task

Individualised therapies are increasingly important for drug development. One reason for individualised therapies is concentration related toxicity and high between-patient variability of exposure levels. In such a case each patient can be monitored for the concentration levels and the drug dose will be adjusted when the individual concentration level exceeds a certain threshold. In this example we demonstrate how Simulx is used to run clinical trial simulations with an adaptive dosing. We will simulate a PK/PD model with a decision rule to adjust the dose based on monitored concentration levels. In addition, we will include a decision rule for the efficacy as well to ensure that desired efficacy level will be maintained.


The PK model

Tne PK model is a one compartment model for oral administration with first order absorption and elimination.

Parameters \(ka\), \(V\) and \(k\) are log-normally distributed

This PK model is implemented in a text file cts1a.txt:

[LONGITUDINAL]
input = {ka, V, k, a1}
EQUATION:
Cc = pkmodel(ka, V, k)
DEFINITION:
y1 = {distribution=lognormal, prediction=Cc, sd=a1}

[INDIVIDUAL]
input={ka_pop, omega_ka, V_pop, omega_V, k_pop, omega_k}
DEFINITION:
ka = {distribution=lognormal, prediction=ka_pop, sd=omega_ka}
V  = {distribution=lognormal, prediction=V_pop,  sd=omega_V}
k  = {distribution=lognormal, prediction=k_pop,  sd=omega_k}


Simulation of PK data

We can use simulx for simulating a clinical trial with 20 patients to whom a dose of 20mg is administrated orally every 12 hours starting at time 0.

We then want to compute for these 20 patients:


adm <- list(time=seq(0,to=440,by=12), amount=20)

ppk <- c(ka_pop=0.4,   V_pop=10,    k_pop=0.05,
         omega_ka=0.3, omega_V=0.5, omega_k=0.1,
         a1=0.05)

g   <- list(size=20, level='individual')

Cc <- list(name='Cc', time=seq(0,to=440,by=1))
y1  <- list(name='y1', time=seq(4,to=440,by=24))

s   <- list(seed=1234)

res1 <- simulx(model     = "model/cts1a.txt",
               parameter = ppk,
               treatment = adm,
               group     = g,
               output    = list(Cc, y1),
               settings  = s)

(we use seed = 1234 for reproducibility reason. Different results will be obtained with different values of the seed.)

We can now plot the predicted and observed concentrations:

library("gridExtra")
library(ggplot2)
theme_set(theme_bw())
plot1 <- ggplot(data=res1$Cc, aes(x=time, y=Cc, colour=id)) +  
  geom_line(size=0.5) + xlab("time (hour)") + ylab("Cc (mg/l)") +
  theme(legend.position="none")
plot2 <- ggplot(data=res1$y1, aes(x=time, y=y1, colour=id)) +  
  geom_line(size=0.5) + geom_point(size=2) +
  xlab("time (hour)") + ylab("measured conc. (mg/l)") + 
  theme(legend.position="none")
grid.arrange(plot1, plot2, ncol=2)


Introducing a single decision rule

Assume now that, for safety reason, the concentration should remain below a given threshold of 5 mg/l.

The exposure, i.e. the ``true’’ concentration \(Cc\), is above this threshold for several patients. Since \(Cc\) cannot be measured in practice, we can then decide to reduce the dose for those patients which measured concentration is above 5 mg/l.

We thus define an adaptive individual design as follows:

This decision rule is implemented as a R list:

r1 <- list(name='y1', 
           time=seq(124,to=440,by=24), 
           condition="y1>5", 
           factor=0.75)

R function titration is a variant of simulx which now takes this rule into account.

Remark: function titration.R is not part of the mlxR package. This is just an example of what can be done: this function can be freely used/modified by anyone.

source("titration.R")
res2 <- titration(model     = "model/cts1a.txt",
                  parameter = ppk,
                  treatment = adm,
                  output    = list(Cc, y1),
                  rule      = r1,
                  group     = g,
                  settings  = s)

Outputs of titration and simulx have similar structure. We can therefore display the predicted and observed concentrations obtained with titration as we did previously with the output of simulx:

plot3 <- ggplot(data=res2$Cc, aes(x=time, y=Cc, colour=id))     +  
  geom_line(size=0.5) + xlab("time (hour)") + ylab("Cc (mg/l)") +
  theme(legend.position="none")  + geom_hline(yintercept=5)
plot4 <- ggplot(data=res2$y1, aes(x=time, y=y1, colour=id))     +  
  geom_line(size=0.5) + geom_point(size=2) + 
  xlab("time (hour)") + ylab("measured conc. (mg/l)")           + 
  theme(legend.position="none") + geom_hline(yintercept=5)
grid.arrange(plot3, plot4, ncol=2)

What was expected is now clearly observed: high concentrations tend to reduce after the first visit at time 124.


Combining several decision rules

We can introduce several decision rules. Imagine for instance that efficacy requires a concentration above 3 mg/l. We can then decide, for example, to multiply the dose by a factor 1.5 when the observed concentration is below this value at a given visit.

We can then implement this new decision rule as a new list:

r2 <- list(name='y1', 
           time=seq(124,to=440,by=24),
           condition="y1<3", 
           factor=1.5)

We will now use titration combining the two decisions rules r1 and r2:

res3 <- titration(model     = "model/cts1a.txt",
                  parameter = ppk,
                  treatment = adm,
                  output    = list(Cc, y1),
                  rule      = list(r1,r2),
                  group     = g,
                  settings  = s)

plot5 <- ggplot(data=res3$Cc, aes(x=time, y=Cc, colour=id))     +  
  geom_line(size=0.5) + xlab("time (hour)") + ylab("Cc (mg/l)") +
  theme(legend.position="none") + geom_hline(yintercept=c(3,5))
plot6 <- ggplot(data=res3$y1, aes(x=time, y=y1, colour=id))     +  
  geom_line(size=0.5) + geom_point(size=2) +
  xlab("time (hour)") + ylab("measured conc. (mg/l)")   + 
  theme(legend.position="none") + geom_hline(yintercept=c(3,5))
grid.arrange(plot5, plot6, ncol=2)


Extension to PKPD data

Assume now that we have PK and PD data. In this example, the PD model is a Emax model.

The joint PKPD model is implemented in a new file cts1b.txt:

[LONGITUDINAL]
input = {ka, V, k, Emax, EC50, a1, a2}
EQUATION:
Cc = pkmodel(ka, V, k)
E  = Emax*Cc/(EC50+Cc)

DEFINITION:
y1 = {distribution=lognormal, prediction=Cc, sd=a1}
y2 = {distribution=normal, prediction=E,  sd=a2}

[INDIVIDUAL]
input={
ka_pop, omega_ka, V_pop, omega_V, k_pop, omega_k,
Emax_pop, omega_Emax, EC50_pop, omega_EC50
}
DEFINITION:
ka   = {distribution=lognormal, prediction=ka_pop,   sd=omega_ka}
V    = {distribution=lognormal, prediction=V_pop,    sd=omega_V}
k    = {distribution=lognormal, prediction=k_pop,    sd=omega_k}
Emax = {distribution=lognormal, prediction=Emax_pop, sd=omega_Emax}
EC50 = {distribution=lognormal, prediction=EC50_pop, sd=omega_EC50}

We need to provide the values of the population PD parameters for simulating PKPD data from this model.

ppd <- c(Emax_pop=100,   EC50_pop=3, 
         omega_Emax=0.1, omega_EC50=0.2, a2=5)

E <-  list(name='E',  time=seq(0,to=440,by=1))
y2 <- list(name='y2', time=seq(16,to=440,by=36))

res4 <- simulx(model     = "model/cts1b.txt",
               parameter = c(ppk,ppd),
               treatment = adm,
               output    = list(Cc, E, y1, y2),
               group     = g,
               settings  = s)

We can now plot the predicted and observed PK and PD data for the 20 patients:

pl1=ggplot(data=res4$Cc, aes(x=time, y=Cc, colour=id))+geom_line(size=0.5)+
  xlab("time (hour)") + ylab("Cc (mg/l)") + theme(legend.position="none")
pl2=ggplot(data=res4$E, aes(x=time, y=E, colour=id))+geom_line(size=0.5)  +
  xlab("time (hour)") + ylab("E") + theme(legend.position="none")
pl3=ggplot(data=res4$y1, aes(x=time, y=y1, colour=id))+geom_line(size=0.5)+ 
  geom_point(size=2) + xlab("time (hour)") + ylab("measured conc (mg/l)") + 
  theme(legend.position="none")
pl4=ggplot(data=res4$y2, aes(x=time, y=y2, colour=id))+geom_line(size=0.5)+ 
  geom_point(size=2) + xlab("time (hour)") + ylab("measured effect")      + 
  theme(legend.position="none")

grid.arrange(pl1, pl2, pl3, pl4, ncol=2)