Estimating ventilation correlation coefficients in the lungs using PREFUL‐MRI in chronic obstructive pulmonary disease patients and healthy adults

Various parameters of regional lung ventilation can be estimated using phase‐resolved functional lung (PREFUL)‐MRI. The parameter “ventilation correlation coefficient (Vent‐CC)” was shown advantageous because it assesses the dynamics of regional air flow. Calculating Vent‐CC depends on a voxel‐wise comparison to a healthy reference flow curve. This work examines the effect of placing a reference region of interest (ROI) in various lung quadrants or in different coronal slices. Furthermore, algorithms for automated ROI selection are presented and compared in terms of test–retest repeatability.


INTRODUCTION
2][3][4] PREFUL-MRI datasets are usually acquired in free tidal breathing using a 2D spoiled gradient echo sequence with high temporal resolution of three to five images per second.A post-processing algorithm detects and separates changes in the signal into ventilation and perfusion components.These are used to calculate parameters of regional lung function which in turn could be used to estimate matched/non-matched ventilation/perfusion defects and calculate pulmonary ventilation and perfusion defect percentages.6][7][8] These methods implement multi-step post-processing with multiple frequency filtering operations on the signal of a dynamic MRI dataset to estimate ventilation and perfusion components at each lung voxel.The post-processing in PREFUL relies more on the calculation of amplitudes in the time domain without a strict low-pass/high-pass filtering; thus, there is no need to define specific frequency thresholds corresponding to ventilation and perfusion components making PREFUL less prone to measurement artifacts in case of respiration variability.PREFUL is also advantageous in the technical feasibility to be performed with conventional spoiled gradient echo sequences allowing for more practicality in the translation to a clinical setting.On the contrary, the balanced SSFP sequences (bSSFP) used in FD-based approaches have limited application at 3T scanners and at 1.5T these might require optimization for lung acquisition and for reduction of motion artifacts, for example, by applying parallel imaging techniques and asymmetric echo sampling. 5Another example is the self-gated non-contrast-enhanced functional lung imaging technique, which requires sequences with self-gating capabilities. 9or detecting early lung disease, previous work showed the importance of assessing the whole dynamics of the lung ventilation cycle instead of looking at fixed time points at end-inspiration and end-expiration. 10Using PREFUL-MRI, flow-volume loops of the ventilation cycle are generated for each voxel in the lung parenchyma.To allow for more practicality in a clinical setting, a quantitative parameter of the whole ventilation cycle based on regional air flow dynamics has been previously presented, that is, the ventilation correlation coefficient (Vent-CC), which was referred to in some previous studies as the cross-correlation metric or as the flow-volume loops correlation coefficient. 11This parameter relies on comparing the shape of the air flow time course of each lung voxel to the mean flow curve of a healthy reference region of interest (ROI).Thus, the use of a suboptimal healthy reference may cause an inaccurate estimation of Vent-CC.Current PREFUL pipeline integrates a fully automated algorithm for the segmentation of the reference ROI.This algorithm selects healthy reference voxels independently in each slice based on regional ventilation (RVent) values.The biggest cluster of voxels having RVent value in the 80th to 90th percentile range is selected as a reference ROI.RVent is calculated based only on specific time points of the ventilation cycle (end expiration, mid-inspiration and end-inspiration) 4 ; thus, some inhomogeneities in ventilation dynamics might still be present in voxels with similar RVent values.Furthermore, in slices with extensive disease and only small areas of healthy ventilation, voxels in the relatively high RVent percentile range might nonetheless be expected to have impaired ventilation, thus, it could be assumed that a significant underestimation of present ventilation disturbance might occur when such non-healthy voxels are used as a reference for the calculation of Vent-CC.A solution could be the use of one ROI for all slices; however, it is expected that various coronal slices could have differences in their ventilation dynamics, for example, due to gravitational effects or due to changes in tidal volume that might occur during the acquisition. 12,13his effect must be also examined to define criteria for selecting the optimal reference ROI.
This work studied whether segmenting the healthy reference ROI at various locations in the lung could affect the estimation of Vent-CC.In addition, we examined whether an identical healthy reference curve could be used for all slices acquired at the same scan session or whether unique ROIs must be segmented separately for each 2D slice.Considering the acquired knowledge obtained from the first two study questions mentioned above, a new algorithm with the goal of a more homogenous ROI was developed.The repeatability of Vent-CC and ventilation defect percentage (VDP) when calculated using the presented algorithms was assessed in a test-retest setting in both healthy participants and in patients with chronic obstructive pulmonary disease (COPD).Inclusion criteria in the CLAIM study were patients aged 40 y or older with a clinical COPD diagnosis and a baseline hyperinflation (residual volume > 135% predicted).Endpoints of the CLAIM study were previously published. 14,15he subgroup of CLAIM participants who were included in this work were patients with stable COPD who received placebo inhalation in the first 2 wk of their participation in the CLAIM study.Few patients had to be excluded from the sub-cohort due to discontinued participation before completing the second baseline MRI scan.Included healthy subjects were the first 32 participants who got prospectively recruited in a previous study. 3Demographic and clinical characteristics of all study participants are summarized in Table 1.

Study cohort and MRI acquisition
Participants in both groups have completed two MRI sessions (including PREFUL in free-breathing) on the same 1.5T scanner (Magnetom Avanto, Siemens Healthcare, Erlangen, Germany) with an eight-channel torso phased array coil.Healthy participants underwent the second scan session after a short break outside the scanner while the COPD group had a 4 wk gap between their two baseline scan sessions.PREFUL datasets were acquired using a 2D spoiled gradient echo sequence with the following scan parameters: TE 0.82 ms, TR 3 ms, flip angle 5 • , matrix size 128 × 96 interpolated to 256 × 256, FOV 50 × 50 cm 2 , slice thickness 15 mm.In each session, a total of three coronal slices of the lung were scanned (the middle slice was always set at the coronal plane of the tracheal bifurcation with the other two slices planned more anteriorly and posteriorly from the middle plane).Two hundred images per slice were acquired with a temporal resolution of 288 ms (scan time was approximately 1 min per slice).

Post-processing
All images were processed with a unified version of an inhouse-developed automated PREFUL pipeline including slice-wise group-oriented registration (GOREG) performed using the advanced normalization toolkit (ANTS). 16,17After sorting the images in a virtual ventilation cycle based on the displacement of the diaphragm, the time series of regional ventilation (RVent (t) ) was calculated for each voxel.Using that, regional flow-volume loops (Figure 1) were calculated as described previously. 10,11In brief, we used the first derivative of RVent (t) as a proxy for flow measurement, while RVent (t) itself provided the volume information.Ventilation correlation maps were computed by comparing the air flow curve of each voxel in lung parenchyma to the mean flow curve of a healthy reference ROI using cross-correlation with a fixed time lag of 0: x n y n with x n ∈ R representing the reference flow component and y n ∈ R the respective regional flow component.The coefficient is normalized according to: .

F I G U R E 1
Mean flow-volume loop from all lung voxels in the chronic obstructive pulmonary disease (COPD) cohort (red line) and in healthy participants (green line).Respective color shaded areas represent the SD.
The Vent-CC parameter was calculated as the mean correlation coefficient from all lung voxels.
Ventilation defect maps were generated by applying a fixed threshold (0.9) to the correlation images. 10VDP was calculated slice-wise as the percentage of defect voxels in the binary ventilation defect map in relation to the voxels count in the corresponding lung segmentation mask.

Data analysis
The following image and statistical analyses were conducted using self-developed scripts written in MATLAB (R2019b, MathWorks, MA, USA).For all analyses in this study, a confidence level of 95% was used.
A) To study the effect of having the reference ROI at various locations in the lung parenchyma or at a different coronal slice, the following analyses were performed on data from the first scan session of the healthy participants.In this cohort, no lung areas with significant function impairment are expected to be present, hence, most voxels in the lung parenchyma could be considered as a healthy reference for the calculation of the Vent-CC.For each coronal slice in that dataset, a trained medical student (J.R.) segmented four free-shape ROIs at different locations in the lung parenchyma (left upper, left lower, right upper, and right lower quadrants), these ROIs will be referred to in the following text with ROI quadrant (Figure 2).During the segmentation, any lung areas with markedly low RVent values (e.g., motion artifacts close to the diaphragm) were avoided.The mean ROI size was 131 voxels with a SD of 48 voxels.For following inter-slice comparisons, slice-wise ROIs (ROI slice ) were generated for each of the three coronal slices by combining the four ROI quadrant from that same slice into one ROI slice .
First, the homogeneity of regional ventilation loops at different ROI locations in the lung parenchyma was evaluated.In each slice, the mean flow curves of the four ROI quadrant were compared to each other using cross-correlation resulting in a 4 × 4 correlation matrix (intra-slice comparison of ROI location).In addition, an inter-slice comparison between the flow curves of each slice-specific ROI slice was also performed using cross-correlation resulting in a 3 × 3 correlation matrix.Correlation values below the threshold of 0.97 could imply regional inhomogeneities in the air flow.This very conservative threshold was set based on knowledge from previous studies and practical experience with the Vent-CC metric. 10The higher the correlation between various ROI quadrant , the less the expected effect of the location of ROI in calculating Vent-CC.The following analysis assessed the variability in regional ventilation between different coronal slices.Mean RVent in lung parenchyma was compared between frontal, middle and dorsal slices using a one-way analysis of variance (ANOVA) test.Although RVent does not assess the dynamics of the whole ventilation cycle, it is still a good basic measure of regional ventilation without the need for a healthy reference.Significant differences in regional ventilation between different coronal slices would suggest the importance of selecting a reference ROI from the same slice.Additional analyses were conducted to look for significant differences in Vent-CC and VDP when calculated using various ROI locations in the same slice or from different coronal slices.For each slice, the Vent-CC image (and hence the VDP) was calculated four times using a different ROI quadrant and three other times using a ROI from different coronal slices (ROI slice ).Mean Vent-CC of entire lung parenchyma and the VDP was compared between the analysis variants using a one-way ANOVA test.

B)
Because the original PREFUL algorithm for selecting a reference ROI depends solely on RVent values, the following analysis was performed to study Representative coronal slices from one participant showing the manually segmented four ROI quadrant .
the homogeneity of ventilation dynamics in relation to RVent.In healthy participants and in the COPD cohort, lung parenchyma voxels from each slice were sorted in 10 groups according to their percentile RVent values (0th-10th percentile, 10th-20th percentile, etc.).A correlation matrix of the group's flow curves was calculated (Pearson's correlation).Group-wise internal mean correlation and SD was examined for all RVent percentile groups.A relatively low mean correlation and high SD would indicate inhomogeneous air flow inside the group's voxels.C) The test-retest repeatability of Vent-CC and VDP was compared between analyses using two different algorithms for automatic selection of the healthy reference ROI.The coefficient of variation was used as a metric for repeatability and was computed using the formula: SD × 100/mean.Both tested algorithms have inherent assumptions that agree with the results of the previously mentioned analyses in section A.
Like with the original PREFUL algorithm, the below presented algorithm assumes that healthy ventilated voxels should have RVent values within normal limits.However, the presented algorithm adds an extra step to ensure that selected reference voxels have homogenous ventilation dynamics.
The algorithms work as follows: 1. Within the framework of PREFUL processing, RVent images are generated, and lung parenchyma is automatically segmented using a previously trained U-NET convolutional neural network (Figure 3A). 182. Lung voxels with a relatively high (normal) RVent (X 1 -X 2 percentile) are selected from each coronal slice (Figure 3B).The parameters X 1 and X 2 were set to 0.8 and 0.9, respectively (in order to exclude outliers with very high RVent values, for example, voxels with motion artifacts or compensating hyperventilated areas).3.Only in the new algorithm (Figure 3C): The air flow time course of each selected voxel is correlated against the flow curves of all other selected voxels from the same slice.The median correlation for each lung voxel is then calculated.Selected voxels are sorted by their median correlation.Voxels with the lowest 50% correlation coefficients are removed from the selection in order to ensure homogenous ventilation dynamics in the healthy ROI. 4. The largest connected voxel cluster in the selection is detected using MATLAB's Image Processing Toolbox (the function "bwareafilt" extracts objects from binary images).Voxels outside that cluster are removed from the selection.The spatial coherence of the selected voxels is assumed to result in a more homogenous healthy reference.This step could be omitted in the new algorithm as the latter inherently excludes voxels with inhomogeneous ventilation dynamics (step 3).
Flowchart of the original and proposed algorithm for the selection of a reference region of interest (ROI).

T A B L E 2
The correlation matrix between mean flow curves of various ROI quadrant which were segmented in the same coronal slice, and the correlation matrix between various ROI slice segmented in different coronal slices.

RESULTS
A) A high correlation between the air flow curves of the four ROI quadrant was present in all examined slices (mean correlation coefficient: 0.994 ± 0.003) with all correlation coefficients higher than the threshold of 0.97 (Table 2).An inter-slice comparison between ROI slice showed lower correlation coefficients of their flow dynamics (mean correlation coefficient: 0.965 ± 0.007) with values below the defined threshold in two of the three test cases (Table 2).
Statistically significant differences in RVent were found in healthy participants between the three coronal slices (ventral, central and dorsal; ANOVA p < 0.001).No significant differences in whole lung Vent-CC or VDP were found between the four analysis variants using four different ROI quadrant from the same slice (ANOVA: p > 0.5).Statistically significant differences in whole lung Vent-CC and VDP were found between the three analysis variants using ROIs from different slices (ANOVA p < 0.001).In all three coronal slices, the use of a reference ROI from the same slice resulted in the lowest VDP in the healthy participant cohort (Table 3).

B)
In healthy participants, all voxel groups in the 40th-100th percentile range of RVent had a comparable high internal correlation with a mean coefficient above 0.97 and a maximum SD of 0.031.Only in the 0-10th percentile range group, an obvious decrease in the correlation coefficient was seen (Figure 4A,C).
On the other hand, COPD patients showed generally lower internal correlation values (best internal correlation was for the 90th-100th percentile group with a mean correlation coefficient of 0.92).The SD was also higher in all voxel groups in this patient cohort compared to healthy participants (Figure 4B,D).Hence, an inhomogeneity of the air flow dynamics was present in COPD patients even in voxels with a relatively high RVent.C) The test-retest coefficient of variation of Vent-CC and VDP was lower for both healthy and COPD cohorts when using the new algorithm as with the original algorithm (Table 4; Figures 5 and 6).Cohen's d effect size was 0.727 for Vent-CC and 0.137 for VDP.Using both algorithms, mean values of Vent-CC and VDP showed a significant difference between COPD patients and healthy participants in the test and re-test datasets (p-values of two-sample t-tests <0.0001 after adjusting for multiple comparisons using Bonferroni correction with n = 8).

DISCUSSION
A non-invasive assessment of regional ventilation dynamics in the lung could be quantified using PREFUL-MRI with the parameter Vent-CC.In a slice-wise calculation of Vent-CC, the anatomical location of the healthy reference ROI (upper or lower, right or left lung zones in the same slice) was shown to be insignificant.Different coronal slices from the same healthy volunteer and the same scan session showed heterogeneities in their air flow dynamics.This might occur due to the presence of physiological anterior-posterior gradient in lung ventilation in supine position which agrees with findings from the literature. 12It is interesting for future work to test whether image acquisition in a prone position would result in a similar heterogeneity gradient of the ventilation dynamics.In addition, other artificial sources for the heterogeneity of the flow curves should be considered.[21][22] Inhomogeneities in ventilation dynamics were present in the same coronal plane especially in COPD patients and even between voxels with high RVent values, implying the necessity to improve the original algorithm for selecting the reference ROI.By excluding those voxels whose air flow dynamics were inhomogeneous to other voxels within the same high RVent range, a more homogeneous healthy reference with exclusion of "non-synchronous" voxels or outliers could be achieved.Using the presented method, Vent-CC and VDP could be calculated with a better test-retest repeatability in COPD and healthy volunteer cohorts.
Our findings might be translated to the 3D PRE-FUL pipeline.While 3D measurements can reduce common 2D artifacts caused by partial-volume effects and through-plane motion, it has been shown that the ventilation dynamic parameters (Vent-CC and VDP) derived by 3D PREFUL-MRI suffer from decreased repeatability. 23,24hus, the gained knowledge from this work may be used to improve the stability of 3D PREFUL-MRI parameters of ventilation dynamics.Further, the improved repeatability of the presented algorithm could be tested on datasets acquired on a 3T scanner as these might suffer more from field inhomogeneities.A previous study on PREFUL found no significant differences in RVent measurements between 1.5T and 3T. 25 However, the Vent-CC parameter was not compared in that study.
Although no ground truth of lung ventilation was available to compare it with the estimated ventilation parameters, the lower variability in test-retest showed improved robustness of the presented algorithm in comparison to the original one.A potential limitation of this study is the inherent assumption in the new algorithm, that most voxels initially included in the percentile range are adequate healthy reference voxels.Thus, in extreme cases where actual healthy voxels represent only a small fraction of voxels in the high RVent percentile range, these healthy voxels might get excluded from the reference in the third

F I G U R E 5
Bland-Altman plots (test-retest) of ventilation correlation coefficient (Vent-CC) and ventilation defect percentage (VDP) in chronic obstructive pulmonary disease (COPD) patients using the original and the new algorithm (with and without selection of the biggest cluster).

F I G U R E 6
Representative slices from a chronic obstructive pulmonary disease (COPD) patient and a healthy participant showing RVent images (top row) and ventilation correlation maps calculated using the original reference region of interest (ROI) algorithm (middle row) and the new algorithm (bottom row).
step of the presented algorithm (removal of voxels with inhomogeneous FVL), which in turn may lead to a higher representation of non-healthy ventilation curves in the reference, and hence, an underestimation of the true disease excess.However, this scenario seems highly unlikely even in end-stage lung disease patients as the Vent-CC parameter is a relative measure of ventilation heterogeneity in the lung parenchyma and the new algorithm uses the highest 80-90 percentile of RVent in each coronal slice.
The calculation of VDP was performed in this work by applying a global threshold to the ventilation correlation maps.This threshold has been previously used and was shown to be sensitive for the detection of ventilation defects in early-stage chronic lung allograft syndrome. 10ith the new algorithm for calculating the ventilation correlation maps, VDP was shown to have a lower test-retest coefficient of variation in healthy subjects and in the COPD cohort, however, the mean test-retest difference in the absolute VDP values was higher than with the original algorithm (Table 4).This could be explained by the use of a static global threshold that was optimized for the original algorithm.A more sophisticated thresholding method to calculate the VDP from the correlation maps, for example using defect distribution index, might further improve the accuracy of PREFUL derived VDP. 26 This needs to be proven in the future by comparison to a reference standard.
To conclude, this study presented a new algorithm for calculating an important parameter of regional ventilation dynamics (i.e., Vent-CC) using the contrast media-free PREFUL-MRI technique.The presented algorithm resulted in better test-retest reproducibility of Vent-CC and VDP in healthy participants and in COPD patients compared to the currently established PREFUL pipeline.
Bar plots of the mean internal correlation of flow curves (A) in healthy participants and (B) in chronic obstructive pulmonary disease (COPD) patients.Each bar represents a group of lung parenchyma voxels with RVent values in a specific percentile range (0th-10th percentile, 10th-20th percentile, and so on).Flow-volume loops of the respective voxel groups are shown (C) in healthy participants and (D) in COPD patients.

n = 28) Healthy (n = 32)
Included datasets are first scan sessions for all healthy patients.Shown values are the mean correlation coefficients with the SD in brackets.Red fields indicate a mean correlation below the threshold of 0.97.
Test-retest measurements of Vent-CC and VDP using the original and the new algorithm (with and without selection of the biggest cluster) in the COPD cohort and in healthy participants.
Abbreviations: COPD, chronic obstructive pulmonary disease; GOLD, global initiative for chronic obstructive lung disease;SE, standard error; ROI, region of interest; VDP, ventilation defect percentage; Vent-CC, ventilation correlation coefficient.