Fits the selected model structure to voxel-wise contrast agent concentration data.
KAT(file = "concatenate.KAT.with.KAT.checkData.RData",
results_file="my_results", method.optimization = "L-BFGS-B",
show.rt.fits = FALSE, param.for.avdt = "Ktrans", range.map = 1.5,
cutoff.map = 0.85, export.matlab = TRUE, export.RData = TRUE,
verbose=FALSE, show.errors = FALSE, try.silent=TRUE, fracGTzero = 0.75,
AIF.shift = "", Force.AIF.peak = FALSE, tlag.Tofts.on = FALSE,
est.per.voxel.tlag = FALSE, ...)Specify the file to be analyzed.
Specify the absolute path to the folder, including file name with no extension (.RData and/or .mat will be automatically added) of the file where results are to be saved.
Optimization method (Nelder-Mead, BFGS, CG, L-BFGS-B, SANN); default=L-BFGS-B.
Shows voxel-wise fits as each ROI is processed (default=FALSE).
Select extended Tofts model param to display using avdt (Ktrans,
kep, ve, vb; default=Ktrans).
Specifies range of color scale relative to the max value within map (typically a value between 1 and 3; default=1.5).
Truncate parametric map values by max value \(\times\) cutoff (typically a value between 0 and 1; default=0.85).
Save results in a matlab file? (default=TRUE). NOTE: matlab files
are intended only for use with MATLAB and are not meant to be read back into R for viewing with the advanced voxel
diagnosis tool.
Save results in an RData file? (default=TRUE).
Should voxel-wise contrast agent curves and other voxel specific
information be printed in the terminal
during the parameter estimation process? (default=FALSE).
Should errors messages be printed while KAT is running? Default is show.errors = FALSE.
Should the silent argument within try functions be TRUE or
FALSE? Default is try.silent=TRUE.
Voxels are excluded from analysis if less than 'fracGTzero' of
contrast agent concentrations in a voxel time series are greater
than zero. (a value between 0 and 1; default is fracGTzero = 0.75).
Specify if your vascular input function is based on arterial or
venous data. Possible values are "ARTERY" or "VEIN" (or
"NONE" if you prefer to not estimate the \(t_{lag}\) parameter). This
argument must be specified when running KAT.
Do you want to force the peak values of the shifted Vascular input
Function to be equal to the peak value of the original, raw VIF that
is read into KAT? The value of this argument is ignored if
AIF.shift="NONE". (default=FALSE).
Do you want to estimate \(t_{lag}\) for the Tofts model
(tlag.Tofts.on = TRUE)? If tlag.Tofts.on = FALSE, \(t_{lag}\) for the Tofts
model will be set equal to the \(t_{lag}\) value estimated for
the extended Tofts model. If tlag.Tofts.on = TRUE a single
\(t_{lag}\) value will be estimated based on the median
contrast agent profile over the slice of interest. Use with caution, since there may be
parameter identifiability issues associated with estimated
\(t_{lag}\) when fitting the Tofts model to data. (default=FALSE).
Do you want to estimate \(t_{lag}\) on a per-voxel basis for
the extended Tofts model
(est.per.voxel.tlag = TRUE)? Use with caution, since there may be
parameter identifiability issues associated with estimated
\(t_{lag}\) on a per-voxel basis. Note that this argument does
not impact per-voxel tlag values within the Tofts model, which will
always use \(t_{lag}\) values based on median contrast agent
profiles. cf. tlag.Tofts.on.
Pass arguments to functions within KAT.
Demo
Run the KAT benchmark test by typing
R> demo(KAT, ask=FALSE)
at the R prompt, to analyze the simulated DCE-MRI data set described in dcemri.data.
Model equations [refs. 1-2]
Tofts Model
$$\frac{dC_t(t)}{dt}=K^{trans}\mbox{VIF}(t \pm t_{lag})-k_{ep}C_t(t)$$
and
$$C(t) = C_t(t).$$
Extended Tofts Model
$$\frac{dC_t(t)}{dt}=K^{trans}\mbox{VIF}(t \pm t_{lag})-k_{ep}C_t(t)$$
and
$$C(t)=v_b\mbox{VIF}(t \pm t_{lag})+C_t(t),$$
where
$$C_{t}=\mbox{tissue/region of interest}$$
and
$$C(t): \mbox{measurement model for ROI corresponding to observed CA conc}.$$
VIF\((t \pm t_{lag})\) represents the Vascular Input Function as VIF\((t + t_{lag})\) if the measured VIF is based on arterial data or VIF\((t - t_{lag})\) if the measured VIF is based on venous data.
Objective Function [ref. 3]
The objective function based on maximum likelihood can be written as
$$\mbox{OF}_{M}=\frac{1}{n_D}\sum_{i=1}^{n_D}\left[ log \left( \frac{1}{\mbox{SD}_i^2} \right) + \left( \frac{(y_i-s(\hat{p},t_i))^2}{\mbox{SD}_i^2} \right) \right],$$
where
$$\mbox{SD}_i = \mbox{standard deviation of each data point } i \mbox{ in the intensity/concentration time curve},$$
$$y_i = \mbox{data point } i \mbox{ in the intensity/concentration time curve},$$
$$s(\hat{p},t_i) = \mbox{simulated data point at parameter vector } \hat{p} \mbox{ and time point } t_i,$$
and
$$n_D = \mbox{number of data points in the intensity/concentration time curve}.$$
If observed data is unweighted, i.e., \(\mbox{SD}_i=1\), \(\mbox{OF}_{M}\) is equal to the mean residual sum of squares, \(\mbox{OF}_R\), or
$$\mbox{OF}_R = \frac{\sum_{i=1}^{n_D} (y_i-s(\hat{p},t_i))^2}{n_D} = \frac{R\!S\!S}{n_D} = \overline{R\!S\!S}.$$
The objective function implemented in KAT is written as
$$\mbox{OF}_{K\!AT} = R\!S\!S = \sum_{i=1}^{n_D} (y_i-s(\hat{p},t_i))^2.$$
Model discrimination [ref. 4]
Fits of the Tofts and Extended Tofts model to the intensity/concentration time curve are compared via the Akaike Information Criterion \(\mbox{(AIC)}\), written generally as
$$\mbox{AIC}=-2 \cdot \mbox{log-likelihood}+n_P.$$
When applying least squares regression, i.e., \(OF_{K\!AT}\), to observed data with Gaussian variability, the \(\mbox{AIC}\) is written as
$$\mbox{AIC}=n_D \cdot log\left( \mbox{OF}_{R} \right) + 2(n_P+1).$$
A correction term for small sample sizes \((n_D/n_P < 40)\) can be derived, yielding
$$\mbox{AIC}_c=n_D \cdot log\left( \mbox{OF}_{R} \right) + 2(n_P+1) + \frac{2(n_P+1)(n_P+2)}{n_D-n_P-2},$$
where
$$n_P=\mbox{number of estimated model parameters}.$$
Coefficients of Variation (CVs) for estimated parameters [refs. 5-6]
The Hessian matrix \(\mathbf{H}(\hat{p})\) is
calculated numerically by R during the
parameter estimation process, so that the covariance matrix
(\(cov\)) and %CVs for model parameters estimated
using the \(\mbox{OF}_{K\!AT}\), i.e., \(R\!S\!S\), objective function may be approximated as
$$cov(\hat{p})=\frac{n_P \cdot \mbox{OF}_{K\!AT}}{n_D-n_P}\mathbf{H}^{-1}(\hat{p})$$
and
$$\mbox{\%CV}(\hat{p})=\frac{\sqrt{\mbox{diag}[cov(\hat{p})]}}{\hat{p}} \times 100\%,$$
where \(\mbox{diag}[cov(\hat{p})]\) is a vector composed of the diagonal elements of \(cov\) and \(\hat{p}\) is the vector of final parameter estimates. Note that this method for calculating %CVs assumes that variability in measured data points follows a Gaussian distribution. Thus, large outlier data points in the concentration/intensity curve, for example those caused by patient motion, may inflate estimated %CVs.
Contents of output file
args: List of all arguments specified when running the
KAT function plus a few additional values generated during the run.
cc: nx x ny x nt array of voxel-wise contrast agent
concentration/intensity-time curves for all voxels within the Field of View.
ccroi: nx x ny x nt array of voxel-wise contrast agent
concentration/intensity-time curves for all voxels within the Region of
Interest.
ccmedian: Median contrast agent-time profile within the Region
of Interest.
maptimes: nt x 1 time vector.
aif: nt x 1 Vascular Input Function.
aifshifted: nt x 1 Time shifted (by \(t_{lag}\))
Vascular Input Function.
maskroi: nx x ny array indicating the region
of interest.
mapKtransxT: nx x ny array of Ktrans values estimated
using the extended Tofts model.
mapKtransxTcv: nx x ny array of %CVs associated with
Ktrans values in mapKtransxT.
mapkepxT: nx x ny array of kep values estimated
using the extended Tofts model.
mapkepxTcv: nx x ny array of %CVs associated with
kep values in mapkepxT.
mapvbxT: nx x ny array of vb values estimated
using the extended Tofts model.
mapvbxTcv: nx x ny array of %CVs associated with
vb values in mapvbxT.
mapvexT: nx x ny array of ve values estimated
using the extended Tofts model. Note that \(v_e\) is not
directly fitted to the concentration/intensity data but is calculated as
\(v_e=K^{trans}/k_{ep}\)
mapOptimValuexT: nx x ny array of objective function values
for voxels where the extended Tofts model has successfully converged (convergence/exit code=0; see
?optim).
mapfitfailuresxT: Map indicating per voxel optimization exit
codes for the extended Tofts model. 0: indicates successful
completion. 1: indicates that the iteration limit maxit had been
reached. 10: indicates degeneracy of the Nelder-Mead simplex. 51:
indicates a warning from the L-BFGS-B method; see component message for
further details. 52: indicates an error from the L-BFGS-B method; see
component message for further details. 99: indicates a try error
(optimization routine crashed). -2: indicates that voxel failed fracGTzero
test. Voxels with exit codes > 0 will appear
white when using the interactive AVDT and are excluded from subsequent
analysis.
paramestmedianxT: List of median values of all extended Tofts
model parameters fitted on a voxel-wise basis to contrast agent curves
within the ROI; includes the percent of total fitted voxels that are
classified as fit failures.
roimedianfittedxTofts: Simulated contrast agent-time profile
generated by fitting the extended Tofts model to median contrast agent
values within the Region of Interest.
paramestwholeroixTofts: Model parameters estimated by fitting the
extended Tofts model to median concentration/intensity data across the
Region of Interest. These parameters are used to generate
roimedianfittedxTofts.
cvwholeroixTofts: %CVs for extended Tofts model parameters
listed in paramestwholeroixTofts.
mapKtransT: nx x ny array of Ktrans values estimated
using the Tofts model.
mapKtransTcv: nx x ny array of %CVs associated with
Ktrans values in mapKtransT.
mapkepT: nx x ny array of kep values estimated
using the Tofts model.
mapkepTcv: nx x ny array of %CVs associated with
kep values in mapkepT.
mapveT: nx x ny array of ve values estimated
using the Tofts model. Note that \(v_e\) is not
directly fitted to the concentration/intensity data but is calculated as
\(v_e=K^{trans}/k_{ep}\).
mapOptimValueT: nx x ny array of objective function values
for voxels where the Tofts model has successfully converged (convergence/exit code=0; see
?optim).
mapfitfailuresT: Map indicating per voxel optimization exit
codes for the Tofts model. 0: indicates successful completion. 1:
indicates that the iteration limit maxit had been reached. 10: indicates
degeneracy of the Nelder-Mead simplex. 51: indicates a warning from the
L-BFGS-B method; see component message for further details. 52:
indicates an error from the L-BFGS-B method; see component message for
further details. 99: indicates a try error (optimization routine
crashed). -2: indicates that voxel failed fracGTzero
test. Voxels with exit codes > 0 will appear
white when using the interactive AVDT and are excluded from subsequent
analysis.
paramestmedianT: List of median values of all Tofts
model parameters fitted on a voxel-wise basis to contrast agent curves
within the ROI; includes the percent of total fitted voxels that are
classified as fit failures.
roimedianfittedTofts: Simulated contrast agent-time profile
generated by fitting the Tofts model to median contrast agent
values within the Region of Interest.
paramestwholeroiTofts: Model parameters estimated by fitting the
Tofts model to median concentration/intensity data across the
Region of Interest. These parameters are used to generate
roimedianfittedTofts.
proctimetotal: Total processing time in minutes.
roiplotparams: Cropping coordinates applied to the FOV and upper
limit of color bar for
visualization of parametric maps via the advanced voxel diagnostic tool
(AVDT).
KATversion: Version of KATforDCEMRI used to generate this output
file.
mapAICxT: nx x ny array of \(\mbox{AIC}_c\) values
for per-voxel fits of the extended Tofts model to
concentration/intensity data.
mapAICT: nx x ny array of \(\mbox{AIC}_c\) values
for per-voxel fits of the Tofts model to
concentration/intensity data.
mapAICcompare: nx x ny nx x ny array that contains
a ``1'' for voxels with a lower (lower=better model) AIC for the
extended Tofts model or a ``2'' for voxels with a lower AIC for the
Tofts model.
nx: x dimension of the FOV.
ny: y dimension of the FOV.
nt: number of elements in the time vector (number of time
points).
ccFittedxT: nx x ny x nt array of extended Tofts model simulations
fitted to voxel-wise contrast agent concentration-time data within the ROI.
ccFittedT: nx x ny x nt array of Tofts model simulations
fitted to voxel-wise contrast agent concentration-time data within the
ROI.
p0T: Initial parameter values for the Tofts model where \(K^{trans}(0)\) and
\(k_{ep}(0)\) are estimated using the numerical deconvolution method
described in the DATforDCEMRI package refs [refs 7-8].
p0xT: Initial values for the extended Tofts model where \(K^{trans}(0)\) and
\(k_{ep}(0)\) are the same as those used for p0T. Initial values for \(v_b\) and
\(t_{lag}\) are set to nominal values (0.05 and 7.2 seconds),
when fitting the
model to median data. Initial values for \(v_b\) and
\(t_{lag}\), when fitting the xTofts model to
per-voxel data, are set to those values estimated based on fits to the median
data.
IRFresults: Contains results of noncompartmental analysis of the
estimated Impulse Response Function (IRF), where AUC is the area
under the curve (AUC) of the IRF and is analogous to the Tofts parameter
\(v_e\), AUCMRT is the AUC of the IRF divided by the Mean
Residence Time (MRT) of the IRF and is analogous to the Tofts parameter
\(K^{trans}\) and AUCMRT.divby.AUC is AUCMRT divided by
AUC (equal to 1/MRT) and is analogous to the Tofts parameter
\(k_{ep}\). Each of these parameters is either corrected for
truncation error and contribution of \(v_b\) using
nominal values (corrnom), as described in the entry for
p0xT, or the actual values estimated based on fitting the
extended Tofts model to median intensity/concentration data
(corr). The method for truncation correction is described in
[ref 7].
mapEF: Map indicating enhancing voxel within the specified region
of interest, where a ``1'' indicates enhancement.
Model equations
[1] Tofts P, Kermode A (1991) 10.1002/mrm.1910170208
[2] Tofts et al. (1999) https://www.ncbi.nlm.nih.gov/pubmed/10508281
[3] Ferl GZ, Port RE (2012) 10.1038/clpt.2012.63
Objective Function
[4] Barrett PHR, Bell BM, Cobelli C, Golde H, Schumitzky A, Vicini P, and Foster DM (1998) 10.1016/S0026-0495(98)90064-6
Model discrimination
[5] Glatting G, Kletting P, Reske SN (2007). 10.1118/1.2794176
Coefficients of Variation (CVs) for estimated parameters
[6] Venables WN, Smith DM, R Core Team (2012) https://CRAN.R-project.org/doc/manuals/R-intro.pdf
R package for numerical deconvolution
[7] Ferl GZ, Xu L, Friesenhahn M, Bernstein LJ, Barboriak DP, Port RE (2010) DOI:10.1002/mrm.22335
[8] Ferl GZ (2011) 10.18637/jss.v044.i03
# NOT RUN {
## Create temporary directory for example code output files
temp_dir <- tempdir(check=FALSE)
##
current_dir <- getwd()
setwd(temp_dir)
##
## Run example code
demo(KAT, ask=FALSE)
##
setwd(current_dir)
##
## ANALYZE DATA FROM A SINGLE DCE-MRI SCAN
##
## Load MATLAB files into R
## R> aif <- readMat("mydatafile-AIF.mat")$aif
## R> ct <- readMat("mydatafile-CT.mat")$ct
## R> roi <- readMat("mydatafile-ROILES.mat")$roi
## R> tvec <- readMat("mydatafile-TVEC.mat")$tvec
##
## Check that the dimensionality of the loaded data is consistent and save
## as a single R object or RData file
## R> dcemri.data <- KAT.checkData(file.name="mydatafile", vector.times=tvec,
## map.CC=ct, mask.ROI=roi, vector.AIF=aif)
##
## Fit the Tofts and extended Tofts model to all ROIs in RData file
## R> KAT(file="mydatafile.RData", results_file= "mydatafile_out",
## show.rt.fits=TRUE, AIFshift="VEIN")
##
## Plot all ROIs in a single figure
## R> KAT.plot(F1="mydatafile_out_slice3.RData",
## F2="mydatafile_out_slice4.RData", F3="mydatafile_out_slice5.RData",
## F4="mydatafile_out_slice6.RData")
##
## Visualize and explore a parametric map for a single ROI
## R> KAT(file="mydatafile_out_slice6.RData")
# }
Run the code above in your browser using DataLab