Robust and Efficient Pharmacokinetic Parameter Estimation: Application to Prostate DCE-MRI

Soudabeh Kargar^{1}, Eric G Stinson^{2}, Eric A Borisch^{2}, Adam T Froemming^{2}, Akira G Kawashima^{3}, Lance A Mynderse^{4}, Joshua D Trzasko^{2}, and Stephen J Riederer^{2}

For each
voxel in the region of interest (ROI), the acquired tissue enhancement curve,
$$$C(t)$$$, is fitted to the perfusion model, $$$P(t)$$$, by minimizing a cost function. The
Tofts^{2} model shown in Eq. 1 and Eq. 2 estimates the tissue
perfusion. In this model, $$$I\left(K_{trans},K_{ep},t\right)$$$ is the tissue
impulse response function, $$$P(t)$$$ is
the estimated tissue perfusion, $$$AIF$$$ is the Arterial Input Function selected
from the L or R iliac artery close to the prostate, and $$$\otimes$$$is the
convolution operator. This model is restated in linear algebraic form, where $$$A$$$ is
a convolution matrix and $$$H\left(K_{ep}\right)$$$ is the discrete-time exponential
part of $$$I$$$. The cost function $$$J$$$ in Eq. 3 is the difference
between the acquired and estimated perfusion, $$$C(t)$$$ and $$$P(t)$$$ respectively, which
is minimized to find the optimum $$$K_{trans}$$$ and $$$K_{ep}$$$. $$$J$$$ depends on both $$$K_{trans}$$$ and $$$K_{ep}$$$, and
thus we
have a 2D optimization problem for each voxel. As shown in Eq. 4 a closed-form
expression for $$$K_{trans}$$$ (dependent on $$$K_{ep}$$$) can be derived. By replacing $$$K_{trans}$$$
from Eq. 4 in Eq. 3, we reduce the dimensionality of $$$J$$$ from 2D to 1D which is
shown in Eq. 5. This technique is called Variable Projection (VARPRO, or VP)^{4}.
The optimum $$$K_{ep}$$$ is efficiently determined via Golden Section search^{5}. The
average runtimes per iteration for one voxel using Matlab 2015b are about 0.53
msec and 1.22 msec for VARPRO and LM method, respectively.

$$I=I\left(K_{trans},K_{ep},t\right)=K_{trans}e^{K_{ep}t}=H\left(K_{ep}\right)K_{trans}\hspace{7cm}\text{(Eq.1)}$$

$$P(t)=AIF\otimes I=AI=AH\left(K_{ep}\right)K_{trans}\hspace{9.1cm}\text{(Eq.2)}$$

$$J\left(K_{trans},K_{ep}\right)=\|AIF\otimes I\left(K_{trans},K_{ep},t\right)-C(t)\|_{2}^{2}=\|AH\left(K_{ep}\right)K_{trans}-C(t)\|_{2}^{2}\hspace{1.2cm}\text{(Eq.3)}$$

$$\nabla_{K_{trans}}(J)=0\Rightarrow K_{trans}=\left(AH\left(K_{ep}\right)\right)^\dagger C(t)=B^\dagger C(t),where, B^\dagger=\left(B^*B\right)^{-1}B^*\hspace{1.7cm}\text{(Eq.4)}$$

$$J\left(K_{ep}\right)=\|BB^\dagger C(t)-C(t)\|^{2}_{2}\equiv-C^{*}BB^\dagger C\hspace{8.8cm}\text{(Eq.5)}$$

We
acquire 55 time frames of the prostate using an accelerated
3D time-resolved DCE-MRI^{6} sequence (256×384×38, 0.86×1.15×3.0 mm^{3},
frame time 6.6 sec). Perfusion parameters from pixels within the
prostate were estimated with the: (i) the above VARPRO, (ii) LM, and (iii) DynaCAD.
The residual defined as root mean square difference between $$$C(t)$$$ and $$$P(t)$$$, is
calculated using a subset of the time frames, (7 to 55). In our preliminary
study, we have performed the VARPRO technique on 15 cases which were matched
with the results of LM and DC.

Fig. 1A shows a plot of $$$J$$$ as a function of $$$K_{ep}$$$ (Eq. 5), for a representative voxel, demonstrating
how the optimum $$$K_{ep}$$$ can be readily found. Fig. 1B shows the fitted perfusion
estimate, $$$P(t)$$$ to the acquired
data, $$$C(t)$$$. Fig. 2A shows the residual, for VP and LM plotted versus
iteration. Fig. 2B shows the residual versus absolute computation time. Fig. 3 shows
the perfusion parameter maps for a representative prostate study (9x4 mm^{2} abnormality in left mid
peripheral zone suggests presence of significant prostate carcinoma;
Brachytherapy performed) using VARPRO, LM,
which are similar to maps produced with DynaCAD. Fig. 4A shows a plot of the normalized
residual versus the number of time frames used as input data. As shown, after about
30 time frames the residual reaches a plateau. The pharmacokinetic maps
according to a few of the subsets are shown in Fig. 4B with $$$K_{trans}$$$
on top and $$$K_{ep}$$$ at the bottom row.

1. American Cancer Society. http://www.cancer.org/cancer/prostatecancer/detailedguide/prostate-cancer-key-statistics

2.Tofts P, et al. Estimating kinetic parameters from dynamic contrast-enhanced T1-weighted MRI of a diffusable tracer standardized, JMRI. 1999; 10: 223-232

3. Marquardt D, et al. An algorithm for least-squares estimation of nonlinear parameters. SIAM J Appl. Math. 1963; 11(2): 431-441

4. Golub G, et al. The differentiation of pseudo-inverse and nonlinear least squares problems whose variables separate SIAM J. Numer. Anal. 1973;10(2): 413-432

5. Trzasko J, et al. Estimating T1 from multichannel variable flip angle SPGR sequences, MRM. 2013; 69(6):1787-1794

6. Froemming A, et al. The Application of sparse reconstruction to high spatio-temporal resolution dynamic contrast enhanced MRI of the prostate: Initial clinical experience with effect on image and parametric perfusion characteristic quality, ISMRM 2015; #1169

** Fig. 2:** Performance comparison for
residual versus **A:** number of
iterations and **B:** absolute computation time for VP and LM methods.

Proc. Intl. Soc. Mag. Reson. Med. 24 (2016)

2463