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)\).


  1. 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,

  2. Compute the least-squares estimate of \(\theta\) (by checking Estimation in the Outputs menu),

  3. 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, Prez-Ruixo JJ, Schropp J. (2014) Modeling of delays in PKPD: classical approaches and a tutorial for delay differential equations. J Pharmacokinet Pharmacodyn. 41(4):291-318.


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