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) }
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
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) \]
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 Estim.
),
Look at the diagnostic plots and the numerical results in the Table Results
,
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