Quantitative Analysis of Dynamic Contrast-Enhanced and Diﬀusion-Weighted Magnetic Resonance Imaging for Oncology in R

The package dcemriS4 provides a complete set of data analysis tools for quantitative assessment of dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI). Image processing is provided for the ANALYZE and NIfTI data formats as input with all parameter estimates being output in NIfTI format. Estimation of T1 relaxation from multiple ﬂip-angle acquisitions, using either constant or spatially-varying ﬂip 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 maximum a posteriori estimation. A non-parametric model, using penalized splines, is also available to characterize the contrast agent concentration time curves. Estimation of the apparent diﬀusion coeﬃcient (ADC) is provided for diﬀusion-weighted imaging. Given the size of multi-dimensional data sets commonly acquired in imaging studies, care has been taken to maximize computational eﬃciency and minimize memory usage. All methods are illustrated using both simulated and real-world medical imaging data available in the public domain.


Introduction
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 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, using standard data formats (such as ANALYZE or NIfTI) that may be visualized and manipulated in R (R Development Core Team 2011) or exported for accessibility across a wide variety of medical image analysis software packages.
The S4 designation in dcemriS4 means that S4 object classes are used throughout to ensure efficient and transparent manipulation of ANALYZE or NIfTI data structures. Data input and output rely upon the oro.nifti package (Whitcher et al. 2011a;Tabelow et al. 2011;Whitcher et al. 2011b). Parameter estimates in dcemriS4 inherit attributes from the incoming ANALYZE/NIfTI objects in order to preserve anatomical and physiological information for appropriate visualization. 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, dcemriS4 supports basic parallel computing by incorporating the multicore package (Urbanek 2011). Computations are easily parallelized with the variable multicore = TRUE in the most computationally expensive functions.

Dynamic contrast-enhanced magnetic resonance imaging
While initial application of DCE-MRI was in the characterization of blood-brain barrier (BBB) integrity (Tofts and Kermode 1984;Larsson et al. 1990;Larsson and Tofts 1992), 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 de novo generation of blood vessels, are key biological processes that supply nutrients to tissue (Collins and Padhani 2004). 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 (Choyke et al. 2003). Thus, the ability to study malignant vasculature through non-invasive MRI methods is advantageous. We will cover a DCE-MRI data analysis pipeline using dcemriS4 that is fully quantitative and produces biologically-relevant parameter estimates from a compartmental model. The interested reader will find several chapters in Jackson et al. (2005) and Parker and Padhani (2003) 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 (Bernstein et al. 2004). 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 seconds before a bolus injection of the low molecular weight contrast agent -a gadolinium chelate -and continue for 8-12 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, 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 dcemriS4 consists of the following steps: 1. Pre-processing of the T1 signal (e.g., motion correction, co-registration, correction of the B1 field) is introduced in Sections 2.1 and 2.2.
2. Estimation of voxel-wise contrast agent concentration time curves is introduced in Section 2.3.
3. Determination of the arterial input function (AIF), either from the literature or by data-driven methods, is introduced in Section 2.4.
4. Parameter estimation for a given compartmental model is introduced in Section 2.5.
5. Statistical inference on kinetic parameters for differences between scans of a single patient or between distinct patients, is discussed in Section 2.6.

Motion correction and co-registration
Basic motion correction within an acquisition, and co-registration between acquired series, is available using the ftm function for fast template matching (Lewis 1995). 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 R is currently provided by the RNiftyReg package (Clayden 2011). Access to the image registration routines from the Insight Toolkit (ITK, Insight Software Consortium 2011) is highly desirable and under investigation.

B1 mapping via the saturated double-angle method
For in vivo MRI at high field (≥ 3 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 (Cunningham et al. 2006).
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 α 1 and I 2 with prescribed angle α 2 = 2α 1 . All other signal-affecting sequence parameters are kept constant. For each voxel, the ratio of magnitude images satisfies I 2 (r) where r represents spatial position and α 1 (r) and α 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 α(r) = arccos I 2 (r) 2I 1 (r) A long repetition time (TR ≥ 5×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 (T1, TR) = f 2 (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 α 2 excitation was used for the preceding acquisition.

Example using phantom data
Using data acquired from a T1 phantom at two flip angles, α 1 = 60 • and α 2 = 120 • , we compute the multiplicative factor relative to the low flip angle using the saturated doubleangle method (Cunningham et al. 2006). Note, repeated acquisitions (five) of each flip angle were obtained and force the additional rowMeans step to average the results from the function dam (double-angle method) in the code below.
R> fpmask <-system.file(file.path("nifti", "t1_phantom_mask.nii.gz"), + package = "dcemriS4") R> t1pmask <-readNIfTI(fpmask) R> pmask <-nifti(array(t1pmask [, , 25], dim(t1pmask))) 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 2. Boxplots summarize the estimated T1 relaxation times, across all voxel in the ROI defined by 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.

T1 relaxation rate and contrast agent concentration
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 (Buckley and Parker 2005). 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 (Buxton 2002). 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 where E 10 = exp(− TR ·R 10 ), relates the observed signal intensity S(·) with the parameters of interest when varying the flip angle θ prior to the injection of the contrast agent. Note, a repetition time of TR ≈ 4 ms is common practice for pulse sequences in clinical applications. The Levenberg-Marquardt algorithm, provided in minpack.lm (Elzhov et al. 2010), is applied to estimate the parameters; see the discussion by Ahearn et al. (2005). 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 where the flip angle θ ∈ [10 • , 30 • ] 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 R1.fast (embedded within CA.fast) rearranges the multi-dimensional structure of the multiple flip-angle acquisitions into a single matrix, to take advantage of internal R functions instead of loops, and calls 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 CA.fast2, where a linear approximation is applied to estimate the parameters (Wang et al. 1987).
The final step of the conversion of the dynamic signal intensities to contrast agent concentration, using the CA.fast function, is performed via 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 (Buckley and Parker 2005). For computational reasons we follow the method of Li et al. (2000).

Arterial input function
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.
where D = 0.1 mmol/kg, a 1 = 3.99 kg/l, a 2 = 4.78 kg/l, m 1 = 0.144 min −1 and m 2 = 0.0111 min −1 (Weinmann et al. 1984;Tofts and Kermode 1984); or D = 1.0 mmol/kg, a 1 = 2.4 kg/l, a 2 = 0.62 kg/l, m 1 = 3.0 and m 2 = 0.016 (Fritz-Hansen et al. 1996). 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 (Parker et al. 2006) and sums of exponentials (Orton et al. 2008). The dcemriS4 package has incorporated the sums-of-exponentials model (Orton et al. 2008), where the unknown parameters β = (A B , µ B , A G , µ G ) are estimated using nonlinear regression. Using the AIF defined in Buckley (2002), we illustrate fitting a parametric model to characterize observed data. The orton.exp.lm function provides this capability using the so-called double-exponential parametric form orton.exp (9).
R> data("buckley") R> aifparams <-with(buckley, orton.exp.lm(time.min, input)) R> fit.aif <-with(aifparams, aif.orton.exp(buckley$time.min, + AB, muB, AG, muG)) Figure 3 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.

Kinetic parameter estimation
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 Collins and Padhani (2004) and Buckley and Parker (2005) for discussions on this point. The standard Kety model (Kety 1960), 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 dcemriS4. Regardless of which parametric model is chosen for the biological system, the contrast agent concentration time curve at each voxel in the ROI is approximated using the convolution of an AIF and the compartmental model; i.e., The K trans parameter represents the volume transfer constant between the plasma and the extravascular extracellular space (EES) per minute, and k ep is the rate constant between EES and blood plasma. The parameter v p , in the so-called "extended" Kety model (11), describes the fraction of contrast agent in the plasma, while 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: 1. dcemri.lm: Non-linear regression using non-linear least squares (Levenberg-Marquardt optimization).

Non-linear regression
Least-square estimates of the kinetic parameters (K trans , k ep ), also for v p for the extended Kety model, are provided in 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 (10) with the literature-based AIF (8) produces a statistical model that is given by where (t) is the observational error at time t, θ 1 = log(K trans ) and θ 2 = log(k ep ). The parametrization (θ 1 , θ 2 ) is used instead of (K trans , k ep ) to ensure positive values for both transfer rates. We assume the expected value of the noise term to be zero; i.e., E( ) = 0. Inference is performed by minimizing the sum of squares of the observational errors min The parameter v e is calculated using (12).
Model parameters are estimated, along with asymptotic standard errors, using the Levenberg-Marquardt algorithm (Moré 1978) in 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 may not be accurate (Schmid et al. 2006).

Bayesian model
A hierarchical Bayesian model can be described in three stages: the data model, the process model and the prior parameters.
1. For the 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 σ 2 , is given by This is the Bayesian analogue to the application of the least-squares fitting method in the non-linear regression approach.
2. For the process model we use the single-compartment model (10) or an extended Kety model (11). Evaluating the convolution for the single-compartment model produces As previously discussed, the kinetic parameters K trans and k ep are transfer rates and must remain positive. Gaussian priors on their logarithmic transforms ensure this constraint is met. In breast tissue, for example, reasonable priors for both parameters should not exceed values of approximately 20 min −1 (Schmid et al. 2006). Hence, we use parameters a(K trans ) = a(k ep ) = 0 and b(K trans ) = b(k ep ) = 1. Thus, the expected value of K trans and k ep is one, and with 99.86% probability a priori neither kinetic parameter will exceed 20 min −1 . For scans covering other tissue types, the hyperparameters a(K trans ), b(K trans ) and a(k ep ), b(k ep ) may be adjusted accordingly when calling 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 a priori expected value 3. For the prior parameter, in this case the variance of the observational error, we apply a flat inverse Gamma prior with default parameters a(σ 2 ) = 1 and b(σ 2 ) = 0.001 that reflect our lack of prior information.
The three stages of the hierarchical model fully specify our a priori knowledge. To combine this with the observed data and produce a posteriori knowledge, we apply Bayes' theorem where h = (K trans , k ep , σ 2 ) denotes the vector of all parameters across all voxel, π(h) the product of the prior PDFs and (C t | h) denotes the (Gaussian) likelihood function of C t (t) from (14). In the Bayesian framework, conclusions are drawn from the joint posterior PDF only. Two functions are provided to exploit the posterior PDF: The function dcemri.map provides voxel-wise maximum a posteriori (MAP) estimators (DeGroot 1970) using the Nelder-Mead algorithm provided in 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.
The function dcemri.bayes provides the posterior median as the summary statistic for (K trans , k ep , v p ), along with the posterior standard deviation for all statistics, by sampling from the posterior using MCMC (Gilks et al. 1996). All samples from the joint posterior distribution may be saved using the option 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 C and linked to 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 multicore package by setting the option multicore = TRUE.

Bayesian penalized splines
An alternative to parametric modeling, the function dcemri.spline may be used to deconvolve and de-noise the contrast agent concentration time curves using an adaptive penalized spline approach (Schmid et al. 2009b). A Bayesian hierarchical model is constructed 1. The data model is the Gaussian observation model (14).
2. For the process model a general approach is used, such that 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 where B is a B-spline design matrix. An adaptive second-difference penalty is used on the spline regression parameters β j ; i.e.
where e j ∼N (0, δ 2 i ). Note, a first-order penalty is also available using the option rw = 1.
3. The prior for the adaptive smoothing parameter δ 2 i is given by with default parameters a(δ) = b(δ) = 10 −5 that provide nearly flat prior information, and an Inverse Gamma prior for the observational error (18).
Making use of Bayes' formula (19), the posterior is assessed using an MCMC algorithm. By default, the function dcemri.spline returns the median of the maximum F p of the response function R(t) per voxel. The median response function (response = TRUE) and the fitted contrast agent concentration time curve (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: 1. Find the minimum time t, for which the contrast agent concentration curve significantly exceeds zero, 2. Compute the gradient of C t at point t, exploiting the derivative of the B-spline,

Compute the onset time as
To provide estimates of the kinetic parameters from a compartmental model, a parametric model may be applied to the estimated response curve (nlr = TRUE). At this point in time a single exponential model ("weinmann") or the adiabatic approximation to the tissue homogeneity model ("AATH") are available. For the AATH model, the response function is given by 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 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 (samples = TRUE). Parallel computing may be accessed using the multicore package (multicore = TRUE).

Estimating the kinetic curves
Using kinetic parameters estimated with one of the methods above, the functions kineticModel or 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 nifti class objects of kinetic parameters can be given to kineticModel to produce voxel-wise estimates of the compartmental model.

Statistical inference
No specific support is provided to perform hypothesis testing on the kinetic parameters in dcemriS4. We recommend one uses built-in functions in 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 (Schmid et al. 2009a), where an implementation is under development in the software project PILFER at 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 K trans and k ep in DCE-MRI using functional data analysis (O'Connor et al. 2010).

Diffusion weighted imaging
Diffusion weighted imaging (DWI) is an imaging biomarker that is rapidly becoming popular and widely applied in oncology (Chenevert et al. 2002;Koh and Padhani 2006). DWI allows one to quantify the diffusion behavior of water by estimating the apparent diffusion coefficient (ADC, Wheeler-Kingshott et al. 2003). 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 in vivo (Moseley et al. 1990;Wheeler-Kingshott et al. 2003). These diffusion measurements reflect the effective displacement of water molecules allowed to migrate for a given time (Le Bihan et al. 1988). Using the Stejskal-Tanner equation one may solve for the unknown diffusion to obtain the apparent diffusion coefficient (ADC) D (Wheeler- Kingshott et al. 2003). For completeness, S 0 is the (unknown) signal intensity without the diffusion weighting, S is the observed signal with the gradient applied, γ is the gyromagnetic ratio, G is the strength of the gradient pulse, δ is the duration of the gradient pulse and ∆ is the time between the two pulses. The micro-parameters (γ, G, ∆, δ) are selected prior to data acquisition and may be combined into a single parameter b = γ 2 G 2 δ 2 (∆ − δ/3), known as the b-value. The functions ADC.fast and 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 s/mm 2 ), at low diffusion weighting (b ≥ 100 s/mm 2 ) and higher diffusion weighting (b ≥ 500 s/mm 2 ). When estimating the ADC value, one should exclude any acquisitions with b ≤ 100 s/mm 2 to minimize the influence of perfusion effects (Padhani et al. 2009).
The diffusivity of water at room temperature without restrictions is approximately 3.0×10 −3 mm 2 /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 (Koh and Padhani 2006;Padhani et al. 2009).

The RIDER Neuro MRI data repository
The National Biomedical Imaging Archive (NBIA; http://cabig.nci.nih.gov/tools/NCIA) is a searchable, national repository integrating 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 DWI data analyzed here were downloaded from the "RIDER Neuro MRI" collection (http://wiki.nci. nih.gov/display/CIP/RIDER).

Dynamic contrast-enhanced MRI
Functions of the 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 CA.fast to estimate the intrinsic relaxation rate R 10 and equilibrium signal intensity m 0 from (3) and the contrast agent concentration curve C t (t) from (7). In order to save computation time and memory, we utilize a binary mask with a very limited region-of-interest (ROI) created in FSLView (Analysis Group, FMRIB, Oxford 2011) and saved in ANALYZE format.

Non-linear regression
A numeric optimization of the least-square criterion, using the Levenberg-Marquardt algorithm, is provided by 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 str2time from 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 rawtimes.txt. The literature-based AIF fritz.hansen is used here for illustrative purposes. Alternatives include estimating values for a parametric AIF (e.g., aif = "orton.exp") from the data and supplying them via the user option in dcemri.lm, or providing an empirical AIF (aif = "empirical") and passing the vector of AIF values via the user option.
R> acqtimes <-str2time(unique(sort(scan("rawtimes.txt", + quiet = TRUE))))$time R> acqtimes <-(acqtimes -acqtimes[9]) / 60 R> conc <-readNIfTI(paste(perf, "gdconc", sep = "_")) R> fit.lm <-dcemri.lm(conc, acqtimes, mask, model = "extended", + aif = "fritz.hansen", control = nls.lm.control(maxiter = 100), + multicore = TRUE, verbose = TRUE) R> writeNIfTI(fit.lm$ktrans, paste(perf, "ktrans", sep = "_")) R> overlay(dynamic, ifelse(mask, fit.lm$ktrans, NA), w = 11, + zlim.x = c(32, 256), zlim.y = c(0, 0.1)) Figure 4 shows the estimated K trans statistical images in the pre-defined ROI overlayed on the dynamic acquisition (for anatomical reference). Two rings of high K trans values are visually apparent from the statistical images, indicating the presence of two tumors. In both cases K trans 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 5 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 k ep 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.

Bayesian maximum a posteriori (MAP) estimation
Caution must be exercised when using non-linear regression algorithms, since the Levenberg-Marquardt algorithm used in dcemri.lm is 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 (Schmid et al. 2006). Two methods for parameter estimation from a Bayesian perspective are implemented in dcemriS4. The function dcemri.bayes uses Markov chain Monte Carlo (MCMC) to explore the posterior PDF (Gilks et al. 1996) and dcemri.map uses numerical optimization to maximize the posterior (DeGroot 1970).
From the non-linear regression analysis of the RIDER Neuro MRI data, it appears that approximately K trans ∈ [0, 0.1] and k ep ∈ [0, 1.25] (Figure 4). Hence, we use a Gaussian distribution with expected value a(K trans ) = log(0.05) and variance b(K trans ) = 1 on log(K trans ), and a Gaussian distribution with expected value a(k ep ) = log(0.7) and variance b(k ep ) = 1 on log(k ep ) as priors. For v p , we use the Beta distribution B(a(v p ) = 1, b(v p ) = 19); i.e., the expected value is given by E(v p ) = 0.05. Parameter estimation via dcemri.map follows a consistent user interface established in dcemri.lm.

Bayesian estimation via Markov chain Monte Carlo
Using the MCMC algorithm provided by 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 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 dcemri.bayes uses the default samples = FALSE and, hence, we are restricted to posterior medians and standard deviations for all parameters in the compartmental model. R> 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)) R> writeNIfTI(fit.bayes$ktrans, paste(perf, "ktrans", "bayes", sep = "_")) Figure 7 displays statistical summaries of K trans (posterior median and standard deviation) provided the default settings of dcemri.bayes. It is clear that the posterior median differs from the MAP estimator (reproduced in Figure 7 for a side-by-side comparison) across the majority of non-tumor voxel in the ROI. However, K trans values around the enhancing rim of the tumor are similar across all three methods: dcemri.lm, dcemri.map and dcemri.bayes. Figure 7 also provides the posterior standard deviation of K trans and an image of the ratio sd(K trans )/median(K trans ). Values of sd(K trans ) are higher in areas of large K trans values, even when one adjusts for the estimated median(K trans ). However, sd(K trans ) is quite low overall, usually less than 0.1 min −1 , in the tumor ROI.

Bayesian penalized splines
An alternative to parametric methods of the biological system is the function dcemri.spline where a non-parametric curve is fit to the data using penalized splines (Eilers and Marx 1996;Marx and Eilers 1998) controlled by two Gamma distributions: a prior for the adaptive smoothness parameters (ab.hyper) and a prior for the variance of the observational error (ab.tauepsilon). Full details on the methodology for Bayesian penalized P -splines for DCE-MRI are provided in Schmid et al. (2009b). For the following application of dcemri.spline default values for the hyperparameters have been selected.
R> mask.spline <-array(FALSE, dim(mask)) R> z <-7 , + mask.spline, model = "weinmann", aif = "fritz.hansen", + multicore = TRUE, nlr = TRUE) R> writeNIfTI(fit.spline$ktrans, paste(perf, "ktrans", "spline", + sep = "_")) R> writeNIfTI(fit.spline$Fp, paste(perf, "Fp", "spline", sep = "_")) As a summary statistic the maximum of the response function may be used. Alternatively, a response model may be derived from the response function (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 dcemri.spline function supports two models, the standard Kety model and the adiabatic approximation of tissue homogeneity (AATH) (St Lawrence and Lee 1998). Figure 8 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 8 also shows the median K trans parameter map estimated from fitting a Kety response model to the estimated response function. Here, compared to the results above, K trans is slightly increased in the top left area of the ROI due to the negligence of the vascular compartment in the standard Kety model.

Diffusion weighted imaging
The RIDER Neuro MRI data repository does not provide DWI acquisitions 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: Horsfield and Jones (2002) and Tofts (2003). 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 s/mm 2 (Barboriak, 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.
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 Koh and Padhani (2006)

Audit trail
The dcemriS4 package supports and enhances the audit.trail functionality of oro.nifti. Hence, from any object that is stored in the nifti class we can trace back all the operations that have been performed on it. Figure 10 displays the XML-based audit trail for the multidimensional array that holds the DWI acquisition (in raw signal intensities). The main blocks of information are the <created>, <saved>, <read> and <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 R being used, version of the 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.

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 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 R. Acquisitions from the MR scanner, assumed to be provided in DICOM format, are converted to the NIfTI format using the oro.dicom and oro.nifti packages. Using dcem-riS4 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 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 dcemriS4 package utilizes the nifti class defined in the 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 oro.nifti. Hence, results may be reproduced in a straightforward manner and errors in the analysis may be identified efficiently.
The dcemriS4 package is available under a BSD license from the Comprehensive R Archive Network at http://CRAN.R-project.org/package=dcemri and also from SourceForge at http://sourceforge.net/projects/dcemri/. The website 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 (http://dcemri.blogspot.com/).