Observations and predictions
Predicted concentration
Parameters
Statistical results
ui.R
#------------------------------------------------------------------------- # 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
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
Estimation
in theOutputs
menu),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