Predictive PK-PD modelling of tumor growth after administration of antiangiogenic and cytotoxic drugs, as single-agents and in combination therapy.
Material extracted from the DDMoRe face-to-face course “Model-informed drug development in Oncology (Advanced)”
Authors of the original exercise
Authors of the online training material
This HTML report was generated with the August 2016 Public Release version of the Interoperability Framework.
Tasks
Clear workspace
rm(list=ls(all=FALSE))
Set working directory
setwd(file.path(.MDLIDE_WORKSPACE_PATH,"OnlineTraining_OncologyAdvanced"))
We will need three helper functions:
getDesignParms() retrieves the original design parameters from PopED outputupdatedDosing2simulx() takes PopED output and a Design Object to create simulx() inputparobj2simulx() creates simulx() input from a Parameter Object.These functions are further documented in helper_functions.R
source("helper_functions.R")
How many model parameters are to be estimated?
The model has 7 parameters to be estimated:
LAMBDA_0, LAMBDA_1, K1, K2, IC50, W0, CV
Among these, only K2 and IC50 are linked to drug effect. K2 is the killing rate constant of NMS-P937, and IC50 is the Bevacizumab concentration that reduces natural growth rate to 50%.
Extract dosing information using functions of the ddmore R package
Set model name for future tasks
filename_noint <- "Rocchetti_combo_no_interaction"
mdlfile_noint <- paste0(filename_noint, ".mdl")
Use ddmore function getDataObjects() to retrieve data object(s) from an existing .mdl file
myDataObj <- getDataObjects(mdlfile_noint)[[1]]
myData <- readDataObj(myDataObj)
Look at the first 6 lines of the data set
head(myData)
## ID TIME DV AMT MDV CMT
## 1 1 0 NA NA 1 1
## 2 1 6 0.1628989 NA 0 1
## 3 1 9 0.2370097 NA 0 1
## 4 1 13 0.4822279 NA 0 1
## 5 1 16 0.8020363 NA 0 1
## 6 1 19 1.0697077 NA 0 1
Extract information on the dosing schedule
unique(myData$ID)
## [1] 1 2 3
myData[!is.na(myData$AMT), ]
## ID TIME DV AMT MDV CMT
## 13 2 8 NA 0.1342282 1 1
## 15 2 12 NA 0.1342282 1 1
## 17 2 16 NA 0.1342282 1 1
## 20 2 20 NA 0.1342282 1 1
## 28 3 9 NA 112.6718245 1 3
## 30 3 10 NA 112.6718245 1 3
## 31 3 11 NA 112.6718245 1 3
## 32 3 13 NA 112.6718245 1 3
## 34 3 14 NA 112.6718245 1 3
## 35 3 15 NA 112.6718245 1 3
## 37 3 17 NA 112.6718245 1 3
## 38 3 18 NA 112.6718245 1 3
## 39 3 19 NA 112.6718245 1 3
ID 1 is control, ID 2 receives 4 doses of Bevacizumab at days 8,12,16,20 and ID 3 receives 9 doses of NMS-P937 at days 9-11,13-15,17-39.
Which of the compounds has a stronger effect?
Plot the data using ggplot (graphs are exported to PDF) and compare the drug effect visually.
Extract only observation records
myEDAData<-myData[myData$MDV==0,]
myEDAData$ID <- as.factor(myEDAData$ID)
levels(myEDAData$ID)<- c("Control", "Only Bevacizumab", "Only NMS-P937")
print(ggplot(data=myEDAData, aes(TIME,DV,col=as.factor(ID))) +
geom_line() +
geom_point() +
xlab("Time (days)") +
ylab("Tumor weight (g)") +
labs(col="group"))
At the given dose level, the antiangiogenic compound Bevacizumab delays tumour growth by approximately 9 days, and the cytotoxic compound NMS-P937 delays growth by approximately 17 days with respect to control.
Tasks
rocchetti_interaction_data.csv.rocchetti_interaction_data.csv).Generate the PharmML code from MDL using the as.PharmML() function
Take a look at the generated code (not a mandatory task). Note that the PharmML code should be directly available for download from the DDMoRe repository, so that this step, in general, is not necessary for models published in the model repository.
model_noint.pharmML <- as.PharmML(mdlfile_noint)
Simulation with simulx()
Store the model parameters in a vector
parObj <- getParameterObjects(mdlfile_noint)[[1]]
parValues <- parobj2simulx(parObj)
Encode the dosing schedule from 'rocchetti_interaction_data.csv'
adm1 <- list(target="Q0_A", time=c(8, 12, 16, 20), amount=0.1340511)
adm2 <- list(target="Q0_B", time=c(9, 10, 11, 13, 14, 15, 17, 18, 19), amount=86.042)
f <- list(name='WTOT', time= seq(0, 70, by=0.1))
y <- list(name='Y', time= seq(0, 70, by=3))
pkA <-list(name='C1_A', time= seq(8,40, by=0.05))
pkB <-list(name='C1_B', time= seq(9,40, by=0.05))
g1 <- list(size=1, level='longitudinal')
g2 <- list(size=1, level='longitudinal', treatment=list(adm1,adm2))
Simulate predictions
resSim_noint <- simulx(model = basename(model_noint.pharmML),
parameter = parValues,
group = list(g1,g2),
output = list(f,pkA,pkB))
Compare the simulated curves with the data from the combination experiment
Read in the interaction data
data_rocchetti_2013 <- read.csv("rocchetti_interaction_data.csv", na.strings = ".")
Reformat the simulx() predictions
mean_noint <- resSim_noint$WTOT$WTOT
sd_noint <- 0.13*resSim_noint$WTOT$WTOT
resSim_noint$WTOT$LL <- qnorm(p = 0.025, mean_noint, sd_noint)
resSim_noint$WTOT$UL <- qnorm(p = 0.975, mean_noint, sd_noint)
levels(resSim_noint$WTOT$group)<-c("Control","Treated")
Plot predictions and data
print(ggplot(data=resSim_noint$WTOT) + geom_line(aes(x = time, y=WTOT, col=group)) +
geom_ribbon(data=resSim_noint$WTOT, aes(x=time,ymin=LL, ymax=UL, fill=group), alpha=0.5) +
geom_point(data=data_rocchetti_2013[data_rocchetti_2013$MDV==0,], aes(x=TIME, y=DV)) +
xlab("Time (days)") +
ylab("Tumor weight (g)") +
ggtitle("Rocchetti TGI no interaction"))
The combined effect of Bevacizumab and NMS-P937 is lower than expected from the model that assumes no interaction. Therefore, there is an antagonistic effect of the combination. We will need to consider an interaction model to capture this behaviour.
Tasks
The interaction model is encoded in the following mdl file:
filename_int <- "Rocchetti_combo_interaction"
mdlfile_int <- paste0(filename_int, ".mdl")
FOCE Estimation with NONMEM
mog_FOCE_int <- createMogObj(dataObj = getDataObjects(mdlfile_int)[[1]],
parObj = getParameterObjects(mdlfile_int)[[1]],
mdlObj = getModelObjects(mdlfile_int)[[1]],
taskObj = getTaskPropertiesObjects(mdlfile_int)[[1]]) #foce
mdlfile_int.FOCE <- paste0(filename_int,"_FOCE.mdl")
writeMogObj(mog_FOCE_int,mdlfile_int.FOCE)
nm_int <- estimate(mdlfile_int.FOCE, target="NONMEM", subfolder="NONMEM")
## -- Thu Aug 18 10:45:28 2016
## New
## Submitted
## Job 9deab2d3-d842-42c0-9b9a-1d50157cac42 progress:
## Running [ ..... ]
## Importing Results
## Copying the result data back to the local machine for job ID 9deab2d3-d842-42c0-9b9a-1d50157cac42...
## From C:\Users\Niklas\AppData\Local\Temp\RtmpglYELE\DDMORE.job15581bb6820 to C:/SEE/RC4_15AUG2016/see-full/MDL_IDE/workspace/OnlineTraining_OncologyAdvanced/NONMEM
## Done.
##
##
## The following main elements were parsed successfully:
## RawResults
## Estimation::PopulationEstimates::MLE
## Estimation::IndividualEstimates::Estimates
## Estimation::Residuals::ResidualTable
## Estimation::Predictions
## Estimation::OFMeasures::Deviance
##
## The following MESSAGEs were raised during the job execution:
## estimation_successful: 1
## covariance_step_run: 0
## rounding_errors: 0
## hessian_reset: 0
## zero_gradients: 0
## final_zero_gradients: 0
## estimate_near_boundary: 0
## s_matrix_singular: 0
## significant_digits: 5.3
## nmoutput2so_version: This SOBlock was created with nmoutput2so version 4.5.27
##
## Completed
## -- Thu Aug 18 10:47:17 2016
Retrieve the parameter estimates
parameterEstimates_nm <- getPopulationParameters(nm_int, what="estimates")$MLE
print(parameterEstimates_nm)
## CV POP_IC50COMBO POP_LAMBDA0 POP_LAMBDA1 POP_K1
## 0.2589510 3.6552100 0.1461500 0.1302580 7.5956700
## POP_K2 POP_IC50 POP_W0
## 0.2341880 1.7596600 0.0612893
Compare model predictions and data
MDL -> PharmML
model_int.pharmML <- as.PharmML(mdlfile_int.FOCE)
Parameters to use in simulation
parValues <- parameterEstimates_nm
Simulate predictions
resSim_int <- simulx(model = basename(model_int.pharmML),
parameter = parValues,
treatment = list(adm1,adm2),
output = list(f,pkA,pkB))
Compare the simulated curves with the data from the combination experiment
mean_int <- resSim_int$WTOT$WTOT
sd_int <- parameterEstimates_nm["CV"]*resSim_int$WTOT$WTOT
resSim_int$WTOT$LL <- qnorm(p = 0.025, mean_int, sd_int)
resSim_int$WTOT$UL <- qnorm(p = 0.975, mean_int, sd_int)
print(ggplot(data=resSim_int$WTOT) + geom_line(aes(x = time, y=WTOT)) +
geom_ribbon(aes(x=time,ymin=LL, ymax=UL), alpha=0.3) +
geom_point(data=data_rocchetti_2013[data_rocchetti_2013$MDV==0,], aes(x=TIME, y=DV)) +
xlab("Time (days)") +
ylab("Tumor weight (g)") +
ggtitle("Rocchetti TGI interaction"))
The interaction model describes the data much better than the non-interaction model.
Tasks
simulx()
(you may use the R functions defined in helper_functions.R to set up simulx input)Load the PopED package for Optimal Design and the gridExtra package for graphical display
library("PopED")
library("gridExtra")
Explore the MDL file “Rocchetti_design.mdl”
Among the DESIGN_PARAMETERS, only the ones for which boundaries are given in the DESIGN_SPACES are optimised over. Therefore, cycle lengths DT1, DT2, DT3 are optimised, each within the permissible interval of 10 to 40 days. In contrast, indDoseA and indDoseB (the individual doses of Bevacizumab and NMS-P937, respecively) are fixed.
PopED maximises an objective function over a parameter space. The file 'tumorcriteria.txt' specifies how this objective function relates to the model prediction which is defined in the mdl file:
model_prediction is used, which is exactly this observation at time 120.Therefore, 'tumorcriteria.txt' together with the mdl file tells PopED to minimise WTOT at time 120.
Set name of .mdl file
mdlfile <- "Rocchetti_design.mdl"
Have a look at the design object
dsgnObj <- getDesignObjects(mdlfile)[[1]]
str(dsgnObj)
Optimise the dosing intervals using PopED
Convert MDL to PharmML
xml_TS <- as.PharmML(mdlfile)
The function as.poped() creates a set of interrelated R objects in the global environment which are needed for PopED estimation:
.oldobjs <- ls()
as.poped(xml_TS)
(poped_objs <- setdiff(ls(), .oldobjs))
## [1] "feps" "ff" "ode_func"
## [4] "poped.db" "proportionalError" "sfg"
## [7] "tumorcriteria"
To learn more about PopED and what the PopED functions do, you can have a look into the package help:
help(package = "PopED")
We will also have a first look into a PopED database, the central element of the PopED package. Since the object is highly nested, we will only show the first level of the structure of the PopED output. Hovering over the object, we can also obtain more information through the MDL-IDE.
str(poped.db, max.level = 1)
## List of 5
## $ settings :List of 83
## $ model :List of 5
## $ design :List of 7
## $ design_space:List of 19
## $ parameters :List of 22
Plot PopED prediction with initial design.
?plot_model_prediction
plot_model_prediction(poped.db, model_minxt = 0, model_maxxt = 120,model_num_points = 1000,
y_lab = "Tumour weight (g)") + ggtitle("Initial design")
Run PopED optimisation
?poped_optim
out_opt <- poped_optim(poped.db, opt_a = TRUE, parallel=TRUE, method = c("LS"), evaluate_fim = FALSE,
out_file = "poped_minimum_tumour_size.txt")
## ===============================================================================
## Initial design evaluation
##
## Initial OFV = -6.09618
##
## Efficiency criterion [usually defined as OFV^(1/npar)] = -6.09618
## *******************************************
## Running Line Search Optimization
## *******************************************
##
## Initial parameters: 25 25 25
## Initial OFV: -6.096182
##
## Searching parameter 3
## Changed from 25 to 35 ; OFV = -4.96735
## Searching parameter 2
## Changed from 25 to 35 ; OFV = -4.05128
## Searching parameter 1
## Changed from 25 to 35 ; OFV = -3.59864
##
## Elapsed time: 228.52 seconds.
##
## Final OFV = -3.598644
## Parameters: 35 35 35
##
## ===============================================================================
## FINAL RESULTS
##
## Optimized Covariates:
## Group 1: 35 : 35 : 35
##
## OFV = -3.59864
##
## Efficiency:
## (ofv_final / ofv_init) = 0.59031
##
## Total running time: 228.88 seconds
A successful PopED optimisation also returns a poped database with updated design parameters
db_opt <- out_opt$poped.db
Extract the design parameters from the PopED database. At present, this is done via a manually defined function. In the future, Standard Output will also be available from PopED and functions of the ddmore R package will allow easy manipulation of these objects in a similar fashion.
(parms_opt <- getDesignParms(db_opt))
## DT_3 DT_2 DT_1
## 35 35 35
Graphical comparison of designs using simulx
Add the dose levels
dosingAB <- c(indDoseA = 0.135, indDoseB = 113)
parms_opt <- c(parms_opt, dosingAB)
For comparison, let us also simulate with the initial design
parms_init <- getDesignParms(poped.db)
parms_init <- c(parms_init, dosingAB)
Extract dosing information from the design object and the parameter vector and create dosing input for simulx(), using a manually defined helper function.
adm_init <- updatedDosing2simulx(dsgnObj, parms_init)
adm_opt <- updatedDosing2simulx(dsgnObj, parms_opt)
Variables to be outputted by simulx
WTOT <- list(name="WTOT", time= seq(0,150, by=0.1))
Define the two groups to be compared:
gr1 <- list(treatment=adm_init)
gr2 <- list(treatment=adm_opt)
Extract the parameters from the MDL file in a simulx-compatible way. This will be also part of the as.simulx() function under development.
parObj <- getParameterObjects(mdlfile)[[1]]
p <- parobj2simulx(parObj)
Call simulx
resSim <- simulx(model = xml_TS, parameter = p, group = list(gr1,gr2), output = WTOT)
Create a meaningful legend for ggplot()
nameLookup <- c("Initial schedule","Minimum final tumour size")
res2 <- resSim$WTOT
res2$Design <- nameLookup[res2$group]
Display the results using simulx()
plot1 <- ggplot(data=res2, aes(x=time, y=WTOT, color=Design)) + geom_line(size=0.05) + ggtitle("Rocchetti TGI combo - Comparison between treatments") + xlab("Time (days)") + ylab("Tumor weight (g)")
print(plot1)
In terms of tumour volume 120 days after start of treatment, it is best to increase the cycle length. This comes at the expense, however, of a higher tumour load earlier on.
Tasks
simulx() (3 designs: initial,
optimised for tumour size and optimised for AUC; 2 outputs: tumour size and AUC)Take a look at the MDL file with appropriate modifications in the solutions folder
mdlfile_AUC <- "Rocchetti_design_AUC.mdl"
An AUC compartment can be interpreted in two ways:
As before, convert MDL -> PharmML -> PopED database
mdlfile_AUC <- "Rocchetti_design_AUC.mdl"
xml_AUC <- as.PharmML(mdlfile_AUC)
Create a PopED database.
as.poped(xml_AUC)
Optimisation in PopED
out_AUC <- poped_optim(poped.db, opt_a = TRUE, parallel=TRUE, method = c("LS"), evaluate_fim=FALSE,
out_file = "poped_minimum_AUC.txt")
## ===============================================================================
## Initial design evaluation
##
## Initial OFV = -285.575
##
## Efficiency criterion [usually defined as OFV^(1/npar)] = -285.575
## *******************************************
## Running Line Search Optimization
## *******************************************
##
## Initial parameters: 25 25 25
## Initial OFV: -285.5747
##
## Searching parameter 1
## Changed from 25 to 22.7551 ; OFV = -285.18
## Searching parameter 3
## Changed from 25 to 32.551 ; OFV = -278.447
## Searching parameter 2
## Changed from 25 to 21.1224 ; OFV = -276.91
##
## Elapsed time: 236.09 seconds.
##
## Final OFV = -276.91
## Parameters: 22.7551 21.12245 32.55102
##
## ===============================================================================
## FINAL RESULTS
##
## Optimized Covariates:
## Group 1: 22.7551 : 21.1224 : 32.551
##
## OFV = -276.91
##
## Efficiency:
## (ofv_final / ofv_init) = 0.96966
##
## Total running time: 236.13 seconds
Extract the design parameters
parms_AUC <- getDesignParms(out_AUC$poped.db)
parms_AUC
## DT_3 DT_2 DT_1
## 22.75510 21.12245 32.55102
parms_opt
## DT_3 DT_2 DT_1 indDoseA indDoseB
## 35.000 35.000 35.000 0.135 113.000
When minimising the final tumour size, it is best to increase cycle length. In contrast, to minimise AUC, the first cycle shoud be longest and subsequent cycles shorter.
Graphical comparison of designs using simulx
Add the dose levels manually
parms_AUC <- c(parms_AUC, dosingAB)
Create dosing input for simulx()
adm_AUC <- updatedDosing2simulx(dsgnObj, parms_AUC)
Variables to be outputted by simulx
WTOT <- list(name="WTOT", time= seq(0,120, by=0.1))
AUC <- list(name="AUC", time=seq(0,120, by=0.1))
Define three groups to be compared:
gr3 <- list(treatment=adm_AUC)
gr_AUC <- list(gr1,gr2,gr3)
Call simulx
res_AUC <- simulx(model = xml_AUC, parameter = p, group = gr_AUC, output = list(WTOT,AUC))
Create a meaningful legend for ggplot()
nameLookup2 <- c("Initial schedule","Minimum final tumour size","Minimum AUC")
res2_WTOT <- res_AUC$WTOT
res2_AUC <- res_AUC$AUC
res2_WTOT$Design <- nameLookup2[res2_WTOT$group]
res2_AUC$Design <- nameLookup2[res2_AUC$group]
Display the results
plot_WTOT <- ggplot(data=res2_WTOT, aes(x=time, y=WTOT, color=Design)) + geom_line(size=0.05) + ggtitle("Modified optimisation criterion") + xlab("Time (days)") + ylab("Tumor weight (g)")
plot_AUC <- ggplot(data=res2_AUC, aes(x=time, y=AUC, color=Design)) + geom_line(size=0.05) + ggtitle("Modified optimisation criterion") + xlab("Time (days)") + ylab("Area under tumour size curve")
grid.arrange(plot_WTOT,plot_AUC)
Optimising AUC at 120 days requires also leads to an improvement in final tumour size, but much less than actually minimising it. At day 120, the schedule minimising final tumour size has a ~20% higher mean tumour volume than the other two schedules.
Tasks
Rocchetti_design.mdl and the same individual dosing amounts of both drugs.simulx()As before, convert MDL -> PharmML -> PopED database
mdlfile_order <- "Rocchetti_design_order.mdl"
xml_order <- as.PharmML(mdlfile_order)
remove the previous poped database before creating a new one
rm(list=poped_objs)
as.poped(xml_order)
Optimisation in PopED
out_order <- poped_optim(poped.db, opt_a = TRUE, parallel=TRUE, method = c("LS"), evaluate_fim=FALSE,
out_file = "poped_minimum_order.txt")
## ===============================================================================
## Initial design evaluation
##
## Initial OFV = -6.62091
##
## Efficiency criterion [usually defined as OFV^(1/npar)] = -6.62091
## *******************************************
## Running Line Search Optimization
## *******************************************
##
## Initial parameters: 0
## Initial OFV: -6.620913
##
## Searching parameter 1
## Changed from 0 to 3.87755 ; OFV = -6.09074
##
## Elapsed time: 81.53 seconds.
##
## Final OFV = -6.090735
## Parameters: 3.877551
##
## ===============================================================================
## FINAL RESULTS
##
## Optimized Covariates:
## Group 1: 3.87755
##
## OFV = -6.09074
##
## Efficiency:
## (ofv_final / ofv_init) = 0.91992
##
## Total running time: 81.55 seconds
Extract the design parameters
parms_order_init <- getDesignParms(poped.db)
(parms_order_opt <- getDesignParms(out_order$poped.db))
## DT
## 3.877551
Show model predictions
model_prediction(poped.db)$PRED
## [1] 6.620913
model_prediction(out_order$poped.db)$PRED
## [1] 6.090735
Graphical comparison of designs using simulx
Add dosing parameters
parms_order_init <- c(parms_order_init, dosingAB)
parms_order_opt <- c(parms_order_opt, dosingAB)
Create a new design object
dsgnObj3 <- getDesignObjects(mdlfile_order)[[1]]
Create dosing input for simulx()
adm_order_init <- updatedDosing2simulx(dsgnObj3, parms_order_init)
adm_order_opt <- updatedDosing2simulx(dsgnObj3, parms_order_opt)
Variables to be outputted by simulx
WTOT <- list(name="WTOT", time= seq(0,120, by=0.1))
Define the two groups to be compared:
gr_order <- list(list(treatment=adm_order_init),
list(treatment=adm_order_opt))
Call simulx
res_order <- simulx(model = xml_order, parameter = p, group = gr_order, output = WTOT)
Create a meaningful legend for ggplot()
res2_order <- res_order$WTOT
res2_order$Design <- nameLookup[res2_order$group]
Display the results using simulx()
plot3 <- ggplot(data=res2_order, aes(x=time, y=WTOT, color=Design)) + geom_line(size=0.05) + ggtitle("Optimisation of combination schedule") + xlab("Time (days)") + ylab("Tumor weight (g)")
print(plot3)
It is best to give first NMS-P937, then Bevacizumab within a cycle. The three day delay used in Exercise 2.1 is almost optimal: the optimum delay is ~3.9 days. However, the combination schedule has a much lower impact on final tumour sizes than cycle length.