Observations and predictions
Predicted concentration
Parameters
Statistical results
ui.R
Download
#-------------------------------------------------------------------------
# This application is governed by the CeCILL-B license.
# You can use, modify and/ or redistribute this code under the terms
# of the CeCILL license: http://www.cecill.info/index.en.html
#
# Marc Lavielle, Inria Saclay
# May 4th, 2015
#-------------------------------------------------------------------------
library(shinydashboard)
sidebar <- dashboardSidebar(
hr(),
sidebarMenu(id="tabs",
menuItem("Plots", icon = icon("line-chart"),
menuSubItem("observations vs time", tabName = "plot1", icon = icon("angle-right")), selected=TRUE,
menuSubItem("observations vs predictions", tabName = "plot2", icon = icon("angle-right")),
menuSubItem("residuals vs time", tabName = "plot3", icon = icon("angle-right")),
menuSubItem("residuals vs predictions", tabName = "plot4", icon = icon("angle-right"))
),
menuItem("Tables", icon = icon("table"), tabName="tabtable",
menuSubItem("parameters", tabName = "tabparam", icon = icon("angle-right")),
menuSubItem("observations and predictions", tabName = "tabobs", icon = icon("angle-right")),
menuSubItem("predicted concentration", tabName = "tabconc", icon = icon("angle-right")),
menuSubItem("statistics", tabName = "tabstat", icon = icon("angle-right"))
),
menuItem("Codes", icon = icon("file-text-o"),
menuSubItem("Mlxtran", tabName = "mlxtran", icon = icon("angle-right")),
menuSubItem("ui.R", tabName = "ui", icon = icon("angle-right")),
menuSubItem("server.R", tabName = "server", icon = icon("angle-right"))
),
menuItem("ReadMe", tabName = "readme", icon = icon("mortar-board")),
menuItem("About", tabName = "about", icon = icon("question"))
),
hr(),
conditionalPanel("input.tabs=='plot1' | input.tabs=='plot2' | input.tabs=='plot3' | input.tabs=='plot4'",
sidebarMenu(
menuItem("Outputs", icon = icon("chevron-circle-right"),
fluidRow(
column(1),
column(10,
checkboxInput("checkIni", "Initialization", value=TRUE),
checkboxInput("checkMLE", "Estimation", value=FALSE),
checkboxInput("rsse", "RSSE", value=TRUE)
)
)
),
menuItem("Scale", icon = icon("chevron-circle-right"),
fluidRow(
column(1),
column(10, radioButtons("log","",c("linear" = "1","log" = "2")))
)
)
),
hr()
)
)
body <- dashboardBody(
tabItems(
tabItem(tabName = "plot1",
fluidRow(
box(width = 12, status = "primary",
plotOutput("plot1", height="400px")
)
)
),
tabItem(tabName = "plot2",
fluidRow(
box(width = 12, status = "primary",
plotOutput("plot2", height="400px")
)
)
),
tabItem(tabName = "plot3",
fluidRow(
box(width = 12, status = "primary",
plotOutput("plot3", height="400px")
)
)
),
tabItem(tabName = "plot4",
fluidRow(
box(width = 12, status = "primary",
plotOutput("plot4", height="400px")
)
)
),
tabItem(tabName = "tabobs",
box( width = NULL, status = "primary", solidHeader = TRUE, title="Observations and predictions",
tableOutput("tabley")
)
),
tabItem(tabName = "tabconc",
box( width = NULL, status = "primary", solidHeader = TRUE, title="Predicted concentration",
tableOutput("tablef")
)
),
tabItem(tabName = "tabparam",
box( width = NULL, status = "primary", solidHeader = TRUE, title="Parameters",
tableOutput("tablep")
)
),
tabItem(tabName = "tabstat",
box( width = NULL, status = "primary", solidHeader = TRUE, title="Statistical results",
tableOutput("tabler")
)
),
tabItem(tabName = "ui",
box( width = NULL, status = "primary", solidHeader = TRUE, title="ui.R",
downloadButton('downloadData2', 'Download'),
br(),br(),
pre(includeText("ui.R"))
)
),
tabItem(tabName = "server",
box( width = NULL, status = "primary", solidHeader = TRUE, title="server.R",
downloadButton('downloadData3', 'Download'),
br(),br(),
pre(includeText("server.R"))
)
),
tabItem(tabName = "mlxtran",
box( width = NULL, status = "primary", solidHeader = TRUE, title="tgiModel.txt",
pre(includeText("tgiModel.txt"))
)
),
tabItem(tabName = "readme",
withMathJax(),
includeMarkdown("readMe.Rmd")
),
tabItem(tabName = "about",
includeMarkdown("../../about/about.Rmd")
)
),
conditionalPanel(condition = "(input.checkIni=='1') & (input.tabs=='plot1'|input.tabs=='plot2'|input.tabs=='plot3'|input.tabs=='plot4') ",
box(width = 4, status = "primary", solidHeader = FALSE,
sliderInput("l0", "l0:", value=0.35, min=0.05, max=0.6, step=0.05),
conditionalPanel(condition = "input.checkMLE == '1'",
tags$style(type='text/css', '#estl0 {color: red; font-size: 14px;}'),
verbatimTextOutput("estl0")
),
sliderInput("l1", "l1:", value=0.4, min=0.05, max=1, step=0.05),
conditionalPanel(condition = "input.checkMLE == '1'",
tags$style(type='text/css', '#estl1 {color: red; font-size: 14px;}'),
verbatimTextOutput("estl1")
)
),
box(width = 4, status = "primary", solidHeader = FALSE,
sliderInput("w0", "w0:", value=0.036, min=0.01, max=0.06, step=0.002),
conditionalPanel(condition = "input.checkMLE == '1'",
tags$style(type='text/css', '#estw0 {color: red; font-size: 14px;}'),
verbatimTextOutput("estw0")
),
sliderInput("kpot", "kpot:", value=0.006, min=0.001, max=0.02, step=0.001),
conditionalPanel(condition = "input.checkMLE == '1'",
tags$style(type='text/css', '#estkpot {color: red; font-size: 14px;}'),
verbatimTextOutput("estkpot")
)
),
box(width = 4, status = "primary", solidHeader = FALSE,
sliderInput("tau", "tau:", value=3, min=0.5, max=10, step=0.5),
conditionalPanel(condition = "input.checkMLE == '1'",
tags$style(type='text/css', '#esttau {color: red; font-size: 14px;}'),
verbatimTextOutput("esttau")
)
)
)
)
dashboardPage(
dashboardHeader(title = "Tumor growth model"),
sidebar,
body
)
server.R
Download
library(mlxR)
shinyServer(function(input, output) {
d <- read.csv(file="tgiData.txt", header=TRUE, sep="\t", quote="\"")
names(d)[3] <- "y"
fc <- list(name='w', time=seq(0,25,length.out=100))
ypred <- list(name='ypred', time=d$time[d$id==1])
adm2 <- list(time=seq(12,16), amount=100)
adm1 <- list(time=12, amount=0)
g <- list(list(treatment=adm1), list(treatment=adm2))
d.param <- 4
n <- length(d$time)
r <- reactive({
p.name <- c('l0','l1','w0','kpot','tau')
p.value <- c(input$l0,input$l1,input$w0,input$kpot,input$tau)
model <- "tgiModel.txt"
parameter <- data.frame(initial=p.value)
row.names(parameter) <- p.name
p <- list(name=p.name, value=p.value)
res <- simulx(model=model,output=list(fc,ypred),group=g,parameter=p)
f <- data.frame(time=res$w$time, id=res$w$id, f_ini=res$w$w)
d$id <- res$ypred$id
d$ypred_ini <- res$ypred$ypred
r=list()
r$sse_ini=sum((d$y-d$ypred_ini)^2)
r$textini <- paste("RSSE =",formatC(r$sse_ini,digits=3))
result <- rsse(r$sse_ini,n,d.param)
row.names(result)="initial"
if (input$checkMLE==TRUE){
p <- list(name=p.name, value=p.value)
rest <- mlxoptim(model=model,
output='w',
group=g,
data=d,
initial=p)
r$sse_est=rest$sse
r$textest <- paste("RSSE =",formatC(r$sse_est,digits=3))
s <- list(load.design=TRUE)
res <- simulx(model=model,output=list(fc,ypred),group=g,parameter=rest$parameter,settings=s)
d$ypred_est <- res$ypred$ypred
f$f_est <- res$w$w
parameter[,2] <-rest$parameter$value
names(parameter)[2]="estimated"
res.est <- rsse(rest$sse,n,d.param)
row.names(res.est) <- "estimate"
result[2,] <- res.est
}
r$parameter <- parameter
r$y <- d
r$f <- f
r$result <- t(result)
return(r)
})
szp <- 3.5
szl <- 1
dxini <- 0.75
dxest <- 0.75
dyini <- 1.07
dyest <- 1
output$plot1 <- renderPlot({
r=r()
pl <- ggplotmlx()
pl <- pl + geom_point(data=r$y, aes(x=time, y=y, colour=id), size=szp)
if (input$checkIni==TRUE){
pl <- pl + geom_line(data=r$f, aes(x=time, y=f_ini, colour=id), size=szl)
}
if (input$checkMLE==TRUE){
pl <- pl + geom_line(data=r$f, aes(x=time, y=f_est, colour=id), size=szl)
}
# pl <- pl + scale_colour_manual(
# values=c("a"="#F8766D", "b"="#00BFC4", "c"="#B79F00",
# "d"="#F564E3", "e"="#00BA38", "f"="#619CFF" ),
# labels=c("a"="first order", "b"="zero order", "c"="alpha order",
# "d"="sequential 0-1", "e"="simultaneous 0-1", "f"="saturated"))
if (input$log==2){
pl <- pl + scale_y_log10()
}
my.ggp.yrange <- ggplot_build(pl)$panel$ranges[[1]]$y.range
my.ggp.xrange <- ggplot_build(pl)$panel$ranges[[1]]$x.range
xini <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxini
xest <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxest
yini <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyini
yest <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyest
if (input$checkIni==TRUE & input$rsse==TRUE){
pl <- pl + annotate("text", label=r$textini, x=xini, y=yini, size=6, colour="black", hjust = 0)
}
if (input$checkMLE==TRUE & input$rsse==TRUE){
pl <- pl + annotate("text", label=r$textest, x=xest, y=yest, size=6, colour="red", hjust = 0)
}
print(pl)
})
output$plot2 <- renderPlot({
r=r()
pl <- ggplotmlx()
if (input$checkIni==TRUE){
pl <- pl + geom_point(data=r$y, aes(x=ypred_ini, y=y), size=szp)
}
if (input$checkMLE==TRUE){
pl <- pl + geom_point(data=r$y, aes(x=ypred_est, y=y), color="red", size=szp)
}
pl <- pl + geom_abline(intercept = 0, slope = 1, color="blue", size=szl) + xlab("y_pred")
my.ggp.yrange <- ggplot_build(pl)$panel$ranges[[1]]$y.range
my.ggp.xrange <- ggplot_build(pl)$panel$ranges[[1]]$x.range
xini <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxini
xest <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxest
yini <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyini
yest <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyest
if (input$checkIni==TRUE & input$rsse==TRUE){
pl <- pl + annotate("text", label=r$textini, x=xini, y=yini, size=6, colour="black", hjust = 0)
}
if (input$checkMLE==TRUE & input$rsse==TRUE){
pl <- pl + annotate("text", label=r$textest, x=xest, y=yest, size=6, colour="red", hjust = 0)
}
print(pl)
})
output$plot3 <- renderPlot({
r=r()
pl <- ggplotmlx()
if (input$checkIni==TRUE){
r$y$e <- r$y$y - r$y$ypred_ini
pl <- pl + geom_point(data=r$y, aes(x=time, y=e), size=szp)
}
if (input$checkMLE==TRUE){
r$y$e <- r$y$y - r$y$ypred_est
pl <- pl + geom_point(data=r$y, aes(x=time, y=e), color="red", size=szp)
}
pl <- pl + geom_hline(aes(yintercept=0), color="blue", size=szl)
my.ggp.yrange <- ggplot_build(pl)$panel$ranges[[1]]$y.range
my.ggp.xrange <- ggplot_build(pl)$panel$ranges[[1]]$x.range
xini <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxini
xest <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxest
yini <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyini
yest <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyest
if (input$checkIni==TRUE & input$rsse==TRUE){
pl <- pl + annotate("text", label=r$textini, x=xini, y=yini, size=6, colour="black", hjust = 0)
}
if (input$checkMLE==TRUE & input$rsse==TRUE){
pl <- pl + annotate("text", label=r$textest, x=xest, y=yest, size=6, colour="red", hjust = 0)
}
print(pl)
})
output$plot4 <- renderPlot({
r=r()
pl <- ggplotmlx()
if (input$checkIni==TRUE){
r$y$e <- r$y$y - r$y$ypred_ini
pl <- pl + geom_point(data=r$y, aes(x=ypred_ini, y=e), size=szp) }
if (input$checkMLE==TRUE){
r$y$e <- r$y$y - r$y$ypred_est
pl <- pl + geom_point(data=r$y, aes(x=ypred_est, y=e), color="red", size=szp) }
pl <- pl + geom_hline(aes(yintercept=0), color="blue", size=szl)
my.ggp.yrange <- ggplot_build(pl)$panel$ranges[[1]]$y.range
my.ggp.xrange <- ggplot_build(pl)$panel$ranges[[1]]$x.range
xini <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxini
xest <- my.ggp.xrange[1] + (my.ggp.xrange[2]-my.ggp.xrange[1])*dxest
yini <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyini
yest <- my.ggp.yrange[1] + (my.ggp.yrange[2]-my.ggp.yrange[1])*dyest
if (input$checkIni==TRUE & input$rsse==TRUE){
pl <- pl + annotate("text", label=r$textini, x=xini, y=yini, size=6, colour="black", hjust = 0)
}
if (input$checkMLE==TRUE & input$rsse==TRUE){
pl <- pl + annotate("text", label=r$textest, x=xest, y=yest, size=6, colour="red", hjust = 0)
}
print(pl)
})
output$tabley <- renderTable({
r()$y
})
output$tablef <- renderTable({
r()$f
})
output$tablep <- renderTable({
r()$parameter
})
output$tabler <- renderTable({
r()$result
})
output$estl0 <- renderText({
r <- round(r()$parameter[1,],3)
pest=NULL
if (length(r)==2)
pest <- paste0("l0_est = ",r[2])
return(pest)
})
output$estl1 <- renderText({
r <- round(r()$parameter[2,],3)
pest=NULL
if (length(r)==2)
pest <- paste0("l1_est = ",r[2])
return(pest)
})
output$estw0 <- renderText({
r <- round(r()$parameter[3,],3)
pest=NULL
if (length(r)==2)
pest <- paste0("w0_est = ",r[2])
return(pest)
})
output$estkpot <- renderText({
r <- round(r()$parameter[4,],3)
pest=NULL
if (length(r)==2)
pest <- paste0("kpot_est = ",r[2])
return(pest)
})
output$esttau <- renderText({
r <- round(r()$parameter[5,],3)
pest=NULL
if (length(r)==2)
pest <- paste0("tau_est = ",r[2])
return(pest)
})
})
rsse <- function(sse,n,d)
{
deviance <- n*(1+ log(sse/n) +log(2*pi))
BIC <- deviance + log(n)*d
AIC <- deviance + 2*d
r <- data.frame(n=n,d=d,sse=sse,deviance=deviance,AIC,BIC)
}
#--------------------------------------
mlxoptim <- function(model=NULL,output=NULL,group=NULL,data=NULL,initial=NULL)
{
errpred <- function(pp,y,dataIn){
dataIn$population_parameters$value=exp(pp)
res = simulx(data=dataIn)
y.pred=res$w$w
e=sum((y.pred-y)^2)
return(e)}
t.obs <- data[data$id==1,2]
y.obs <- data[,3]
dd <- simulx(model=model,
parameter=initial,
output=list(name=output,time=t.obs),
group=group,
settings=list(data.in=TRUE,load.design=TRUE))
dd$population_parameters=dd$individual_parameters
dd$population_parameters$value=dd$population_parameters$value[1,]
dd$individual_parameters=NULL
res = simulx(data=dd)
l0=log(initial$value)
r <- optim(l0,errpred,y=y.obs,dataIn=dd,control=list(maxit=500))
n <- length(t.obs)
pest=initial
pest$value=exp(r$par)
final = list(parameter=pest,sse=r$value)
return(final)
}
tgiModel.txt
[LONGITUDINAL]
input = {l0,l1,w0,kpot,tau}
PK:
depot(target=Ad)
EQUATION:
ka = 20
V = 2.8
ke = 2.4
t0 = 0
Ac_0 = 0
p_0 = w0
d_0 = 0
Cc = Ac/V
ddt_Ad = -ka*Ad
ddt_Ac = ka*Ad - ke*Ac
ddt_p = (2*l0*l1*p^2)/((l1+2*l0*p)*(p+d)) - kpot*Cc*p
ddt_d = kpot*Cc*p - kpot*delay(Ac,tau)/V*delay(p,tau)
w = p+d
ypred = w
Objective
Fit the same tumor growth model to tumor size data from two individuals who received two different treatments.
The data
Simulated tumor sizes from two patients are displayed:
patient 1 (red dots) didn't receive any treatment
- patient 2 (blue dots) received a dose of 100 units every month between month 12 and month 16
The model
Simeoni et al. presented a model for tumor growth and anti-cancer effects where a transit compartment model describes the apoptosis of tumor cells attacked by a drug. Following Koch et al. (2014), this model can be rewritten with a lifespan model
\[ \begin{aligned} \dot{Ad}(t) &= -k_a \, A_d(t) \\ \dot{Ac}(t) &= k_a \, A_d(t) - k_e \, A_c(t) \\ \dot{p}(t) &= \frac{2 l_0 l_1 p(t)^2}{(l_1+2 l_0 p(t))(p(t)+d(t))} - k_{\rm pot} C_c(t) p(t) \\ \dot{d}(t) &= k_{\rm pot} C_c(t) p(t) - k_{\rm pot} p(t-\tau) \frac{A_c(t-\tau)}{V} \end{aligned} \] where \(C_c(t) = A_c(t)/V\) and assuming the following initial conditions for \(t\leq 0\):
\[ \begin{aligned} A_c(t) &= 0 \\ p(t) &= w_0 \\ d(t) &= 0 \end{aligned} \]
The predicted tumor size given by this model is \[ w(t) = p(t) + d(t) \]
Tasks
The objective of this exercice is to fit the same tumor growth model to this data, using the same parameters for both patients.
The PK parameters \((k_a, V, k_e)\) are supposed to be known. We must therefore estimate \(\theta = (l_0,l_1,w_0,k_{\rm pot},\tau)\).
Check the checkbox
Init. Then, using the sliders, try to find some values of the components of \(\theta\) that fit well the observed data for both patients,Compute the least-squares estimate of \(\theta\) (by checking
Estimationin theOutputsmenu),Look at the diagnostic plots and the numerical results in the Table
Results,
References
Simeoni M, Magni P, Cammia C, De Nicolao G, Croci V, Pesenti E, Germani M, Poggesi I, Rocchetti M (2004) Predictive pharmacokinetic-pharmacodynamic modeling of tumor growth kinetics in xenograft models after administration of anticancer agents. Cancer Res 64(3):1094-1101.
Koch G, Krzyzanski W, P
This application is governed by the CeCILL-B license.
You can use, modify and/or redistribute these codes under the terms of the
CeCILL license.
This Shiny application requires the mlxR package.
Marc Lavielle
Inria Saclay, Popix team
April 29th, 2015