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