Individual fitting - TG model
















[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
shinyUI(fluidPage(
  titlePanel(
    list(HTML('<p style="color:#4C0B5F; font-size:24px" fontsize=14>Individual fitting - TG model</p>' )),
    windowTitle="TG - indiv. fit."),
  tabsetPanel(
    tabPanel("Plot",
             fluidRow(
               column(2,
                      br(),br(),br(),
                      sliderInput("l0", "l0:", value=0.3, min=0.05, max=0.5, step=0.05),
                      br(),
                      sliderInput("l1", "l1:", value=0.3, min=0.05, max=0.5, step=0.05),
                      br(),
                      sliderInput("w0", "w0:", value=0.02, min=0.002, max=0.04, step=0.002),
                      br(),
                      sliderInput("kpot", "kpot:", value=0.01, min=0.001, max=0.02, step=0.001),
                      br(),
                      sliderInput("tau", "tau:", value=1, min=0.5, max=10, step=0.5),
                      br()
               ),
               column(10,
                      column(10, 
                             tabsetPanel( 
                               tabPanel("Fit", plotOutput("plot1")),
                               tabPanel("Obs vs Pred", plotOutput("plot2")), 
                               tabPanel("Residuals", plotOutput("plot3")) 
                             )),
                      column(2, 
                             br(),br(),br(),br(),br(),br(),
                             checkboxInput("checkObs", "Obs.", value=TRUE),
                             checkboxInput("checkIni", "Init.", value=TRUE),
                             checkboxInput("checkMLE", "Estim.", value=FALSE),
                             br(),
                             radioButtons("log", "", c("linear scale" = FALSE,"log scale" = TRUE))                    
                      )
               ))),
    tabPanel("Tables",
             tabsetPanel( 
               tabPanel("Parameters", tableOutput("tablep")),      
               tabPanel("Observations", tableOutput("tabley")),
               tabPanel("Concentration", tableOutput("tablef")),
               tabPanel("Results", tableOutput("tabler"))
             )),
    tabPanel("Mlxtran", pre(includeText("tgi_model.txt"))),
    tabPanel("ui.R", pre(includeText("ui.R"))),
    tabPanel("server.R", pre(includeText("server.R"))),
    tabPanel("ReadMe", withMathJax(), includeMarkdown("ReadMe.Rmd"))
  )
))
my.dir=getwd()
source("../../initMlxR.R")
setwd(my.dir)

shinyServer(function(input, output) {
    
  d <- read.csv(file="tgi_data.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])
  
  adm1 <- list(time=seq(12,16), amount=100) 
  adm2 <- list(time=12, amount=0) 
  g    <- list(list(treatment=adm1),list())
  
  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 <- "tgi_model.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[[1]][,2])
    f$id <- res[[1]][,1]
    f$f_ini <- res[[1]][,3]
    d$id <- res[[2]][,1]
    d$ypred_ini <- res[[2]][,3]
    
    r=list()
    r$sse_ini=sum((d$y-d$ypred_ini)^2)
    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
      s <- list(load.design=TRUE)
      res <- simulx(model=model,output=list(fc,ypred),group=g,parameter=rest$parameter,settings=s)
      d$ypred_est <- res[[2]][,3]
      f$f_est <- res[[1]][,3]
      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)
  })
  
  
  output$plot1 <- renderPlot({
    withProgress(message = 'Creating plot', detail='Wait...', value = 0.1, expr={
      r=r()
      pl <- ggplot() 
      if (input$checkObs==TRUE) {
        pl <- pl + geom_point(data=r$y, aes(x=time, y=y, colour=id),size=3) 
      }      
      if (input$checkIni==TRUE){
        pl <- pl + geom_line(data=r$f, aes(x=time, y=f_ini, colour=id), size=0.75) # +  
      }      
      if (input$checkMLE==TRUE){
        pl <- pl + geom_line(data=r$f, aes(x=time, y=f_est, colour=id), size=0.75) # +
      }    
      if (input$log==TRUE){
        pl <- pl + scale_y_log10()
      }
      print(pl)
    })
  })
  
  
  output$plot2 <- renderPlot({
    if (input$checkIni==TRUE | input$checkMLE==TRUE)  {
      r=r()
      pl <- ggplot() 
      
      if (input$checkIni==TRUE){
        pl <- pl + geom_point(data=r$y, aes(x=ypred_ini, y=y, colour=id),size=3) 
      }
      
      if (input$checkMLE==TRUE){
        pl <- pl + geom_point(data=r$y, aes(x=ypred_est, y=y, colour=id),size=3) 
      }
      pl <- pl + geom_abline(intercept = 0, slope = 1, color="blue") + xlab("y_pred")
      print(pl)
    }
  })
  
  
  output$plot3 <- renderPlot({
    if (input$checkIni==TRUE | input$checkMLE==TRUE){
      r=r()
      pl <- ggplot() 
      
      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, colour=id),size=2) 
      }
      
      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, colour=id),size=2)  
      }
      pl <- pl + geom_hline(aes(yintercept=0), color="blue")
      print(pl)
    }
  })
  
  output$tabley <- renderTable({ 
    r()$y
  })
  
  output$tablef <- renderTable({ 
    r()$f
  })
  
  output$tablep <- renderTable({
    r()$parameter
  })
  
  output$tabler <- renderTable({
    r()$result
  })
})

#--------------------------------------
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)
}

#--------------------------------------
list <- structure(NA,class="result")
"[<-.result" <- function(x,...,value) {
  args <- as.list(match.call())
  args <- args[-c(1:2,length(args))]
  length(value) <- length(args)
  for(i in seq(along=args)) {
    a <- args[[i]]
    if(!missing(a)) eval.parent(substitute(a <- v,list(a=a,v=value[[i]])))
  }
  x
}

#--------------------------------------
mlxoptim <- function(model=NULL,output=NULL,group=NULL,data=NULL,initial=NULL)
{
  errpred <- function(pp,y,dataIn,s){
    dataIn$population_parameters$value=exp(pp)
    res = simulx(data=dataIn,settings=s)
    #    y.pred=res[[1]][,3]
    y.pred=res$w$w
    e=sum((y.pred-y)^2)
    return(e)}
  
  t.obs <- data[data$id==1,2]
  y.obs <- data[,3]
  list[res,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)
  s <- list(load.design=FALSE)
  l0=log(initial$value)
  r <- optim(l0,errpred,y=y.obs,dataIn=dd,s=s,control=list(maxit=150))
  n <- length(t.obs)
  pest=initial
  pest$value=exp(r$par)
  final = list(parameter=pest,sse=r$value)
  return(final)
}
The data

Simulated tumor sizes from two patients are displayed:

  • patient 1 (red dots) received a dose of 100 units every month between month 12 and month 16

  • patient 2 (blue dots) didn't receive any treatment

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

Task

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 Estim.),

  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.