3D nonrigid motion correction for quantitative assessment of hepatic lesions in DCE-MRI.

PURPOSE: To provide nonrigid respiratory motion-corrected DCE-MRI images with isotropic resolution of 1.5 mm, full coverage of abdomen, and covering the entire uptake curve with a temporal resolution of 6 seconds, for the quantitative assessment of hepatic lesions.
METHODS: 3D DCE-MRI data were acquired at 3 T during free breathing for 5 minutes using a 3D T1 -weighted golden-angle radial phase-encoding sequence. Nonrigid respiratory motion information was extracted and used in motion-corrected image reconstruction to obtain high-quality DCE-MRI images with temporal resolution of 6 seconds and isotropic resolution of 1.5 mm. An extended Tofts model was fitted to the dynamic data sets, yielding quantitative parametric maps of endothelial permeability using the hepatic artery as input function. The proposed approach was evaluated in 11 patients (52 ± 17 years, 5 men) with and without known hepatic lesions, undergoing DCE-MRI.
RESULTS: Respiratory motion produced artifacts and misalignment between dynamic volumes (lesion average motion amplitude of 3.82 ± 1.11 mm). Motion correction minimized artifacts and improved average contrast-to-noise ratio of hepatic lesions in late phase by 47% (p < .01). Quantitative endothelial permeability maps of motion-corrected data demonstrated enhanced visibility of different pathologies (e.g., metastases, hemangiomas, cysts, necrotic tumor substructure) and showed improved contrast-to-noise ratio by 62% (p < .01) compared with uncorrected data.
CONCLUSION: 3D nonrigid motion correction in DCE-MRI improves both visual and quantitative assessment of hepatic lesions by ensuring accurate alignment between 3D DCE images and reducing motion blurring. This approach does not require breath-holds and minimizes scan planning by using a large FOV with isotropic resolution.

INTRODUCTION

Dynamic contrast‐enhanced MRI (DCE‐MRI) of the liver is an approach yielding structural and functional information regarding tissue microvasculature. This is especially important in oncology, as it allows one to differentiate between healthy liver parenchyma as well as benign and malignant lesions.1 For DCE‐MRI in clinical practice, a series of T1‐weighted MR images are acquired at different timepoints (TA) during intravenous injection of a paramagnetic contrast agent. In addition to assessing the contrast uptake visually, this technique allows one to calculate quantitative maps of pharmacokinetic parameters related to perfusion,2 such as the volume transfer constant (ktrans), which is associated with endothelial permeability, the volume fraction of blood plasma (ν), or the extravascular extracellular fractional volume (ν).3, 4 These DCE maps are commonly obtained by fitting a pharmacokinetic model on a voxel‐by‐voxel basis.5, 6

To date, organ motion is one of the main challenges for the acquisition of accurate DCE images of the liver. In particular, respiratory motion can lead to substantial motion artifacts within each T1‐weighted DCE‐MR image and to misalignment between different images. This misalignment is especially problematic for the calculation of quantitative DCE maps, as the voxel‐by‐voxel fitting might be carried out over different voxels for each TA. To minimize respiratory motion artifacts, breath‐holding techniques have been used commonly in clinical setting for 3D DCE‐MRI. Often patients are asked to hold their breath at multiple instances for each TA for a predefined duration (10‐15 seconds).7, 8, 9 Although breath‐holding can ensure high image quality for each TA, it can still suffer from misalignment between DCE images, as each breath‐hold might be carried out at a slightly different breathing position. Different postprocessing techniques have been suggested to minimize this misalignment by registering each volume in the time series to a reference volume. This approach, however, has limitations when applied to DCE, as the signal in each voxel varies over time not only due to breathing motion, but also due to contrast uptake,10, 11, 12 therefore introducing a bias that can be misinterpreted by registration algorithms. Additionally, breath‐holding‐based techniques are especially limited in patients with limited respiratory reserve (e.g., elderly, postoperative patients, or patients with compromised cardiopulmonary functions), and the total number of acquired images is restricted.13, 14 Hence, the number of timepoints that can be sampled along the contrast uptake curves is also limited, providing a poor temporal resolution for the fitting of the pharmacokinetic models. Free breathing DCE‐MRI techniques have been proposed primarily for patients who have difficulties in complying with breathing instructions (such as pediatric patients), using soft‐gating or radial imaging and compressed sensing to minimize respiratory motion artifacts.13, 14, 15, 16, 17 In general, these approaches suffer from residual motion artifacts and are associated with limitations in terms of through‐plane resolution and SNR.

Recently, Feng et al proposed an approach in which DCE‐MRI data are acquired continuously during free breathing using a golden‐angle stack‐of‐stars scheme and then retrospectively sorted into different respiratory motion states.15, 18 Each motion state is then separated into multiple DCE‐MR images at different TAs. Motion and TA‐resolved 3D images are reconstructed directly from these bins using iterative reconstruction schemes and exploiting sparsity along the motion dimension. This method minimizes respiratory motion artifacts while providing accurate DCE‐MR information. Other approaches using nonrigid motion correction of each dynamic phase have also been proposed.19, 20 Nevertheless, so far these techniques have only been demonstrated for a few timepoints TA with an intermediate slice resolution. In addition, the calculation of quantitative DCE maps has not been evaluated based on these methods.

In this study, we propose a novel approach that provides motion‐corrected DCE‐MRI images with full abdominal coverage, isotropic spatial resolution of 1.5 mm, and covering the entire contrast uptake curve with a temporal resolution of 6 seconds. All data are acquired continuously during free breathing by using a golden‐angle radial phase‐encoding (GRPE) acquisition scheme. In the proposed technique, DCE images and deformable nonrigid respiratory motion information are obtained directly from the same data, without the need for an additional scan. Motion correction is carried out iteratively during image reconstruction. As a result, quantitative DCE maps can be calculated directly from the motion‐corrected DCE images.

METHODS

The approach investigated in this study can be divided into 3 main steps, as shown in Figure 1. First, continuous data acquisition during free breathing is carried out with a 3D‐GRPE sampling trajectory. Second, the acquired data are retrospectively binned into multiple respiratory motion states (Nresp) along the breathing cycle, which are predicted according to a respiratory self‐navigator extrapolated from central k‐space lines. The retrospective binning carried out in this approach ensures that the same amount of data is used to form each motion state and that it is taken evenly throughout contrast uptake. This step averages the intensity variations induced by contrast uptake evenly throughout in each motion state, which would otherwise introduce a bias in the image registration procedures. 3D images describing different respiratory motion states are reconstructed with a non‐Cartesian iterative SENSE reconstruction with spatial and temporal regularization. 3D nonrigid respiratory motion of each state with respect to a reference state in the breathing cycle (end‐expiration) is then estimated using a nonrigid spline‐based image registration with normalized mutual information metric. This provides Nresp nonrigid transformations or motion fields (MFs), which are then directly incorporated into the final image reconstruction of 3D dynamic motion‐corrected (MC) images. Data sets without motion correction (NMC) are also reconstructed for comparison.

Figure 1

Magnetic resonance continuous acquisition during free breathing using a 3D golden‐angle radial phase‐encoding (GRPE) k‐space sampling trajectory (“Acquisition”). Data are retrospectively binned in Nresp respiratory motion states that characterize the patient's breathing cycle, according to a respiratory self‐navigator. After reconstruction of respiratory‐resolved images, nonrigid motion fields (MFs) are obtained using spline‐based image registration (“Respiratory Data Binning and Motion Estimation”). The MFs are incorporated directly into the final image reconstruction of the dynamic motion‐corrected 3D images (“Dynamic Motion Correction”)

Data acquisition

A 3D fat‐suppressed T1‐weighted gradient‐echo acquisition is carried out continuously for 4 minutes during free breathing with a GRPE sampling trajectory.21 This trajectory combines Cartesian frequency encoding with a radial‐like sampling in the phase‐encoding plane, where the individual phase‐encoding steps are sampled along radial lines. Successively acquired lines are separated by the golden angle (111.25º). This strategy ensures a homogenous sampling of the phase‐encoding plane through time, flexibility regarding temporal retrospective data binning, and the possibility to derive a respiratory self‐navigator directly from the acquired data.22, 23 Characteristic imaging parameters are TR/TE = 3.3 ms/1.36 ms, flip angle = 12°, FOV = 288 × (288 – 345) × (288 – 345) mm3, isotropic spatial resolution = 1.5 mm, and partial Fourier factor = 5/8. A total of 120 k‐space points along each GRPE line are acquired with a total number of 640 GRPE lines. Magnetic resonance contrast injection always takes place during the acquisition, approximately 1 minute after the start of the sequence.

A 2D fat‐suppressed T2‐weighted turbo spin‐echo sequence is also acquired for further characterization of lesions (TR/TE = 6427 ms/100 ms).

Retrospective respiratory data binning and motion estimation

After the GRPE sequence is acquired, a full 1‐dimensional projection of the signal along the foot–head direction is reconstructed from the central k‐space spokes. These are acquired once every 120 TRs, yielding 640 samples over the entire acquisition and therefore yielding a temporal resolution for the navigator of 0.5 seconds. After 1‐dimensional Fourier transform, these projections showed organ motion along the foot–head direction. This procedure is repeated for each receiver coil separately. These signal changes are caused primarily by respiratory motion but are also due to the injection of the MR contrast agent in the imaging volume, resulting in changes of the main signal. However, respiratory motion and contrast agent dynamics show very different temporal behavior and hence can be separated with a principal component analysis. Nevertheless, depending on the area of interest (i.e., depending on the receiver coil), either respiratory or contrast changes are the dominating components; hence, simply always taking the first or second phase contrast (PC) is not sufficient to reliably extract a respiratory self‐navigator. To ensure that the correct PC is extracted, we use the following approach in all patients:

Reconstruction of 1‐dimensional projection from central k‐space spokes for each receiver coil;

Application of principal component analysis to each coil separately, to obtain the PC of signal changes;

Transformation of each PC into the frequency domain (frequency range: ±1.1 Hz) by applying a Fourier transform;

Identification of PC with the highest respiratory contribution (PCresp), achieved by looking at the frequency spectrum of each component and selecting the one with the greatest area under the curve in the frequency range of 0.1‐0.5 Hz24; and

Application of principal component analysis along the PCresp of each coil, to extract 1 respiratory self‐navigator.

The GRPE data are then retrospectively separated into Nresp bins according to the amplitude of the respiratory self‐navigator. This is achieved in a sliding window fashion such that each bin contains the same amount of k‐space data.

Respiratory motion estimation

Once the GRPE data have been reassigned to the Nresp different motion states, a non‐Cartesian iterative sensitivity‐encoding approach25, 26, 27 is used to reconstruct for each respiratory motion state (including the reference) an image , as follows:

where the terms TV = ||∇ 1 and TV = ||∇ 1 are the spatial total variation within each motion state and temporal total variation across different but confining motion states; ∇ and ∇ represent the spatial and temporal gradients, respectively (these terms minimize aliasing artifacts, which in turn allow for a more precise registration procedure to the reference image); λ and λ represent the spatial and temporal regularization parameters, respectively, and are used to tune regularization strength to reach the desired equilibrium between data consistency and minimization of undersampling artefacts; K is the k‐space data present in the bin of interest; and E is the encoding operator that applies the Fourier transform and gridding operations while encompassing information on the coil sensitivities. To minimize reconstruction times, coil compression to 6 channels was carried out prior to image reconstruction. The regularization parameters λ and λ were manually optimized once and then kept constant for all patients. Equation 1 was solved using a nonlinear conjugate gradient approach with a backtracking line search.25

Next, these images are registered through a 3D nonrigid spline‐based image registration algorithm based on normalized mutual information.28 The MFs resulting from this procedure are then used in the final image reconstructions to transform all data to the reference motion state (i.e., end‐expiration).

Motion‐corrected DCE‐MRI reconstruction

A motion‐corrected non‐Cartesian iterative kt‐SENSE image reconstruction algorithm29, 30 is used to obtain the final 3D dynamic DCE‐MR images. The kt‐SENSE reconstruction uses signal similarities along the temporal dimension to minimize undersampling artifacts. Radial phase‐encoding undersampling artifacts show high incoherence over time, as each dynamic phase is reconstructed with a different radial phase‐encoding pattern. The contrast agent dynamics, however, lead to slower intensity variations; therefore, the 2 signals can be separated in xf‐space. This separation is achieved using previous information obtained from low‐spatial‐resolution training data (M), which have the same temporal resolution as the final dynamic DCE‐MRI. For GRPE, M can be obtained by applying a low‐pass Gaussian filter in the phase‐encoding plane and reconstructing a set of dynamic DCE‐MR images with a standard iterative SENSE reconstruction.31 In this study we have set sigma of the Gaussian filter to 0.1, leading to a spatial resolution of approximately 8 mm for M. The training data are then used to regularize the iterative image reconstruction in xf‐space, which reconstruct the dynamic image series () at once in the following manner:

where represents the encoding operator with an additional temporal Fourier transform , K represents the k‐space data, and is the regularization parameter. In this equation, is given in xf‐space, and the final diagnostic DCE images in xt‐space are obtained by applying to the result of Equation 2. The retrospective binning of the data to obtain is carried out in a sliding window fashion to yield 48 DCE phases, each formed by 3134 GRPE points. The sliding window approach was chosen to have an overlap of 50% (i.e., 1567 GRPE points), leading to a temporal spacing between neighboring DCE phases of 6 seconds. The motion correction is then introduced directly in Equation 2 through the encoding operator E as follows:

where is the motion field relative to the motion state j; is the coil sensitivity map for each coil c; is the Fourier transform; and is a logical operator that collects only k‐space data acquired in the motion state j. This means that the data for each dynamic time frame are split into all motion states and transformed to the reference motion state. For comparison, kt‐SENSE reconstruction was also carried out without respiratory motion correction.

Similar to the reconstruction of the respiratory‐resolved images, coil compression using 6 channels was applied. The strength of the regularization parameter was determined manually to achieve a good tradeoff between artifact suppression and temporal accuracy of the DCE signals. Image reconstruction was carried out with a preconditioned conjugate gradient approach.

The reconstruction algorithm was written in MATLAB (R2016b, The MathWorks, Natick, MA) and Python (Python Software Foundation). It was run in parallel on a 24‐core workstation equipped with 256 GB of memory. On average, the MC reconstruction of a 48‐phase DCE time series took 5 hours to complete versus the 40 minutes needed for the reconstruction without motion correction.

Patient population

The study was prospectively approved by and registered with the local institutional review board. All data were acquired on a 3T Biograph mMR hybrid scanner (Siemens Healthcare, Erlangen, Germany). A total of 11 oncological patients undergoing whole‐body PET‐MR examination were included and provided written, informed consent. Five patients who were suspected or known to have hepatic lesions received 0.01 mmol/kg of gadoxeate disodium, while the remaining 6 patients were examined for reasons primarily not related to liver imaging of metastases and were injected with 0.1 mmol/kg of gadobutrol. The former contrast agent is hepato‐specific and used solely in MR liver imaging of focal lesions, whereas the latter is extracellular and used in a wide range of MRI applications including hepatic imaging. All contrast injections were administered manually during the scan.

Quantitative DCE maps

Signal intensity curves (S(t)) are extracted from the dynamic 3D‐DCE images into contrast uptake curves (C(t)) using a conversion approach suggested by the Quantitative Imaging Biomarkers Alliance,32 which was further developed by Medved et al33 for heavily T1‐weighted sequences as follows:

where is the time‐dependent concentration of contrast agent expressed in mM; S(t) and S(0) are the time‐dependent and precontrast signal intensities, respectively; r is the relaxivity of the contrast agents (5 mM‐1 s‐1 for gadobutrol and 6.2 mM‐1 s‐1 for gadoxetate disodium) at a field strength of 3 T34; T1ref is the native precontrast T1 value of a reference tissue (809 ms was chosen for the liver parenchyma)35; and S(0)ref is the precontrast signal intensity value of the same reference tissue calculated in a region of interest (ROI).

After the conversion from signal intensity to contrast agent concentration, a standard model such as the extended Toft model was applied,5, 6 to derive voxel‐wise maps of the volume transfer constant (ktrans) as follows:

where is the temporal evolution of the contrast uptake curve for any given voxel in the tissue of interest; represents the blood plasma contrast uptake curve as obtained from the arterial input function (AIF) measured in the hepatic artery; is the volume fraction of blood plasma; is the extravascular extracellular fractional volume; ΔT is the time delay of the tissue enhancement relative to the dynamic contrast agent concentration in arterial plasma ; and * is the convolution operator. The hepatic artery was used as source for the input function to fit the parametric maps.

Quantitative evaluation

To assess the effect of respiratory motion correction on the visualization of hepatic lesions, contrast‐to‐noise ratio (CNR) at the border between the inspected lesions and their surrounding tissue for each data set was determined as follows:

where represents the mean intensity value measured inside an ROI enclosing the lesion; and is the mean signal intensity value in an ROI that directly borders this lesion, whereas represents its SD. The mean CNR () and SD () of the motion‐corrected and uncorrected data set were then calculated. Statistical significance in the difference of values was assessed through a paired t‐test after ensuring the normality assumption of distribution of values through the 1‐sample Kolmogorov‐Smirnov test. Differences were considered to be significant for p < .05 in all tests. The CNR evaluation was carried out on the latest phase of the reconstructed dynamic data sets (approximately 4 minutes after contrast injection) and on the ktrans parametric maps fitted from the pharmacokinetic model. In addition, the motion of each lesion examined in the former evaluation was calculated as the RMS displacement (RMSD), derived from the amplitude of the MF that characterizes the transformation between the end‐inspiration and end‐expiration state.

RESULTS

Evaluation of motion correction

Figure 2 shows a coronal slice of the reference motion state together with 4 different respiratory‐resolved motion states (Figure 2A). The difference images between the motion states and the reference state (Figure 2B), and between the transformed images of each motion state to the reference state and the reference state itself (Figure 2C), show that the obtained MFs accurately transform to and minimize anatomical differences between the motion states. The value of shows strong differences seen especially around small features such as edges, the heart, and vessels, which increase from I1 to I7 as the states being examined become more distant from the reference end‐expiration state. For image differences can be seen primarily around the heart due to uncompensated cardiac motion. Respiratory movement concerning the liver in general and the central part of the abdominal area appears to be effectively corrected for.

Figure 2

Depiction of 4 respiratory‐resolved motion states and reference state (A) and their absolute difference (B). C, Absolute difference between the motion states transformed to the reference state by using the estimated MFs

Supporting Information Figure S1 and Figure 3 depict the NMC and MC coronal sections of the liver, taken from a patient injected with gadobutrol and a patient injected with gadoxetate disodium, respectively, showing dynamic contrast uptake. The patient in Supporting Information Figure S1 is exemplary of a healthy liver and shows how the contrast uptake evolves in time in healthy tissue and in the vasculature. Figure 3 is an example of a liver affected by neuroendocrine tumors, which can be clearly identified (arrows) in the late stages of the contrast uptake. In the uncorrected images, respiratory motion strongly impairs the visibility of the vasculature and neuroendocrine tumor lesions. Motion blurring of the kidneys and at the edge of the liver can also be seen. The proposed respiratory MC approach compensates for this effect and ensures accurate depiction of the anatomy and pathologies leading to sharp edges and clear visibility of vasculature and neuroendocrine tumor lesions, which is maintained throughout the dynamic sequence.

Figure 3

Contrast dynamics of 8 temporal timepoints out of the 48 DCE volumes for patients injected with gadoxetate disodium. For visibility, only a part of the FOV is shown. Respiratory motion correction (MC) reduces motion‐induced blurring (red arrows) of the anatomy and pathology. In addition, some of the small lesions only become distinguishable from healthy liver with MC (blue arrows). Abbreviations: NMC, no motion correction

Temporal profiles in 6 patients of tumors (gadoxetate disodium patients, Figure 4A,B) and vessels (gadobutrol patients, Figure 4C) are shown in Figure 4 for uncorrected and corrected images. In each dynamic DCE image, a line along the foot–head direction was selected and plotted over time, showing the contrast dynamics. Blurring induced by respiratory movement reduces the detectability of small features. This can lead to the apparent disappearance of features for the entire (Figure 4A,C) dynamic image sequence. In addition, small features in close proximity cannot be distinguished anymore (Figure 4B,C). Motion‐corrected images generally show well‐defined feature boundaries and contrast uptake in both space and time.

Figure 4

Temporal dynamics of signal intensity of lesions (gadoxetate disodium) and vessels (gadobutrol). For the temporal plots, a line of voxels was sampled along the foot–head direction (i.e., along the main respiratory motion direction) in 6 different patients. A,C, Blurring induced by respiratory movement led to the apparent disappearance of features for the entire dynamic image sequence. B,C, Small features in close proximity cannot be distinguished anymore. Motion‐corrected images generally show well‐defined feature boundaries and contrast uptake in both space and time

The quantitative CNR evaluation of the MC and NMC data sets was carried out on 19 lesions in the latest dynamic phase that were chosen across all patients affected by hepatic metastases (Table 1). The CNR of the MC data set showed a statistically significant increase, yielding a of 1.46 ± 0.51 against of 1.09 ± 0.33, therefore scoring an overall average percentage increase of 47% ( ). Table 1 summarizes the RMSD of each investigated lesion and a mean RMSD over all lesions () of 3.82 ± 1.11 mm.

Table 1

Contrast‐to‐noise ratio of the last DCE phase and displacement statistics

Lesion No.	CNR_NMC	CNR_MC	RMSD(mm)
1	1.69	2.39	3.01
2	1.32	1.90	2.06
3	1.27	1.72	3.01
4	1.48	2.03	2.34
5	0.45	2.08	5.79
6	1.15	1.44	4.72
7	2.25	3.59	5.49
8	1.76	3.06	4.19
9	2.32	3.15	2.95
10	2.52	2.09	1.92
11	0.76	2.28	4.41
12	2.14	3.67	3.05
13	1.33	1.86	4.16
14	1.18	2.11	4.36
15	1.00	1.49	3.83
16	1.90	2.79	3.88
17	1.94	1.91	5.52
18	2.34	2.62	4.66
19	1.07	1.71	3.44
	CNR¯_NMC±σ_NMC	CNR¯_MC±σ_MC	RMSD¯±σ_RMSD(mm)
	1.57 ± 0.59	2.31 ± 0.67	3.82 ± 1.11

These results are from the contrast‐to‐noise ratio (CNR) analysis of the MC and NMC () dynamic data sets carried out together with the lesion RMS displacement (RMSD) evaluation. The value of () represents the mean CNR (± SD) calculated across all lesions of the MC data set; () refers to the mean CNR (± SD) of the NMC data set.

Dynamic contrast‐enhanced quantitative improvement assessment

Figure 5 shows the anatomical images from the NMC (Figure 5A) and MC (Figure 5B) data sets of a liver affected by metastases, together with the contrast uptake curves measured in an ROI encircling a lesion and a sample of healthy tissue (Figure 5C). The uptake curves regarding the lesion show how the MC and NMC dynamic data sets differ, with a focus on the estimation of the enhancement peak (1.5 minutes) and of the later dynamic phases (3‐5 minutes). The MC contrast uptake curve measured in a lesion shows the characteristic rapid wash‐in followed by rapid washout, which can be clearly distinguished from healthy tissue. In the NMC images, the signal from the lesion is strongly blurred and the measured uptake curve appears to be very similar to healthy tissue.

Figure 5

Transversal slice acquired at the last sampled phase (4 minutes after contrast injection) and extracted from the uncorrected (NMC) (A) and MC (B) dynamic data sets. A region of interest (ROI) is drawn in a hepatic lesion in the MC (blue) and NMC (red) data sets, together with an ROI of healthy tissue (green) in the MC data set. C, Respiratory motion leads to a blurring of signal from the lesion and surrounding healthy tissue, making the measured uptake curve in the lesion difficult to distinguish from healthy tissue. Uptake curves from the MC image show the characteristic uptake behavior of hepatic lesions. Transversal slice at the offspring of the hepatic artery (HA) also show the aorta (AO) at an arterial phase (15 seconds after contrast injection) for the uncorrected (D) and motion‐corrected (E) data sets. Respiratory motion has little effect on the uptake curve from the aorta (F), but can strongly impair the uptake curves in small vessels such as HA (G)

The achieved spatial resolution of 1.5 mm3 is sufficiently high to sample the input functions (Figure 5D,E) necessary for the underlying pharmacokinetic model, in both the aorta (Figure 5F) and hepatic artery (Figure 5G). The input functions sampled in the aorta are very similar due to the position and small size of the sampling ROI, which is hardly affected by respiratory movement. Input functions sampled in the hepatic artery, in contrast, differ strongly between NMC and MC, especially during the enhancement peak, which appears underestimated in the NMC data set.

Figure 6 shows the quantitative ktrans maps calculated on the data sets of patients affected by hepatic lesions injected with gadoxetate disodium. In the first patient (Figure 6A), larger lesions have diminished contrast in the maps derived from NMC dynamic data sets, whereas small lesions disappear completely. Both large and small lesions can be distinguished clearly from healthy liver tissue using the proposed MC scheme. In Figure 6B, the contrast provided in the uncorrected ktrans map is not sufficient to distinguish the highlighted lesion from the background tissue, whereas it is clearly visible using MC. Figure 7 displays the ktrans maps calculated using 2 patients without known or suspected haptic malignant lesions injected with gadobutrol. The first patient (Figure 7A) is affected by a liver hemangioma, which is typically less perfused than healthy tissue. The hemangioma can be better distinguished in the MC ktrans map compared with the uncorrected ktrans map. The second patient (Figure 7B) shows a cyst that can only be properly distinguished from the surrounding healthy tissue using MC. In general, the NMC maps obtained from both gadoxetate disodium and gadobutrol contrast agents show an overall higher range of values for the healthy liver parenchyma. This can reduce the contrast between healthy tissue and the features of interest. Figure 8A,B shows 2 examples of lesions characterized by a central necrotic core (blue arrows), identified on the corresponding ktrans maps as areas of low permeability. In the T2‐weighted images, the necrotic core is confirmed as a hyperintense region.36 In the MC maps, the necrotic cores can be better identified compared with the NMC ktrans maps.

Figure 6

The ktrans maps of 2 patients (A,B) injected with gadoxetate disodium suffering from hepatic malign lesions. Motion correction improves the visibility especially of small lesions in the ktrans maps

Figure 7

The ktrans maps of 2 patients injected with gadobutrol. A, Visibility of a hemangioma is improved with MC. B, A hepatic cyst can be seen clearly in MC, but can hardly be distinguished from the surrounding healthy tissue in NMC

Figure 8

The ktrans maps of a patient injected with gadoxetate disodium and suffering from necrotic hepatic lesions (A,B). Identification and delineation of the necrotic substructure is more accurate in the MC map compared with the NMC map. The necrotic core is confirmed by the hyperintense region in the turbo spin‐echo (TSE) image. It is important to note that the TSE images were acquired during a breath‐hold with a low slice resolution of 6 mm. Therefore, our ktrans maps and the TSE images do not show the exact same anatomy

The quantitative CNR evaluation was carried out on a set of 10 lesions visible in the MC and NMC permeability maps (Table 2). Similar to the CNR evaluation of the late phase, a statistically significant increase ( = 4.12 ± 1.01 against = 2.55 ± 0.80) was obtained, yielding an overall percentage increase averaged over all examined lesions of 62% ( ).

Table 2

Contrast‐to‐noise ratio of parametric ktrans

Lesion No.	CNR_NMC	CNR_MC
1	2.85	3.88
2	3.30	4.91
3	2.23	3.53
4	3.74	4.59
5	3.43	6.32
6	2.16	4.53
7	1.51	3.40
8	2.10	3.93
9	1.42	2.95
10	2.77	3.13
	CNR¯_NMC±σ_NMC	CNR¯_MC±σ_MC
	2.55 ± 0.80	4.12 ± 1.01

These results are from the CNR analysis of the MC ) and NMC () parametric ktrans maps. The value of () represents the mean CNR (± SD) calculated across all lesions of the MC‐derived map; () refers to the mean CNR (± SD) of the NMC map.

DISCUSSION

Abdominal organ motion caused by respiration is a major challenge in the accurate acquisition of the DCE‐MRI time series and subsequent derivation of parametric maps, as it introduces motion artifacts and misalignment between the images at different timepoints. This study demonstrates that based on the proposed MC algorithm, breathing artifacts can be minimized, substantially improving the image quality in DCE‐MRI and increasing the reliability in the derivation of quantitative ktrans maps obtained during free breathing.

In this study, respiratory MFs were estimated directly from the DCE‐MR data. To obtain MFs that accurately describe respiratory motion, it is necessary that image registration of the different motion states is not affected by intensity changes due the contrast agent uptake. That was achieved with the retrospective data‐binning procedure of our approach, which used the same amount of data in each motion state, taken evenly from the entire contrast uptake curve according to the respiratory self‐navigator. Therefore, temporal intensity changes due to contrast agent were averaged out for all motion states, which can be seen in Figure 2.

The measured across all investigated lesions reveals an average absolute displacement of approximately 4 mm, ranging from 2 mm to 6 mm. These measures agree with recent studies that reported an average translational rigid liver motion comprised between 1.2 mm and 19.9 mm, depending on the examined location.37, 38 However, the observed increase in CNR values is not only directly related to the motion amplitude of the lesions before MC is applied. The size of the lesion (with respect to the motion amplitude) and the patient's breathing pattern are also important factors that affect CNR estimation.

An approach to minimize rather than correct for respiratory motion of DCE‐MRI of the liver has been proposed by Feng et al.20 In this XD‐GRASP approach, the DCE‐MRI data acquired with a golden‐angle stack‐of‐stars sequence are retrospectively sorted into 2 additional temporal dimensions (contrast‐enhancement phases and respiratory phases) and reconstructed iteratively, exploiting data similarity between neighboring timepoints along these dimensions. The golden‐angle stack‐of‐stars trajectory obtains a high number of radial frequency–encoded lines, ensuring sufficient data in each contrast and motion‐resolved bin. Golden‐angle radial phase‐encoding, on the other hand, obtains much fewer radial spokes in the ky‐kz plane. Therefore, each bin would not have a sufficient amount of data for image reconstruction. Furthermore, 8 respiratory bins and 48 dynamic time frames would require the simultaneous reconstruction of 384 3D volumes, leading to very high memory demands.

Recently, Johansson et al39 presented a deformable MC approach for DCE‐MRI of the liver. In contrast to our work, motion is not applied during image reconstruction, but a simplified back‐projection deformation approach is used, which carries out MC as a preprocessing step before image reconstruction.40 In addition, their voxel size is approximately a factor of 4 larger than in our study, potentially making the assessment of small features such as necrotic cores challenging.

The parametric ktrans maps show that the proposed MC scheme significantly improves contrast between features of interest (e.g., metastases, cysts, hemangiomas) and healthy tissue and can more accurately resolve tumor substructures that are closely related to perfusion, such as necrotic cores. In addition, ktrans values of healthy liver tissue are more reproducible among patients. This is due to 2 main factors that have a direct effect on the fitting of the pharmacokinetic model. First, minimizing motion artifacts and image misalignment in the NMC data sets leads to better‐defined features and ensures that data fitting is carried out over well‐aligned voxels. Second, motion artifacts can also impair the AIF signal, especially in the case of the relatively small hepatic artery. These artifacts primarily lead to an underestimation of contrast uptake during the enhancement peak (Figure 5C) and therefore to an overestimation of ktrans. Following MC, the detection and delineation of hepatic lesions is substantially improved in both the dynamic images and parametric ktrans maps. This assumption is confirmed by the CNR evaluations conducted on the latest dynamic phase recorded and the ktrans maps, achieving an average increase of 47% and 62%, respectively. These results suggest that successful MC can improve the visualization of various features within each patient, but at the same it could also provide more accurate comparisons concerning group studies.

Further improvements in ktrans estimation could be achieved by using native T1 maps32 for signal conversion and automatic contrast injection. Both NMC and MC ktrans maps are affected by these parameters; therefore, the results of the comparison between NMC and MC carried out in this study would not change.

Here, we obtained the AIF from the hepatic artery to demonstrate the effect that respiratory motion can have on the calculation of the AIF, but we did not use an uptake model that is specifically designed for this AIF. Nevertheless, the improvements of ktrans achieved with our approach in this study are also expected to apply to advanced pharmacokinetic models, which specifically require information from the hepatic artery such as the dual‐input single‐compartment and 2‐compartment models.41, 42

CONCLUSIONS

In this study, we introduced a novel MC method allowing for free breathing 3D abdominal DCE‐MRI with isotropic resolution of 1.5 mm covering the entire contrast uptake curve with a temporal resolution of 6 seconds. Nonrigid respiratory MC minimized motion‐induced blurring and artifacts, leading to an improved visualization and delineation of lesions. In addition, it ensured accurately aligned dynamic images, improving the calculation of quantitative 3D ktrans permeability maps. This approach does not require any breath‐holds, and the large FOV with isotropic resolution minimizes scan planning, making this technique straight forward to apply in a clinical setting.

CONFLICT OF INTEREST

Dr. Mathias Lukas is an employee of Siemens Healthcare.

FIGURE S1 Contrast dynamics of 8 temporal timepoints out of the 36 DCE volumes for a patient injected with gadobutrol. Note: For visibility, only a part of the FOV is shown. Respiratory MC reduces motion‐induced blurring (red arrows) of the anatomy

Click here for additional data file.