inst/doc/dtat.R

In DTAT: Dose Titration Algorithm Tuning

## ----setup, include=FALSE---------------------------------------------------------------
.First <- function() {
    library(knitr)
    library(rms)
    library(tidyr)
    library(latticeExtra)
    library(RColorBrewer)
    library(pomp)
}
.First()
# set global chunk options
options(replace.assign=TRUE,width=90)
options(datadist="dd")
options(contrasts=c("contr.treatment","contr.treatment"))
options(tinytex.engine_args = '-shell-escape')
opts_knit$set(eval.after="fig.cap")
opts_chunk$set(include=TRUE, tidy=FALSE, tidy.opts=list(width.cutoff=70), cache=FALSE, autodep=TRUE)
set.seed(2016)

## ----Population, results='asis'---------------------------------------------------------
N <- 25
dtx.mm <- 0.808 # molar mass (mg/µM) of docetaxel
pop <- data.frame(id=1:N
                  # From Friberg et al 2002 (Table 4, row 1), taking sdlog ~= CV
                  ,Circ0=rlnorm(N, meanlog=log(5050), sdlog=0.42) # units=cells/mm^3
                  ,MTT=rlnorm(N, meanlog=log(89.3), sdlog=0.16)  # mean transit time
                  ,gamma=rlnorm(N, meanlog=log(0.163), sdlog=0.039) # feedback factor
                  ,Emax=rlnorm(N, meanlog=log(83.9), sdlog=0.33)
                  ,EC50=rlnorm(N, meanlog=log(7.17*dtx.mm), sdlog=0.50)
                  # PK params from 2-compartment docetaxel model of Onoue et al (2016)
                  ,CL=rlnorm(N, meanlog=log(32.6), sdlog=0.295)
                  ,Q =rlnorm(N, meanlog=log(5.34), sdlog=0.551)
                  ,Vc=rlnorm(N, meanlog=log(5.77), sdlog=0.1) # Onoue gives no CV% for V1
                  ,Vp=rlnorm(N, meanlog=log(11.0), sdlog=0.598) # Called 'V2' in Onoue
)

pop <- upData(pop
             ,kTR = 4/MTT
             ,units = c(Circ0="cells/mm^3"
                       ,MTT="hours"
                       ,kTR="1/hour"
                       ,CL="L/h"
                       ,Q="L/h"
                       ,Vc="L"
                       ,Vp="L"
                       )
             ,print=FALSE
             )

latex(describe(pop), file="")

## ----pomp-PKPD--------------------------------------------------------------------------
pkpd.skel <- "
  double c2p = Q*( Cc - Cp ); // central-to-peripheral flux
  DCc = (dose/duration)*(t < duration ? 1.0 : 0.0)/Vc - (CL/Vc)*Cc - c2p/Vc;
  DCp = c2p/Vp;
  // Myelosuppression model (Emax model, then dynamics per eqs 3-7 from Friberg et al 2002
  double Edrug = Emax*Cc/(Cc + EC50); // classic 'Emax model'
  DProl = (1-Edrug) * Prol * kTR * pow((Circ0 / Circ), gamma)  -  kTR * Prol;
  DTx_1 = kTR * (Prol - Tx_1);
  DTx_2 = kTR * (Tx_1 - Tx_2);
  DTx_3 = kTR * (Tx_2 - Tx_3);
  DCirc = kTR * (Tx_3 - Circ);
"
pkpd.txform <- "
  T_Circ0 = log(Circ0);
  T_kTR = log(kTR);
  T_Emax = log(Emax);
  T_EC50 = log(EC50);
  T_CL = log(CL);
  T_Q = log(Q);
  T_Vc = log(Vc);
  T_Vp = log(Vp);
  T_sigma = log(sigma);
  T_dose = log(dose);
  T_duration = log(duration);
"
pkpd.txback <- "
  Circ0 = exp(T_Circ0);
  kTR = exp(T_kTR);
  Emax = exp(T_Emax);
  EC50 = exp(T_EC50);
  CL = exp(T_CL);
  Q = exp(T_Q);
  Vc = exp(T_Vc);
  Vp = exp(T_Vp);
  sigma = exp(T_sigma);
  dose = exp(T_dose);
  duration = exp(T_duration);
"

Tmax <- 21*24 # solve for full 21 days of 3-week cycle
df <- data.frame(time=c(seq(0.0, 1.95, 0.05) # q3min for 2h, 
                       ,seq(2.0, Tmax, 1.0)) # then hourly until Tmax
                ,y=NA)
pkpd <- pomp(data = df
            , times="time", t0=0
            , skeleton=vectorfield(Csnippet(pkpd.skel))
            , statenames=c("Cc", "Cp", "Prol", "Tx.1", "Tx.2", "Tx.3", "Circ")
            , paramnames=c("Circ0","kTR","gamma","Emax","EC50","CL","Q","Vc","Vp"
                          ,"sigma","dose","duration"
                          ,"Cc.0","Cp.0","Prol.0","Tx.1.0","Tx.2.0","Tx.3.0","Circ.0")
           , partrans=parameter_trans(
               toEst = Csnippet(pkpd.txform)
              ,fromEst = Csnippet(pkpd.txback)
             )
)

## ----JustForTesting, echo=FALSE, eval=FALSE---------------------------------------------
#  id <- 7
#  traj <- trajectory(pkpd,
#                     params=c(pop[pop$id==id, -which(names(pop) %in% c('id','MTT'))]
#                              , sigma=0.05
#                              , dose=50
#                              , duration=1
#                              , Cc.0 = 0.0
#                              , Cp.0 = 0.0
#                              , Prol.0 = Circ0
#                              , Tx.1.0 = Circ0
#                              , Tx.2.0 = Circ0
#                              , Tx.3.0 = Circ0
#                              , Circ.0 = Circ0)) |> as.data.frame()

## ----FirstDose, tidy=FALSE--------------------------------------------------------------
# initial conditions and output times
dose1 <- 100 # mg
Tinfusion <- 1 # 1-hour infusion
allout <- data.frame() # accumulator for nrow(pop) individual ODE solutions
parms <- function(id, dose, duration){
  parms <- unlist(pop[id,c('Circ0','kTR','gamma','Emax','EC50','CL','Q','Vc','Vp')])
  parms['sigma'] <- 0.05
  parms['dose'] <- dose
  parms['duration'] <- duration
  parms['Cc.0'] <- parms$Cp.0 <- 0.0
  parms[c('Prol.0','Tx.1.0','Tx.2.0','Tx.3.0','Circ.0')] <- parms$Circ0
  parms
}
for (id in 1:nrow(pop)) {
  Circ0 <- pop$Circ0[id]
  out <- trajectory(pkpd,
                    params=c(pop[pop$id==id, -which(names(pop) %in% c('id','MTT'))]
                             , sigma=0.05
                             , dose=dose1
                             , duration=Tinfusion
                             , Cc.0 = 0.0
                             , Cp.0 = 0.0
                             , Prol.0 = Circ0
                             , Tx.1.0 = Circ0
                             , Tx.2.0 = Circ0
                             , Tx.3.0 = Circ0
                             , Circ.0 = Circ0)) |> as.data.frame()
  out <- out[,c('time','Cc','Cp','Prol','Tx.1','Tx.2','Tx.3','Circ')]
  out$id <- paste("id",id,sep="")
  allout <- rbind(allout, out)
}

library(data.table)
allout <- as.data.table(allout)

## When a function y(x) is sampled around a minimum at equispaced (x-dx, x, x+dx)
## with corresponding values (y-dy1, y, y+dy2) such that dy1 < 0 < dy2 (i.e., with
## the middle value y being lowest value), then a quadratic interpolation yields
## estimated minimum at (x_, y_) given by:
##   x_ = x - dx/2 * (dy1 + dy2)/(dy2 - dy1)
##   y_ = y - 1/8 [(dy1 + dy2)/(dy2 - dy1)]^2.

for(.id in unique(allout$id)) {
  with(subset(allout, id == .id), {
    nadirIdx <- which.min(Circ)
    Dy1Dy2 <- diff(Circ[nadirIdx + (-1:1)])
    SdyDdy <- sum(Dy1Dy2)/diff(Dy1Dy2)
    allout[id == .id, CircMin := Circ[nadirIdx] - (1/8)*SdyDdy^2]
    allout[id == .id, tNadir := time[nadirIdx] - 0.5*diff(time[nadirIdx+(0:1)])*SdyDdy]
  })
}

allout <- upData(allout
      ,id = ordered(id, levels=paste("id",1:N,sep="")) # Note that we may
      ,units=c(Prol="cells/mm^3" # treat the non-circulating compartments
              ,Tx.1="cells/mm^3" # on a 'circulating-equivalent basis';
              ,Tx.2="cells/mm^3" # thus we attach 'cells/mm^3' units to
              ,Tx.3="cells/mm^3" # these quantities as well.
              ,Circ="cells/mm^3"
              ,Cc="ng/mL"
              ,Cp="ng/mL"
              ,time="hours")
      ,print=FALSE
       )

## ----PKplot, fig.cap="Two-compartment pharmacokinetics of the drug."--------------------
cout <- gather(allout, key="Series", value="Concentration"
               , Cc, Cp
               , factor_key = TRUE)

label(cout$Concentration) <- "Drug Concentration"

xYplot(Concentration ~ time | id, group=Series
       , data=cout, subset=time<6
       , layout=c(5,5)
       , type='l', as.table=TRUE
       , label.curves=FALSE
       , par.settings=list(superpose.line=list(lwd=2,col=brewer.pal(4,"PRGn")[c(1,4)]))
       , scales=list(y=list(log=TRUE, lim=c(10^-3,10^1)))
       , auto.key=list(points=FALSE, lines=TRUE, columns=2))

## ----MyelosuppressionPlot, fig.cap=paste("Myelosuppression in the 5-compartment model of Friberg *et al.*, with an infused dose of", dose1, "mg. **Prol:** proliferative compartment; **Tx.n:** transit compartments; **Circ:** circulating cells.")----
mout <- gather(allout, key="Series", value="ANC"
               , Prol, Tx.1, Tx.2, Tx.3, Circ
               , factor_key = TRUE)

mout <- upData(mout
               , time = time/24
               , units = c(time="days")
               , print = FALSE)

xYplot(ANC ~ time | id, group=Series
       , data=mout
       , layout=c(5,5)
       , type='l', as.table=TRUE
       , label.curves=FALSE
       , par.settings=list(superpose.line=list(lwd=2,col=brewer.pal(11,"RdYlBu")[c(1,3,4,8,10)]))
       , scales=list(y=list(log=TRUE, lim=c(100,15000), at=c(200, 500, 1000, 2000, 5000, 10000)))
       , auto.key=list(points=FALSE, lines=TRUE, columns=5))

## ----defseq, echo=FALSE-----------------------------------------------------------------
# Define a seq.function method supporting custom-scaled plot axes.
seq.function <- function(scalefun, from, to, length.out, digits=NULL){
  x <- seq(from=scalefun(from), to=scalefun(to), length.out=length.out)
  y <- numeric(length.out)
  for(i in seq(length(y)))
    y[i] <- uniroot(function(y) scaled(y)-x[i], c(from, to))$root
  if(!is.null(digits))
    y <- round(y, digits=digits)
  y
}

## ----ANCnadirs--------------------------------------------------------------------------
scaled <- function(dose, a=4.0) dose^(1/a)
dose.scale <- list(lim=scaled(c(40,550))
                   , at=scaled(c(50, 100, 250, 500))
                   , lab=c('50','100','250','500')) # TODO: Redo limits & labels when right dosing is found
doChemo <- function(draw.days=NULL, Tcyc=3*7*24, adapt.dosing=c('Newton'), omega=0.75) {
  # Find the ANC nadirs of all 20 IDs, checking ANCs on (integer-vector) draw.days
  # We will accumulate data about each course of treatment into this data frame.
  # TODO: Consider doing away entirely with columns Cc..Circ, since these inits are available in 'out'
  hourly <- which(abs(time(pkpd) - round(time(pkpd))) < .Machine$double.eps^0.5)
  anc.ts <- data.frame() # This will be used to collect an hourly 'Circ' time series
  course <- expand.grid(cycle=1:10, id=1:nrow(pop), Cc=0.0, Cp=0.0
                        , Prol=NA, Tx.1=NA, Tx.2=NA, Tx.3=NA, Circ=NA
                        , dose=NA, ANC.nadir=NA, time.nadir=NA, scaled.dose=NA
                        , ANC.nadir.est=NA, time.nadir.est=NA)
  for (day in draw.days) {
    newcolumn <- paste("ANC", day, sep="_d")
    course[,newcolumn] <- NA
    units(course[,newcolumn]) <- "cells/mm^3"
    label(course[,newcolumn]) <- paste("Day-",day," ANC", sep="")
  }
  course$dose <- seq(scaled, from=50, to=500, length.out=max(course$cycle), digits=0)[course$cycle]
  statevector <- c('Cc','Cp','Prol','Tx.1','Tx.2','Tx.3','Circ')
  course[,statevector[-(1:2)]] <- pop$Circ0[course$id] # Prol(0)=Tx.1(0)=Tx.2(0)=Tx.3(0)=Circ(0):=Circ0
  for (id in 1:nrow(pop)) { # outer loop over IDs permits state cycling
    params <- unlist(pop[id,c('Circ0','gamma','Emax','EC50','CL','Q','Vc','Vp','kTR')])
    params['sigma'] <- 0.05
    params['duration'] <- Tinfusion
    params[c('Cc.0','Cp.0')] <- 0.0
    params[c('Prol.0','Tx.1.0','Tx.2.0','Tx.3.0','Circ.0')] <- params['Circ0']
    for (cycle in 1:max(course$cycle)) {
      idx <- which(course$cycle==cycle & course$id==id)
      if (cycle>1) {
        lag_1 <- which(course$cycle==(cycle-1) & course$id==id)
        if (adapt.dosing=='Newton') { # Override preconfigured dose
          params[paste(statevector,'0',sep='.')] <- traj[nrow(traj),statevector] # recycle end-state
          if (cycle==2) {
            slope <- -2.0
          } else { # cycle >= 3 so we also have lag -2 to look back at
            lag_2 <- which(course$cycle==(cycle-2) & course$id==id)
            dY <- log(course[lag_1,'ANC.nadir'] / course[lag_2,'ANC.nadir'])
            dX <- scaled(course[lag_1,'dose']) - scaled(course[lag_2,'dose'])
            slope <- dY/dX
            if (slope >= 0) slope <- NA # to avoid instability from dY/dX>=0 due to hysteresis
          }
          delta.scaleddose <- ifelse(is.na(slope), 0, log(500 / course[lag_1,'ANC.nadir']) / slope)
          # For safety's sake, we (asymmetrically) apply relaxation factor 'omega' to any proposed dose increase:
          delta.safer <- ifelse(delta.scaleddose > 0
                                , omega*delta.scaleddose
                                , delta.scaleddose)
          new.scaleddose <- scaled(course[lag_1,'dose']) + delta.safer
          course$dose[idx] <- uniroot(function(y) scaled(y)-new.scaleddose, c(0,100000))$root
        } 
      }
      params['dose'] <- course$dose[idx]
      traj <- trajectory(pkpd, params=params) |> as.data.frame()
      to.add <- data.frame(id=rep(id,length(hourly))
                           , time=traj$time[hourly]+(cycle-1)*Tmax
                           , ANC=traj$Circ[hourly])
      anc.ts <- rbind(anc.ts, to.add)
      course[idx,c('ANC.nadir','time.nadir')] <- traj[which.min(traj$Circ),c('Circ','time')]
      course[idx,statevector] <- traj[which.max(traj$time),statevector]
      for (day in draw.days) {
        day.idx <- which(traj$time==day*24)
        course[idx,paste("ANC", day, sep="_d")] <- traj[day.idx,'Circ']
      }
      if (length(draw.days)) {
        # TODO: Here, estimate nadir level and timing based on fitted spline
        ys <- course[idx,paste("ANC", draw.days, sep="_d")]
        fit <- spline(draw.days, log(ys), n=200)
        course[idx,'ANC.nadir.est'] <- exp(min(fit$y))
        course[idx,'time.nadir.est'] <- fit$x[which.min(fit$y)]
      }
    }
  }
  
  course <- upData(course[order(course$cycle),]
                   , id = ordered(paste("id",id,sep="")
                          ,levels=paste("id",1:N,sep=""))
                   , time.nadir = time.nadir/24
                   , scaled.dose = scaled(dose)
                   , labels=c(ANC.nadir="ANC nadir"
                              ,time.nadir="Time of ANC nadir"
                              ,dose="Drug Dose"
                              ,scaled.dose="Drug Dose (rescaled)")
                   , units=c(ANC.nadir="cells/mm^3"
                             ,time.nadir="days"
                             ,dose="mg"
                             ,scaled.dose="mg")
                   , print=FALSE
  )

  anc.ts <- upData(anc.ts
                   , id = ordered(paste("id",id,sep="")
                          ,levels=paste("id",1:N,sep=""))
                   , time = time/(24*7)
                   , labels=c(ANC="ANC")
                   , units=c(ANC="cells/mm^3"
                             ,time="weeks")
                   , print=FALSE
  )
  
  list(course=course, anc.ts=anc.ts)
} # end of function

## ----Linearize--------------------------------------------------------------------------
chemo <- doChemo(draw.days=c(5,6,7,8,11), adapt.dosing=FALSE)
course <- chemo$course
anc.ts <- chemo$anc.ts

## ----ANCnadirsPlot, fig.cap="ANC nadir vs cytotoxic dose. Note the dose axis is scaled to achieve near-linearity."----
xYplot(ANC.nadir ~ scaled.dose | id, data=course, as.table=TRUE
      , layout=c(5,5)
      , scales=list(y=list(log=TRUE, lim=c(100,10000), at=c(200, 500, 1000, 2000, 5000)),
                    x=dose.scale)
)

## ----DoseSlopes, fig.cap="Slopes of $\\log(\\textit{ANC}_{nadir})$ vs $\\sqrt[4]{\\textit{dose}}$ for 25 simulated study subjects."----
slopeForID <- function(x){
  fit <- lm(log(ANC.nadir) ~ scaled.dose, data=subset(course, id==x))
  slope <- fit$coefficients['scaled.dose']
  unname(slope)
}

slope <- sapply(levels(course$id), slopeForID)
densityplot(~slope, plot.points="rug")

## ----NadirTime, fig.cap="Timing of ANC nadir vs initial cytotoxic dose. Note the dose axis is scaled to achieve near-linearity."----
xYplot(time.nadir ~ scaled.dose | id, data=course, as.table=TRUE
      , layout=c(5,5)
      , scales=list(x=dose.scale)
)

## ----NadirTimePaths, fig.cap="Timing and level of ANC nadir vs initial cytotoxic dose. The doses are equally spaced on a power-law scale that yields a roughly linear parametrization of the nadir coordinate paths."----
##, subset=(dose %in% c(500,1000,1500,2500,4000))
xYplot(ANC.nadir ~ time.nadir | id, group=dose, data=course
      , as.table=TRUE
      , layout=c(5,5)
      , auto.key=list(columns=5, title="Dose (mg)", lines=FALSE)
      , par.settings=list(superpose.symbol=list(col=gray.colors(10, start=0.7, end=0.0)))
      , scales=list(y=list(log=TRUE, lim=c(10,6000), equispaced.log=FALSE))
       )


## ----Nadir_vs_Day7, fig.cap="\\label{fig:Nadir_vs_Day7}Day-7 ANC does not serve as suitable proxy for ANC nadir across the whole modeled population. Note that for some individuals (e.g., id3, id11, id20), the day-7 ANC may dangerously underestimate the actual nadir."----
xYplot(ANC.nadir ~ ANC_d7, data=course, group=id, type='l', as.table=TRUE
       , aspect=1
       , label.curves=list(method="offset")
       , scales=list(x=list(log=TRUE, lim=c(40, 6000), equispaced.log=FALSE),
                     y=list(log=TRUE, lim=c(40, 6000), equispaced.log=FALSE))
       , par.settings=list(superpose.line=list(col="black"))
)

## ----Newtonize, fig.cap="Dose titration by Newton's method, to a goal ANC nadir of 500 cells/mm$^3$."----
chemo <- doChemo(adapt.dosing='Newton')
newton <- chemo$course
new.ts <- chemo$anc.ts
anc.tics <- c(200,500,1500,4000,10000)
right <- xYplot(ANC ~ time | id, data=new.ts
                , as.table=TRUE, type="l"
                , layout=c(5,5)
                , scales=list(y=list(log=TRUE, lim=c(100,1.5e4)
                                    , at=anc.tics
                                    , lab=as.character(anc.tics)),
                              x=list(at=seq(0,30,3)))
                )
newton <- upData(newton
                 , time = 3*(cycle-1)
                 , labels = c(time="Time")
                 , units = c(time="weeks")
                 , print = FALSE)
dose.tics <- c(50, 200, 600, 1500, 3000)
left <- xYplot(scaled.dose ~ time | id, data=newton
               , as.table=TRUE, type='p', pch='+', cex=1.5
               , layout=c(5,5)
               , scales=list(y=list(lim=scaled(c(30,3200))
                                    , at=scaled(dose.tics)
                                    , lab=as.character(dose.tics)),
                             x=list(lim=c(-1,31)
                                    , at=seq(0,30,3)
                                    , lab=c('0','','6','','12','','18','','24','','30')))
               )
update(doubleYScale(left, right, add.ylab2=TRUE)
       , par.settings = simpleTheme(col=brewer.pal(4,"PRGn")[c(4,1)])
       )

## ----SessionInfo, echo=TRUE-------------------------------------------------------------
sessionInfo()