inst/doc/custom.R

In tci: Target Controlled Infusion (TCI)

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  echo = TRUE, 
  message = FALSE, 
  warning = FALSE, 
  fig.align="center",
  fig.height= 4, 
  fig.width = 6
)

## ----mrgsolve-model, eval = TRUE, echo = TRUE---------------------------------
library(tci)
library(mrgsolve) # implement ODE equation
library(xtable)   # printing tables
library(ggplot2)  # plotting results
library(reshape2) # melt function

form <- '
$PARAM V1 = 7.88, V2=23.9, V3=13.8, CL1=5, CL2=0.828, CL3=0.0784,
k10 = 0.172, k12=0.373, k21=0.103, k13=0.0367, k31=0.0124
$CMT A1 A2 A3
$ODE
dxdt_A1 = k21*A2 + k31*A3 - (k12+k13+k10)*A1 - CL1/V1*A1;
dxdt_A2 = k12*A1 - k21*A2 - CL2/V2*A2;
dxdt_A3 = k13*A1 - k31*A3 - CL3/V3*A3;
'
mrg_mod_remif <- mcode("remifentanil", form)

## ----remif-wrapper------------------------------------------------------------
pk_remif <- function(tm, kR, pars, init = c(0,0,0)){

  # allow lowercase names
  names(pars) <- toupper(names(pars))
  # store volume
  vols <- pars[c("V1","V2","V3")]
  A0 <- init*vols # initial amounts
  names(A0) <- c("A1","A2","A3") # names required by mrgsolve
  
  # pass parameters as list
  pars <- sapply(pars, as.list)

  # update parameters and initial values (as amounts)
  mrg_mod_remif <- update(mrg_mod_remif, param = pars, init = A0)

  # dosing regimen - mrgsolve function in terms of amount infused
  event <- ev(amt =  kR*max(tm), time = 0, tinf = max(tm))

  # simulate responses (skip tm=0 unless specified)
  dat <- mrgsim_q(x = mrg_mod_remif, # pk model
                  data = event, # dosing event
                  stime = tm) # evaluation times

  # skip tm=0 unless specified in tm
  dat <- dat@data[-1,]

  # return concentrations with compartments in rows and times in columns
  cons <- t(dat[,c("A1","A2","A3")]) / vols

  rownames(cons) <- colnames(cons) <- NULL
  return(cons)
}

## ----Cascone-table, results='asis', echo = FALSE------------------------------
tab_cascone <- data.frame(Parameter1 = c("$V_1$","$V_2$","$V_3$","$CL_1$",
                                         "$CL_2$","$CL_3$"),
                          Optimized_val1 = c("7.88 L","23.9 L","13.8 L",
                                             "2.08 L/min","0.828 L/min","0.0784 L/min"),
                          Parameter2 = c("$k_{10}$","$k_{12}$","$k_{21}$","$k_{13}$","$k_{31}$",""),
                          Optimized_val2 = c("0.172/min","0.373/min","0.103/min","0.0367/min","0.0124/min",""))

names(tab_cascone) <- c("Parameter", "Optimized value", "Parameter", "Optimized value")
tab_cascone_df <- xtable(tab_cascone,
                         caption = "Reproduction of Table 1 from Cascone et al. (2013): Values and dimensions of the three-compartmental model parameters.",
                         label = "tab:Cascone-tab1")
print(tab_cascone_df,
      include.rownames = FALSE,
      caption.placement = "top",
      timestamp = NULL,
      booktabs = TRUE,
      hline.after = c(-1,0,6),
      sanitize.text.function = function(x){x},
      include.colnames=TRUE)

## -----------------------------------------------------------------------------
dose_remi <- inf_manual(inf_tms = 0, inf_rate = 60, duration = 20)
pars_remif <- c(V1 = 7.88, V2=23.9, V3=13.8, CL1=5, CL2=0.828, CL3=0.0784,
                k10 = 0.172, k12=0.373, k21=0.103, k13=0.0367, k31=0.0124)
mod_remif <- pkmod(pkfn = pk_remif, pars_pk = pars_remif)
p1 <- predict(mod_remif, inf = dose_remi, tms = 0:80)
ggplot(melt(data.frame(time = 0:80, p1), id = "time"), 
       aes(x = time, y = value, color = variable)) +
  geom_line()

## ----alternate-effect-site-alg, echo = TRUE-----------------------------------
tci_plasma_lim <- function(Ct, pkmod, dtm = 1/6, maxrt = 1200,
                          lim_amt = 0.5, ecmpt = NULL, tmax_search = 20, 
                          cetol = 0.05, cptol = 0.1, ...){

  pkmod <- update(pkmod,...)
  
  # if effect-site concentration is close to target,
  # switch to plasma targeting
  if(with(pkmod,(Ct - init[ecmpt]) / Ct <  cetol &
     (Ct - init[pcmpt])/Ct <= cptol))
    return(tci_plasma(Ct, pkmod = pkmod, dtm = dtm, maxrt = maxrt))
  
  # maximum tolerable plasma concentration
  Cp_max <- Ct + lim_amt
  
  # infusion required to reach Cp_max
  pinf <- tci_plasma(Ct = Cp_max, pkmod = pkmod, dtm = dtm, maxrt = maxrt)
  
  # Administer dtm-minute infusion
  unit_inf <- inf_manual(inf_tms = 0, inf_rate = pinf, duration = dtm)
  
  # Calculate maximum effect-site concentration
  CeP <- function(tm) predict(pkmod, inf = unit_inf, tms = tm)[,pkmod$ecmpt]
  Ce_max <- optimize(CeP, c(0,20), maximum = TRUE)$objective
  
  # if max Ce < Ct administer infusion to reach maximum target
  if(Ce_max <= Ct + cetol*Ct) 
    infrt <- pinf
  else
    infrt <- tci_effect_only(Ct, pkmod, dtm, maxrt = maxrt)

  return(infrt)
}

## ----base-custom-alg, echo = TRUE---------------------------------------------
mod3ecpt <- pkmod(pars_pk = c(cl = 10, q2 = 2, q3 =20, v = 15, v2 = 30, v3 = 50, ke0 = 1.2))
tci_plasma_lim(Ct = 2, pkmod = mod3ecpt, lim_amt = 0.25)

## ----tci-custom-alg, echo = TRUE, fig.cap= "Evaluation of user-defined TCI effect-site algorithm. \\label{fig:custom-tci-alg}"----
# tci target concentrations
tci_targets <- cbind(value = c(1,2,2.5,2.5), time = c(0,3,7,10))
# calculate infusion schedule using plasma-limiting algorithm
plim_inf <- inf_tci(target_vals = c(1,2,2.5,2.5),
                    target_tms = c(0,3,7,10),
                    pkmod = mod3ecpt, 
                    custom_alg = tci_plasma_lim, 
                    lim_amt = 0.25)
head(plim_inf)

## -----------------------------------------------------------------------------
# effect-site targeting
eff_inf <- inf_tci(target_vals = c(1,2,2.5,2.5),
                    target_tms = c(0,3,7,10), 
                   pkmod = mod3ecpt, 
                   type = "effect")

## -----------------------------------------------------------------------------
# predict responses
tms_pred <- seq(0,10,0.1)
plim_pred <- predict(mod3ecpt, plim_inf, tms_pred)
eff_pred <- predict(mod3ecpt, eff_inf, tms_pred)

# plot results
dat <- data.frame(time = tms_pred, 
                  `plasma (custom)` = plim_pred[,"c1"], 
                  `effect (custom)` = plim_pred[,"c4"],
                  `plasma (effect)` = eff_pred[,"c1"], 
                  `effect (effect)` = eff_pred[,"c4"],
                  check.names = FALSE)
datm <- melt(dat, id = "time")
datm$algorithm <- ifelse(datm$variable %in% c("plasma (custom)","effect (custom)"),
                         "Plasma-limiting", "Effect-site")
ggplot(datm, aes(x = time, y = value, color = variable, linetype = algorithm)) + 
  geom_line() + 
  xlab("Minutes") + 
  ylab("Concentration (mg/L)") +
  ggtitle(label = "Plasma-limiting effect-site TCI algorithm")