Effect of intravenously injected gadolinium‐based contrast agents on functional lung parameters derived by PREFUL MRI

To evaluate the influence of intravenously administered gadolinium‐based contrast agents on functional ventilation and perfusion parameters derived by phase‐resolved functional lung (PREFUL) MRI.


| INTRODUCTION
The current clinical standard for pulmonary disease monitoring and diagnosis is spirometry. Because of the lack of regional information, detection of early disease remains challenging. Therefore, robust and sensitive markers that reliably quantify pulmonary function on a regional level are of tremendous interest. As an alternative to hyperpolarized multi-nuclear MRI, several proton MRI techniques have been introduced to assess pulmonary ventilation and perfusion on a regional scale without the usage of a contrast agent during free breathing. 1,2 Fourier decomposition (FD)-based methods use cardiopulmonary frequency-dependent signal changes to obtain either cardiacor respiration-derived signal alteration and thereby enable the calculation of functional parameters. [3][4][5] Recently, techniques including phase-resolved functional lung MRI (PREFUL) 6 with increased temporal resolution were introduced to extract full cardiac and respiratory cycles assessing the dynamic component. 7,8 Functional lung parameters derived by PREFUL MRI and other FD-based methods have been shown to be biomarkers for detection and monitoring of many pulmonary diseases [9][10][11] and have served as endpoints in clinical trials. 12 The benefit of radiation-free acquisition in combination with patient-friendly examination during free breathing using the widely available and stable spoiled gradient echo sequence makes PREFUL a potential candidate for functional lung assessment in daily patient care. However, further validation and development of new biomarkers and pipelines that enable fast processing will be necessary to establish this method in clinical practice. 13 Despite its known risks of allergic reactions, renal failure and nephrogenic systemic fibrosis (NSF), 14 intravenous gadolinium-based contrast agents (GBCAs) are still used for various clinical indications. [15][16][17][18] Moreover, dynamic contrastenhanced (DCE) MR-perfusion is considered the current clinical standard for perfusion quantification and therefore serves as the primary endpoint in many studies. Because the total scan time is limited and optional sequences are mostly included at the end of the protocol, the influence of i.v.administered GBCAs on PREFUL and similar techniques is highly relevant. Additionally, an influence of injected GBCAs on PREFUL parameters would provide further insight regarding the level of hematocrit or blood oxygenation, because they also influence T 1 relaxation times. [19][20][21] Therefore, the purpose of this study was to test the effect of i.v.-administered GBCAs on perfusion and ventilation parameters derived by PREFUL MRI.

| Participants
Data sets of 17 participants consisting of 9 men and 8 women were processed in this study. This data was acquired specifically for the presented work either from clinically indicated lung MRIs (CTEPH population) and 1 research study (COPD population). The MR imaging protocols were modified with additional FD series after pulmonary MR angiographies to generate the data sets. Patients with known chronic obstructive pulmonary disease (COPD) or chronic thromboembolic pulmonary hypertension (CTEPH) and >18 y of age were included. Subjects with a general inability to undergo MRI (e.g. due to claustrophobia, paramagnetic implants, hypersensitivity to i.v. GBCA's, pregnancy or breast feeding) and subjects with a previous lung surgery, pneumonia or exacerbation within the last 4 weeks were excluded. The study was approved by the local ethics committee. Furthermore, written informed consent was obtained from all participants before the MR examination.

| Imaging technique
All MR scans were performed at 1.5T (Magnetom Aera or Avanto, Siemens Healthcare, Erlangen, Germany). For PREFUL MRI, a spoiled gradient echo sequence with following parameters was used: field of view = 50 × 50 cm 2 , matrix size = 128 × 128, slice thickness = 15 mm, TE = 0.82 ms, TR = 3 ms, flip angle = 5°, pixel bandwidth = 1500 Hz/pixel and parallel imaging with GRAPPA with acceleration factor 2. The participants were scanned in head first supine position. In ~60 s, 250 images of a single coronal slice at the level of the carina were acquired during free breathing. Before reconstruction, all images were interpolated to the final in-plane resolution of 256 × 256 pixels. Furthermore, the non-uniform intensity profiles of the applied receiver coils were corrected using the sensitivity profiles of the surface and body coil before PREFUL analysis.
The 1st PREFUL measurement was obtained 33:46 min (SD = 6:20 min) before an MR-angiography with i.v. application of (0.1 mL/kg*body mass) gadobutrol. 22 Subsequently, 2 further PREFUL measurements were acquired to assess the effect of the applied contrast agent on the PREFUL technique. The starting points of the 2nd and the 3rd PREFUL measurements were 43 s (SD = 1.9 s) and 91 s (SD = 1.9 s) after the beginning of the MR-angiography.

| Postprocessing
The first 20 images of every PREFUL data set were discarded to ensure a steady state of the stationary tissue. Subsequently, all 3 data sets were combined into 1 and linear and diffeomorphic image registration using advanced normalization tools (ANTs) with a previously described group-oriented registration scheme (GOREG) was performed to achieve a uniform respiratory position for each image and each PREFUL acquisition. 23,24 Further postprocessing steps were performed separately on each data set. Image guided filtering was applied to all images of every data set using an edge-preserving filter. 25

| Full respiratory cycle
For the calculation of ventilation-weighted images, a lowpass filter at 0.8 Hz was used to remove perfusion effects from the registered images. Based on the averaged signal amplitude of a region of interest (ROI) positioned at the lungdiaphragm, images were sorted according to their respiratory phase. Subsequently, the sorted images were interpolated to a uniform time grid. Finally, ventilation including phase information was calculated. 6 Ventilation was quantified using regional ventilation (RV) according to Klimeš et al 26 S Exp and S Insp represent the averaged signals of the endexpiratory and end-inspiratory images. S Mid represents the mean signal of the images in intermediate lung inflation level, which is also used as the reference volume during registration.

| Full cardiac cycle
A high pass filter at 0.75 Hz was applied to the registered images to exclude signal changes because of respiration. To estimate the cardiac phase of each time point, similar to ventilation analysis, the averaged signal of a manually segmented ROI of the aorta was fitted to a piecewise cosine function. Subsequently, images were sorted according to their phase and interpolated to an equidistant time grid. Perfusion(Q) was calculated by taking the difference between the signal at t Max and t Min representing the time with the highest and the lowest signal of the perfusion weighted time series of the parenchyma ROI. Furthermore, Q was quantified (Q Quant ) in mL/min/100 mL by normalizing the signal of every voxel to the signal of a completely blood-filled voxel (S Blood ) according to Kjørstad et al 27 Q represents the cardiac frequency dependent signal difference ∆ of the parenchymal signal, S Blood represents the corresponding ∆ of the aorta signal between t Max and t Min in the aorta ROI and t exp represents the time between 2 heartbeats.

| Assessment of ventilation and perfusion maps
The lung parenchyma was segmented manually by individual thresholding. Visual criteria were the inclusion of the vast majority of the lung parenchyma and the exclusion of large central vessels up to the segmental level. Insufficiently segmented areas were corrected manually to receive optimal segmentations. As all 3 series of each patient were coregistered, only 1 segmentation for all series was obtained on native images. Mean of all maximal and minimal values in the time series of the segmented lung parenchyma and aorta ROI were calculated. To calculate ventilation defect percentage (VDP), voxels were rated as ventilation defect, if RV was below the 75th percentile of all RV values multiplied by a factor of 0.7. 28 Similarly, perfusion defect percentage (QDP) maps were calculated with a threshold defined as the 75th percentile of all Q values multiplied by a factor of 0.6. 29 These factors were determined empirically to achieve maximal accordance between FD and DCE. Furthermore, V/Q match, which equals V/Q Quant match, was calculated.

| Statistics
Statistical analysis was performed using JMP Pro 13 (SAS Institute, Cary, NC) and MATLAB R2018b (The MathWorks, Natick, MA). Shapiro-Wilk test revealed no normal distribution of all functional parameters. Hence, non-parametric Wilcoxon signed-rank test was used to assess significant differences of functional parameters (RV, Q, Q Quant , QDP, VDP, VQ-match percentages, and mean maximal and minimal values of the lung parenchyma ROI and the aorta ROI) between the 3 PREFUL acquisitions. Because of multiple testing, the Bonferroni-Holm correction of significant P-values was included. Furthermore, Sørensen-Dice coefficients were calculated to compare the spatial overlap of QDP und VDP between each series. As QDP and VDP maps (Supporting Information Figures S1 and S2) contain only binary information, the resulting DICE maps are also binary and present ether matching or mismatching areas of 2 series. Results with P < 0.05 were considered as significant. A retrospective power calculation for same day repeatability without leaving the MR scanner in-between scans was performed based on the calculated SD of QDP and VDP in the Bland-Altman plots.

| Patient characteristics
Three out of 17 processed data sets had to be excluded. Exclusion criteria was incomparability with the other data sets and scanning procedures. All 3 exclusions showed individual but no systematic errors. The first patient was excluded as the 3 data sets could not be registered toward 1 image. Analysis of the 3 data sets revealed a significant movement between the series that led to repetitive abortion of the registration procedure. As we desired a voxel-by-voxel comparability between the series, the registration to 1 image was obligatory.
The second patient was excluded because a wrong scaling factor in the reconstruction procedure of the scanner prevented comparability of the signal intensities of the individual data sets. This problem was detected after the raw data was deleted from the scanner and the reconstruction could not be repeated once again. The third patient was excluded because an interruption of the examination led to an exceptional pause between the angiogram and the second series. Whereas the second series of all other patients followed precisely 43 s after the angiogram, the interruption in this examination caused a delay of 146 s. Hence, the gadolinium concentration was significantly lower and comparability to the other data sets could not be achieved. Table S1 presents demographic, physiological, and technical data. VQ maps revealed no significant differences (see Table 1). Representative binary VDP, QDP, and VQ maps of an 81-year-old male patient with COPD are displayed in Figure 2. Sørensen-Dice coefficients, depicting the spatial overlap of QDP and VDP maps between the different series, varied up to ±9% (see Table 2 and Supporting Information Figures S1   and S2). The Bland-Altman analysis revealed a SD for VDP of ≤4.03% and for QDP of ≤10.1%. Using the SD of these parameters for same day repeatability without leaving the MR scanner in-between scans, a retrospective power analysis showed that a sample size of 14 patients can detect a VDP difference of 3.3% and a QDP difference of 8.2%.

Functional parameters
Wilcoxon signed-rank test

| Quantified perfusion (Q Quant )
As indicated in Figure 3 and Table 1, quantified perfusion (Q Quant ) increased significantly after the application of gadobutrol between the 1st and the 2nd series and between the 1st and the 3rd series. Analysis of signal intensities of systole and diastole revealed a pronounced decrease of signal range in the aorta (see Table 3). More specifically, all mean values increased after the application of gadobutrol, but the systolic and diastolic signal range was affected differently for the lung parenchyma and the aorta. Heart frequency, which has a linear influence on quantified perfusion, indicated no significant difference between each series (see Table 1).

| DISCUSSION
This study is the first to assess the influence of i.v. applied GBCA on ventilation and perfusion parameters derived by PREFUL MRI and related techniques, including Fourier decomposition. The major 2 findings of this work are as follows. First, PREFUL MRI-derived regional ventilation (RV) and perfusion (Q) maps, as well as secondary parameters including ventilation defect percentage (VDP), perfusion defect percentage (QDP), and ventilation/perfusion match (VQM) were not influenced significantly by the administration of gadobutrol. RV was not significantly changed by gadobutrol, because in a linear correlation between lung volume and ventilation values-as it was shown by Zapke et al 1 -the increase of signal is proportional between expiration and inspiration and no alteration is observed. Similarly to RV, Q values were not perturbed significantly by gadobutrol. Assuming that there is only a small difference in flow velocity between systole and diastole in the lung parenchyma (∆ parenchyma signal), the influence of shortening of T 1 relaxation time after the application of gadobutrol-as it is dependent on velocity-would be similar and no significant change of Q values would be expected.
Second, quantified perfusion maps (Q Quant ) derived by the method as described by Kjørstad et al 27 delivered significantly elevated values after the application of gadobutrol. The significant difference of Q Quant can be explained by the flowrelated enhancement and its dependency on T 1 relaxation time. During the systolic phase, spins in the aorta provide higher signal, because these fast-flowing spins experienced less saturation. In this phase, the signal is almost independent of the T 1 relaxation time (see Supporting Information Figure S3, t < 500 ms). Whereas during diastole, the process of saturation is distinctly advanced through very slow flow and the influence of different T 1 relaxation time becomes noticeable (see Supporting Information Figure S3, t > 500 ms).

T A B L E 2 Mean Sørensen-Dice coefficients
Although the decrease of the T 1 relaxation time by gadobutrol leads to a signal increase in steady-state regime (described by Equation 3), the effective signal enhancement is a function of flow as demonstrated in Supporting Information Figure S3 in a Bloch-Simulation. 30 K is a constant and [H] refers to spin density. This effect is pronounced in arterial vessels close to the heart where high differences of velocities occur. As the systolic flow velocity is very high in the aorta, only diastolic signal increases noticeably, and the ∆ of the aorta signal shrinks. The signal difference of the aorta (∆ aorta signal) decreased because of a higher gain of signal during diastole compared to systole, whereas in the lung parenchyma, the systolic, and diastolic signal increased equally and ∆ parenchyma signal stayed constant (see Figure 4). Regarding the Kjørstad method for quantified perfusion (Equation 2), results increase as ∆ aorta signal (S Blood ) decreases.
In this context, the formula used for Q Quant (Equation 2) should be discussed critically. Besides the mentioned influence of variable T 1 relaxation time on perfusion (1 − cos ) e −TR/T1 × e −TR/T2 * .

Series 1 Series 2 Series 3 P-value 1 vs. 2 P-value 2 vs. 3 P-value 1 vs. 3
Systole aorta (a.u.) 1724. 11  Mean values of all subjects ± SD. All values increase after application of gadobutrol. Highest gain is observed in the diastolic signal of the aorta, which results in a significantly decreased aortic signal range ∆. As expiratory and inspiratory values increase almost proportionally, the percentage of parenchymal signal change does not alter significantly.

F I G U R E 4 Noticeable decrease of
aortic signal range because of a stronger gain of diastolic signal (pink line) compared to systolic signal (red line) after contrast enhancement. The parenchymal signal range remains almost unchanged parameters, a major disadvantage of this method is the mandatory requirement of a central lung vessel as a ROI for the completely blood-filled voxel. This required central lung vessel isn't depicted orthogonally on all slice orientations. Specifically, flow-related enhancement depends on slice thickness, slice orientation, and the direction of vessels. Furthermore, different flow velocities in these ROIs lead to distinct values of Q Quant . In our case-as we chose the aorta as the reference ROI for the completely blood-filled voxelpatients with aortic valve stenosis would present lower values of Q Quant as expected because of an altered and more turbulent flow velocity during systole. Whether this second observation is present under all clinical circumstances (e.g., decreased perfusion and different pathophysiologic constitutions of vessel structure in COPD) 31 cannot be answered by this study. Further limitations of our study are the limited data points after the application of gadobutrol. The temporal resolution of measured data is not sufficient to allow for the detection of complex changes over a time course, which may depend on cardiac output, diffusion parameters, and renal function. The limited and asymmetrically distributed time points also mask the critical role of time delay in regard to reproducibility of VDP and QDP maps. Although the technical reproducibility of ventilation-and perfusionweighted maps was already proven, 4 Voskrebenzev et al 32 showed that tidal volume variability is the main limiting factor to FV reproducibility. As perfusion is strongly connected to ventilation via hypoxic pulmonary vasoconstriction (Euler-Liljestrand mechanism), it is not surprising to see both ventilation and perfusion mismatch between 2 scans. The highest Dice-coefficients can be found between the 2nd and the 3rd series, which were scanned directly after another (see Supporting Information Table S1). Accordingly, the Bland-Altman plots of QDP show a SD of only 3.62% between these connected series, whereas the SD between the 1st and the 2nd series is 10.10% and 7.97% between the 1st and the 3rd series (see Supporting Information Figure S4). Although Figure 2 presents a good agreement of defect ventilation and perfusion areas, not all defects of the 1st series are present in the 2nd, and new defects occurred. Supporting Information Figure S5 reveals the different patient groups of this study. The crowded data points with minimal mean delta belong to the patient group with pulmonary hypertension where many patients had a VDP of 0%. Patients with COPD show more spread of the VDP values. Nevertheless, comparing the SD of ≤4.03% for VDP same-day reproducibility without leaving the scanner in this study to a SD of 10.14% for VDP of COPD patients in a previous study with a time delay of 2 weeks between the scans, 12 one can estimate the major influences of time delay and type of disease on reproducibility by local changes of ventilation.
Further minor causes for mismatched areas might be technical inconsistencies like misregistrations or local field inhomogeneities. Even with inhomogeneous lung parenchymal enhancement, perfusion and ventilation parameters remain unchanged by GBCA, because both diastolic and systolic values are enhanced proportionally in each voxel.
Despite the mentioned limitation regarding quantified perfusion, perfusion-weighted maps, regional ventilation metrics, and defect percentage maps derived by PREFUL and related techniques have been shown to be promising parameters capable of monitoring different lung diseases 9,29,30 with good agreement to DCE-MRI perfusion parameters 33 and spirometry. 29

| CONCLUSION
In conclusion, all ventilation and perfusion parameters derived by PREFUL MRI, with exception of quantified perfusion using the Kjørstad method, 27 are not influenced significantly by the administration of gadobutrol, which is important for the design of MRI protocols in future studies.

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of the article.

FIGURE S1
Example of Dice maps showing the spatial overlap for all voxels between the QDP maps of the 1st, the 2nd, and the 3rd series of an 81-year-old male patient with COPD FIGURE S2 Example of Dice maps showing the spatial overlap for all voxels between the VDP maps of the 1st, the 2nd, and the 3rd series of a 68-year-old female patient with pulmonary hypertension FIGURE S3 Bloch-simulation of transverse magnetization toward steady state. Whereas at t (0) the maximal value is equal, steady-state signal indicates a clear difference between variable T 1 relaxation times (red line T 1 = 1400 ms; blue line T 1 = 1100 ms). Accordingly, fast flowing spins with short length of stay in slice experience less saturation and their signal is almost independent of T 1 time (t < 500 ms), whereas slow flowing spins shift toward steady state and T 1 effect becomes noticeable (t > 500 ms) FIGURE S4 Bland-Altman plots comparing QDP of the 3 series. The red line shows mean delta between the 2 measurements. The dashed lines indicate the upper and lower limits of the 95% confidence interval. Reduced reproducibility in the plots comparing series with longer time delay (1 vs. 2 and 1 vs. 3). Mean SD of all plots is 7.23% FIGURE S5 Bland-Altman plots comparing VDP of the 3 series. The red line shows the mean delta between the 2 measurements. The dashed lines indicate the upper and lower limits of the 95% confidence interval. Two different groups can be distinguished. One crowded point cloud close to 0 average and 0 delta and 1 dispersed point cloud with increased variability. Mean SD of all plots is 3.68% TABLE S1 Demographic and technical data. Mean values with standard deviation in parenthesis. Absolute numbers of the underlying diseases chronic obstructive pulmonary disease (COPD) and pulmonary hypertension. An additional reduced amount of gadobutrol was applied prior in a TWIST sequence to determine the optimal trigger time for the MR angiography