---
title: "Code and figures for 'DTAT should supersede MTD' v1"
author: "David C. Norris <david@precisionmethods.guru>"
date: "March 2023"
output:
  pdf_document:
    keep_tex: true
    fig_width: 6.5
    fig_height: 6
vignette: <
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Code and figures for 'DTAT should supersede MTD' v1}
header-includes:
- \usepackage{setspace,relsize}
- \usepackage{tikz}
#- \usepackage[autosize,dot]{dot2texi}
---

```{r 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)
```

The following code reproduces all analyses and figures presented in article `Dose Titration Algorithm Tuning (DTAT) should supersede the Maximum Tolerated Dose (MTD) concept in oncology dose-finding trials' (v1) submitted to *F1000Research*. Several analyses and figures excluded from the article are included here; numbering of figures does not necessarily correspond to that in the article.

## Pharmacokinetic model

\newcommand{\Vc}{\mathrm{V_c}}
\newcommand{\Vp}{\mathrm{V_p}}
\newcommand{\Q}{\mathrm{Q}}
\newcommand{\CL}{\mathrm{CL}}

The simulations presented in the article are based on a notional cytotoxic drug, modeled after docetaxel. The pharmacokinetics are supposed to follow a standard 2-compartment pharmacokinetic model with *c*entral and *p*eripheral compartments having volumes $\Vc$ and $\Vp$, respectively, and drug concentrations $C_c(t)$ and $C_p(t)$. Denoting the intercompartmental and peripheral-compartment clearances $\Q$ and $\CL$, respectively, and (time-dependent) cumulative dose by $D(t)$, this model is a system of two ordinary differential equations (ODEs):

\begin{align}
\dot{C_c} &= \frac{\dot{D}}{\Vc} - \frac{\mathrm{CL}}{\Vc} C_c - \frac{\Q}{\Vc} ( C_c - C_p )
\\
\dot{C_p} &=  \frac{\Q}{\Vp} ( C_c - C_p ).
\end{align}

## Myelosuppression model

Chemotherapy-induced neutropenia (CIN) is supposed to occur according to the semimechanistic model of Friberg *et al.* (2002).
<!-- See Figure \ref{fig:Friberg-diagram}. -->
The model may be expressed by the following system of ODEs:

\newcommand{\kTR}{k_\textit{tr}}
\newcommand{\Prol}{\textit{Prol}}
\newcommand{\Tx}{\textit{Tx}}
\newcommand{\Circ}{\textit{Circ}}

\begin{align}
  \dot{\Prol} &= \kTR \cdot \Prol \cdot (1 - E_\textit{drug}) \left(\frac{\Circ_0}{\Circ}\right)^\gamma - \kTR \cdot \Prol
  \\
  \dot{\textit{Tx}_1} &= \kTR \, (\Prol - \textit{Tx}_1)
  \\
  \dot{\textit{Tx}_2} &= \kTR \, (\textit{Tx}_1 - \textit{Tx}_2)
  \\
  \dot{\Tx_3} &= \kTR \, (\Tx_2 - \Tx_3)
  \\
  \dot{\Circ} &= \kTR \, (\Tx_3 - \Circ)
\end{align}

```{dot Friberg-diagram, eval=FALSE, echo=FALSE, fig.cap="\\label{fig:Friberg-diagram}Semimechanistic model of chemotherapy-induced myelosuppression, due to Friberg *et al.* (2002). ***Prol:*** proliferative compartment; $\\boldsymbol{Transit_n}$: transit compartments; ***Circ:*** circulating cells; $\\boldsymbol{k_{tr}}$: transition rate constant; $\\boldsymbol{k_{prol}}$: proliferation rate constant; $\\boldsymbol{\\gamma > 0}$: feedback exponent; $\\boldsymbol{Circ_0}$: circulating neutrophil concentration at baseline."}
digraph g {
  node [fontname="Helvetica"]
  Prol [label="Prol"]
  Tx1 [label=<Transit<SUB>1</SUB>>]
  Tx2 [label=<Transit<SUB>2</SUB>>]
  Tx3 [label=<Transit<SUB>3</SUB>>]
  Circ [label="Circ"]
  Exit [style=invis]
  Prol -> Tx1 -> Tx2 -> Tx3 -> Circ [label=< k<SUB>tr</SUB>>]
  Prol -> Prol [label=< k<SUB>prol</SUB>=k<SUB>tr</SUB> (Circ0/Circ(t))<SUP>&gamma;</SUP>>]
  Circ -> Exit [label=< k<SUB>circ</SUB>&equiv;k<SUB>tr</SUB>>]
}
```

\clearpage

## Simulated population (inter-individual heterogeneity)

A population of 25 individuals is simulated, using population pharmacokinetic parameters reported for docetaxel in Onoue *et al.* (2016) and parameters estimated in Friberg *et al.* (2002) for their myelosuppression model.

```{r 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="")
```

\clearpage

## PK/PD simulation model

To simulate the pharmacokinetics and (myelosuppressive) pharmacodynamics of our notional cytotoxic drug, we define a **pomp** model as follows:

```{r 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)
             )
)
```

```{r 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()
```

In each subject, we now infuse the same initial dose over 1 hour, and integrate the PK and myelosuppression ODEs to obtain time series for plotting:

```{r 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
       )
```

\clearpage

```{r 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))
```

\clearpage

```{r 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))
```

```{r 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
}
```

\clearpage

## Adaptive dosing based on a simple Newton-Raphson heuristic

The `doChemo` function defined below implements multiple, 3-week courses of chemotherapy with optional adaptive dosing based on the heuristic of Newton-Raphson root-finding. When `adapt.dosing==FALSE`, the infusion doses are stepped from 50 up to 500 mg, in increments that are evenly spaced on the scale of $\sqrt[4]{\textit{dose}}$. 

```{r 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 dynamics

We now demonstrate graphically an approximate linearization of ANC nadir level and timing, under logarithmic transformation of ANC and fourth-root transformation of dose.

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

```{r 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)
)
```

```{r 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")
```

\clearpage

```{r 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)
)
```

\clearpage

```{r 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))
       )

```

\clearpage

## Would day-7 ANC suffice as a proxy for nadir?

Evidently not; see Figure \ref{fig:Nadir_vs_Day7}. Thus, CIN monitoring constitutes a nontrivial aspect of any practical implementation of DTAT, and requires explicit modeling in follow-on work. No doubt, the implementation must take the form of an adaptive process designed to balance (a) the burden of multiple blood draws against (b) the need for early warning of an impending severely neutropenic nadir that would indicate prophylactic administration of colony-stimulating factor.

```{r 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"))
)
```

\clearpage

## Simulated dose titration in 25 study subjects

```{r 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)])
       )
```

\clearpage

## SessionInfo

This document was produced using:
```{r SessionInfo, echo=TRUE}
sessionInfo()
```