\documentclass[article,nojss,shortnames]{jss}
%\VignetteIndexEntry{Using dcemriS4 for the analysis of DCE-MRI and DWI}
%\VignettePackage{dcemriS4}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% declarations for jss.cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\usepackage{amsmath}
%\definecolor{BrickRed}{rgb}{1,0,0}

\newcommand{\ktrans}{K^\text{trans}}
\newcommand{\iaugc}[1]{\text{IAUGC}_{#1}}
\newcommand{\kep}{k_\text{ep}}
\newcommand{\kepi}{k_\text{ep,i}}
\newcommand{\vp}{v_p}
\newcommand{\ve}{v_e}
\newcommand{\TR}{\mathop{\text{TR}}\nolimits}
\newcommand{\bh}{\mathbf{h}}

%% almost as usual
\author{Brandon Whitcher\\Mango Solutions \And 
        Volker J. Schmid\\Ludwig-Maximilians Universit\"at M\"unchen}
\title{Quantitative Analysis of Dynamic Contrast-Enhanced and Diffusion-Weighted Magnetic Resonance Imaging for Oncology in \proglang{R}}

%% for pretty printing and a nice hypersummary also set:
\Plainauthor{Brandon Whitcher, Volker Schmid} %% comma-separated
\Plaintitle{Quantitative Analysis of Dynamic Contrast-Enhanced and Diffusion-Weighted Magnetic Resonance Imaging for Oncology in R} %% without formatting
\Shorttitle{\pkg{dcemriS4}: Analysis of DCE-MRI and DW-MRI in \proglang{R}} %% a short title (if necessary)

%% an abstract and keywords
\Abstract{

The package \pkg{dcemriS4} provides a complete set of data analysis
tools for the quantitative assessment of dynamic contrast-enhanced
MRI.
% and diffusion-weighted MRI.  
Image processing is provided for the ANALYZE and NIfTI data formats as
inputs with all parameter estimates being output in NIfTI format.
Estimation of T1 relaxation from multiple flip-angle acquisitions,
using either constant or spatially-varying flip angles, is performed
via nonlinear regression.  Both literature-based and data-driven
arterial input functions are available and may be combined with a
variety of compartmental models.  Kinetic parameters are obtained from
nonlinear regression, Bayesian estimation via Markov Chain Monte Carlo
or Bayesian \emph{maximum a posteriori} estimation.  A non-parametric
model, using penalized splines, is also available to characterize the
contrast agent concentration time curves.  Methods for quantification
of diffusion-weighted MRI acquisitions common in oncology are also
provided.  Given the size of multi-dimensional data sets commonly
acquired in imaging studies, care has been taken to maximize
computational efficiency and minimize memory usage.  All methods are
illustrated using both simulated and real-world medical imaging data
available in the public domain.

}

\Keywords{contrast, \pkg{dcemriS4}, diffusion, dynamic, enhanced,
  imaging, magnetic, resonance}

\Plainkeywords{contrast, dcemriS4, diffusion, dynamic, enhanced,
  imaging, magnetic, resonance} %% without formatting
%% at least one keyword must be supplied

%% publication information
%% NOTE: Typically, this can be left commented and will be filled out by the technical editor
%% \Volume{13}
%% \Issue{9}
%% \Month{September}
%% \Year{2004}
%% \Submitdate{2004-09-29}
%% \Acceptdate{2004-09-29}

%% The address of (at least) one author should be given
%% in the following format:
\Address{
  Brandon Whitcher\\
  Mango Solutions\\
  Office 202, Second Floor\\
  14 Greville Street\\
  London EC1N 8SB, United Kingdom\\
  E-mail: \email{bwhitcher@mango-solutions.com}\\
  URL: \url{http://www2.imperial.ac.uk/~bwhitche}, \url{http://www.dcemri.org}\\
  
  Volker J. Schmid\\
  Bioimaging group\\
  Department of Statistics\\
  Ludwig-Maximilians Universit\"at M\"unchen, Germany\\
  E-mail: \email{volker.schmid@lmu.de}\\
  URL: \url{http://volkerschmid.de}
}
%% It is also possible to add a telephone and fax number
%% before the e-mail in the following format:
%% Telephone: +43/1/31336-5053
%% Fax: +43/1/31336-734

%% for those who use Sweave please include the following line (with % symbols):
%% need no \usepackage{Sweave.sty}

%% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

<<preliminaries,echo=FALSE,results=hide>>=
require("oro.dicom")
require("oro.nifti")
require("dcemriS4")
require("bitops")
require("minpack.lm")
require("splines")
options(prompt = "R> ", continue = "+  ", width = 70, useFancyQuotes = FALSE)
@

\begin{document}

%% include your article here, just as usual
%% Note that you should use the \pkg{}, \proglang{} and \code{} commands.

\section{Introduction}
%% Note: If there is markup in \(sub)section, then it has to be escape as above.

Quantitative analysis of tissue perfusion using dynamic
contrast-enhanced magnetic resonance imaging (DCE-MRI) is achieved
through a series of processing steps, starting with the raw data
acquired from the MRI scanner, and involves a combination of physics,
mathematics, engineering and statistics to produce a set of
statistical images based on parameter estimates from a compartmental
model.  The purpose of the \pkg{dcemriS4} package is to provide the
user with a collection of functions and subroutines that move
experimental data through all steps of this data analysis pipeline
\citep{whi-sch:JSS}, using standard data formats (such as ANALYZE or
NIfTI) that may be visualized and manipulated in \proglang{R}
\citep{R} or exported for accessibility across a wide variety of
medical image analysis software packages.

The \proglang{S}4 designation in \pkg{dcemriS4} means that
\proglang{S}4 object classes are used throughout to ensure efficient
and transparent manipulation of ANALYZE or NIfTI data structures.
Data input and output rely upon the \pkg{oro.nifti} package
\citep{oro.nifti,neuroimageing_with_R}.  Parameter estimates in
\pkg{dcemriS4} inherit attributes from the incoming ANALYZE/NIfTI
objects in order to preserve anatomical and physiological information
for appropriate visualization in \proglang{R} or other software
packages.  All functions for parameter estimation may also be applied
to aggregated data; i.e., a mean curve across an anatomical region of
interest.  However, the methodology presented here is intended to be
applied on a voxel-by-voxel basis to the ANALYZE or NIfTI objects and
all statistical summaries are output as valid NIfTI objects to
facilitate interoperability.  As voxel-wise quantitative analysis can
be time consuming, \pkg{dcemriS4} supports basic parallel computing by
incorporating the \pkg{parallel} package.  Computations are easily
parallelized with the variable \code{multicore=TRUE} in the most
computationally expensive functions.

\section{Dynamic contrast-enhanced magnetic resonance imaging}

While initial application of DCE-MRI was in the characterization of
blood-brain barrier (BBB) integrity
\citep{tof-ker:measurement,lar-etal:quantitation,lar-tof:measurement},
we will instead focus on the application of DCE-MRI to the
quantitative analysis of tumor perfusion in oncology.  Angiogenesis,
the creation of new capillaries from existing blood vessels, and
vasculogenesis, the \textit{de novo} generation of blood vessels, are
key biological processes that supply nutrients to tissue
\citep{col-pad:ieee}.  New drug therapies, such as antiangiogenic or
vaccine therapies, are expected to be cytostatic (i.e., inhibiting or
suppressing cellular growth) and may not produce changes in tumor size
as measured by traditional structural imaging techniques
\citep{cho-dwy-kno:functional}.  Thus, the ability to study malignant
vasculature through non-invasive MRI methods is advantageous.  We will
cover a DCE-MRI data analysis pipeline using \pkg{dcemriS4} that is
fully quantitative and produces biologically-relevant parameter
estimates from a compartmental model.  The interested reader will find
several chapters in \citet{jac-buc-par:dce-mri} and
\citet{par-pad:DCE-MRI} pertaining to the scientific background,
methodology and application of DCE-MRI.

The acquisition protocol for a DCE-MRI scan involves several steps.
Conversion of signal intensity to contrast agent concentration
requires the preliminary step of estimating the T1 relaxation value of
the tissue; e.g., using multiple flip-angles or inversion recovery
\citep{ber-kin-zho:handbook}.  For higher field strengths one may also
consider characterizing the inhomogeneity of the magnetic field in the
scanner.  The dynamic set of T1-weighted images are started
approximately $30-60\;\text{seconds}$ before a bolus injection of the
low molecular weight contrast agent -- a gadolinium chelate -- and
continue for $8-12\;\text{minutes}$ in oncology studies.  The length
of the dynamic acquisition may easily exceed one hour when
characterizing the integrity of the blood brain barrier.  Temporal
sampling rates depend on multiple factors: type of cancer, anatomy,
spatial resolution, spatial coverage, field strength of the scanner,
breathing method, etc. and must be carefully considered on a
study-by-study basis.  DCE-MRI relies on the reduction in T1
relaxation time, corresponding to positive signal enhancement, caused
by the presence of the contrast agent.  The quantitative analysis of
DCE-MRI data in \pkg{dcemriS4} consists of the following steps
\begin{enumerate}
  \item Pre-processing of the T1 signal (e.g., motion correction,
    co-registration, correction of the B1 field) is introduced in
    Sections~\ref{sec:motion} and~\ref{sec:b1},
  \item Estimation of voxel-wise contrast agent concentration time
    curves is introduced in Section~\ref{sec:t1},
  \item Determination of the arterial input function (AIF), either
    from the literature or by data-driven methods, is introduced in
    Section~\ref{sec:aif},
  \item Parameter estimation for a given compartmental model is
    introduced in Section~\ref{sec:kinetic}, and
  \item Statistical inference on kinetic parameters for differences
    between scans of a single patient or between distinct patients,
    is discussed in Section~\ref{sec:statistics}.
\end{enumerate}

\subsection{Motion correction and co-registration}
\label{sec:motion}

Basic motion correction within an acquisition, and co-registration
between acquired series, is available using the \code{ftm} function
for fast template matching \citep{lew:template-matching}.  A reference
volume must be pre-specified where a mask has been applied to remove
all voxel that should not be included in the algorithm.  Note, only
three-dimensional translations are allowed and no interpolation is
used (i.e., only whole-voxel translations) at this time.  More
sophisticated image registration methodology in \proglang{R} is being
developed in the \pkg{RNiftyReg} package \citep{RNiftyReg}.  Access to
the image registration routines from the Insight Toolkit
(\proglang{ITK}, \url{http://www.itk.org}) is highly desirable and
currently under investigation.

\subsection{B1 mapping via the saturated double-angle method}
\label{sec:b1}

For \emph{in vivo} MRI at high field ($\geq{3}\text{Tesla}$) it is
essential to consider the homogeneity of the active B1 field (B1+).
The B1 field is the magnetic field created by the radio frequency coil
passing an alternating current at the Larmor frequency.  The B1+ field
is the transverse, circularly polarized component of B1 that is
rotating in the same sense as the magnetization.  When exciting large
collections of spins, non-uniformity in B1+ results in nonuniform
treatment of spins.  This leads to a spatially varying signal and
contrast, and to difficulty in image interpretation and quantification
\citep{cun-pau-nay:saturated}.

The proposed method uses an adaptation of the double angle method
(DAM).  Such methods allow calculation of a flip-angle map, which is
an indirect measure of the B1+ field.  Two images are acquired: $I_1$
with prescribed angle $\alpha_1$ and $I_2$ with prescribed angle
$\alpha_2=2\alpha_1$.  All other signal-affecting sequence parameters
are kept constant.  For each voxel, the ratio of magnitude images
satisfies
\begin{equation}
  \frac{I_2(r)}{I_1(r)} = \frac{\sin\alpha_2(r)f_2(\text{T1},\TR)}
       {\sin\alpha_1(r)f_1(\text{T1},\TR)}
\end{equation}
where $r$ represents spatial position and $\alpha_1(r)$ and
$\alpha_2(r)$ are flip angles that vary with the spatially varying B1+
field.  If the effects of T1 and T2 relaxation can be neglected,
then the actual flip angles as a function of spatial position satisfy
\begin{equation}
  \alpha(r) = \text{arccos}\left(\left|\frac{I_2(r)}{2I_1(r)}\right|\right)
\end{equation}
A long repetition time ($\TR\geq{5{\cdot}\text{T1}}$) is typically
used with double-angle methods so that there is no T1 dependence in
either $I_1$ or $I_2$; i.e.,
$f_1(\text{T1},\TR)=f_2(\text{T1},\TR)=1$.  Instead, the proposed
method includes a magnetization-reset sequence after each data
acquisition with the goal of putting the spin population in the same
state regardless of whether the or $\alpha_2$ excitation was used for
the preceding acquisition.

\subsubsection{Example using phantom data}

Using data acquired from a T1 phantom at two flip angles,
$\alpha_1=60^\circ$ and $\alpha_2=120^\circ$, we compute the
multiplicative factor relative to the low flip angle using the
saturated double-angle method \citep{cun-pau-nay:saturated}.  Note,
repeated acquisitions (five) of each flip angle were obtained and force
the additional \code{rowMeans} step to average the results from the
function \code{doubleAngleMethod} in the code below.
  
<<doubleanglemethod>>=
f60 <- system.file(file.path("nifti", "SDAM_ep2d_60deg_26slc.nii.gz"), 
                   package="dcemriS4")
sdam60 <- readNIfTI(f60)
f120 <- system.file(file.path("nifti", "SDAM_ep2d_120deg_26slc.nii.gz"),
                    package="dcemriS4")
sdam120 <- readNIfTI(f120)
sdam.image <- rowMeans(doubleAngleMethod(sdam60, sdam120, 60), dims=3)
mask <- (rowSums(sdam60, dims=3) > 64)
@ 
<<dam+png,echo=FALSE,results=hide>>=
png("sdam.png", width=480, height=480)
@ 
<<dam+figure>>=
# A smooth version of "sdam.image"
fsmooth <- system.file(file.path("nifti", "SDAM_smooth.nii.gz"), 
                       package="dcemriS4")
SDAM <- readNIfTI(fsmooth)
overlay(sdam120, ifelse(mask, SDAM, NA), z=13, zlim.x=range(sdam120), 
        zlim.y=c(0.5,1.5), plot.type="single")
par(cex=4,col="white")
@ 
<<dam+dev.off,echo=FALSE,results=hide>>=
dev.off()
@ 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[tbp]
  \begin{center}
    \begin{tabular}{cc}
      \includegraphics*[width=0.4\textwidth]{sdam.png} &
      \includegraphics*[width=0.4\textwidth]{t1_phantom.png}\\
      \textbf{(a)} & \textbf{(b)} 
    \end{tabular}
  \end{center}
  \caption{\textbf{(a)} Estimated B1+ field (with isotropic Gaussian
    smoothing) using the saturated double-angle method.  The colors
    correspond to a multiplicative factor relative to the true flip
    angle ($60^\circ$).  \textbf{(b)} Estimated T1 relaxation rates
    for the phantom data acquisition.  The colormap covers the range
    $[0,2.5]\;\text{sec}$.}
  \label{fig:sdam+t1-phantom}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Three-dimensional isotropic smoothing should be applied before using
this field to modify flip angles associated with additional
acquisitions; e.g., in the \pkg{AnalyzeFMRI} package
\citep{AnalyzeFMRI,neuroimageing_with_R}.
Figure~\ref{fig:sdam+t1-phantom}a is the estimated B1+ field (with
isotropic Gaussian smoothing) for a gel-based phantom containing a
variety of T1 relaxation times.  The center of the phantom exhibits a
flip angle $>60^\circ$ while the flip angle rapidly becomes
$<60^\circ$ when moving away from the center in either the $x$, $y$ or
$z$ dimensions.  The function \code{overlay} is part of the
\pkg{oro.nifti} package, additional functions for the visualization of
ANALYZE/NIfTI data are \code{image} (overloaded for classes
\code{nifti} and \code{anlz}) and \code{orthographic}.
  
Assuming the smoothed version of the B1+ field has been computed (in
the \code{SDAM} object here), multiple flip-angle acquisitions can be
used to estimate the T1 relaxation rate from the subject (or phantom).
The multiplicative factor, derived from the saturated double-angle
method, is used to produce a spatially-varying flip-angle map and
input into the function \code{R1.fast}.

<<t1estimation>>=
alpha <- c(5,10,20,25,15)
nangles <- length(alpha)
fnames <- file.path("nifti", paste("fl3d_vibe-", alpha, "deg.nii.gz", sep=""))
X <- Y <- 64
Z <- 36
flip <- fangles <- array(0, c(X,Y,Z,nangles))
for (w in 1:nangles) {
  vibe <- readNIfTI(system.file(fnames[w], package="dcemriS4"))
  flip[,,1:nsli(vibe),w] <- vibe
  fangles[,,,w] <- array(alpha[w], c(X,Y,Z))
}
TR <- 4.22 / 1000 # in seconds
fanglesB1 <- fangles * array(SDAM, c(X,Y,Z,nangles))
zi <- 13
maskzi <- mask
maskzi[,,(! 1:Z %in% zi)] <- FALSE
R1 <- R1.fast(flip, maskzi, fanglesB1, TR, verbose=TRUE)
@ 
<<t1estimation+png,echo=FALSE,results=hide>>=
png("t1_phantom.png", width=480, height=480)
@ 
<<t1estimation+figure>>=
overlay(vibe, 1/R1$R10[,,1:nsli(vibe)], z=13, zlim.x=c(0,1024), 
        zlim.y=c(0,2.5), plot.type="single")
@ 
<<t1estimation+dev.off,echo=FALSE,results=hide>>=
dev.off()
@

Figure~\ref{fig:sdam+t1-phantom}b displays the quantitative T1 map for
a gel-based phantom using information from the estimated B1+ field.
The horizontal tubes embedded within the phantom cover a range of
$\text{T1}\in[350,1543]\;\text{ms}$, where shorter T1 relaxation times
are darker and longer relaxation times are brighter.  By defining
regions of interest (ROIs) in FSLView one may construct a mask that
separates voxels belonging to the 10 unique gels.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[tbp]
  \begin{center}
    \includegraphics*[width=0.55\textwidth]{boxplots.pdf}
  \end{center}
  \caption{Boxplots of the estimated T1 values for the gel-based
    phantom, grouped by user-specified regions of interest.  True
    T1 values are displayed as colored circles for each distinct ROI.}
  \label{fig:t1-phantom-boxplot}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

<<FSLmask>>=
fpmask <- system.file(file.path("nifti", "t1_phantom_mask.nii.gz"), 
                      package="dcemriS4")
t1pmask <- readNIfTI(fpmask)
pmask <- nifti(array(t1pmask[,,25], dim(t1pmask))) # repeat slice 25
@ 
<<t1estimation+boxplots,echo=FALSE,results=hide>>=
pdf(file="boxplots.pdf", width=6, height=6)
T1 <- c(0.484,0.350,1.07,0.648,0.456,1.07,0.660,1.543,1.543,0.353)
par(mfrow=c(1,1), mar=c(5,4,4,2)+.1)
boxplot(split(1/drop(R1$R10), as.factor(drop(pmask)))[-1], 
        ylim=c(0,2.5), xlab="Region of Interest", ylab="T1 (seconds)")
points(1:10, T1, col=rainbow(10), pch=16, cex=2)
dev.off()
@ 

We compare the ``true'' T1 values for each ROI with those obtained
from acquiring multiple flip angles with the application of B1+
mapping in Figure~\ref{fig:t1-phantom-boxplot}.  Boxplots summarize
the estimated T1 relaxation times, across all voxel in the ROI defined
by \code{pmask}, with the true T1 values (large circles).  The first
seven ROIs correspond to the cylinders that run around and through the
phantom, clockwise starting from approximately one o'clock.  The
eighth and ninth ROIs are taken from the main compartment in the gel
phantom; ROI~\#8 is drawn in the middle of the phantom while ROI~\#9
is drawn from the outside of the phantom.  The final ROI is taken from
the central cylinder embedded in the phantom.

\subsection{T1 relaxation rate and contrast agent concentration}
\label{sec:t1}

Estimation of the longitudinal relaxation time T1 is the first step in
converting signal intensity, obtained in the dynamic acquisition of
the DCE-MRI protocol, to contrast agent concentration
\citep{buc-par:measuring}.  Note, the longitudinal (or spin-lattice)
relaxation time T1 is the decay constant of the recovery of the $z$
component of the nuclear spin magnetization towards its thermal
equilibrium value \citep{bux:introduction}.  Multiple flip-angle
acquisitions are commonly used to estimate the intrinsic relaxation
rate maps $\{m_0,R_{10}\}$ of the tissue, where $m_0$ is the
equilibrium signal intensity and $R_{10}$ is the pre-injection
longitudinal relaxation rate.  The non-linear equation
\begin{equation}
  S(\theta) = \frac{m_0 \sin(\theta) (1 - E_{10})} {1 - \cos(\theta)
    E_{10}},\label{eqn:signal}
\end{equation}
where $E_{10}=\exp(-\TR{\cdot}R_{10})$, relates the observed signal
intensity $S(\cdot)$ with the parameters of interest when varying the
flip angle $\theta$ prior to the injection of the contrast agent.
Note, a repetition time of $\TR\approx4\;\text{ms}$ is common practice
for pulse sequences in clinical applications.  The Levenberg-Marquardt
algorithm, provided in \pkg{minpack.lm} \citep{minpack.lm}, is applied
to estimate the parameters; see the discussion by
\cite{ahe-etal:levenburg-marquardt}.  It is worthwhile to consult
known T1 values $(T_{10}=1/R_{10})$ for different tissue types (e.g.,
muscle, grey matter, white matter) to ensure the parameter estimates
obtained are sensible.

Estimation of the post-injection longitudinal relaxation rate $R_1(t)$
using the time-varying signal intensity $S(t)$ from pre- and
post-contrast acquisitions is performed via
\begin{eqnarray}
  A(t) & = & \frac{S(t) - S(0)}{m_0\sin(\theta)}\\
  B & = & \frac{1 - E_{10}}{1 - \cos(\theta) E_{10}}\\
  R_1(t) & = & -\frac{1}{\TR} \cdot \ln\left\{\frac{1 - [A(t)+B]}{1
    - \cos(\theta)[A(t) + B]}\right\},
\end{eqnarray}
where the flip angle $\theta\in[10^\circ,30^\circ]$ is common for the
dynamic acquisitions, but will depend on both the field strength of
the magnet and the anatomical region of interest.

The function \code{R1.fast} (embedded within \code{CA.fast})
rearranges the multi-dimensional structure of the multiple flip-angle
acquisitions into a single matrix, to take advantage of internal
\proglang{R} functions instead of loops, and calls \code{E10.lm} to
perform the non-linear regression using the Levenberg-Marquardt
algorithm.  If only two flip angles have been acquired it is possible
to use the function \code{CA.fast2}, where a linear approximation is
applied to estimate the parameters \citep{wan-rie-lee:optimizing}.

The final step of the conversion of the dynamic signal intensities to
contrast agent concentration, using the \code{CA.fast} function, is
performed via
\begin{equation}\label{eqn:conc}
  C_t(t) = \frac{1}{r_1}\left(\frac{1}{T_1} - \frac{1}{T_{10}}\right),
\end{equation}
where $r_1$ is the spin-lattice relaxivity constant (depends on the
gadolinium chelate and magnet field strength) and $T_{10}=1/R_{10}$ is
the spin-lattice relaxation time in the absence of contrast media
\citep{buc-par:measuring}.  For computational reasons we follow the
method of \cite{li-etal:improved}.

\subsection{Arterial input function}
\label{sec:aif}

Whereas quantitative PET (positron emission tomography) studies
routinely perform arterial cannulation on the subject in order to
characterize the arterial input function (AIF) directly, it has been
common to use literature-based AIFs in the quantitative analysis of
DCE-MRI.  Examples include
\begin{equation}\label{eqn:aif}
  C_p(t) = D \left(a_1 e^{-m_1t} + a_2 e^{-m_2t}\right),
\end{equation}
where $D=0.1\,\text{mmol/kg}$, $a_1=3.99\,\text{kg/l}$,
$a_2=4.78\,\text{kg/l}$, $m_1=0.144\,\text{min}^{-1}$ and
$m_2=0.0111\,\text{min}^{-1}$
\citep{wei-lan-mut:pharmacokinetics,tof-ker:measurement}; or
$D=1.0\,\text{mmol/kg}$, $a_1=2.4\,\text{kg/l}$,
$a_2=0.62\,\text{kg/l}$, $m_1=3.0$ and $m_2=0.016$
\citep{fri-etal:measurement}.  There has been progress in measuring
the AIF using the dynamic acquisition and fitting a parametric model
to the observed contrast agent concentration time curves.  Recent
models include a mixture of Gaussians \citep{par-etal:derived} and
sums of exponentials \citep{ort-etal:efficient}.  The \pkg{dcemriS4}
package has incorporated the sums-of-exponentials model
\begin{equation}\label{eqn:orton}
  C_p(t) = A_B t e^{-\mu_Bt} + A_G \left(e^{-\mu_Gt} + e^{-\mu_Bt}\right)
\end{equation}
\citep{ort-etal:efficient}, where the unknown parameters
$\beta=(A_B,\mu_B,A_G,\mu_G)$ are estimated using nonlinear
regression.  Using the AIF defined in \cite{buc:uncertainty}, we
illustrate fitting a parametric model to characterize observed data.
The \code{orton.exp.lm} function provides this capability using the
so-called double-exponential parametric form \code{orton.exp}
\eqref{eqn:orton}.

<<buckley.aif>>=
data("buckley")
aifparams <- with(buckley, orton.exp.lm(time.min, input))
fit.aif <- with(aifparams, 
                aif.orton.exp(buckley$time.min, AB, muB, AG, muG))
@
<<buckley.aif+figure,echo=FALSE,results=hide>>=
pdf(file="buckley_aif.pdf", width=6, height=6)
with(buckley, plot(time.min, input, type="l", lwd=2, xlab="Time (minutes)", 
                   ylab=""))
with(buckley, lines(time.min, fit.aif, lwd=2, col=2))
legend("topright", c("Simulated AIF", "Estimated AIF"), lwd=2, col=1:2, 
       bty="n")
dev.off()
@ 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[tbp]
  \begin{center}
    \includegraphics*[width=0.55\textwidth]{buckley_aif.pdf}
  \end{center}
  \caption{Simulated arterial input function (AIF) from
    \cite{buc:uncertainty} and the best parametric fit using the
    sums-of-exponentials model in \cite{ort-etal:efficient}.}
  \label{fig:fitted-aif}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Figure~\ref{fig:fitted-aif} shows both the true AIF and the best
parametric description using a least-squares fitting criterion.  It is
apparent from the figure that the sums-of-exponentials model cannot
match the underlying AIF from the simulated data.  This illustrates an
inherent deficiency in parametric models regardless of their
application -- the fact that it may not be appropriate to describe the
true process.

\subsection{Kinetic parameter estimation}
\label{sec:kinetic}

The focus in this section is fully quantitative pharmacokinetic
modeling of tissue perfusion and assumes the raw scanner data has been
converted to contrast agent concentration.  Please see
\cite{col-pad:ieee} and \cite{buc-par:measuring} for discussions on
this point.

The standard Kety model \citep{ket:blood-tissue}, a single-compartment
model, and the extended Kety model, the standard Kety model with an
extra ``vascular'' term, form the collection of basic parametric
models one can apply using \pkg{dcemriS4}.  Regardless of which
parametric model is chosen for the biological system, the contrast
agent concentration time curve at each voxel in the region of interest
(ROI) is approximated using the convolution of an AIF and the
compartmental model; i.e.,
\begin{eqnarray}
  C_t(t) &=& \ktrans \left[ C_p(t) \otimes e^{-\kep t} \right], 
  \label{eqn:kety}\\
  C_t(t) &=& \vp C_p(t) + \ktrans \left[ C_p(t) \otimes e^{-\kep t}
    \right].  \label{eqn:extendedkety}
\end{eqnarray}
The $\ktrans$ parameter represents the volume transfer constant
between the plasma and the extravascular extracellular space (EES) per
minute, and $\kep$ is the rate constant between EES and blood plasma.
The parameter $\vp$, in the so-called ``extended'' Kety model
\eqref{eqn:extendedkety}, describes the fraction of contrast agent in
the plasma, while
\begin{equation}\label{eqn:vie}
  \ve = \frac{\ktrans}{\kep}
\end{equation}
is the fraction of the contrast agent in the EES.

Parameter estimation may be performed using one of four options in
the current version of this software:
\begin{enumerate}
\item \code{dcemri.lm}: Non-linear regression using non-linear least
  squares (Levenberg-Marquardt optimization),
\item \code{dcemri.map}: Bayesian \emph{maximum a posteriori} (MAP)
  estimation (Nelder-Mead algorithm)
\item \code{dcemri.bayes}: Fully Bayesian inference using Markov chain
  Monte Carlo (MCMC) \citep{sch-etal:TMI},
\item \code{dcemri.spline}: Deconvolution via non-parametric curve
  fitting using Bayesian penalized splines (with MCMC)
  \citep{sch-etal:TMI09}.
\end{enumerate}

\subsubsection{Non-linear regression}

Least-square estimates of the kinetic parameters $(\ktrans,\kep)$,
also for $\vp$ for the extended Kety model, are provided in
\code{dcemri.lm}.  In each voxel a nonlinear regression model is
applied to the contrast agent concentration time curves.  All
convolutions between compartmental models and AIFs are evaluated
analytically to increase computational efficiency.  For example,
the convolution in \eqref{eqn:kety} with the literature-based AIF
\eqref{eqn:aif} produces a statistical model that is given by
\begin{equation}\label{eqn:analytical-convolution}
  C_t(t) = D \exp(\theta_1) \sum_{i=1}^2 \frac{a_i
    \{\exp(-m_it) - \exp[-\exp(\theta_2)t]\}}{\exp(\theta_2)-m_i} +
  \epsilon(t),
\end{equation}
where $\epsilon(t)$ is the observational error at time $t$,
$\theta_1=\log(\ktrans)$ and $\theta_2=\log(\kep)$.  The
parametrization $(\theta_1,\theta_2)$ is used instead of
$(\ktrans,\kep)$ to ensure positive values for both transfer rates.
We assume the expected value of the noise term to be zero; i.e.,
$E(\epsilon)=0$.  Inference is performed by minimizing the sum of
squares of the observational errors
$\min_\theta\{\sum_t[\epsilon(t)]^2\}$.  The parameter $\ve$ is
calculated using \eqref{eqn:vie}.

Model parameters are estimated, along with asymptotic standard errors,
using the Levenberg-Marquardt algorithm \citep{more78} in
\pkg{minpack.lm}.  Note, for the typical number of time points used in
DCE-MRI, the estimation procedure is not well-behaved asymptotically
and, thus, the asymptotic standard errors are not accurate
\citep{sch-etal:TMI}.

\subsubsection{Bayesian model}

A hierarchical Bayesian model can be described in three stages: the
data model, the process model and the prior parameters.

\begin{enumerate}

\item For the \emph{data model} we assume a signal-plus-noise model;
  such that the observed concentration of contrast agent $C_t(t)$ in a
  single voxel at time point $t$ with additive Gaussian error variance
  $\sigma^2$, is given by
  \begin{equation}\label{datamodel}
    C_t(t) \sim N\left(f(\ktrans,\kep,t),\sigma^2\right).
  \end{equation}
  This is the Bayesian analogue to the application of the
  least-squares fitting method in the non-linear regression
  approach.

\item For the \emph{process model} we use the single-compartment
  model \eqref{eqn:kety} or an extended Kety model
  \eqref{eqn:extendedkety}.  Evaluating the convolution for
  the single-compartment model produces
  \begin{equation}\label{processmodel}
    f\left(\ktrans,\kep,t\right) = D \ktrans \sum_{i=1}^2
    \frac{a_i [\exp(-m_it) - \exp(-\kep{t})]}{\kep-m_i}.
  \end{equation}
  As previously discussed, the kinetic parameters $\ktrans$ and $\kep$
  are transfer rates and must remain positive.  Gaussian priors on
  their logarithmic transforms
  \begin{eqnarray}
    \log(\ktrans) &\sim& N(a(\ktrans),b(\ktrans)), \label{priorktrans}\\
    \log(\kep) &\sim& N(a(\kep),b(\kep)), \label{priorkep}
  \end{eqnarray}
  ensure this constraint is met.  In breast tissue, for example,
  reasonable priors for both parameters should not exceed values of
  approximately $20\,\text{min}^{-1}$ \citep{sch-etal:TMI}.  Hence, we
  use parameters $a(\ktrans)=a(\kep)=0$ and $b(\ktrans)=b(\kep)=1$.
  Thus, the expected value of $\ktrans$ and $\kep$ is one, and with
  $99.86\%$ probability \emph{a priori} neither kinetic parameter will
  exceed $20\,\text{min}^{-1}$.  For scans covering other tissue
  types, the hyperparameters $a(\ktrans),b(\ktrans)$ and
  $a(\kep),b(\kep)$ may be adjusted accordingly when calling
  \code{dcemri.bayes}.  In case of the ``extended'' Kety model, a Beta
  prior with parameters $a(v_p)$ and $b(v_p)$ is used for the vascular
  fraction $v_p$, with \emph{a priori} expected value
  $E(v_p)=a(v_p)/[a(v_p)+b(v_p)]$.
  
\item For the \emph{prior parameter}, in this case the variance of the
  observational error, we apply a flat inverse Gamma prior
  \begin{equation}\label{priorsigma}
    \sigma^2 \sim IG\left(a(\sigma^2),b(\sigma^2)\right),
  \end{equation}
  with default parameters $a(\sigma^2)=1$ and $b(\sigma^2)=0.001$
  that reflect our lack of prior information.

\end{enumerate}

The three stages of the hierarchical model fully specify our \emph{a
  priori} knowledge.  To combine this with the observed data and
produce \emph{a posteriori} knowledge, we apply Bayes' theorem
\begin{equation}\label{bayestheorem}
  \pi(\bh\,|\,C_t) = \frac{\pi(\bh) ~
    \ell(C_t\,|\,\bh)}{\int\pi(\bh^{\ast}) ~
    \ell(C_t\,|\,\bh^{\ast})},
\end{equation}
where $\bh=(\ktrans,\kep,\sigma^2)$ denotes the vector of all
parameters across all voxel, $\pi(\bh)$ the product of the prior PDFs
and $\ell(C_t\,|\,\bh)$ denotes the (Gaussian) likelihood function of
$C_t(t)$ from \eqref{datamodel}.  In the Bayesian framework,
conclusions are drawn from the joint posterior PDF only.  Two
functions are provided to exploit the posterior PDF:

\begin{itemize}

\item The function \code{dcemri.map} provides voxel-wise
  \emph{maximum a posteriori} (MAP) estimators \citep{degroot70}
  using the Nelder-Mead algorithm provided in \code{optim}.  Note, the
  posterior may be multi-modal and, hence, a global optimization may
  not be appropriate and/or feasible.  No standard errors are provided
  with this method.

\item The function \code{dcemri.bayes} provides the posterior median
  as the summary statistic for $(\ktrans,\kep,\vp)$, along with the
  posterior standard deviation for all statistics, by sampling from
  the posterior using MCMC \citep{gil-ric-spi:MCMC}.  All samples from
  the joint posterior distribution may be saved using the option
  \code{samples=TRUE}, allowing one to interrogate the posterior
  probability density function (PDF) of all parameter estimates.  To
  increase computational efficiency sampling from the posterior
  distribution is implemented in \proglang{C} and linked to
  \proglang{R}.  It is useful to retain all samples from the joint
  posterior when one wants to construct, for example, voxel-wise
  credible intervals on the kinetic parameters.  The algorithm is
  computationally expensive and parallel computation has been enabled
  with the \pkg{parallel} package by setting the option
  \code{multicore=TRUE}.

\end{itemize}

\subsubsection{Bayesian penalized splines}

An alternative to parametric modeling, the function
\code{dcemri.spline} may be used to de-convolve and de-noise the
contrast agent concentration time curves using an adaptive penalized
spline approach \citep{sch-etal:TMI09}.  A Bayesian hierarchical model
is constructed

\begin{enumerate}

\item The \emph{data model} is the Gaussian observation model
  \eqref{datamodel}.

\item For the \emph{process model} a general approach is used, such
  that
  \begin{equation}
    f(t) = C_p(t) \otimes R(t),
  \end{equation}
  where $R(t)$ is the response function in the tissue.  The
  convolution is derived through the discretization of $C_p(t)$ and
  $R(t)$ (e.g., at the observation time points), allowing one to use
  the observed AIF instead of a parametric model.  The response is
  assumed to be a smooth function, modeled as an adaptive penalized
  spline
  \begin{equation}
    R(t) = \sum_{j=1}^p \beta_j B_j(t),
  \end{equation}
  where $\mathbf{B}$ is a $B$-spline design matrix.  An adaptive
  second-difference penalty is used on the spline regression
  parameters $\beta_j$; i.e.
  \begin{equation}
    \beta_j=2\beta_{j-1}-\beta_{j-2}+e_j \quad j=3,\ldots,p,
  \end{equation}
  where $e_j{\sim}N(0,\delta_i^2)$. Note, a first-order penalty is
  also available using the option \code{rw=1}.

\item The \emph{prior} for the adaptive smoothing parameter
  $\delta_i^2$ is given by
  \begin{equation}
    \delta_i^2 \sim IG(a(\delta),b(\delta))
  \end{equation}
  with default parameters $a(\delta)=b(\delta)=10^{-5}$ that
  provide nearly flat prior information, and an Inverse Gamma prior
  for the observational error \eqref{priorsigma}.

\end{enumerate}

Making use of Bayes' formula \eqref{bayestheorem}, the posterior is
assessed using an MCMC algorithm.  By default, the function
\code{dcemri.spline} returns the median of the maximum $F_p$ of the
response function $R(t)$ per voxel.  The median response function
(\code{response=TRUE}) and the fitted contrast agent concentration
time curve (\code{fitted=TRUE}) may also be provided.

An automated method for estimation of the onset time of contrast agent
(from a bolus injection) has been implemented.  The algorithm follows
these steps:

\begin{enumerate}

\item Find the minimum time $t$, for which the contrast agent
  concentration curve significantly exceeds zero,

\item Compute the gradient of $C_t$ at point $t$, exploiting the
  derivative of the $B$-spline,

\item Compute the onset time as
  \begin{equation}
    t_0 = t - \frac{C_t(t)}{dC_t(t)/dt}.
  \end{equation}
  
\end{enumerate}

To provide estimates of the kinetic parameters from a compartmental
model, a parametric model may be applied to the estimated response
curve (\code{nlr=TRUE}).  At this point in time a single exponential
model (\code{"weinmann"}) or the adiabatic approximation to the tissue
homogeneity model (\code{"AATH"}) are available.  For the AATH model,
the response function is given by
\begin{equation}
  R(t) = F_p \cdot \left\{
  \begin{array}{ll}
    E \exp\left[-\left(t-T_c\right)EF_p/\ve\right] & 
    \text{for} \quad t \geq T_c,\\
    0 & \text{for} \quad t < 0,\\
    1 & \text{for} \quad 0 \leq t < T_c,
  \end{array} \right.
\end{equation}
where $T_c$ is the transit time through the capillary, $E$ is the
extraction fraction and $F_p$ is the mean plasma flow. These
parameters may be re-expressed as the kinetic parameters from the
(extended) Kety model via
\begin{eqnarray}
  \ktrans &=& EF_p,\\
  \kep &=& EF_p/\ve,\\
  \vp &=& T_cF_p.
\end{eqnarray}
The response model is applied to each sample of the estimated response
curve, and the median and standard error of the kinetic parameters are
provided.  Samples from the posterior density for the kinetic
parameters, the maximum response $F_p$, the onset time $t_0$, the
response functions, and the fitted curves are also available
(\code{samples=TRUE}).  Parallel computing may be accessed using the
\pkg{parallel} package (\code{multicore=TRUE}).

\subsubsection{Estimating the kinetic curves}

Using kinetic parameters estimated with one of the methods above, the
functions \code{kineticModel} or \code{orton.exp.lm} may be used to
compute the estimated contrast agent concentration time curves for the
given parametric models.  A list of arrays or \code{nifti} class
objects of kinetic parameters can be given to \code{kineticModel} to
produce voxel-wise estimates of the compartmental model.

\subsection{Statistical inference}
\label{sec:statistics}

No specific support is provided for hypothesis testing in
\pkg{dcemriS4}.  We recommend one uses built-in functions in
\proglang{R} to perform ANOVA (analysis of variance) or mixed-effects
models based on statistical summaries of the kinetic parameters over a
given ROI on a per subject per visit basis.  An alternative to this
traditional approach is to analyze an entire study using a Bayesian
hierarchical model \citep{sch-etal:hierarchical}, where an
implementation is under development in the software project
\pkg{PILFER} (\url{http://pilfer.sourceforge.net}).  One may also
question the rationale for hypothesis testing in only one kinetic
parameter.  Preliminary work has been performed in looking at the
joint response to treatment of both $\ktrans$ and $\kep$ in DCE-MRI
using functional data analysis \citep{oco-etal:fPCA}.

\section{Diffusion weighted imaging}

Diffusion weighted imaging (DWI) is an imaging bio\-marker that is
rapidly becoming popular and widely applied in oncology
\citep{che-etal:diffusion,koh-pad:DW-MRI}.  DWI allows one to quantify
the diffusion behavior of water by estimating the apparent diffusion
coefficient (ADC) \citep{whe-etal:D}.  That is, assuming completely
unrestricted motion of water, how is the motion of water impeded by
the biological structure of tissue?  The reduction in the ability of
water to diffuse in tissue has been used to infer
biologically-relevant information in oncology; e.g., in tumor
detection, disease progression and the evaluation of treatment
response.

DWI is an MR technique that provides a unique insight into tissue
structure through MRI diffusion measurements \emph{in vivo}
\citep{mos-etal:diffusion,whe-etal:D}.  These diffusion measurements
reflect the effective displacement of water molecules allowed to
migrate for a given time \citep{leb-del-lev:temperature}.  Using the
Stejskal--Tanner equation
\begin{equation}
  S = S_0 \exp\left(-\gamma^2G^2\delta^2(\Delta-\delta/3)D\right) =
  S_0 \exp\left(-bD\right),
\end{equation}
one may solve for the unknown diffusion to obtain the \emph{apparent
  diffusion coefficient} (ADC) $D$ \citep{whe-etal:D}.  For
completeness, $S_0$ is the (unknown) signal intensity without the
diffusion weighting, $S$ is the observed signal with the gradient
applied, $\gamma$ is the gyromagnetic ratio, $G$ is the strength of
the gradient pulse, $\delta$ is the duration of the gradient pulse and
$\Delta$ is the time between the two pulses.  The micro-parameters
$(\gamma,G,\Delta,\delta)$ are selected prior to data acquisition and
may be combined into a single parameter
$b=\gamma^2G^2\delta^2(\Delta-\delta/3)$, known as the $b$-value.  The
functions \code{ADC.fast} and \code{adc.lm} perform parameter
estimation using a similar interface to kinetic parameter estimation
previously introduced for DCE-MRI with the Levenberg-Marquardt
algorithm.

Acquisition protocols typically involve obtaining a volume without
diffusion weighting ($b=0\;\text{s/mm}^2$), at low diffusion weighting
($b\geq100\;\text{s/mm}^2$) and higher diffusion weighting
($b\geq500\;\text{s/mm}^2$).  When estimating the ADC value, one
should exclude any acquisitions with $b\leq100\;\text{s/mm}^2$ to
minimize the influence of perfusion effects
\citep{pad-etal:neoplasia}.

The diffusivity of water at room temperature without restrictions is
approximately $3.0{\times}10^{-3}$ $\text{mm}^2/\text{s}$.  Once the
ADC is estimated in the tissue of interest at baseline, treatment
response may be assessed at subsequent time points.  The most
appropriate timings depend on both the type of tumor and treatment
regime.  Observing an decrease in diffusivity, via a decrease in the
ADC values post-treatment, may be a result of cell swelling after
initial chemotherapy or radiotherapy followed by an increase in
diffusivity, via an increase in the ADC values, from cell necrosis and
lysis.  A decrease in ADC values may be observed directly through
tumor apoptosis after treatment
\citep{koh-pad:DW-MRI,pad-etal:neoplasia}.

\section{The RIDER Neuro MRI data repository}

The National Biomedical Imaging Archive (NBIA;
\url{http://cabig.nci.nih.gov/tools/NCIA}) is a searchable, national
repository integrating \emph{in vivo} cancer images with clinical and
genomic data.  The NBIA provides the scientific community with public
access to DICOM images, image markup, annotations, and rich metadata.
The DCE-MRI and DW-MRI data analyzed here were downloaded from the
``RIDER Neuro MRI'' collection
(\url{http://wiki.nci.nih.gov/display/CIP/RIDER}).

\subsection{Dynamic contrast-enhanced MRI}

Functions of the \pkg{oro.nifti} package are utilized to read the
signal intensity files, in ANALYZE or NIfTI format, obtained from the
MRI scanner (after conversion from DICOM).  In this example
pre-contrast multiple flip angle acquisitions are available for
estimation of contrast agent concentration.  We use \code{CA.fast} to
estimate the intrinsic relaxation rate $R_{10}$ and equilibrium signal
intensity $m_0$ from \eqref{eqn:signal} and the contrast agent
concentration curve $C_t(t)$ from \eqref{eqn:conc}.  In order to save
computation time and memory, we utilize a binary mask with a very
limited region-of-interest (ROI) created in FSLView
(\url{http://fsl.fmrib.ox.ac.uk/fsl/fslview/}) and saved in ANALYZE
format.

<<RiderNeuroMRI+pre0,eval=FALSE>>=
perf <- paste("281949", "19040721", "perfusion.nii.gz", sep="_")
fmask <- system.file(file.path("nifti", sub(".nii", "_mask.hdr", perf)),
                     package="dcemriS4")
mask <- readANALYZE(fmask)
mask <- ifelse(mask > 0, TRUE, FALSE)
dynamic <- readNIfTI(perf)
@ 
<<RiderNeuroMRI+pre1,eval=FALSE>>=
TR <- 4.43 / 1000 # taken from CSV file
dangle <- 25      # taken from CSV file
(fflip <- list.files(pattern="ax[0-9]"))
(fangles <- as.numeric(sub(".*ax([0-9]+).*", "\\1", fflip)))
flip <- array(NA, c(dim(mask), length(fangles)))
for (fa in 1:length(fangles)) {
  flip[,,,fa] <- readNIfTI(fflip[fa])
}
@ 
<<RiderNeuroMRI+pre2,eval=FALSE>>=
ca <- CA.fast(dynamic, mask, dangle, flip, fangles, TR)
writeNIfTI(ca$M0, paste(perf, "m0", sep="_"))
writeNIfTI(ca$R10, paste(perf, "r10", sep="_"))
writeNIfTI(ca$conc, paste(perf, "gdconc", sep="_"))
@ 

Note, we have used information in the file names to provide the flip
angles (\code{fangles}) for input into \code{CA.fast} ensuring that
the flip angles match the flip-angle array (\code{flip}).  After
estimating the contrast agent concentration time curve in each voxel
we fit a compartmental model to obtain estimates of the kinetic
parameters that describe the simplified biological process of
perfusion.  Here, the ``extended'' Kety model is used, which includes
a vascular compartment.  An AIF must be defined in order to complete
the compartmental model in \eqref{eqn:extendedkety}.  It is relatively
straightforward to estimate such an AIF from contrast agent
concentration time curves from an appropriate voxel or collection of
voxels.  However, for simplicity we select a literature-based AIF, the
sum of two exponentials with values taken from
\cite{fri-etal:measurement}, that is available for all compartmental
model fitting procedures.

\subsubsection{Non-linear regression}

A numeric optimization of the least-square criterion, using the
Levenberg-Marquardt algorithm, is provided by \code{dcemri.lm} and
illustrated below.  Note that the acquisition times for the dynamic
series are read in from a pre-existing text file, converted from
seconds to minutes (using \code{str2time} from \pkg{oro.dicom}) and
offset by the bolus injection time (at the ninth acquisition).  This
information was obtained from the original DICOM data and saved in
\code{rawtimes.txt}.  The literature-based AIF \code{fritz.hansen} is
used here for illustrative purposes.  Alternatives include estimating
values for a parametric AIF (e.g., \code{aif="orton.exp"}) from the
data and supplying them via the \code{user} option in
\code{dcemri.lm}, or providing an empirical AIF
(\code{aif="empirical"}) and passing the vector of AIF values via the
\code{user} option.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[tbp]
  \begin{center}
    \includegraphics*[width=0.9\textwidth]{281949_19040721_perfusion_ktrans.png}
  \end{center}
  \caption{Statistical images of $\ktrans$ overlayed on the dynamic
    acquisition for the RIDER Neuro MRI data.  Two potential tumors
    are visible, in the region-of-interest, by enhanced rims of high
    $\ktrans$ and central cores of low $\ktrans$.  The values of
    $\ktrans$ are $[0,0.1]\,\text{min}^{-1}$.}
  \label{fig:perfusion-ktrans}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

<<RiderNeuroMRI+lm0,eval=FALSE>>=
acqtimes <- str2time(unique(sort(scan("rawtimes.txt", quiet=TRUE))))$time
acqtimes <- (acqtimes - acqtimes[9]) / 60 # minutes
conc <- readNIfTI(paste(perf, "gdconc", sep="_"))
@ 
<<RiderNeuroMRI+lm1,eval=FALSE>>=
fit.lm <- dcemri.lm(conc, acqtimes, mask, model="extended", 
                    aif="fritz.hansen", control=nls.lm.control(maxiter=100), 
                    multicore=TRUE, verbose=TRUE)
writeNIfTI(fit.lm$ktrans, paste(perf, "ktrans", sep="_"))
@ 
<<RiderNeuroMRI+lm2,eval=FALSE,echo=FALSE>>=
writeNIfTI(fit.lm$kep, paste(perf, "kep", sep="_"))
writeNIfTI(fit.lm$vp, paste(perf, "vp", sep="_"))
writeNIfTI(fit.lm$ve, paste(perf, "ve", sep="_"))
writeNIfTI(fit.lm$sse, paste(perf, "sse", sep="_"))
rm(fit.lm)
@ 
<<RiderNeuroMRI+lm3,echo=FALSE,eval=FALSE>>=
fit.lm <- list(ktrans=readNIfTI(system.file(file.path("nifti", sub(".nii", "_ktrans.nii", perf)), package="dcemriS4")))
x <- 41:220
y <- 21:220
@ 
<<RiderNeuroMRI+lm4,echo=FALSE,eval=FALSE>>=
png(file=paste(paste(perf, "ktrans", sep="_"), "png", sep="."), 
    width=480, height=480)
@ 
<<RiderNeuroMRI+lm5,eval=FALSE>>=
overlay(dynamic, ifelse(mask, fit.lm$ktrans, NA), w=11, zlim.x=c(32,256), 
        zlim.y=c(0,0.1))
@ 
<<RiderNeuroMRI+lm6,echo=FALSE,eval=FALSE,results=hide>>=
dev.off()
fit.lm$kep <- readNIfTI(paste(perf, "kep", sep="_"))
fit.lm$vp <- readNIfTI(paste(perf, "vp", sep="_"))
fit.lm$ve <- readNIfTI(paste(perf, "ve", sep="_"))
fit.lm$sse <- readNIfTI(paste(perf, "sse", sep="_"))
zrx <- c(32,256)
png(file=paste(paste(perf, "ktrans7", sep="_"), "png", sep="."), 
    width=480, height=480)
overlay(as(dynamic[x,y,,], "nifti"), 
        ifelse(mask[x,y,], fit.lm$ktrans[x,y,], NA), z=7, w=11, 
        zlim.x=zrx, zlim.y=c(0,0.1), plot.type="single")
dev.off()
png(file=paste(paste(perf, "kep", sep="_"), "png", sep="."), 
    width=480, height=480)
overlay(as(dynamic[x,y,,], "nifti"), 
        ifelse(mask[x,y,], fit.lm$kep[x,y,], NA), z=7, w=11, 
        zlim.x=zrx, zlim.y=c(0,1.25), plot.type="single")
dev.off()
png(file=paste(paste(perf, "vp", sep="_"), "png", sep="."), 
    width=480, height=480)
overlay(as(dynamic[x,y,,], "nifti"), ifelse(mask[x,y,], fit.lm$vp[x,y,], NA),
        z=7, w=11, zlim.x=zrx, zlim.y=c(0,0.03), plot.type="single")
dev.off()
png(file=paste(paste(perf, "ve", sep="_"), "png", sep="."), 
    width=480, height=480)
overlay(as(dynamic[x,y,,], "nifti"), ifelse(mask[x,y,], fit.lm$ve[x,y,], NA),
        z=7, w=11, zlim.x=zrx, zlim.y=c(0,0.3), plot.type="single")
dev.off()
png(file=paste(paste(perf, "sse", sep="_"), "png", sep="."), 
    width=480, height=480)
overlay(as(dynamic[x,y,,], "nifti"), ifelse(mask[x,y,], fit.lm$sse[x,y,], NA),
        z=7, w=11, zlim.x=zrx, zlim.y=c(0,0.05), plot.type="single")
dev.off()
@ 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[tbp]
  \begin{center}
  \begin{tabular}{cc}
    $\kep$ & $v_e$\\
    \includegraphics*[width=0.4\textwidth]{281949_19040721_perfusion_kep.png} &
    \includegraphics*[width=0.4\textwidth]{281949_19040721_perfusion_ve.png}\\
    $v_p$ & $\text{SSE}$\\    
    \includegraphics*[width=0.4\textwidth]{281949_19040721_perfusion_vp.png} &
    \includegraphics*[width=0.4\textwidth]{281949_19040721_perfusion_sse.png}
    \end{tabular}
  \end{center}
  \caption{Statistical images of the RIDER Neuro MRI data using
    non-linear regression for slice $z=7$.  The parameter $\kep$ is
    the rate constant between EES and blood plasma (units are
    $[0,1.25]\,\text{min}^{-1}$), $v_p$ is the vascular space fraction
    of plasma (units are $[0,30]\,\text{\%}$) and $v_e$ is the EES
    space fraction (units $[0,30]\,\text{\%}$).  The sums-of-squared
    error (SSE) measures the quality of fit over the tumor
    region-of-interest (units are $[0,0.05]$).}
  \label{fig:perfusion-parametermaps}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
Figure~\ref{fig:perfusion-ktrans} shows the estimated $\ktrans$
statistical images in the pre-defined ROI overlayed on the dynamic
acquisition (for anatomical reference).  Two rings of high $\ktrans$
values are visually apparent from the statistical images, indicating
the presence of two tumors.  In both cases $\ktrans$ is nearly zero in
the center of both rings, indicating that no blood is being supplied
to the core of the tumors (most likely caused by necrotic or apoptotic
processes).  Figure~\ref{fig:perfusion-parametermaps} provides
statistical images that summarize the entire model-fitting procedure
for slice $z=7$.  For both tumors the fraction of contrast agent in the
extravascular extracellular space (EES) $v_e$ is high at the tumor
rim, a common feature in a tumor that is hypervascular compared to the
surrounding tissue.  The perfusion/permeability is vastly diminished
in the core of the tumor, exhibited by low values of $\kep$ and $v_p$.
Goodness-of-fit for the compartmental model may be assessed using the
sums-of-squared error (SSE).  The SSE over the given ROI covers a
variety of tissue types; e.g., white matter, gray matter,
cerebrospinal fluid (CSF), skull and air.  The SSE is high for tissue
types in which the compartmental model is not appropriate.  In
contrast, the SSE is nearly spatially invariant across the healthy
brain tissue and tumor.

\subsubsection{Bayesian maximum a posteriori (MAP) estimation}

Caution must be exercised when using non-linear regression algorithms,
since the Levenberg-Marquardt algorithm used in \code{dcemri.lm} is
\emph{not} guaranteed to converge and is susceptible to noise.  More
robust results may be achieved by using biologically-relevant prior
information in a Bayesian framework \citep{sch-etal:TMI}.  Two methods
for parameter estimation from a Bayesian perspective are implemented
in \pkg{dcemriS4}.  The function \code{dcemri.bayes} uses Markov Chain
Monte Carlo (MCMC) to explore the posterior PDF
\citep{gil-ric-spi:MCMC} and \code{dcemri.map} uses numerical
optimization to maximize the posterior \citep{degroot70}.

From the non-linear regression analysis of the RIDER Neuro MRI data,
it appears that approximately $\ktrans\in[0,0.1]$ and
$\kep\in[0,1.25]$ (Figure~\ref{fig:perfusion-ktrans}).  Hence, we use
a Gaussian distribution with expected value $a(\ktrans)=\log(0.05)$
and variance $b(\ktrans)=1$ on $\log(\ktrans)$, and a Gaussian
distribution with expected value $a(\kep)=\log(0.7)$ and variance
$b(\kep)=1$ on $\log(\kep)$ as priors.  For $\vp$, we use the Beta
distribution $B(a(v_p)=1,b(v_p)=19)$; i.e., the expected value is
given by $E(\vp)=0.05$.  Parameter estimation via \code{dcemri.map}
follows a consistent user interface established in \code{dcemri.lm}.

<<RiderNeuroMRI+map1,eval=FALSE>>=
fit.map <- dcemri.map(conc, acqtimes, mask, model="extended", 
                      aif="fritz.hansen", ab.ktrans=c(log(0.05),1),
                      ab.kep=c(log(0.7),1), ab.vp=c(1,19),
                      multicore=TRUE)
writeNIfTI(fit.map$ktrans, paste(perf, "ktrans", "map", sep="_"))
@ 
<<RiderNeuroMRI+map2,eval=FALSE,echo=FALSE>>=
writeNIfTI(fit.map$kep, paste(perf, "kep", "map", sep="_"))
writeNIfTI(fit.map$ve, paste(perf, "ve", "map", sep="_"))
writeNIfTI(fit.map$vp, paste(perf, "vp", "map", sep="_"))
writeNIfTI(fit.map$sigma2, paste(perf, "sigma2", "map", sep="_"))
@ 
<<RiderNeuroMRI+map3,echo=FALSE,eval=FALSE>>=
fit.map <- list(ktrans=readNIfTI(system.file(file.path("nifti", sub(".nii", "_ktrans_map.nii", perf)), package="dcemriS4")))
@ 
<<RiderNeuroMRI+map4,echo=FALSE,eval=FALSE,results=hide>>=
png(file=paste(paste(perf, "ktrans", "map", sep="_"), "png", sep="."), 
    width=480, height=480)
overlay(as(dynamic[x,y,,], "nifti"), 
        ifelse(mask[x,y,], fit.map$ktrans[x,y,], NA),
        z=7, w=11, zlim.x=zrx, zlim.y=c(0,0.1), plot.type="single")
dev.off()
pdf(file=paste(paste(perf, "ktrans", "compare", sep="_"), "pdf", sep="."), 
    width=8, height=8)
plot(fit.lm$ktrans, fit.map$ktrans, xlim=c(0,0.3), ylim=c(0,0.3),
     xlab=expression(paste(K^{trans}, " (Levenberg-Marquardt)")), 
     ylab=expression(paste(K^{trans}, " (MAP)")), 
     pch=19)
abline(0, 1, col="red", lwd=2)
dev.off()
@ 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[tbp]
  \begin{center}
    \begin{tabular}{cc}
      Levenberg-Marquardt & MAP\\
      \includegraphics*[width=0.4\textwidth]{281949_19040721_perfusion_ktrans7.png}  &
      \includegraphics*[width=0.4\textwidth]{281949_19040721_perfusion_ktrans_map.png}
    \end{tabular}
  \end{center}
  \caption{Statistical images of $\ktrans$ for the RIDER Neuro MRI
    data (slice $z=7$).  Two methods for parameter estimation are
    displayed: non-linear regression using the Levenberg-Marquardt
    algorithm (left) and \emph{maximum a posteriori} (MAP) estimation
    (right).  The values of $\ktrans$ are $[0,0.1]\,\text{min}^{-1}$
    for both images.}
  \label{fig:perfusion-map}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Figure~\ref{fig:perfusion-map} shows the $\ktrans$ statistical image
for the slice $z=7$ obtained using the Bayesian MAP estimator.
Estimated values are similar to those obtained using non-linear
regression (reproduced in Figure~\ref{fig:perfusion-map} to facilitate
a side-by-side comparison).  The estimation is similar, but subtly
different.  For example, by using a biological prior on the kinetic
parameters the number of voxels where the MAP estimator does not
converge is essentially eliminated when compared with non-linear
regression.

<<RiderNeuroMRI+lm-versus-map,echo=FALSE,eval=FALSE>>=
sum.lm <- sum(is.na(fit.lm$ktrans[mask]))
sum.map <- sum(is.na(fit.map$ktrans[mask]))
100 * c("LM"=sum.lm, "MAP"=sum.map) / sum(mask > 0)
@ 

Because we have avoided a computationally expensive procedure, the
computing times for \code{dcemri.lm} and \code{dcemri.map} are roughly
the same.  However, the MAP estimator does not make use of the
complete joint posterior PDF.

\subsubsection{Bayesian estimation via Markov Chain Monte Carlo}

Using the MCMC algorithm provided by \code{dcemri.bayes} is
computationally expensive when compared with all previous estimation
procedures.  However, the MCMC algorithm explores the complete
posterior PDF.  Statistical summaries of the marginal posteriors,
associated with all parameters of interest, are provided by default
and all samples from the joint posterior may be obtained using the
option \code{samples=TRUE} (internal memory may become an issue when
using this option).  Using samples from the joint posterior,
additional statistics may be derived from the model-fitting procedure;
e.g., the reliability of the estimated parameters using credible
intervals.  The following application of \code{dcemri.bayes} uses the
default \code{samples=FALSE} and, hence, we are restricted to
posterior medians and standard deviations for all parameters in the
compartmental model.

<<RiderNeuroMRI+bayes1,eval=FALSE>>=
fit.bayes <- dcemri.bayes(conc, acqtimes, mask, model="extended", 
                          aif="fritz.hansen", ab.ktrans=c(log(0.05),1), 
                          ab.kep=c(log(0.7),1), ab.vp=c(1,19))
writeNIfTI(fit.bayes$ktrans, paste(perf, "ktrans", "bayes", sep="_"))
@ 
<<RiderNeuroMRI+bayes2,echo=FALSE,eval=FALSE>>=
writeNIfTI(fit.bayes$ktranserror, paste(perf, "ktrans", "bayes", "sd", sep="_"))
writeNIfTI(fit.bayes$kep, paste(perf, "kep", "bayes", sep="_"))
writeNIfTI(fit.bayes$keperror, paste(perf, "kep","bayes", "sd", sep="_"))
writeNIfTI(fit.bayes$ve, paste(perf, "ve", "bayes", sep="_"))
writeNIfTI(fit.bayes$vp, paste(perf, "vp", "bayes", sep="_"))
writeNIfTI(fit.bayes$vperror, paste(perf, "vp","bayes", "sd", sep="_"))
rm(fit.bayes)
@ 
<<RiderNeuroMRI+bayes3,echo=FALSE,eval=FALSE,results=hide>>=
fit.bayes <- list(ktrans=readNIfTI(paste(perf, "ktrans","bayes", sep="_")))
png(file=paste(paste(perf, "ktrans", "bayes", sep="_"), "png", sep="."), 
    width=480, height=480)
overlay(as(dynamic[x,y,,], "nifti"), 
        ifelse(mask[x,y,], fit.bayes$ktrans[x,y,], NA), z=7, w=11, 
        zlim.x=zrx, zlim.y=c(0,0.1), plot.type="single")
dev.off()
fit.bayes$ktranserror <- readNIfTI(paste(perf, "ktrans","bayes","sd", sep="_"))
png(file=paste(paste(perf, "ktrans", "bayes", "sd", sep="_"), "png", sep="."), 
    width=480, height=480)
overlay(as(dynamic[x,y,,], "nifti"), 
        ifelse(mask[x,y,], fit.bayes$ktranserror[x,y,], NA), z=7, w=11, 
        zlim.x=zrx, zlim.y=c(0,0.0075), plot.type="single")
dev.off()
png(file=paste(paste(perf, "ktrans", "bayes", "cv", sep="_"), "png", sep="."), 
    width=480, height=480)
overlay(as(dynamic[x,y,,], "nifti"), 
        ifelse(mask[x,y,], (fit.bayes$ktranserror/fit.bayes$ktrans)[x,y,], NA),
        z=7, w=11, zlim.x=zrx, zlim.y=c(0,0.2), plot.type="single")
dev.off()
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[tbp]
  \begin{center}
    \begin{tabular}{cc}
      MAP & $\text{median}(\ktrans)$\\
      \includegraphics*[width=0.4\textwidth]{281949_19040721_perfusion_ktrans_map.png}  &
      \includegraphics*[width=0.4\textwidth]{281949_19040721_perfusion_ktrans_bayes.png}\\
      $\text{sd}(\ktrans)$ & $\text{sd}(\ktrans)/\text{median}(\ktrans)$\\
      \includegraphics*[width=0.4\textwidth]{281949_19040721_perfusion_ktrans_bayes_sd.png} &
      \includegraphics*[width=0.4\textwidth]{281949_19040721_perfusion_ktrans_bayes_cv.png}
    \end{tabular}
  \end{center}
  \caption{Statistical images of $\ktrans$ for the RIDER Neuro MRI
    data (slice $z=7$) using Bayesian estimation.  The posterior
    median $\ktrans$, posterior standard deviation of $\ktrans$ and
    the ratio of these two statistics have been obtained from a fully
    Bayesian MCMC algorithm.  The MAP estimator of $\ktrans$ is also
    displayed for comparison.  The units of $\ktrans$ are
    $[0,0.1]\,\text{min}^{-1}$ for both estimation techniques.}
  \label{fig:perfusion-bayes}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Figure~\ref{fig:perfusion-bayes} displays statistical summaries of
$\ktrans$ (posterior median and standard deviation) provided the
default settings of \code{dcemri.bayes}.  It is clear that the
posterior median differs from the MAP estimator (reproduced in
Figure~\ref{fig:perfusion-bayes} for a side-by-side comparison) across
the majority of non-tumor voxel in the ROI.  However, $\ktrans$ values
around the enhancing rim of the tumor are similar across all three
methods: \code{dcemri.lm}, \code{dcemri.map} and \code{dcemri.bayes}.
Figure~\ref{fig:perfusion-bayes} also provides the posterior standard
deviation of $\ktrans$ and an image of the ratio
$\text{sd}(\ktrans)/\text{median}(\ktrans)$.  Values of
$\text{sd}(\ktrans)$ are higher in areas of large $\ktrans$ values,
even when one adjusts for the estimated $\text{median}(\ktrans)$.
However, $\text{sd}(\ktrans)$ is quite low overall, usually less than
$0.1$, in the tumor ROI.

\subsubsection{Bayesian penalized splines}

An alternative to parametric methods of the biological system is the
function \code{dcemri.spline} where a non-parametric curve is fit to
the data using penalized splines \citep{Eilers1996,Marx1998}.
Smoothness of the curve and goodness-of-fit to the data are controlled
by two Gamma distributions: a prior for the adaptive smoothness
parameters (\code{ab.hyper}) and a prior for the variance of the
observational error (\code{ab.tauepsilon}).  Full details on the
methodology for Bayesian penalized $P$-splines for DCE-MRI are
provided in \cite{sch-etal:TMI09}.  For the following application of
\code{dcemri.spline} default values for the hyperparameters have been
selected.

<<RiderNeuroMRI+spline1,eval=FALSE>>=
mask.spline <- array(FALSE, dim(mask))
z <- 7
mask.spline[,,z] <- mask[,,z]
fit.spline <- dcemri.spline(conc[,,,-(1:8)], acqtimes[-(1:8)], mask.spline,
                            model="weinmann", aif="fritz.hansen", 
                            multicore=TRUE, nlr=TRUE)
writeNIfTI(fit.spline$ktrans, paste(perf, "ktrans","spline", sep="_"))
writeNIfTI(fit.spline$Fp, paste(perf, "Fp","spline", sep="_"))
@ 
<<RiderNeuroMRI+spline2,echo=FALSE,eval=FALSE,results=hide>>=
fit.spline <- list(ktrans=readNIfTI(paste(perf, "ktrans","spline", sep="_")))
png(file=paste(paste(perf, "ktrans", "spline", sep="_"), "png", sep="."), 
    width=480, height=480)
overlay(as(dynamic[x,y,,], "nifti"), 
        ifelse(mask[x,y,], fit.spline$ktrans[x,y,], NA), z=7, w=11, 
        zlim.x=zrx, zlim.y=c(0,0.1), plot.type="single")
dev.off()
png(file=paste(paste(perf, "Fp", "spline", sep="_"), "png", sep="."), 
    width=480, height=480)
fit.spline$Fp <- readNIfTI(paste(perf, "Fp","spline", sep="_"))
overlay(as(dynamic[x,y,,], "nifti"), 
        ifelse(mask[x,y,], fit.spline$Fp[x,y,], NA), z=7, w=11, 
        zlim.x=zrx, zlim.y=c(0,0.2), plot.type="single")
dev.off()
@

As a summary statistic the maximum of the response function may be
used.  Alternatively, a response model may be derived from the
response function (\code{nlr=TRUE}).  Please note that a sample of the
posterior PDF is given for the response function and, hence, a
non-linear fit to the response model is performed for each response
function in the sample.  The \code{dcemri.spline} function supports
two models, the standard Kety model and the adiabatic approximation of
tissue homogeneity (AATH) \citep{stl-lee:aath}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[tbp]
  \begin{center}
    \begin{tabular}{cc}
      $F_p$ & $\ktrans$\\
      \includegraphics*[width=0.4\textwidth]{281949_19040721_perfusion_Fp_spline.png}&
      \includegraphics*[width=0.4\textwidth]{281949_19040721_perfusion_ktrans_spline.png}
    \end{tabular}
  \end{center}
  \caption{Statistical images of $\ktrans$ for the RIDER Neuro MRI
    data (slice $z=7$) using Bayesian penalized splines.  The
    parameter $F_p$ (left) is given by the maximum of the response
    function after deconvolution using an adaptive spline fitting
    procedure, while the parameter $\ktrans$ (right) is obtained from
    fitting a response model.  The units of $F_p$ are
    $[0,0.2]\;\text{min}^{-1}$ and the units of $\ktrans$ are
    $[0,0.1]\;\text{min}^{-1}$.}
  \label{fig:perfusion-spline}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Figure~\ref{fig:perfusion-spline} depicts the maximum perfusion $F_p$
parameter map for the central tumor slice.  Increased perfusion is
visible in this image, but overall the quality of the statistical
image is poor.  This is most likely due to the fact that the
acquisition protocol was not optimized for the AATH model, where high
temporal resolution is required for accurate parameter estimation.
Figure~\ref{fig:perfusion-spline} also shows the median $\ktrans$
parameter map estimated from fitting a Kety response model to the
estimated response function.  Here, compared to the results above,
$\ktrans$ is slightly increased in the top left area of the ROI due to
the negligence of the vascular compartment in the standard Kety model.

\subsection{Diffusion weighted imaging}

The RIDER Neuro MRI data repository does not provide DWI acquisitions
\emph{per se} but a diffusion tensor imaging (DTI) acquisition was
performed at each visit.  The analysis of DTI data is beyond the scope
of this article, but the interested reader is pointed to the following
references: \citet{Horsfield2002,tof:book}.  The methodology behind
DWI and DTI are virtually identical, so we will ignore the extra
information provided by a DTI acquisition and analyze the
non-directional aspects of the diffusion process here.  We acknowledge
the fact that the ADC values derived from this DTI acquisition may
differ slightly from those estimated using a more common DWI sequence.

There are 13 data volumes in the DWI acquisition: a single T2-weighted
image without diffusion weighting ($b=0$) and 12 volumes with
different gradient encodings but the same diffusion weighting.  The
$b$-value for this acquisition is $b=1000\;\text{s/mm}^2$ (Barboriak,
\emph{personal communication}), a common value in clinical practice.
As previously noted this acquisition protocol has not been optimized
for ADC estimation, and will include both perfusion and diffusion
effects, but is adequate for the purpose of illustration.

<<RIDERNeuroMRI+dwi1,eval=FALSE>>=
tensor <- system.file(file.path("nifti", sub("perfusion", "axtensor", perf)),
                      package="dcemriS4")
(dwi <- readNIfTI(tensor))
tmask <- readANALYZE(sub(".nii", "_mask.hdr", tensor))
tmask <- ifelse(tmask > 0, TRUE, FALSE)
b <- c(0, rep(1000, ntim(dwi)-1)) # from Daniel Barboriak!
fit.adc <- ADC.fast(dwi, b, tmask)
writeNIfTI(fit.adc$S0, paste(tensor, "S0", sep="_"))
writeNIfTI(fit.adc$D, paste(tensor, "D", sep="_"))
@ 
<<RiderNeuroMRI+dwi2,echo=FALSE,eval=FALSE,results=hide>>=
fit.adc <- list(D=readNIfTI(paste(tensor, "D", sep="_")))
png(file=paste(paste(tensor, "D", sep="_"), "png", sep="."), 
    width=480, height=480)
overlay(as(dwi[,,17:20,], "nifti"),
        ifelse(tmask[,,17:20], fit.adc$D[,,17:20], NA), zlim.x=c(32,1024), 
        zlim.y=c(0.0005,0.003))
dev.off()
@ 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[tbp]
  \begin{center}
    \includegraphics*[width=0.75\textwidth]{281949_19040721_axtensor_D.png}
  \end{center}
  \caption{Statistical images of the apparent diffusion coefficient
    (ADC) values for the RIDER Neuro MRI data.  The range of displayed
    values of the ADC is $[0.0005,0.003]\,\text{mm$^2$/s}$.}
  \label{fig:tensor}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Given the larger voxel dimensions for the DTI acquisition
($5\;\text{mm}$ slice thickness), four axial slices with a generous
ROI were selected for ADC estimation and are displayed in
Figure~\ref{fig:tensor}.  The range of physical units for the ADC
values is $[0.0005,0.003]\,\text{mm$^2$/s}$, where high ADC values
correspond to the high mobility of water molecules in tissue.  Thus,
bright areas in the ROI may be found, for example, in the ventricles
or major blood vessels and to a lesser extent the tumor(s).  There
appear to be areas of high diffusion in the core of each tumor when
compared to either the rim of the tumor or ``normal'' brain tissue
(white or grey matter).  This may indicate sparse cell density in the
tumor ``cores'' due to necrotic or apoptotic processes and the
subsequent removal of cells in the tissue.

The application of DWI to oncology is still a relatively immature
field and caution should be used when interpreting any results because
of the indirect nature of MRI data acquisition; the pharmcodynamic
effects of treatment are being measured not directly in tissue, but
via the diffusivity of water molecules in the tissue.  However,
numerous authors have published on this topic and the interested
reader is encouraged to look at
\cite{koh-pad:DW-MRI,yan-etal:integration,pad-etal:neoplasia}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[tbp]
\small
% <<AuditTrail>>=
% audit.trail(dwi)
% @ 
\begin{verbatim}
R> audit.trail(dwi)
<audit-trail xmlns="http://www.dcemri.org/namespaces/audit-trail/1.0">
  <created>
    <call>oro.nifti::nifti(img = img, datatype = datatype)</call>
    <system>
      <r-version.version.string>R version 2.11.0 (2010-04-22)
        </r-version.version.string>
      <date>Thu May 27 08:40:18 PM 2010 BST</date>
      <user.LOGNAME>brandon</user.LOGNAME>
      <package-version.Version>0.1.4</package-version.Version>
    </system>
  </created>
  <saved>
    <workingDirectory>/home/guest/rider</workingDirectory>
    <filename>281949_19040721_axtensor</filename>
    <call>writeNIfTI(nim = uid.nifti, filename = fname)</call>
    <system>
      <r-version.version.string>R version 2.11.0 (2010-04-22)
        </r-version.version.string>
      <date>Thu May 27 08:40:23 PM 2010 BST</date>
      <user.LOGNAME>brandon</user.LOGNAME>
      <package-version.Version>0.1.4</package-version.Version>
    </system>
  </saved>
  <read>
    <workingDirectory>/home/guest/rider</workingDirectory>
    <filename>281949_19040721_axtensor.nii.gz</filename>
    <call>readNIfTI(fname = tensor)</call>
    <system>
      <r-version.version.string>R version 2.14.1 (2011-12-22)
        </r-version.version.string>
      <date>Fri Dec 30 10:01:15 2011 GMT</date>
      <user>bwhitcher</user>
      <package-version.Version>0.3.1</package-version.Version>
    </system>
  </read>
  <event>
    <type>processing</type>
    <call>ADC.fast(dwi, b, tmask)</call>
    <date>Fri Dec 30 10:01:20 2011 GMT</date>
    <user>bwhitcher</user>
  </event>
  <event>
    <type>completed</type>
    <call>ADC.fast(dwi, b, tmask)</call>
    <date>Fri Dec 30 10:01:23 2011 GMT</date>
    <user>bwhitcher</user>
  </event>
</audit-trail> 
\end{verbatim}
\normalsize
\caption{\proglang{XML}-based ``audit trail'' for the DWI acquisition
  of the RIDER Neuro MRI data.}
\label{fig:audit-trail}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsubsection{Audit trail}

The \pkg{dcemriS4} package supports and enhances the
\code{audit.trail} functionality of \pkg{oro.nifti}.  Hence, from any
object that is stored in the \code{nifti} class we can trace back all
the operations that have been performed on it.
Figure~\ref{fig:audit-trail} displays the \proglang{XML}-based audit
trail for the multi-dimensional array that holds the DWI acquisition
(in raw signal intensities).  The main blocks of information are the
\code{<created>}, \code{<saved>}, \code{<read>} and \code{<event>}
tags.  The first two tags occurred during the initial DICOM-to-NIfTI
conversion process and the last three tags were performed during
compilation of this document.  In each block pertinent information has
been recorded; such as the function call, version of \proglang{R}
being used, version of the \pkg{oro.nifti} package being used,
user ID, date, etc.  Notice that some of these properties have changed
over time, allowing one to accurately reproduce the data processing
stream (if necessary).

\section{Conclusions}

Quantitative analysis of dynamic contrast-enhanced MRI (DCE-MRI) and
diffusion-weighted imaging (DWI) data requires a series of processing
steps, including pre-processing of the MR signal, voxel-wise curve
fitting, and post-processing (e.g., statistical analysis of kinetic
parameters from a series of scans).  The \pkg{dcemriS4} package
provides a comprehensive set of functions for pre-processing and
parametric models for quantifying DCE-MRI and DWI data.

A (nearly) complete pipeline for the analysis of DCE-MRI and DWI data
has been established in \proglang{R}.  Acquisitions from the MR
scanner, assumed to be provided in DICOM format, are converted to the
NIfTI format using the \pkg{oro.dicom} and \pkg{oro.nifti} packages.
Using \pkg{dcemriS4} dynamic T1-weighted acquisitions are converted
into contrast agent concentration time curves on a voxel-by-voxel
basis.  A variety of compartmental models for the tissue kinetics, and
models for the arterial input function (AIF), are available.  Point
estimates for kinetic parameters are provided in a fast and robust
manner using either least-squares or \emph{maximum a posteriori}
techniques, and information about the uncertainty in these parameter
estimates may be obtained from the Bayesian MCMC (Markov Chain Monte
Carlo) algorithm; e.g., by looking at standard errors, credible
intervals or the entire posterior distribution.

The \pkg{dcemriS4} package utilizes the \code{nifti} class defined in
the \pkg{oro.nifti} package.  This allows one to retain metadata
information stored in the original DICOM data (e.g., patient ID or the
scan date) when performing an analysis.  In addition, each step in the
data analysis pipeline are recorded using the audit trail capability
provided by \pkg{oro.nifti}.  Hence, results may be reproduced in a
straightforward manner and errors in the analysis may be identified
efficiently.

The \pkg{dcemriS4} package is available from CRAN
(\url{http://CRAN.R-project.org}) and also from SourceForge
(\url{http://sourceforge.net/projects/dcemri}) under a BSD license.
The website \url{http://www.dcemri.org} has been established as a
convenient front end to the software development project and news
items are regularly provided on the blog
(\url{http://dcemri.blogspot.com}).

\section*{Acknowledgments}

The authors would like to thank the National Biomedical Imaging
Archive (NBIA), the National Cancer Institute (NCI), the National
Institute of Health (NIH) and all institutions that have contributed
medical imaging data to the public domain.  VS is supported by the 
German Research Council (DFG SCHM 2747/1-1).

\bibliography{dcemriS4}

\end{document}

