Accelerated Brain DCE-MRI Using Iterative Reconstruction With Total Generalized Variation Penalty for Quantitative Pharmacokinetic Analysis: A Feasibility Study

Abstract

Purpose:

To investigate the feasibility of using undersampled k-space data and an iterative image reconstruction method with total generalized variation penalty in the quantitative pharmacokinetic analysis for clinical brain dynamic contrast-enhanced magnetic resonance imaging.

Methods:

Eight brain dynamic contrast-enhanced magnetic resonance imaging scans were retrospectively studied. Two k-space sparse sampling strategies were designed to achieve a simulated image acquisition acceleration factor of 4. They are (1) a golden ratio–optimized 32-ray radial sampling profile and (2) a Cartesian-based random sampling profile with spatiotemporal-regularized sampling density constraints. The undersampled data were reconstructed to yield images using the investigated reconstruction technique. In quantitative pharmacokinetic analysis on a voxel-by-voxel basis, the rate constant K ^trans in the extended Tofts model and blood flow F _B and blood volume V _B from the 2-compartment exchange model were analyzed. Finally, the quantitative pharmacokinetic parameters calculated from the undersampled data were compared with the corresponding calculated values from the fully sampled data. To quantify each parameter’s accuracy calculated using the undersampled data, error in volume mean, total relative error, and cross-correlation were calculated.

Results:

The pharmacokinetic parameter maps generated from the undersampled data appeared comparable to the ones generated from the original full sampling data. Within the region of interest, most derived error in volume mean values in the region of interest was about 5% or lower, and the average error in volume mean of all parameter maps generated through either sampling strategy was about 3.54%. The average total relative error value of all parameter maps in region of interest was about 0.115, and the average cross-correlation of all parameter maps in region of interest was about 0.962. All investigated pharmacokinetic parameters had no significant differences between the result from original data and the reduced sampling data.

Conclusion:

With sparsely sampled k-space data in simulation of accelerated acquisition by a factor of 4, the investigated dynamic contrast-enhanced magnetic resonance imaging pharmacokinetic parameters can accurately estimate the total generalized variation-based iterative image reconstruction method for reliable clinical application.

Introduction

Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is an advanced MR technique that can be used as a noninvasive imaging tool for the measurement of microvascular properties. Based on the rapid acquisition of a series of images before, during, and after the intravenous administration of a low-molecular-weight contrast agent (CA), DCE-MRI is sensitive to differences in blood flow, blood volume, and vascular permeability that can be associated with tumor angiogenesis. ¹ By fitting the DCE-MRI data to certain pharmacokinetic (PK) models, volumetric maps of quantitative physiological parameters can be obtained in PK model analysis. Compared with computed tomography in probing microvascular characteristics, DCE-MRI has a higher sensitivity to CA with no use of ionizing radiation. ² In addition, as a functional imaging method against the functional nuclear medicine imaging methodology including positron emission tomography and single-photon emission computed tomography, DCE-MRI has the distinguishing feature of combining morphological and functional information acquirement in a single imaging session with potentially improved spatial resolution. ³ Because of these merits, DCE-MRI is a promising method for various clinical oncologic tasks and has been investigated in early diagnosis, ^4,5 tumor staging, ⁶ treatment planning, ^7,8 and treatment response assessment. ^{9 –13}

Of all possible factors that may influence the accuracy in either semiquantitative or quantitative DCE-MRI analysis, temporal resolution (Δt) is well documented as an important one. ^14,15 Because of the dynamic nature of DCE-MRI, high temporal resolution is desirable. In particular, the arterial input function (AIF), which depicts the CA dynamics within the feeding arterial blood, has to be determined via imaging a major arterial structure that exists in the imaging volume with high temporal resolution. ^16,17 To ensure the accuracy of PK analysis, the reliable AIF information derivation demands 1 second or faster temporal resolution in image acquisition to accurately capture the rapid washin and washout processes of CA in arterial blood. ^18,19 However, owing to the limited volume acquisition time or lack of arterial structure within the image volume, the individualized AIF measurement may not always be feasible. ²⁰ By assuming the similar CA dynamics of different subjects, population-derived mathematical AIF models have been proposed. ^{21 –23} However, these simplified models ignore the potentially large physiological variability among subjects and may introduce additional errors to the PK analysis. ²⁴ Additionally, insufficient temporal resolution may also affect the accuracy of PK analysis by introducing intensity uncertainties in the CA concentration measurement of tissue of interest. It has been noted by both imaging and simulation studies that when accurate AIF information was available, the error of CA rate constants increases as the temporal resolution degrades. ^25,26 As a result, improving temporal resolution is of particular interest in DCE-MRI studies toward the clinical application. ¹²

To accelerate DCE-MRI acquisition, various fast imaging methods using reduced k-space data sampling, such as sliding window (SW), 3D keyhole, ^27,28 and k-t Broad-use linear acquisition speed-up technique, ²⁹ have been proposed and investigated. In recent years, the theory of compress sensing (CS) proved that at the violation of classic Nyquist theorem, a small sampling fraction in certain patterns is sufficient to reconstruct compressible signals using a linear reconstruction subject to a specified penalty. ^{30 –32} The sparsity implicit in MR images under Fourier transformation suggests the great potential of CS in the acceleration of MRI studies. To date, CS application in dynamic MR has been successfully demonstrated in many studies. ^{33 –36} As for DCE-MRI, several studies have been conducted in the past to investigate CA rate constants using the reconstructed images with undersampled patient and animal data. ^{37 –41} In these studies, the concept of total variation (TV) based on the first-order derivative calculation was commonly adopted in the constrained image reconstruction to minimize the gradient of the reconstructed image. ⁴² In spite of its popularity, the concept of TV relies on the assumption that MR image consists of piecewise constant regions, yet this assumption may not always be valid in human MR imaging with the potential inhomogeneity of exciting field. ⁴³ With the existence of staircasing artifacts, the accuracy of TV-based reconstructed image may be compromised. The staircasing artifacts may further affect the accuracy of PK parameter map calculation in DCE-MRI. Another problem in these studies is that the AIF was approximated by published models or previously reported measurements, and these approximations may lead to errors in the evaluation of undersampling effect. Furthermore, few studies have been conducted to explore the feasibility of using undersampled k-space data for quantitative blood perfusion analysis.

In this study, we explored the clinical feasibility of using undersampled k-space data for quantitative analysis with 2 commonly adopted DCE-MRI PK models. Clinical brain DCE-MRI scans were used, and 2 k-space data sampling strategies were proposed to achieve an acceleration factor of 4. In the iterative image reconstruction, a penalty term based on the recently developed mathematical concept of second-order total generalized variation (TGV) was adopted. ⁴⁴ In the PK analysis, individualized AIF information was measured, and the accuracy of our approach was investigated via the quantitative evaluation of the derived PK parametric maps.

Materials and Methods

Reference Data

Eight sets of brain DCE-MRI scans from 6 patients with recurrent malignant glioma were selected from an institutional review board-approved clinical study as the reference images. The scans were acquired in the axial plane using a bird-cage quadrature head coil and a 3-dimensional (3D) spoiled-gradient recalled echo sequence on a clinical 1.5T scanner (General Electric, Fairfield, Connecticut). Each set of scan consists of preinjection calibration volumes at various flip angles and 60 postinjection volumes, and the temporal resolution was about 5.25 seconds of each volume (image matrix size: 256 × 256; acquisition matrix size: 256 × 128). The k-space data of the clinical scans were saved as the full sampling data in this study, and the undersampled k-space data were generated via sampling the full data with 2 specific undersampling strategies discussed subsequently.

Radial-Based Sampling Scheme

Because of a good balance of central/peripheral sampling weights with good low-frequency coverage and potentially improved point spread function, ⁴⁵ an optimized radial-based undersampling profile with multiple spokes was investigated to demonstrate the mathematical accuracy of this TGV-based reconstruction algorithm. Specifically, in this work, the undersampled k-space data were generated by sampling each slice of the full sampling data with a 32 non-Cartesian-ray radial sampling profile defined on a Cartesian grid. ³⁹ Related to the golden ratio, the sampling profile was generated by a constant azimuthal increment of Δβ = 111.25° (Figure 1), ⁴⁶ and radial rays were evenly spaced with time. ⁴⁷ This sampling pattern acquired approximately 11.5% of full k-space area and simulated an accelerated image acquisition by a factor of 4. Compared to the even-distributed sampling profile, the sampling efficiency of this golden ratio–based profile has no difference. ⁴⁸ Radial rays were evenly spaced with time, and only 3 different azimuthal gaps occurred in comparison with the evenly distributed rays. Thus, this sampling profile allows a posteriori adjustment of the temporal resolution. ⁴⁹

Figure 1.

Generation of the 32-ray radial sampling pattern at a golden ratio angular increment.

Cartesian-Based Sampling Scheme

Although radial sampling profile has its advantages for limited k-space acquisition, the higher demand for hardware implementation inhabits its application on most clinical task-oriented scanners. Besides, image reconstruction from radially sampled data requires certain regridding approaches that are often far from trivial. ⁵⁰ In contrast, sampling profiles defined in Cartesian coordinate system are straightforward and are suitable for easy implementation. However, because of evenly distributed sampling density on the frequency-encoding direction, more readout lines are required to achieve an adequate k-space coverage for backprojection. Early efforts were made to integrate SW technique into Cartesian undersampling profile sequences. ⁵¹ Although the feasibility was demonstrated, the achieved acceleration factor was limited without proper indication of possible application in DCE-MRI.

Based on above issues, a novel Cartesian-based sampling profile sequence with integrated SW technique was designed as the second sampling scheme of this work. First, a variable Poisson ellipsoid spatial sampling density function was built as demonstrated in Figure 2.

Figure 2.

Illustration of a 2-dimensional (2D) sampling profile with variable sampling density. The central red region is fully sampled, and the rest in gray level is sampled with varying sampling density.

Specifically, for each 2-dimensional (2D) full profile covered in the original acquisition of the clinical scan, the central 15% of k-space along the phase-encoding direction was fully sampled (indicated by red color). Then, the rest k-space region was treated as high-frequency (HF) sampling region, and each side of the region was evenly divided into 4 subregions along the phase-encoding direction. Each subregion was assigned a uniform sampling density ρ _i , which decreases exponentially from k-space center (indicated by gray levels with color bar):

${\text{ρ}}_i= {\text{ρ}}_0{\mathrm{e}}^{-r_i/\tau }$

1

where ρ₀ is the scaling factor to ensure the probability density normalization. The parameter r_i indicates the relative length (in phase-encoding direction) if the ith HF subregion in reference to the total HF region length, and the dimensionless decay constant τ was chosen as 1. In this work, HF data in 10% of the original k-space size was acquired in the designated HF region. Thus, an acceleration of 4 was achieved by sampling each 2D reference data with 25% data size (15% in central region and 10% in HF region).

For a series of 2D k-space data at a certain location in the dynamic scan, the SW technique was applied by inserting HF data acquired in 2 adjacent frames into the HF region of the k-space to be constructed. Figure 3 shows the schematic of SW in the HF region. In the first row, at each time point, k-space data were acquired using the variable density profile in Figure 2. For the k-space to be reconstructed at time point t₁ in the second row, the HF region was filled by HF data from t₀ , t₁ , and t₂ points, while central k-space region was filled with data only from t₁ central k-space data. To achieve a higher k-space coverage ratio, temporal constraints were added to Equation 1 to ensure that sampled HF regions have no overlap within each k-space triplet. The synthesized k-space slices were then used for image reconstruction.

Figure 3.

A schematic of sliding window (SW) with variable density Cartesian sampling profiles. Each black line represents a sampled readout line.

Image Reconstruction

Using each data sampling scheme, the undersampled data were then used for Cartesian MR reconstruction in an inverse problem with a penalty term of TGV theory:

${\text{argmin}}_u{\left|\left|F\left(u\right)-k\right|\right|}_2^2+ {\beta \text{TGV}}^2\left(u\right),$

2

where u represents the reconstructed image, F represents Fourier transform, k represents the undersampled k-space data, and β denotes the regularization parameter for the TGV penalty term. The TGV penalty term can be organized as ⁴⁴ :

${\text{TGV}}^2\left(u\right)= {\text{min}}_v \left\{\underset{\Omega }{{\displaystyle \int }}\left|\nabla u-\text{v}\right|\mathrm{dx}+ \alpha \underset{\Omega }{{\displaystyle \int }}\frac{1}{2}\left|\nabla v+\nabla {\text{v}}^T\right|\mathrm{dx}\right\},$

3

where the minimum is taken over all vector fields on Ω, and v is an introduced test function in the same order as the first-order derivative of the reconstructed image (∇u). At the neighborhood of edges on the image, the test function v approximates the first-order derivative of the image, and the first integration term in Equation 3 preserves the edge information as the optimal outcome. On the other hand, the symmetric derivative of v in the second integration term of Equation 3 improves the smoothness of the image based on the observation that second-order derivative (∇²u) is locally small within the smooth region of the image. The parameter α is used to keep a balance between the first- and second-order derivatives of the image. Thus, the TGV penalty term regularizes the features on both smooth and sharp region of the image. The nonlinear inversion in Equation 2 can be approximated using the iterative-regularized Gaussian-Newton method. ^{52 –54} At each step n for the given uⁿ , update the δu of the minimization problem:

${\text{min}}_{\text{δ}}{\left|\left|F^{\prime }\left(u^n\right)\text{δ}u+F\left(u^n\right) -k\right|\right|}^2+2{\beta }_k{\text{TGV}}^2\left(u^n+\text{δ}u\right),$

4

which can be estimated efficiently using the first-order primal dual algorithm. ^{55 –57} As in the reconstruction, no reference image was used. In Equation 2, α was set to 2 based on our experience for best brain imaging accuracy ⁴³ ; the initial value of β _k was set as 1, and its attenuation factor along iterations was empirically chosen as 5. The maximum iteration step was set to 8 in the consideration of a balance between image accuracy and computation cost.

Pharmacokinetics Analysis

Pharmacokinetics analysis with same procedures was conducted using the reference clinical image series as well as the image series reconstructed from undersampled k-space data with 2 difference sampling schemes, respectively. In order to quantify the CA concentration C(t), a linear relationship between the 2 quantities is assumed:

$\frac{1}{T_{\text{1}}\left(t\right)}= r_{1}C\left(t\right)+\frac{1}{T_{\text{10}}},$

5

where T ₁₀ is the preinjection T ₁ value, T ₁(t) is the changed T ₁ value at each time point after the injection, and r ₁ is the longitudinal relaxivity of CA at 1.5T magnetic field strength. Both T ₁₀ and T ₁(t) were calculated voxel-by-voxel using the standard dual flip angle method, ⁵⁸ and flip angles 5° and 15° were used. The CA concentration distribution at each time point was then calculated by Equation 5 and was processed by a denoising method. ⁵⁹ The CA concentration–time curve in the arterial blood C _p(t) (AIF) was measured in the way as described below: first, an arterial structure within the intracranial imaging volume was carefully contoured in 3D fashion by a clinical professional. Then a histogram of each arterial structure voxel’s time-to-peak (T _max) was generated. ⁶⁰ Next, the T _max of the C _p(t) was determined by the peak position of the generated histogram, and the C _p(t) was calculated by averaging the CA concentration of the voxels that constituted the histogram’s peak.

Two commonly used PK models were investigated in this study. First, the extended Tofts model was analyzed to quantify K ^trans, which has served as the primary quantitative imaging biomarker in many oncologic DCE-MRI studies. ^14,61 This model assumes that CA exchanges between the blood plasma and extravascular extracellular space (EES) and the CA concentration in tissue of interest C _t (t) can be expressed as ⁶² :

$C_t\left(t\right)=K^{\text{trans}}\underset{0}{\overset{t}{{\displaystyle \int }}}{C}_{\text{p}}\left(u\right)\cdot {\mathrm{e}}^{-k_{\text{ep}}\left(t-u\right)}\mathrm{du}+v_{\text{p}}\cdot C_{\text{p}}\left(t\right),$

6

where K ^trans describes CA extravasation rate from blood vessel to EES, k_ep = K^trans/v_e represents the rate constant of CA from EES back to blood plasma, and v _e and v _p denote the volume fractions of EES and blood plasma, respectively. The parameters were determined by a linear least squares-based curve fitting method. ⁶³

The 2-compartment exchange model was used to quantify the blood flow and blood volume. ^64,65 Following the conservation of CA mass, C _t (t) is expressed as:

$C_t\left(t\right)=F_{\text{B}}\underset{0}{\overset{t}{{\displaystyle \int }}}{C}_{\text{p}}\left(u\right)\cdot \left[{\mathrm{e}}^{-K_{+}\left(t-u\right)}+E\left({\mathrm{e}}^{-K_{-}\left(t-u\right)}-{\mathrm{e}}^{-K_{+}\left(t-u\right)}\right)\right]\mathrm{du},$

7

where F _B stands for blood flow, K _± were related to the reciprocal CA exchange rate constants between blood plasma and tissue of interest, and E is a parameter associated with K _± and circulation transit time T_T. Fitting Equation 7 with C_p(t) and C _t (t) by Levenberg-Marquardt algorithm with zero initial conditions, ⁶⁶ T_T can be calculated as:

$T_T{ }^{-1}= K_{+}-E\left(K_{+}-K_{-}\right),$

8

and blood volume V _B can be estimated as:

$V_{\text{B}}= F_{\text{B}}\cdot T_T.$

9

Quantitative Accuracy Analysis

The investigated PK parameters (K ^trans, F _B, and V _B) were calculated in a voxel-by-voxel pattern, and 2 sets of parameter maps were generated from the original clinical scan image (full sampling set) and the reconstructed images using undersampled data (reduced sampling set), respectively. For the purpose of quantitative accuracy evaluation, a 3D region of interest (ROI) was defined using gross tumor volume from the radiotherapy treatment planning system. As the most commonly reported metrics in clinical studies, the means of K ^trans, F _B, and V _B across the ROI were compared for relative differences. For each parameter, we used the “error in volume mean (EIVM)” for the PK parameters estimation accuracy ⁴⁰ :

$\text{EIVM}\left(P\right)= \frac{\stackrel{\_}{P_{\text{r}}\left(x,y,z\right)}- \stackrel{\_}{P_{\text{f}}\left(x,y,z\right)}}{\stackrel{\_}{P_{\text{f}}\left(x,y,z\right)}} \times 100\%,$

10

where the P _r(x,y,z) and P _f(x,y,z) denote the parameter’s distributions within the ROI from the reduced sampling set and the full sampling set, respectively, and the overhead bar represents the average across ROI. For voxel-based analysis, a histogram depicting the parameter’s calculation difference using reduced sampled data was generated. The parameter’s difference map within the ROI was also evaluated. The accuracy was further evaluated using the 2 quantitative metrics: total relative error (TRE) and cross-correlation (CC). The TRE is defined as follows ⁶⁷ :

$\text{TRE}= \frac{\sqrt{{{\displaystyle \sum }}_{x,y}{\left[P_r\left(x,y,z\right)-P_f\left(x,y,z\right)\right]}^2}}{\sqrt{{{\displaystyle \sum }}_{x,y}{\left[P_f\left(x,y,z\right)\right]}^2}}\text{ ,}$

11

where P _r(x,y,z) and P _f(x,y,z) follow the same definitions in Equation 10. The CC was calculated to measure the similarity of a certain parameter’s distributions from the 2 parameter maps sets:

$\text{CC}= \frac{{{\displaystyle \sum }}_i\left(P_{\text{r},i}- \overline{P_{\text{r}}}\right)\left(P_{\text{f},i}- \overline{P_{\text{f}}}\right)}{{\left({{\displaystyle \sum }}_i{\left(P_{\text{r},i}- \overline{P_{\text{r}}}\right)}^2\right)}^{1/2}{\left({{\displaystyle \sum }}_i{\left(P_{\text{f},i}- \overline{P_{\text{f}}}\right)}^2\right)}^{1/2}}\text{ ,}$

12

where the summation was calculated within the 3D ROI, and the overhead bar represents the mean value. To determine the significance levels of the differences between the PK parameter results using full sampling data and reduced sampling data, Wilcoxon signed rank test was adopted. Statistical significance was defined as α < .05 with multiple comparison correction.

All image reconstruction works and PK analysis were carried out in a software developed in-house using MATLAB (R2012b; MathWorks Inc, Natick, Massachusetts) on a personal desktop with 16 GB RAM and 3.4 GHz central processing unit. As for efficiency evaluation, computation time was measured using functions “tic” and “toc” in MATLAB.

Results

Radial-Based Sampling Scheme Results

Figure 4 shows a selected patient’s 2D comparison of the original clinical DCE image (A) and the reconstructed image using the undersampled data from the radial-based sampling scheme (B). The images were from a postinjection volume, and the red contour line indicates the ROI for quantitative analysis. As can be seen, (B) preserved most structural details in (A), though some details were lost at high-intensity region. The average TRE value of all patients’ reconstructed images using the radial-based sampling scheme was about 0.113. Figure 5 presents the comparison of measured AIF information, and C _p,f(t) and C _p,r(t) denote the AIF measured on the CA concentration maps derived from original clinical image set (blue) and the reconstructed image set (red), respectively. As can be seen, the peak positions of the 2 curves were identical, and the peak enhancements of C _p,f(t) and C _p,r(t) were 0.350 and 0.366 mmol/L, respectively. The shape features of C _p,f(t) were generally preserved by C _p,r(t), though C _p,r(t) had an overestimated intensity both before and after the first pass peak.

Figure 4.

A comparison of the original postinjection image (A) and the reconstructed image using the radial-based undersampled data (B). The intensity difference between (A) and (B) is shown in (C). The red contour indicates the region of interest (ROI) that contains the tumor.

Figure 5.

The comparison of the measured C _p(t) from the contrast agent (CA) concentration maps calculated from original clinical images (blue curve) and the images reconstructed from radial-based undersampled data (red curve).

Figures 6 to 8 demonstrate the comparison of PK parameter maps (K ^trans, F _B, and V _B) that correspond to the same imaging position shown in Figure 4. In each figure, the parametric maps calculated using the full sampling data and the reduced sampling data are presented in (A) and (B), respectively. The histogram (C) presents the distribution of the relative error of the 3D ROI parameter map in (B) versus (A). The PK parameters’ percentage difference map within the ROI is shown in (D). The results of K ^trans are presented in Figure 6. As can be observed, the 2 K ^trans maps were morphologically similar, and a K ^trans intensity-elevated region within the ROI (indicated by white arrows) was clearly identified on both maps. The heterogeneity of K ^trans intensity within the identified region in Figure 6A was successfully reproduced in Figure 6B. The histogram in (C) suggested an overall acceptable 3D K ^trans estimation accuracy using the reduced sampling data. The median value of absolute percentage difference in (C) was about 6.70%. The difference map in Figure 6D also suggested a good voxel-based quantitative calculation accuracy of K ^trans at the center of ROI, though minor discrepancies of K ^trans underestimation in Figure 6B were found at ROI peripheral region.

Figure 6.

A comparison of K ^trans map calculated using the original full sampling data (A) and using the radial-based undersampled data (B). The histogram about relative error of (B) versus (A) of all 3-dimensional (3D) region of interest (ROI) voxels (C). The histogram about relative error of (B) versus (A) of all 3D region of interest (ROI) voxels. (D). The percentage difference map of the 2D ROI between (A) and (B)

Figure 7.

A comparison of F _B map calculated using the original full sampling data (A) and using the radial-based undersampled data (B). The histogram about relative error of (B) versus (A) of all 3-dimensional (3D) region of interest (ROI) voxels (C). The histogram about relative error of (B) versus (A) of all 3D region of interest (ROI) voxels. (D). The percentage difference map of the 2D ROI between (A) and (B).

Figure 8.

A comparison of V _B map calculated using the original full sampling data (A) and using the radial-based undersampled data (B). The histogram about relative error of (B) versus (A) of all 3-dimensional (3D) region of interest (ROI) voxels (C). The histogram about relative error of (B) versus (A) of all 3D region of interest (ROI) voxels. (D). The percentage difference map of the 2D ROI between (A) and (B).

Figure 7 demonstrates the blood flow F _B results. On both Figure 7A and B, a circular region with elevated F _B was identified approximately at the K ^trans hot region shown in Figure 6A and B. Compared to the F _B-elevated region in Figure 7A, the identified region in Figure 7B had a minor difference at its right posterior boundary. The difference map in Figure 7D shows a good accuracy of F _B estimation in Figure 7B, and the left anterior region indicated a minor F _B underestimation in Figure 5A. The results of blood volume V _B are illustrated in Figure 8. As is the pattern indicated by K ^trans and F _B results, a region with elevated V _B values was indicated in Figure 8A and B. The V _B difference map in Figure 8D also exhibited the good V _B calculation accuracy using the reduced sampling data set.

As for the quantitative accuracy evaluation within the ROI in 3D fashion, Table 1 summarizes each PK parameter’s EIVM, TRE, and CC within the defined ROI of all 8 scans using the radial-based undersampling profile. The EIVM values demonstrated in Table 1 were all around 5% or less except one V _B EIVM value that indicated an 8.7% error. These results suggest that just as the commonly reported metrics in clinical studies, the mean values of the investigated PK parameters can be accurately calculated using the undersampled data. The TRE values varied from 0.08 to 0.15 with an average value of 0.11. As a description of the morphological similarity of 2 images, most CC values were higher than 0.95, and the observed smallest CC was 0.83. In the Wilcoxon signed rank tests, the statistical significance level of differences in K ^trans, F _B, and V _B were 0.528, 0.889, and 0.123, respectively. These results suggested that there is no significant differences between the PK maps calculated from the full sampling set and the one calculated from the reduced sampling set.

Table 1.

The Results of Error in Volume Mean (EIVM), Total Relative Error (TRE), and Cross-Correlation (CC) Within the 3D ROI Using Radial-Based Undersampled Data.

Patient No.	K trans			F B		V B
	EIVM	TRE	CC	EIVM	TRE	CC	EIVM	TRE	CC
1	−0.65%	0.108	0.951	0.78%	0.086	0.981	−0.76%	0.100	0.987
2	5.17%	0.110	0.977	−1.29%	0.098	0.971	−4.33%	0.108	0.941
3a	−1.53%	0.104	0.990	−3.00%	0.109	0.950	−0.12%	0.096	0.989
3b	0.85%	0.129	0.979	−0.59%	0.128	0.955	6.49%	0.101	0.978
4	2.66%	0.102	0.960	1.54%	0.075	0.966	6.72%	0.121	0.996
4b	−1.76%	0.129	0.926	−0.28%	0.131	0.967	−2.09%	0.134	0.981
5	5.90%	0.103	0.941	0.66%	0.097	0.827	−2.23%	0.105	0.943
6	<0.01%	0.133	0.944	−1.66%	0.086	0.944	−8.73%	0.146	0.943

Abbreviations: F _B, blood flow; V _B, blood volume; ROI, region of interest.

^aImages are selected from this scan.

^bSecond scan of this patient.

Cartesian-Based Sampling Scheme Results

Figure 9 shows the comparison of reference DCE image (A) and the reconstructed image using undersampled data from the Cartesian-based undersampling scheme (B) in the same organization as Figure 4. The DCE image in (B) preserves most structural details in (A), and the difference map in (C) indicates the discrepancies at boundaries. The average TRE values of all reconstructed images using the Cartesian-based undersampling scheme was about 0.109, and this number was very close to the corresponding number (0.113) when using the aforementioned radial-based sampling scheme. Figure 10 presents the comparison of measured AIF information from original clinical image (blue) and the reconstructed image (red), respectively. Similar to the results in Figure 5, the red curve has same peak position and fluctuation patterns as of the blue curve. The peak enhancement of the red curve was 0.382 mmol/L, which is more deviant from the 0.350 mmol/L value of the blue curve than the red curve in Figure 5 (0.366 mmol/L) with radial-based sampling scheme.

Figure 9.

A comparison of the original postinjection image (A) and the reconstructed image using the Cartesian-based undersampled data (B). The intensity difference between (A) and (B) is shown in (C). The red contour indicates the location of gross tumor volume (GTV) as region of interest (ROI).

Figure 10.

The comparison of the measured C _P(t) from the contrast agent (CA) concentration maps calculated from original clinical images (blue curve) and the images reconstructed from Cartesian-based undersampled data (red curve).

Figures S1, S2, and S3 in the supplementary material document summarize the K ^trans, F _B, and V _B results in the same organization as in Figures 6 to 8. For a short conclusion, all the 3 PK parametric maps that were reconstructed using the Cartesian-based undersampled data were very similar to the reference maps. For quantitative evaluations, each PK parameter’s EIVM, TRE, and CC of all the 8 scans are summarized as Table 2. The EIVM results in Table 2 were generally around 5% or less, but in some cases, the results were worse than the corresponding values in Table 1. The largest EIVM value in Table 2 is 16.67%. All observed CC values in Table 2 were higher than 0.9, and this result indicates the good morphological match of PK parametric maps using undersampled data. The average TRE values of all calculated PK parametric maps was about 0.122, and this number is higher than the corresponding result (0.111) using the radial sampling scheme. Additionally, no statistically significant results were found between the mean K ^trans (P = .950), F _B (P = .401), and V _B (P = .888) using reference data and using the undersampled data with Cartesian-based scheme.

Table 2.

The Results of Error in Volume Mean (EIVM), Total Relative Error (TRE), and Cross-Correlation (CC) Within the 3D ROI Using Cartesian-Based Undersampled Data.

Patient No.	K trans				F B		V B
	EIVM	TRE	CC	EIVM	TRE	CC	EIVM	TRE	CC
1	6.92%	0.135	0.982	5.11%	0.119	0.961	4.40%	0.106	0.993
2	1.40%	0.135	0.967	6.81%	0.103	0.949	1.98%	0.135	0.994
3a	2.01%	0.082	0.994	5.02%	0.110	0.977	0.80%	0.084	0.993
3b	4.55%	0.147	0.935	−3.98%	0.110	0.968	−6.49%	0.121	0.972
4	−1.45%	0.137	0.980	1.40%	0.081	0.923	−5.94%	0.092	0.995
4b	−4.41%	0.139	0.938	6.22%	0.111	0.994	−1.27%	0.152	0.983
5	−7.89%	0.108	0.990	−6.26%	0.162	0.763	−2.35%	0.149	0.965
6	2.81%	0.111	0.988	−3.96%	0.096	0.990	16.67%	0.155	0.995

Abbreviations: 3D, 3-dimensional; F _B, blood flow; V _B, blood volume; ROI, region of interest.

^aImages are selected from this scan.

^bSecond scan of this patient.

Discussion

In this pilot study toward potential clinical application, we investigated the feasibility of using undersampled k-space data for quantitative PK analysis in clinical brain DCE-MRI study. Particularly, we adopted a TGV penalty term in the iterative image reconstruction, and we analyzed the accuracy of PK parameters from 2 commonly used models using the undersampled data. Our results demonstrate that the investigated PK parameters can be accurately estimated using 2 data undersampling strategies at an acceleration factor of 4. In this part of the article, we further discuss several technical issues regarding the proposed methods in this work.

To ensure the accurate DCE image reconstruction, central k-space region needs to be sampled properly. Because of the superior central k-space sampling scheme, the radial sampling profile could use fewer rays than the Cartesian sampling profile to reach an acceptable central/peripheral k-space coverage balance. The golden ratio–optimized sampling profile is noteworthy. Mathematical deduction proved that when N = 32 rays were used, only 3 different azimuthal gaps occurred in comparison with the evenly distributed rays. ⁴⁶ Thus, this sampling profile allows a posteriori adjustment of the temporal resolution. ⁴⁹ In terms of sampling efficiency, this golden ratio–based profile has no difference with the uniform profile. ⁴⁸ On the other hand, to achieve a same acceleration factor using Cartesian grid sampling, SW technique becomes important to guarantee the adequate k-space coverage. The comparison of the numerical TRE and EIVM results in Tables 1 and 2 suggest that the accuracy of PK parameters using Cartesian-based undersampled data was generally and yet slightly worse than using the radial-based undersampled data. In this work, the radial sampling profile was simulated at the Cartesian grid. Using this approach, the potential error from data gridding/regridding was not included, and the difference in reconstructed images to the reference came from the image reconstruction process only. As a feasibility investigation of the TGV algorithm, this is a good demonstration of the algorithm’s mathematical accuracy. For potential real clinical application, however, the regridding process must be included and could be an important source error. In addition, radial sampling may have gradient delay, off-resonance effects, and strong streaks, which require additional physics works for imaging implementation. In contrast, the proposed Cartesian-based undersampling with spatiotemporal sampling density constraints is relatively easier for implementation and thus more promising for clinical application. The use of 2D sampling strategy is a good start as a feasibility study demonstrating the TGV application in DCE-MRI acceleration. As a natural extension, the 3D radial/Cartesian sampling strategy may further improve the acceleration factor with the k-space undersampling on 2 phase-encoding directions. ^68,69 Since the second-order derivative calculation involved in TGV would become computation intensive, the adoption of Graphics Processing Unit (GPU)-based reconstruction implementation may become necessary for future DCE-MRI acceleration in clinical use. ⁷⁰

As discussed previously, the TGV term can maintain edge information while enforcing the smoothness in low-contrast region. In Figure 4B, the sharp contrast between the major blood vessels and nearby soft tissue of (A) was preserved, and the homogeneity of soft tissue within ROI in (A) was successfully reproduced. At the same time, the blurring effect can be found at the proximity of some high-intensity region. Similar results can be found in Figure 9 with the Cartesian-based strategy. For comparison purpose, Figure 11 provides an example of DCE image comparison with and without the implementation of the investigated TGV reconstruction algorithm. As can be seen, the image quality improvement after adopting the investigated TGV reconstruction algorithm is obvious. Since SW was adopted in the Cartesian-based undersampling strategy to increase k-space coverage, the quality improvement from (E) to (D) may less obvious than the improvement from (C) to (B). As a quantitative image quality evaluator, the TRE of (C) and (E) was about 0.332 and 0.205, respectively. These 2 values were higher than the TRE values reported in Figures 4 and 9. Since the overall image quality in (C) and (E) are bad with unidentifiable key structures, the PK analysis based on the direct inverse Fourier transform image series becomes meaningless. Although the simulating acceleration factor of 4 in this study seems to be conservative in view of the recently reported studies with up to ×36 acceleration factors, ^69,70 the rigorous PK parameter accuracy evaluation in Tables 1 and 2 indicates the good acceptability of the current simulated acceleration of 4; on the hand, in spite of the appealing technical capabilities, more comprehensive evaluation of PK parameter calculation accuracy following the clinical standard is in demand before the further investigation of ultrafast DCE acceleration.

Figure 11.

The comparison of undersampled dynamic contrast-enhanced (DCE) image with and without total generalized variation (TGV) reconstruction. A, Reference full-sampled image; (B) TGV-reconstructed image using radial-based undersampling strategy; (C) direct inverse Fourier transform from undersampled data using radial-based undersampling strategy; (D) TGV-reconstructed image using Cartesian-based undersampling strategy; (E) direct inverse Fourier transform from undersampled data using Cartesian-based undersampling strategy.

In a balance of calculation accuracy and efficiency, we limited the maximum iteration step to 8 for DCE image reconstruction in this work, and this choice can be justified by the curves shown in Figure S4 in the supplementary material document. As the average results of 100 reconstructed images from different scans, the blue curve shows the time (in seconds) consumption of each iteration step, and the red curve summarizes the residue (normalized to 100) of the reconstructed images in logarithmic scale at each step. After the eighth step, the blue computation time curve reached a stable plateau, whereas the red residue curve indicated no significant reduction. The total reconstructed time of an image was about 90 seconds.

To improve the PK model fitting quality, an improved texture preserving variational denoising method was applied to the CA concentration maps to reduce the spatial noise effect. The same implementation was adopted for both fully sampled data sets and reduced sampled data sets. In the DCE-MRI acquisition, the image quality may be compromised as the acquisition speed increased. So far, the noise effect in CA concentration maps has not been fully characterized yet. As a typical approach of pixel-by-pixel analysis, the spatial averaging strategy has been reported in various clinical studies for PK analysis. ⁷¹ In addition, for certain PK parameters, the spatial filtering was needed to remove the noise effect. ⁶⁴ In the denoising method used in this study, the noise reduction was optimized with a spatially varying local power constraint. ^59,72 Compared to the classic variational denoising technique and moving average kernel, this method has been reported to improve the image quality in terms of signal-to-noise ratio (SNR) at no change in spatial patterns. ^59,73 Defined by the general shape and spatial distribution of a certain biomarker’s map, texture features are independent of the biomarker’s absolute value and are less sensitive to temporal and spatial resolution effects. ^{74 –76} It has been reported that texture features may have great potentials in monitoring the treatment response as a robust biomarker. ⁷⁷ As a future work, it would be interesting to investigate whether texture feature can be accurately reproduced as a robust biomarker with highly accelerated DCE-MRI studies.

In this study, we used the individualized AIF measurement for the PK analysis, and such measurement requires the existence of a major arterial structure within the imaging volume. Due to current technical limit, a limited volume of interest is covered with an adequate spatial resolution at a reasonable temporal resolution, and the lesion has to be close to a major intracranial arterial structure. Unfortunately, this requirement was hard to be satisfied, and thus, only a small number of patients were qualified for this pilot study. We acknowledge that 8 DCE-MRI data sets may not broad enough for a population study, but we think that this feasibility study with these 8 scans would be a good pilot study, and the current work is an early demonstration of the proposed method. Certainly, future works with a much larger number of data sets on various tumors are essential toward the clinical application of the proposed DCE-MRI acceleration method.

The inclusion of individual AIF measurement is an important factor of this work. It has been claimed that characterization of an individual AIF in preclinical studies is very difficult and could be misleading. ³⁷ For brain analysis, the AIF extraction based on time-to-peak information has shown its reproducibility in the longitudinal study. ⁷⁸ To maintain AIF’s enhancement feature, some studies used only a certain number or percentage of voxels based on the peak enhancement ratio in AIF generation. ^60,78 However, such a method is sensitive to the absolute CA concentration value, and the data reproducibility may be reduced when intensity-based image reconstructed methods were used. In this work, we selected a group of voxels based on the time-to-peak distribution and average the resultant time course to create AIF. To ensure the accuracy, this approach requires a large volume of arterial structure to generate a sufficient number of qualified voxels. In practice, however, the requirement was difficult to satisfy. The noise content of the CA concentration curve has been quantified as contrast-to-noise ratio (CNR) that is the ratio of maximum peak to the standard deviation of Gaussian random noise. ⁷⁹ The reported CNR value of CA concentration curves ranged from 4 to 20. ^64,80 In earlier studies, polynomial filters have been used in the model-free DCE-MRI analysis. ⁸¹ For the quantitative PK analysis, the similar filtering strategies may further improve analysis accuracy when using undersampled data, but additional future works are in demand to demonstrate the feasibility of such filtering strategy.

Conclusion

The feasibility of accelerated brain DCE-MRI has been successfully demonstrated in the quantitative PK parameters. As shown in this study, the quantitative parameters from 2 commonly adopted PK models can be accurately calculated from the images reconstructed with 2 data undersampling strategies at a simulated acceleration factor of 4. These results suggested that the investigated acceleration method can be used for reliable clinical image acquisition. Future works with a large population are in demand toward the exploration of this method before routine clinical implementation.