Impact of temporal resolution and motion correction for dynamic contrast-enhanced MRI of the liver using an accelerated golden-angle radial sequence

This study was based on free-breathing acquisitions from 11 patients which were reconstructed using the XD-GRASP technique with different number of motion states followed by several analyses (figure 1). This prospective study was approved by our institutional review board and consent was obtained from all patients. Eleven patients, 6 men and 5 women (70 years ±11 [standard deviation]) diagnosed with HCC, with tumors exceeding 10 mm, underwent DCE-MRI examinations between September 2018 and December 2018.

Zoom In Zoom Out Reset image size

Patients were scanned on a 3.0T clinical MRI system (Ingenia, Philips Healthcare, Best, The Netherlands) using the combination of a 16-channel anterior array and a 12-channel posterior array for signal reception.

To limit peristalsis, patients were asked to fast for four hours before the examination. T1 mapping was performed prior to contrast injection using the variable flip-angle method with fat-saturated Cartesian 3D gradient-echo acquisitions with flip angles of 4°, 8° and 12°. These acquisitions were performed in breath-hold using the same spatial resolution and field-of-view as the subsequent DCE acquisition.

The DCE acquisition was performed in free-breathing using a 3D stack-of-stars golden-angle gradient-echo radial acquisition with fat-saturation (Pahwa et al 2018). Thirty seconds after the beginning of the DCE acquisition, gadobenate dimeglumine (MultiHance®, Bracco Imaging SpA, Milan, Italy) was injected with an adjusted dose according to patient weight (0.1 mmol kg⁻¹ with a maximum dose of 20 ml) at 1.5 ml s⁻¹ followed by a 15 ml saline flush at 1.5 ml s⁻¹. Data was acquired continuously for 6 min. In total 50 slices were acquired with 5 mm spacing and then the acquired data was reconstructed to obtain 100 slices with 2.5 mm. The FOV was 250 mm. Imaging parameters for the DCE sequence are listed in table 1.

XD-GRASP is a novel method for dynamic image reconstruction based on a continuous scan during free breathing developed by Feng et al (2016). It is based on reconstructing the extra respiratory motion dimension from the acquired k-space data (figure 2). In the context of 3D liver imaging, this reconstruction is performed by solving the following optimization equation using compressed sensing:

$$\begin{equation}\hat x\,{\rm{ = }}\,{\rm{argmin||}}y\, - \,{\rm{FS}}x{\rm{||}}_2^2\, + \,{\lambda _1}{\left| {\left| {{C_1}x} \right|} \right|_1}\, + \,{\lambda _2}{\left| {\left| {{C_2}x} \right|} \right|_1}.\end{equation} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

F is the Non-Uniform Fast Fourier Transform (NUFFT) operator defined on the radial sampling pattern of each frame of x, where x is the 3D dynamic image-series and y is the corresponding multicoil radial k-space data sorted into several contrast phases. Furthermore, S represents the coil sensitivity maps estimated using Walsh method (Grimm 2015) from the multicoil reference image given by the coil-by-coil NUFFT (Fessler 2019) reconstruction of k-space data by averaging all the acquired spokes. Here, C₁ and C₂ are the sparsifying transforms with weighting parameters ${{\lambda }_1}$ and${{\,\,\lambda }_2\,\,}$applied in the contrast-enhancement dimension and along the extra respiratory state dimension, respectively. The same sparsifying transform 'Total Variation' constraint was used for both the respiratory dimension and the contrast enhanced phase but with different weighting parameters. The sparsifying transform is based on removing unwanted details while preserving important details such as edges, while the duration of each contrast-enhancement phase was 6 s. This transform plays an important role in removing streaking artifacts or the interspokes k-space information in frequency domain (Block et al 2007). The total number of spokes acquired was 691.

The coil elements used for signal reception may have different noise levels. Therefore the k-space data is obtained after a noise whitening step which decorrelates and normalizes the channel and thus suppresses the influence of the defective coil elements.

The MR pulse sequence used in this work represents Cartesian encoding over k_z and a stack of stars trajectory with radial sampling over k_x–k_y (Usman et al 2013) (figure 3). The radial sampling angle is the golden angle ∼111.25°. The sampling motion artifacts are averaged out by oversampling the center of k-space in the radial sampling. This oversampling may lead to an overweighting of low frequencies so a density compensation is applied where it serves the same aim as a ramp filter (Fessler 2019):

$$\begin{equation}d\left( {\emptyset ,{\rm{r}}} \right) = {\left| {\rm{r}} \right|^*}{1 \over {{N_a}}},\end{equation} \tag{ 2 }$$

where r is the k-space data and N_a the number of radial samplings.

Zoom In Zoom Out Reset image size

Finally, the coil sensitivity is an important parameter to be estimated in the pre-processing steps. The Walsh approach (Lustig et al 2007), using an adaptive matched filter for the multi-coil data to reduce the noise in the background and dark areas of the images, was adapted for these purposes. It is based on calculating the covariance matrices of the noise and signal at each image pixel where the eigenvector of these matrices provides an almost optimal coil sensitivity estimation. The individual coil images used for the sensitivity estimation are extracted from the radial data itself using a NUFFT.

The XD-GRASP reconstructions for each contrast enhancement phase was performed by sorting the acquired data into 4, 3, 2 and 1 motion states. A total 100 consecutive spokes with footprint of 25 spokes were used for each motion state and for 1 motion state (GRASP) where, single dynamic dimension represents the contrast enhancement disrupted by motion. For comparison purposes, reconstructions using 100, 75, 50 and 25 spokes for 4, 3, 2 and 1 motion states, respectively, were also generated. An effective temporal resolution of 6 s/volume, 12 s/volume, 18 s/volume and 24 s/volume, was associated to 1, 2, 3 and 4 motion states, respectively while the nominal temporal resolution was 6 s/volume. These reconstruction were compared to a baseline which is based on an approach without applying the compressed sensing constraint using Non-Uniform Fast Fourrier transform (NUFFT) reconstructions. It was implemented by means of convolutions with a Kaiser Bessel kernel. All images were reconstructed in MATLAB R2018a (The MathWorks, Inc. Natick, MA, USA) using the framework implementation by Feng et al.¹

For quantitative analysis, the contrast-enhanced liver MR perfusion image series was registered with MoCo (Corstem, Montreal, QC, Canada), a non-rigid motion compensation software using subpixel deformation maps. This software used non-rigid estimation of displacement based on an optical flow formulation adapted to dynamics of perfusion series (Brox and Malik 2011). This implementation allows to register MRI perfusion image series with varying breathing patterns (Jacobs et al 2016, Benovoy et al 2017). All the images were then registered to the first image used as the initial target.

In the T1 mapping, to handle the varying breath-hold, the same non rigid registration using the MoCO technique was performed between images with different flip angles.

The reconstructed 5D images obtained using XD-GRASP are composed of several respiratory motion states for each phase of contrast enhancement phase. To completely suppress these movements and their effect, even after their reconstruction, the motion compensation software was used to register all the motion states to an intermediate one for each phase of contrast enhancement. And finally, these motion states were averaged to obtain 4D (x-y-z-t) images.

The reference frame was chosen as one of the intermediate motion state to avoid registration to an extreme phase (end inhale or end exhale). In the case of four motion states, we chose the third one as reference and in case of three motions states, we chose the second one as the reference.

The optimal solution for registration is obtained by minimizing the following cost function:

$$\begin{equation}E\left( w \right) = {E_{\mathrm{int}}}\left( w \right) + \gamma {E_{\mathrm{grad}}}\left( w \right) + \alpha {E_{\mathrm{smooth}}}\left( w \right) + \beta {E_{\mathrm{match}}}\left( {w,{w_1}} \right) + {E_{\mathrm{desc}}}\left( {{w_1}} \right)\end{equation} \tag{ 3 }$$

where E_int indicates the pixel intensity energy term, E_grad is the gradient constraint, E_smooth is the regularization term, E_match is histogram oriented gradient (HOG) based landmark correspondence energy term and E_desc is a continuous formulation for discrete descriptor matching.

This equation is controlled with three parameters: β the weight of the HOG descriptor matching term, α flow fields smoothness and γ the gradient consistency term weight.

We did not perform any registration over time so the contrast enhancement was not affected. The registration was just applied between motion states.

For each set of images, a manual selection of four regions of interest (ROI) was performed during the vascular phases providing the best visualization of the portal vein, the abdominal aorta, the liver tumors and the normal-appearing liver. The abdominal aorta was identified after contrast arrival when it had the highest signal intensity whereas the portal vein was most clearly visible during the portal venous phase.

For the aorta, portal vein and liver, only one slice was selected. For the liver, the slice including the complete liver was chosen and for the aorta and the vein, the slice which best showed these structures according to a fellowship-trained abdominal radiologist was selected.

There was at least one tumor (HCC) in each patient confirmed by an experienced abdominal radiologist (AT) using previous diagnostic imaging. If a patient had multiple tumors, the segmentation was performed for each one. Each tumor was segmented on multiples slices which includes all slices where the tumor appears, to minimize sampling errors. Segmentations were performed by a third-year radiology resident (LB, with two years of experience), using MATLAB R2018a (The MathWorks, Inc. Natick, MA, USA).

To measure the average signal intensity versus time curves, each selected ROI was then extended to all motion compensated DCE-MRI time points.

In order to assess the image quality, we first calculated the signal-to-noise ratio (SNR) as a measure of image quality (Welvaert and Rosseel 2013). It is defined as the ratio of the average normal liver signal and the background standard deviation:

$$\begin{equation}\mathrm{SNR} = \frac{{{\mu _{\mathrm{liver}}}}}{{{\sigma _{\mathrm{background}}}}}.\end{equation} \tag{ 4 }$$

Then, the contrast-to-noise ratio (CNR) was calculated. It is defined as the ratio of the absolute difference between average normal liver signal and average tumor signal to the background standard deviation:

$$\begin{equation}CNR = \frac{{\left| {{\mu _{\mathrm{liver}}} - {\mu _{\mathrm{tumor}}}} \right|}}{{{\sigma _{\mathrm{background}}}}}.\end{equation} \tag{ 5 }$$

The background region was selected just outside the subject at the closet position to the liver. This region was selected in same position for all reconstructions of each patient. From patient to another, slight differences depending on the liver position may have occurred.

A visual scoring of image quality based on five point Likert scoring (used as standard in the hepatic imaging field) which was performed by an abdominal fellowship radiologist with five years of experience. The scoring scale use quality index: 1 (Extremely poor), 2 (Poor), 3 (Fair), 4 (Good) and 5 (Excellent), where index 1 and 2 are considered clinically unacceptable and 3 to 5 clinically acceptable (Martinez et al 2020).

We quantified the motion of the tumors between extreme phases (full inspiration and full expiration) of the XD-GRASP reconstruction with four respiratory motion states, before the motion compensation, through deformable registration using NiftyReg software (Modat et al 2010). A cubic B-splines-based transformation model with a pyramidal multi-resolution strategy at three different levels was adopted. The level of alignment between source and target volumes was quantified by Normalized Mutual Information. A regularizer term based on the bending-energy as proposed in (Rueckert et al 1999) was used.

In addition, we applied the same deformable registration after the motion compensation step using the third phase as the reference to estimate residual tumor displacement.

To measure the lesion motion range we look only at the tumor motion where dense deformation vectors were multiplied by a segmentation mask in order to consider only voxels belonging to the tumor. Then we calculated the three-dimensional Euclidean distance as follows:

$$\begin{equation}d = \sqrt {{{\left( {{\Delta _X}} \right)}^2} + {{\left( {{\Delta _Y}} \right)}^2} + {{\left( {{\Delta _Z}} \right)}^2}} \end{equation} \tag{ 6 }$$

Finally, we reported the lesion mean motion value over all temporal points. For visualisation purposes, the resulted images of this registration were computed with an open-source ITK module.

A linear interpolation to a 1 s time stamp was performed on the signal-time curves followed by a non-parametric analysis that directly calculated simple perfusion parameters correlating with the physiology of tissue. The parameters defined in table 2 were subsequently used to perform a direct comparison between reconstructed images with increasing number of motion states. These parameters are illustrated in figure 4 where TTP is the time that contrast agent need to reach its maximum value (TTP = T_maxC–T₀) and T₀ is the onset time representing the contrast agent arrival in tissue.

Table 2. Definition of perfusion measurements.

Measurement	Definition
TTP	Time that contrast agent needs to reach its maximum value
T0	Contrast agent arrival time in tissue
ΔS	Smax−S0, difference between maximum and baseline intensities
nMITR	$\frac{{\Delta \mathit{S}}}{{{\rm{TTP*}}{\mathit{S}_0}{\rm{}}}}$, normalized maximum intensity time ratio between TTP and peak enhancement
PER	$\frac{{\Delta {S}}}{{{{S}_0}}}$, peak enhancement ratio
Wash-out slope	$\frac{{{{S}_{{\text{max}}}} - {S_f}}}{{{{T}_{{\text{max}}}} - {\text{TTP}}}}$, slope between $S_{\text{max}}$ and $S_f$ points (final signal intensity)
Wash-in slope	$\frac{{\Delta {S}}}{{{\text{TTP}} - {{T}_0}}}$, slope between $S_{0}$ and $S_{\text{max}}$ points

TTP = time to peak, T₀ = onset time

Zoom In Zoom Out Reset image size

T1 mapping linear conversion was used to convert the signal intensity for each previously selected ROI at different times into a contrast agent concentration. The linear relation between 1/T1 and the contrast agent concentration is considered in tissue and blood using the following equation:

$$\begin{equation}\frac{1}{{{\text{T}}1{\text{post}}}} = \frac{1}{{{\text{T}}1{\text{pre}}}} + {r}1\left[ {{\text{Contrast Agent}}} \right].\end{equation} \tag{ 7 }$$

where r1 is the relaxation of T1.

Patient-specific B1 shimming (two transmit channels, phase and amplitude) was performed, but no additional correction based on the acquisition of a B1 map was performed, primarily for acquisition time constraints. B1 inhomogeneities could add to inter-lesion variability in the DCE parameters, but would most likely not affect the direct comparison between the different reconstructions performed from the same raw data and thus under the same B1 conditions.

The dual-input single-compartment model, a direct approach normally used for parametric analysis, was adjusted to experimental tumor data and the liver tissue surrounding it. It needs an accurate assessment of the arterial input function (AIF) and the portal venous input function (PIF).

After linear interpolation, the ROI contrast agent concentration in the portal vein was used to estimate the portal venous input function (PIF). And the arterial supply extracted from the drawn ROI on the abdominal aorta (celiac trunk) and its contrast agent concentration were used for the estimation of the arterial input function (AIF) (Hormuth et al 2014).

From this model, we had the ability to extract the transfer constant from the liver tissue to the central vein K₂ and the transfer constants from the arterial plasma to the surrounding extravascular space k_a and the transfer constants from the portal venous plasma to the surrounding tissue k_p (Breuer et al 2005). Using the relations developed by Materne et al (2002), we estimated the tumor arterial fraction (ART), where the arterial blood flow was calculated and expressed in ml · min⁻¹ · 100 ml⁻¹ and then as a fraction (%):

$$\begin{equation}{\text{ART}}\left( {\text{% }} \right) = 100{\,{{*}}}\,\frac{{{{{k}}_{{a}}}}}{{{{{k}}_{{a}}} + {{k}_{p}}}}\end{equation} \tag{ 8 }$$

and the distribution volume (DV) (Materne et al 2002), using:

$$\begin{equation}{\text{DV}}\left( \% \right) = 100{\,*\,}\frac{{{k_a} + {k_p}}}{{{k_2}}},\end{equation} \tag{ 9 }$$

which estimates the ratio of arterial and portal plasma over the central vein blood flow.

For this study, statistical analysis were performed in Excel, using ANOVA test between XD-GRASP with different number of motion states (1, 2, 3 and 4), and using paired t-test between XD-GRASP with 1 and 4 motion states. In both tests, p < 0.05 was considered as statistically significant.

Patients and Tumors Characteristics: table 3 summarized demographic descriptive statistics (sex and age) and clinical data, including the tumor characteristics and the etiology of liver disease. Variables of categories were represented as percentages and numbers and the rest as minimum, maximum and mean ±standard deviation (SD).

Parametric and Nonparametric Analyses: comparisons of mean values for perfusion parameters states as well as the comparison of the p-value of contrast-to-noise ratio (CNR) and signal-to-noise ratio (SNR) was made for XD-GRASP with different numbers of respiratory motion. The tumors were then divided into two categories, the first representing the tumor displacement greater than 1 cm while the second representing a displacement less than 1 cm. Then, the mean values of perfusion data parameters for each category were calculated and a comparison was made between 1 motion state and 4 motion states.

In this pilot study, eleven patients totaling 22 HCCs with tumor diameters between 11 and 52 mm were included. Figure 5 shows the patient selection flowchart and table 3 represents the patients and tumors characteristics. 169 patients diagnosed with HCC and awaiting their clinical MRI between September 2018 and December 2018 were eligible. 11 patients were selected for additional MRI acquisition using the proposed protocol. The most common reasons for excluding patients were related to MRI scheduling (timeline too short to schedule MRI, MRI already done, postponed MRI or impossible to book MRI), distance from the hospital and fasting too long. Thus, a total of 11 patients were included.

Zoom In Zoom Out Reset image size

In order to analyse the performance of the technique and to assess the image quality from the reconstructed images, we calculated SNR and CNR.

XD-GRASP reconstructed with increased motion states with different number of spokes improves the contrast-to-noise ratio and the signal-to-noise ratio (p < 0.05) but reduces temporal resolution (0.04 volume/s vs 0.17 volume/s for one motion state) (p < 0.05) (figure 6). XD-GRASP with 3 and 4 motion states had higher overall image quality scores compared to XD-GRASP with 1 and 2 motion states as it showed reduced streak artefact (figure 8). Results for image reconstructed with different number of spokes show that on average, the reconstruction with 4 phases received the higher scores (4.6) while the reconstruction with one phase received the lower one (1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the case of XD-GRASP reconstructed with 100 spokes, there was no significant difference in SNR and CNR and the visual scoring show an average score of 4 for all cases (4, 3, 2 and 1 motion state). Since the reconstruction had a footprint of 25 spokes, the nominal temporal resolution was the same 6 s/volume while the effective one was as discussed previously.

Compared to the baseline reconstruction, we can clearly see the importance of XD-GRASP with regards to quality image as shown in figure 7. This reference reconstruction was evaluated using the same ROI and same length of time points to generate the corresponding perfusion curves and compared it to GRASP and XD-GRASP. There is a significant difference between the baseline and the images reconstructed with XD-GRASP, particularly concerning image quality, and shows higher peak of aorta concentration, but is more noisy and unstable.

Zoom In Zoom Out Reset image size

Figure 6 shows that increasing the number of respiratory motion states improves the SNR and CNR but reduce the temporal resolution thus influencing quantitative parameters that capture rapid signal changes. This is confirmed with figure 8, which clearly shows the improvements of the image quality by increasing the number of motion states for several patients.

Figure 9(a) compares the diaphragm boundary deformation between the 4 motion states of the XD-GRASP reconstructed image. The diaphragm boundaries are marked by the white line. The results show that there is a clear motion between different phases which was confirmed by the deformation vector field between the first and last motion state represented in figure 9(b).

Zoom In Zoom Out Reset image size

The mean Euclidean distance calculated before registration for each tumor varies between 1.3 mm to 31 mm (mean 11.9 ± 8 mm) depending on the respiration of each patient as shown by the blue bars in figure 9(c). Following motion compensation, tumor target displacements varied between 0.2 and 4 mm (mean 0.9 ± 1 mm), as shown by the orange bars in figure 9(c).

Figure 10 represents global enhancement curves for the portal vein, the liver, the abdominal aorta and the tumor from data reconstructed with increasing number of motion states (1–4) and same number of spokes. The HCC enhancement clearly occurs before the liver's one in the case of 1 and 2 motion states.

Zoom In Zoom Out Reset image size

Comparison of all perfusion parameters between different numbers of phases is summarized in table 4 as mean ± SD. Time to peak (TTP) did not significantly change with decreasing number of motion states. The normalized maximum intensity time ratio (nMITR), the peak enhancement ratio (PER) (p < 0.05), and the tumor arterial fraction increased with decreasing number of motion states (p < 0.05) while the transfer constant from the portal venous plasma to the surrounding tissue (K_p) significantly decreased (p < 0.05).

Table 5 as well summarized the perfusion parameters between 4 motion states and 1 motion state depending on the tumor displacement. For displacement greater than 1 cm, the normalized maximum intensity time ratio (nMITR), the peak enhancement ratio (PER), the tumor arterial fraction (ART) and the transfer constant from the portal venous plasma to the surrounding tissue (K_p) significantly change (p < 0.05). On the other hand, in case of tumors displacement less than 1 cm, it was clear that there was no significant changes in all parameters except for the tumor arterial fraction.