Nothing
Code and figures for 'DTAT should supersede MTD' v1
In DTAT: Dose Titration Algorithm Tuning
.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 central and peripheral 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).
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=1>]
Tx2 [label=2>]
Tx3 [label=3>]
Circ [label="Circ"]
Exit [style=invis]
Prol -> Tx1 -> Tx2 -> Tx3 -> Circ [label=< ktr>]
Prol -> Prol [label=< kprol=ktr (Circ0/Circ(t))γ>]
Circ -> Exit [label=< kcirc≡ktr>]
}
\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
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:
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)
)
)
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:
# 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
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
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))
# 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}}$.
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.
chemo <- doChemo(draw.days=c(5,6,7,8,11), adapt.dosing=FALSE)
course <- chemo$course
anc.ts <- chemo$anc.ts
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_{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
xYplot(time.nadir ~ scaled.dose | id, data=course, as.table=TRUE
, layout=c(5,5)
, scales=list(x=dose.scale)
)
\clearpage
##, 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.
```rDay-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
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:
sessionInfo()
Try the DTAT package in your browser
Any scripts or data that you put into this service are public.
DTAT documentation built on March 31, 2023, 10:26 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
.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 central and peripheral 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). 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=
\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 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:
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) ) )
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:
# 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
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
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))
# 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}}$.
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.
chemo <- doChemo(draw.days=c(5,6,7,8,11), adapt.dosing=FALSE) course <- chemo$course anc.ts <- chemo$anc.ts
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_{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
xYplot(time.nadir ~ scaled.dose | id, data=course, as.table=TRUE
, layout=c(5,5)
, scales=list(x=dose.scale)
)
\clearpage
##, 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.
```rDay-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 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:
sessionInfo()
Try the DTAT package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.