
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:

• the predicted concentration $$Cc$$ every hour between times 0h and 440h
• the observed concentration $$y_1$$ every 24 hours between time 4h and 440h

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,
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:

• Visits are scheduled every 24 hours, starting at time 124h.
• If, at a given visit, the observed concentration $$y_1$$ of a patient is greater than 5, then the next doses administrated to this patient will be reduced, by multiplying the current amount by a factor 0.75.

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,
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,
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),
settings  = s)
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)  +
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)+
grid.arrange(pl1, pl2, pl3, pl4, ncol=2)