Pulmonary hypertension (PHT) is a complex pathophysiological disorder characterized by an elevation in pulmonary arterial pressure (PAP). This diagnosis is given when mean PAP exceeds 25 mmHg at rest using invasive right heart catheterization [1]. In early stages of the disease, elevations in PAP and architectural remodeling typically induce dilation and overload-induced hypertrophy of the right ventricle. Over time, this results in impairment of cardiac function, eventual right heart failure, and secondary pathophysiology such as right atrial dilation [2]. Although two-dimensional (2D) Doppler echocardiography is commonly used to provide evidence of the disease hallmarks in the right ventricle, and right ventricular catherization is needed for a conclusive diagnosis. Thus, other non-invasive diagnostic imaging biomarkers would be beneficial in clinical settings [3]. Magnetic resonance imaging (MRI) and angiography (MRA) could address this unmet need. Relative to 2D Doppler echocardiography, MRI possesses two advantages.

**Citation:** Hong, Z.M.; Garcia, J. Pulmonary Artery Remodeling and Advanced Hemodynamics: Magnetic Resonance Imaging Biomarkers of Pulmonary Hypertension. *Appl. Sci.* **2022**, *12*, 3518. https://doi.org/ 10.3390/app12073518

Academic Editor: Vladislav Toronov

Received: 23 February 2022 Accepted: 29 March 2022 Published: 30 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

First, it can perform multiplanar analysis allowing for accurate volumetric measurement of the vasculature. Second, it can evaluate hemodynamic changes in the pulmonary pathway via 4D-flow MRI. In both cases, volume and flow measurements demonstrate low variability and high test–re-test reliability [3,4]. Although MRI has improved the characterization of right ventricular pathophysiology, it has also identified potential biomarkers of the disease in the pulmonary artery (PA). Patients diagnosed with PHT often exhibit hemodynamic changes such as decreased vorticity, peak velocity, and wall shear stress (WSS) [5–7]. Furthermore, anatomic remodeling, including PA dilation and medial thickening, can be simultaneously observed [8]. Although promising, these novel biomarkers and their relationship must be investigated.

The purpose of this study was to test the hypothesis that hemodynamic changes within the pulmonary artery are associated with remodeling of the surrounding vascular architecture. In addition, we aimed to further characterize PHT biomarkers in the PA by investigating novel 3D flow biomarkers and bifurcation distances.

#### **2. Materials and Methods**

#### *2.1. Study Population*

A total of 51 subjects were recruited retrospectively, including 17 PHT patients and 33 controls. All subjects were registered within the Cardiovascular Imaging Registry of Calgary (CIROC). The study was approved by the University of Calgary Research Ethics Board and all subjects provided written informed consent. All research activities were performed in accordance with the Declaration of Helsinki. The study was coordinated by commercial software (cardioDITM, Cohesic Inc., Calgary, AB, Canada) for the routine capture of patient informed consent, for health questionnaires, and for standardized collection of MRI-related variables. PHT was clinically defined as a mean PA pressure ≥25 mmHg previously assessed by clinical invasive catheterization or echocardiography. Patients were recruited at the time of visit, with exclusions occurring for patients with prior surgical interventions in the pulmonary valve and/or artery; complex congenital heart disease; and/or contraindications to MRI. Control subjects (≥18 years of age) underwent similar workflow and confirmed no prior history of cardiovascular disease with a certified nurse. No catherization was performed in healthy controls. Prior to scanning, demographic measurements including age, sex, height, weight, and heart rate were obtained. Volume and mass measurements were normalized to body surface area using the Mosteller formula.

#### *2.2. Cardiac Magnetic Resonance Imaging Protocol*

Cardiac imaging examination was performed using 3T MRI scanners (Skyra and Prima, Siemens, Erlangen, Germany). Indication-based protocolling ensured consistent imaging procedures for all subjects, and cardiac imaging was performed in accordance with published recommendations [9]. Standard routine retrospective electrocardiographic gating, and timeresolved balanced steady-state free precession (SSFP) cine imaging in four-chamber, threechamber, two-chamber, and short-axis views of LV at end-expiration, was performed. Contrast usage of gadolinium contrast volume of 0.2 mmol/kg (Gadovist®, Bayer Inc., Mississauga, ON, Canada) was administrated to acquire a contrast-enhanced 3D magnetic resonance angiogram (CE MRA) of the cardiovascular structure. Time-resolved three-dimensional phase-contrast MRI with three-directional velocity encoding and retrospective ECG-gating (4D-flow, Siemens WIP 785A) was performed for 5–10 min, following contrast administration to measure in-vivo blood flow velocities within the whole heart. We have previously reported this whole-heart protocol [10–12]. Briefly, 4D-flow data was acquired during free breathing using navigator gating of diaphragmatic motion; sequence parameters were as follows: flip angle = 15 degrees, spatial resolution = 2.0−3.5 × 2.0−3.5 × 2.5−3.5 mm; temporal resolution = 39–48 ms; and velocity sensitivity = 150−250 cm/s. Total acquisition time varies between 5–10 min, depending on heart rate and respiratory navigator efficiency. The number of phases was adjusted to 25.

#### *2.3. Standard Cardiac Imaging Analysis*

Standard cardiac images were analyzed by a blinded reader to the study the same day of the acquisition using dedicated software cvi42 version 5.11.5 (Circle Cardiovascular Imaging Inc, Calgary, AB, Canada) to determine left and right end-diastolic volume (LVEDV; RVEDV), LV and RV end-systolic volume (LVSEV; RVESV), and LV and RV ejection fraction (LVEF; RVEF). Furthermore, the contrast-enhanced angiogram (CE MRA) was used to measure pulmonary artery diameters (LPA, MPA, and RPA), bifurcation distances, and bifurcation angles (Figure 1D).

**Figure 1.** Data processing and analysis workflow. Panel (**A**) shows an example of corrected 4D-flow images. Panel (**B**) illustrates a time-average phase-contrast (PC) magnetic resonance angiogram (MRA) with an optimal threshold for visualizing the anatomy. The pulmonary artery (PA) was segmented by isolating the corresponding volume and identifying the main PA (MPA), right PA (RPA), and left PA (LPA). Panel (**C**) shows the 2D analysis planes for basic flow quantification in the right ventricle outflow track (RVOT), MPA, RPA, and LPA. Panel (**D**) shows contrast-enhanced MRA (CE MRA) maximum intensity projection (MIP) and the corresponding 3D volume of a sample case. Multi-planar reconstruction (MPR) was used for planning of angle (in blue), distance (d in red), and diameter (location in pink, diameters in golden). Pink and golden lines in View 1 represent the projection planes corresponding to View 2 and 3, respectively. Distance samples (orange lines in View 2) were defined for MPA the closer point to the pulmonary valve to the PA bifurcation; right and left PA distances were defined from the PA bifurcation to the closer visible in-plane ramification. Panel (**E**) shows examples of sub-volume regions of MPA, RPA, and LPA; wall shear stress (WSS); energy loss (EL); kinetic energy (KE); and vortex core size (red arrows point to vortex cores identified by λ2).

#### *2.4. 4D-Flow Data Analysis*

All 4D-flow MRI data was pre-processed using in-house program developed in MAT-LAB 2020b (Mathworks, Natick, MA, USA), and the following tasks were performed: corrections for Maxwell terms, eddy currents, and aliasing, Figure 1A. After pre-processing, a 3D phase-contrast (PC) angiogram (PC MRA) was generated, Figure 1B. This angiogram was used to segment selected anatomical regions of the pulmonary artery, Figure 1B red inset, using in-house MATLAB-based tool "4D-Flow Analysis Tool" [13,14]. Analysis planes were created at the right ventricular outflow tract (RVOT), main pulmonary artery (MPA), the left pulmonary artery (LPA), and right pulmonary artery (RPA) using a specialized visualization software (Ensight 10.2, CEI Inc, Research Triangle Park, North Carolina, USA) to measure peak velocity, net flow, retrograde flow, regurgitation fraction, and forward flow, Figure 1C. Furthermore, the PC MRA was used to measure pulmonary artery diameters (LPA, MPA and RPA), bifurcation distances, and bifurcation angles (Figure 1D) in cvi42 when standard CE MRA was not available in controls (*n* = 2). Volume sub-regions were created for local advanced 4D-flow analysis. These volume subregions were used to obtain the maximum velocities from the LPA, PA, and RPA using the approach proposed by Rose et al. [15].

Regional WSS was calculated at peak systole as previously described [16,17]. Energy loss (EL) was calculated as proposed by Barker et al. [18], which provides the rate of EL in a volume of interest (i.e., power in Watts) at a given time of the cardiac cycle, and peak systole was used in this study. Kinetic energy (KE) was calculated as reported by Geeraert et al. [13] in each sub-region at peak systole. Vortex size analysis at peak systole was performed as introduced by Garcia et al. [10] using the Lambda2 (λ2) method, which is an accepted method to identify 3D vortices. Valid vortex networks were defined as those >75 mm3 (considering a voxel resolution of 2.5 × 2.5 × 3.0 mm) with a shared connectivity neighborhood ≥4 voxels. This rule reduces the impact of possible incorrect derived vortices due to noise and/or partial derivative errors, which are directly impacted by 4D-flow spatial resolution. Examples for these advanced 4D-flow parameters are presented in Figure 1E.

#### *2.5. Statistical Analysis*

Statistical analysis was performed using SPSS 25 (SPSS, Chicago, IL, USA). Normality was assessed using normal plots and the Shapiro–Wilk test. Between-group comparisons for demographic, baseline cardiac MRI, hemodynamics, and anatomical measurements were performed using an independent samples t-test incorporating Levene's test for equality of variances. To assess the relationship between variables, a series of linear regression tests was performed using Pearson R-Values. If a significant relationship was found, a univariate linear model was generated using variables demonstrating statistical significance during the linear regression tests. Pearson correlation coefficients were also generated in relation to LVEF and RVEF. The significance threshold was <0.05.

#### **3. Results**

#### *3.1. Patient Characteristics*

Patient and control group characteristics inclusive of demographic and cardiac function data are summarized in Table 1. There were differences in age (55 ± 17 years vs. 39 ± 15 years, *p* = 0.001), sex (76% vs. 31% female, *p* = 0.002), height (1.64 ± 0.14 m vs. 1.73 ± 0.10 m, *p* = 0.023), and heart rate (76 ± 16 bpm vs. 66 ± 11 bpm, *p* = 0.016) between PHT patients and controls. BSA (1.88 ± 0.28 m2 vs. 1.92 ± 0.34 m2, *p* = 0.643) and weight (78 ± 18 kg vs. 79 ± 21 kg, *p* = 0.791) did not differ. Of the cardiac function parameters, PHT patients demonstrated lower LVEDV (127 ± 48 mL vs. 160 ± 35 mL, *p* = 0.011). However, other measurements, including LVEF (59 ± 13% vs. 61 ± 6%, *p* = 0.752), LVESV (58 ± 42 mL vs. 66 ± 23 mL, *p* = 0.381), RVEF (56 ± 7% vs. 55 ± 30%, *p* = 0.878), RVEDV (190 ± 79 mL vs. 173 ± 48 mL, *p* = 0.456), and RVESV (106 ± 72 mL vs. 77 ± 28 mL, *p* = 0.152), did not differ.

**Table 1.** Demographic, cardiac function baseline, and medications.

Data are expressed as mean ± SD. SD: standard deviation. LVEF: left ventricle ejection fraction. LVEDV: LV end-diastolic volume. LVESV: LV end-systolic volume. RVEF: right ventricle ejection fraction. RVEDV: RV end-diastolic volume. RVESV: RV end-systolic volume.

#### *3.2. Anatomical Remodeling*

Anatomical measurements of the PA are detailed in Table 2. Of the regions measured, only the MPA (AP: 32 ± 6 mm vs. 30 ± 4 mm, *p* = 0.050; RL: 32 ± 6 mm vs. 28 ± 3 mm, *p* = 0.005) and the RPA-AP (21 ± 4 mm vs. 19 ± 3 mm, *p* = 0.050) demonstrated dilatation in PHT patients. The LPA (AP: 20 ± 4 mm vs. 19 ± 4 mm, *p* = 0.370; CC: 21 ± 4 mm vs. 19 ± 3 mm, *p* = 0.226) and the RPA-CC (22 ± 4 mm vs. 21 ± 3 mm, *p* = 0.558) did not differ. The bifurcation distances of PHT patients were greater in the branches (LPA: 34 ± 6 mm vs. 29 ± 8 mm, *p* = 0.043; RPA: 34 ± 5 mm vs. 27 ± 6 mm, *p* < 0.001), but not the MPA (33 ± 5 vs. 30 ± 8 mm, *p* = 0.302) itself. Bifurcation angles did not differ (93 ± 9 degrees vs. 98 ± 15 degrees, *p* = 0.208).

**Table 2.** Anatomical measurements.

Data are expressed as mean ± SD. SD: standard deviation. LPA: left pulmonary artery. MPA: main PA. RPA: right PA. deg: degrees. AP: anterior–posterior. CC: cranio-caudal.

#### *3.3. Hemodynamic Characteristics*

4D flow data is reported in Table 3, and sample datasets are illustrated in Figure 2. Maximum EL was lower in PHT patients in all regions of the PA (LPA: 1.06 ± 0.75 μW vs. 2.12 ± 1.22 μW, *p* = 0.002; MPA: 2.76 ± 1.20 μW vs. 4.30 ± 1.93 μW, *p* = 0.001; PA: 3.94 ± 1.71 μW vs. 7.27 ± 3.24 μW, *p* < 0.001; RPA: 0.88 ± 0.70 μW vs. 2.06 ± 1.22 μW, *p* = 0.001) relative to controls. Although mean EL did not differ in the MPA (1.24 ± 0.56 μW vs. 1.27 ± 0.45 μW, *p* = 0.800) or the entirety of the PA (1.76 ± 0.83 μW vs. 2.00 ± 1.09 μW, *p* = 0.439), it was increased in PHT patients in the LPA (1.76 ± 0.39 μW vs. 0.49 ± 0.51 μW, *p* = 0.022) and decreased in the RPA (0.06 ± 0.24 μW vs. 0.46 ± 0.56 μW, *p* = 0.001).

Data are expressed as mean ± SD. PA: pulmonary artery. LPA: left pulmonary artery. MPA: main PA. RPA: right PA. EL: energy loss. WSS: wall shear stress. KE: kinetic energy.

**Figure 2.** Sample dataset of advanced 4D-flow based measurements. Top raw shows a healthy control, man, 36-year-old. Bottom raw shows a pulmonary hypertension (PHT) patient, man, 58-year-old. WSS: wall shear stress; EL: energy loss; KE: kinetic energy.

Maximum WSS for PHT patients was lower in all regions (LPA: 0.75 ± 0.29 Pa vs. 1.17 ± 0.26 Pa, *p* < 0.001; PA: 0.98 ± 0.28 Pa vs. 1.28 ± 0.28 Pa, *p* = 0.001; RPA: 0.64 ± 0.26 Pa vs. 1.08 ± 0.30 Pa, *p* < 0.001), except the MPA (0.97 ± 0.28 Pa vs. 1.06 ± 0.18 Pa, *p* = 0.224). Mean WSS (LPA: 0.14 ± 0.05 Pa vs. 0.19 ± 0.05 Pa, *p* = 0.001; MPA: 0.14 ± 0.03 Pa vs. 0.18 ± 0.04 Pa, *p* = 0.022; PA: 0.14 ± 0.03 Pa vs. 0.18 ± 0.04 Pa, *p* = 0.002; RPA: 0.12 ± 0.04 Pa vs. 0.17 ± 0.05 Pa, *p* < 0.001) and the mean of the highest 5% WSS values (LPA: 0.48 ± 0.20 Pa vs. 0.76 ± 0.15 Pa, *p* < 0.001; MPA: 0.52 ± 0.13 Pa vs. 0.68 ± 0.14 Pa, *p* < 0.001; PA: 0.52 ± 0.14 Pa vs. 0.74 ± 0.14 Pa, *p* < 0.001; RPA: 0.40 ± 0.16 Pa vs. 0.70 ± 0.18 Pa, *p* < 0.001) were lower in all regions of the PA. The mean of all WSS values above 95% was decreased in the LPA (0.73 ± 0.28 Pa vs. 1.14 ± 0.26 Pa, *p* < 0.001) and RPA (0.63 ± 0.25 Pa vs. 1.03 ± 0.32 Pa, *p* < 0.001), but not the MPA (0.94 ± 0.27 Pa vs. 1.03 ± 0.17 Pa, *p* = 0.234) or the entirety of the PA (0.96 ± 0.27 Pa vs. 1.51 ± 1.44 Pa, *p* = 0.123) relative to controls. Maximum vortex volume was only increased in the MPA (27.08 ± 10.66 mm3 vs. 20.33 ± 7.09 mm3, *p* = 0.01) but otherwise did not differ (LPA: 5.60 ± 3.01 mm3 vs. 4.44 ± 1.85 mm3, *p* = 0.090; PA: 33.99 ± 12.14 mm3 vs. 27.47 ± 7.08 mm3, *p* = 0.053; RPA: 4.11 ± 2.22 mm3 vs. 3.70 ± 1.89 mm3, *p* = 0.493). Mean vortex volume was increased in all regions (LPA: 3.55 ± 2.04 mm3 vs. 2.14 ± 1.11 mm3, *p* = 0.014; MPA: 19.00 ± 8.71 mm3 vs. 11.97 ± 4.76 mm3, *p* = 0.005; PA: 24.55 ± 10.60 mm3 vs. 15.58 ± 5.41 mm3, *p* = 0.004) except for the RPA (2.61 ± 1.41 mm3 vs. 2.09 ± 1.19 mm3, *p* = 0.175). Maximum KE attained was also lower in PHT patients (LPA: 1.1 ± 0.7 mJ vs. 2.2 ± 1.2 mJ, *p* < 0.001; MPA: 2.8 ± 1.1 mJ vs. 4.2 ± 1.9 mJ, *p* = 0.002; PA: 4.0 ± 1.7 mJ vs. 7.3 ± 3.3 mJ, *p* < 0.001; RPA: 0.7 ± 0.5 mJ vs. 2.0 ± 1.2 mJ, *p* < 0.001).

Peak velocity did not differ in any region of the PA (RVOT: 0.71 ± 0.26 m/s vs. 0.82 ± 0.17 m/s, *p* = 0.145; LPA: 0.71 ± 0.27 m/s vs. 0.81 ± 0.25 m/s, *p* = 0.210; MPA: 0.72 ± 0.26 m/s vs. 0.84 ± 0.22 m/s, *p* = 0.092; RPA: 0.69 ± 0.31 m/s vs. 0.81 ± 0.18 m/s, *p* = 0.173). Forward flow was lower in the RPA (26.08 ± 13.12 mL vs. 33.82 ± 12.08 mL, *p* = 0.042) of PHT patients but otherwise did not differ (RVOT: 74.86 ± 37.25 mL vs. 83.11 ± 23.11 mL, *p* = 0.413; LPA: 43.10 ± 27.79 mL vs. 43.15 ± 13.29 mL, *p* = 0.995; MPA: 80.52 ± 40.39 mL vs. 89.27 ± 19.41 mL, *p* = 0.408). Similarly, net flow was also lower in the RPA (25.84 ± 13.35 mL vs. 33.28 ± 11.83 mL, *p* = 0.049) of PHT patients but otherwise did not differ (RVOT: 72.38 ± 36.23 mL vs. 81.06 ± 23.71 mL, *p* = 0.380; LPA: 43.04 ± 27.84 mL vs. 42.58 ± 13.52 mL, *p* = 0.949; MPA: 79.70 ± 40.75 mL vs. 88.25 ± 20.00 mL, *p* = 0.424). Retrograde flow (0.06 ± 0.14 mL vs. 0.57 ± 1.15 mL, *p* = 0.017) was lower in the LPA of PHT patients but otherwise did not differ (RVOT: 2.48 ± 5.89 mL vs. 2.05 ± 2.48 mL, *p* = 0.715; MPA: 0.81 ± 1.01 mL vs. 1.02 ± 1.65 mL, *p* = 0.637; RPA: 0.24 ± 0.65 mL vs. 0.55 ± 0.79 mL, *p* = 0.162). Furthermore, regurgitation fractions were decreased in the LPA (0.33 ± 0.83% vs. 1.57 ± 3.18% *p* = 0.043) and RPA (0.23 ± 0.54% vs. 1.56 ± 1.89%, *p* = 0.001) but did not differ in the RVOT (3.48 ± 5.60% vs. 2.92 ± 4.20%, *p* = 0.695) or MPA (1.78 ± 3.78% vs. 1.32 ± 2.21%, *p* = 0.588).

#### *3.4. Univariate Analysis*

Correlational analyses were reported in Table 4. LVESV correlated with LVEF (R = −0.610, *p* < 0.001). RVEDV correlated with LVEF and RVEF (R = −0.454, *p* = 0.005; R = −0.610, *p* < 0.001). A similar correlation was observed for RVESV (R = −0.677, *p* < 0.001;R= −0.864, *p* < 0.001). MPA-AP (R = 0.385, *p* = 0.017) and MPA-RL (R = 0.418, *p* = 0.009) are significantly correlated to RVEF but not LVEF (MPA-AP: R = 0.097, *p* = 0.520; MPA-RL: R = 0.079; *p* = 0.603). All peak velocity measurements were significantly correlated with LVEF (RVOT: R = 0.451, *p* = 0.002; MPA: R = 0.489, *p* = 0.001; LPA: R = 0.454, *p* = 0.002; RPA: R = 0.527, *p* < 0.001) and RVEF (RVOT: R = 0.416, *p* = 0.009; MPA: R = 0.373, *p* = 0.021; LPA: R = 0.401, *p* = 0.013; RPA: R = 0.420, *p* = 0.009). The MPA RF (R = −0.447, *p* = 0.002), MPA max WSS (R = 0.519, *p* < 0.001) and MPA mean WSS (R = 0.393, *p* = 0.007) correlated with LVEF. LPA mean WSS also correlated with LVEF (R = 0.306, *p* = 0.038). MPA mean Vortex (R = −0.335, *p* = 0.023) and MPA max KE (R = 0.332, *p* = 0.024) also showed a correlation with LVEF.

**Table 4.** Univariate analysis for ejection fraction.

Data are expressed as mean ± SD. BSA: body surface area. SBP: systolic blood pressure. DBP: diastolic blood pressure. LVEDV: left ventricle end-diastolic volume. LVESV: left ventricle end-systolic volume. LVEF: left ventricle ejection fraction. RVEDV: right ventricle end-diastolic volume. RVESV: right ventricle end-systolic volume. RVEF: right ventricle ejection fraction. RVOT: right ventricular outflow track. LPA: left pulmonary artery. MPA: main PA. RPA: right PA. AP: anterior–posterior. CC: cranio-caudal. RL: right Left. RF: regurgitant fraction. EL: energy loss. WSS: wall shear stress. KE: kinetic energy.

#### *3.5. Sex Matching for Women*

Sex matching analysis is reported in Table 5. When matching a women subset of controls and patients (*n* = 10), significant differences (*p* < 0.05) were identified for weight; for MPA-AP, MPA-RL, RPA-AP diameters; and for RPA distance. Similarly, advanced parameters showed differences for mean EL MPA (*p* = 0.025), RPA max WSS (*p* = 0.003), LPA max WSS (*p* = 0.019), MPA max vortex (*p* = 0.003), MPA mean vortex (*p* = 0.001), RPA mean vortex (*p* = 0.037), LPA max vortex (*p* = 0.017), LPA mean vortex (*p* = 0.001), and RPA max KE (*p* = 0.004).

**Table 5.** Matching analysis for women.

Data are expressed as mean ± SD. BSA: body surface area. SBP: systolic blood pressure. DBP: diastolic blood pressure. LVEDV: left ventricle end-diastolic volume. LVESV: left ventricle end-systolic volume. LVEF: left ventricle ejection fraction. RVEDV: right ventricle end-diastolic volume. RVESV: right ventricle end-systolic volume. RVEF: right ventricle ejection fraction. RVOT: right ventricular outflow track. LPA: left pulmonary artery. MPA: main PA. RPA: right PA. AP: anterior–posterior. CC: cranio-caudal. RL: right Left. RF: regurgitant fraction. EL: energy loss. WSS: wall shear stress. KE: kinetic energy.

#### **4. Discussion**

In this study, we demonstrate that PHT is associated with both hemodynamic and anatomic changes throughout the PA. Patients with PHT demonstrated multiple changes across a variety of hemodynamic parameters, including EL, WSS, vorticity, flow direction and flow velocity. In addition, pulmonary artery dilation was observed alongside changes in bifurcation geometry.  Taken together, these findings provide additional evidence that hemodynamic and anatomic biomarkers possess clinical relevance for diagnosis of the disease.

Although hemodynamic changes are emerging biomarkers of the disease, they remain poorly understood. Previous research suggests patients with PHT demonstrate lower wall shear stress [5,7]. Although these studies were conducted using an adult PHT population, the findings possess relevance to the broader definition of PHT. For vorticity, a recent study conducted by Kroeger et al. found that vortices are observed in healthy controls and patients with PHT [19]. Similarly, Schafer et al. identified reduced vorticity but increased helicity in patients with PHT [6]. For EL, Han et al. demonstrated greater total energy losses over the cardiac cycle throughout the entirety of the PA [20] in PHT patients. PHT patients may also demonstrate greater retrograde flow and lower peak velocity [21,22]. In contrast, anatomical changes are more well-studied using CT, identifying dilation of the main pulmonary artery in various stages of the disease [23,24]. These findings have been replicated using MRI and appear to be prevalent in the later stages [25].

Our findings generally corroborate results of previous studies and provide insight into novel biomarkers. We observed lower WSS, maximum energy attained, and regurgitation fractions in the LPA and RPA in patients with PHT. Maximum EL at peak systole remained lower in patients than in controls in all PA sections. Mean EL tended to be also lower, except for the LPA. EL reported in our study was quantified at peak systole, which provides further support in addition of previous findings reported over the cardiac cycle. Contrary to other studies, we found that vorticity tended to increase in PHT patients. For anatomical changes, we observed dilation in the MPA and part of the RPA. Of note, we are the first study to our knowledge to identify greater bifurcation distances associated with PHT. We speculate this could be an additional form of anatomical remodeling associated with disease progression. However, the exact significance is unknown. Further studies investigating changes in bifurcation geometry at various disease stages are needed to better understand their significance. The difference between the RPA and LPA characterization may be influenced by the patient's anatomy and the hemodynamic resistance from the lungs. However, the latter requires additional investigation. Furthermore, sex-matching analysis confirmed the impact of advanced hemodynamic parameters in PHT patients and highlighted the importance of subject matching in future studies.

This study has a few limitations. First, our study was limited by the spatiotemporal resolution constraints of MRI. For example, while we could visualize larger vortex networks, smaller vortices occurring at faster speeds exceeding the available resolution could be missed. Second, as a single centre study, the addition of multiple centres would improve the generalization of results. Furthermore, due to its retrospective nature, the sample size and power is limited. Third, our control patients were not age and sex-normalized to the patient cohort, which suggests normal cardiac aging and sex differences could have impacted results. Fourth, PHT is a progressive multi-factorial condition with diverse underlying etiology. Clinical protocols did not mandate further characterization of the disease into smaller subgroups such as pulmonary arterial hypertension, or patients at varying stages of the disease, limiting the characterization of disease progression. The relatively small cohort size did not allow one to explore the effect of medical therapy in 4D-flow derived parameters. The latter remains an important aspect to assess in future studies. Finally, we did not explore the right ventriculo-pulmonary arterial coupling in the present cohort. This coupling analysis may produce a better understanding of the bidirectional effects between the RV and the pulmonary artery tree.

#### **5. Conclusions**

Adult PHT patients demonstrate hemodynamic and anatomic differences throughout the PA. In addition to providing further support for established disease phenotypes such as reduced WSS and PA dilatation, these findings suggest that novel biomarkers such as energy changes and changes in bifurcation geometry may possess clinical relevance for PHT patients. Future studies are needed to investigate these phenomena at varying stages of the disease and determine how they may be affected by underlying etiology.

**Author Contributions:** Conceptualization, Z.M.H. and J.G.; methodology, J.G.; software, J.G.; validation, Z.M.H. and J.G.; formal analysis, Z.M.H. and J.G.; investigation, Z.M.H. and J.G.; resources, J.G.; data curation, Z.M.H. and J.G.; writing—original draft preparation, Z.M.H. and J.G.; writing—review and editing, J.G.; visualization, Z.M.H. and J.G.; supervision, J.G.; project administration, J.G.; funding acquisition, J.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by The University of Calgary, URGC SEM #1054341; J.G. start-up funding. We acknowledge the support of the Natural Science and Engineering Research Council of Canada/Conseil de recherche en science naturelles et en génie du Canada, RGPIN-2020-04549 and DGECR-2020-00204.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Conjoint Health Research Ethics Board of University of Calgary (REB13-0902 approved on 6/18/2014 and currently active.

**Informed Consent Statement:** Written informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The anonymized data presented in this study are available on request from the corresponding author. The data are not publicity available due to privacy and ethical restrictions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Impact of Respiratory Gating on Hemodynamic Parameters from 4D Flow MRI**

**Esteban Denecken 1,2,3,4, Julio Sotelo 1,3,4,5, Cristobal Arrieta 1,3, Marcelo E. Andia 1,3,4,6 and Sergio Uribe 1,3,4,6,\***

	- Santiago 782 0436, Chile

**Abstract:** The hemodynamic parameters from 4D flow datasets have shown promising diagnostic value in different cardiovascular pathologies. However, the behavior of these parameters can be affected when the 4D flow data are corrupted by respiratory motion. The purpose of this work was to perform a quantitative comparison between hemodynamic parameters computed from 4D flow cardiac MRI both with and without respiratory self-gating. We considered 4D flow MRI data from 15 healthy volunteers (10 men and 5 women, 30.40 ± 6.23 years of age) that were acquired at 3T. Using a semiautomatic segmentation process of the aorta, we obtained the hemodynamic parameters from the 4D flow MRI, with and without respiratory self-gating. A statistical analysis, using the Wilcoxon signed-rank test and Bland–Altman, was performed to compare the hemodynamic parameters from both acquisitions. We found that the calculations of the hemodynamic parameters from 4D flow data that were acquired without respiratory self-gating showed underestimated values in the aortic arch, and the descending and diaphragmatic aorta. We also found a significant variability of the hemodynamic parameters in the ascending aorta of healthy volunteers when comparing both methods. The 4D flow MRI requires respiratory compensation to provide reliable calculations of hemodynamic parameters.

**Keywords:** hemodynamics; 4D flow; self-gating; aorta

#### **1. Introduction**

Several strategies have been proposed to suppress respiratory motion in cardiac MRI. Most of these techniques are used to track the position of the diaphragm and can be combined with a respiratory-based k-space reordering scheme [1–6]. These techniques include pencil beam navigators [1,2], cross-pair excitation to acquire a column of pixels across the lung-to-liver interface, self-gating strategies based on k0 points [3], k0 profiles [4] or 2D image navigators [5].

Newer approaches have achieved respiratory resolved cardiac images with a near 100% respiratory navigator efficiency [7,8]. These approaches have been used to obtain singleor dual-phase 3D whole heart balanced steady state free precession (b-SSFP) images [2], 2D or 3D b-SSFP cine images [9], and 4D flow MRI data [4]. The use of cine b-SSFP respiratory compensated data has remarkably improved the image quality, as well as the anatomical and functional accuracy [7,8]. Nevertheless, there is a debate whether the use of respiratory motion compensation techniques improves the quantification of velocity-derived parameters from the 4D flow MRI [10].

**Citation:** Denecken, E.; Sotelo, J.; Arrieta, C.; Andia, M.E.; Uribe, S. Impact of Respiratory Gating on Hemodynamic Parameters from 4D Flow MRI. *Appl. Sci.* **2022**, *12*, 2943. https://doi.org/10.3390/ app12062943

Academic Editor: Julio Garcia Flores

Received: 15 December 2021 Accepted: 26 January 2022 Published: 14 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Applying respiratory compensation to 4D flow MRI increases the already long scanning time of this sequence. Moreover, its impact on the accuracy of the hemodynamic parameters remains unclear. Some authors propose that the use of respiratory self-gating techniques has a clear impact on 4D flow measurements by looking at the stroke volume [4,11]. Other studies propose that aorta flow volume quantification does not need respiratory compensation [12]. Other approaches that used profile reordering for respiratory compensation in 4D flow MRI [13,14] suggest that the use of a motion compensation technique increases the accuracy of 4D flow measurements. Most of this discussion has been based on the net flow going through a vessel of interest. However, the net flow is a parameter that averages the velocities in a region of interest, which may hide the distortion of respiratory motion on the velocity field.

Over recent years, several parameters from 4D flow datasets have been derived, such as Wall Shear Stress (WSS), vorticity, helicity, kinetic energy, and viscous dissipation, among others [10,15]. These parameters have been applied in several groups of patients and have shown promising diagnostic value in different cardiovascular pathologies. Nevertheless, these parameters need to be carefully used and compared since their value depends on the resolution and the level of noise of the acquisitions [16]. We also know little about the behavior of these novel parameters when the 4D flow data are corrupted by respiratory motion.

Some of these new hemodynamic parameters are directly proportional to velocity, but others are related to spatial derivatives. Therefore, any rapid perturbation of the velocity field that is introduced by noise or motion could impact the precision and accuracy of the computations. We hypothesized that there are significant differences in the hemodynamic parameters obtained from 4D flow MRI that are acquired with and without respiratory gating in healthy volunteers. The purpose of this work was to perform a quantitative comparison of hemodynamic parameters obtained from 4D flow MRI data that were acquired both with and without respiratory self-gating. We compared the differences between the following hemodynamic parameters, computed in the thoracic aorta: velocity; WSS; WSS axial (WSSA) and circumferential (WSSC) components; vorticity; viscous dissipation; kinetic energy; and energy loss. We also considered the effect of the segmentation process, comparing the segmentations from respiratory self-gated and non-gated data using the DICE as a global index and computing mesh distances to study the local differences.

#### **2. Materials and Methods**

#### *2.1. Data Acquisition*

This study included thoracic aorta 4D flow MRI data from 15 healthy volunteers (10 men and 5 women, 30.40 ± 6.23 years of age), which had been previously acquired. The study was conducted according to the guidelines of the Declaration of Helsinki and was approved by the Ethics Committee at Pontificia Universidad Católica (protocol code 171,127,001, 20 June 2019). A 3T Philips MRI scanner was used to acquire the data, with and without respiratory self-gating. The respiratory self-gated strategy consists of the continuous acquisition of the k0 profile and the cross-correlation to determine the respiratory positions [4]. A window of 7 mm was used to collect the data during expiration [4]. A summary of the acquisition parameters is shown in Table 1. The non-respiratory gating scan was acquired with two repetitions (NSA, NEX = 2) in order to have a similar acquisition time to the respiratory self-gated scan.

#### *2.2. Segmentation and Data Processing*

Both 4D Flow MRI datasets (with and without respiratory self-gating) were postprocessed using an in-house MATLAB toolbox (MathWorks, Natick, MA, USA) [17]. The image processing steps involved the segmentation of the thoracic aorta and the generation of a tetrahedral finite element mesh (characteristic length 1 mm), as described in [18]. The velocity, WSS, WSSA, WSSC, vorticity, viscous dissipation, kinetic energy, and energy loss were obtained from the thoracic aorta 4D flow images, using a least square projection

method based on the finite elements described in [18–20]. All parameters were obtained at a representative peak systolic cardiac phase. These parameters have shown potential value in diagnostics [21]. Angiographic 4D flow images were segmented using a semiautomatic process based on the thresholding, labeling, and manual separation of the vessels. The analysis was performed at the peak systolic phase.

The computations of the hemodynamic parameters were performed using the in-house MATLAB toolbox developed by Sotelo et al. [17] and were exhaustively validated [18–20]. The velocity values were transferred to each node of the mesh from the 3D PC-MRI datasets using a cubic interpolation [18].

**Table 1.** The MRI acquisition parameters.

#### *2.3. Statistical Analysis*

We assessed the differences between the hemodynamic parameters that were obtained with and without respiratory self-gating over four different regions of the aorta: region I = ascending aorta (between the sinus of Valsalva and before the brachiocephalic trunk); region II = aortic arch (before the brachiocephalic trunk and isthmus level); region III = superior descending aorta (between the isthmus level and the Valsalva level); and region IV = posterior descending aorta (between the Valsalva level and the diaphragmatic aorta). We computed the mean of the measurements in each region to perform the statistical analysis that is presented herein.

Furthermore, we calculated the DICE coefficient between the segmentations of the data that were obtained with and without respiratory self-gating and the local differences between the mesh surfaces of each section. The local differences between both meshes were calculated using the Euclidean distance, as reported in [22].

We analyzed the mean, standard deviation, and coefficient of variation (CV) of each hemodynamic parameter. The Wilcoxon signed-rank test and Bland–Altman plots were also performed between the respiratory self-gated and non-gated data, using a statistical significance of 0.05. We computed the Wilcoxon signed-rank test using the SPSS statistical analysis software and the Bland–Altman plots using MATLAB.

#### **3. Results**

Figure 1 shows a representative case of the hemodynamic parameters in the aorta that was obtained from respiratory and non-respiratory gated 4D flow MRI data. By visual inspection, the images look similar. Nevertheless, differences can be observed in the WSS, viscous dissipation, energy loss, and kinetic energy.

Table 2 shows the mean and CV obtained for each hemodynamic parameter at each section of the aorta. The mean of the hemodynamic parameters was similar in section I, but it showed differences in the other sections. Further, the CV values were greater in the 4D flow data acquired without respiratory self-gating compared to the self-gated data. The CV values were greater in section I.

**Figure 1.** The hemodynamic parameters obtained from 4D flow MRI, both with and without selfgating. In this figure, we show visual representations of the hemodynamic parameters for one volunteer, which were measured at a representative peak systolic cardiac phase.

**Table 2.** The hemodynamic parameter means and the coefficient of variation (CV) of the mean values at each section of the vessel for the self-gated and non-gated data, which were measured at a representative peak systolic cardiac phase.

Table 3 shows the results obtained from the statistical analysis. The Wilcoxon signedrank test showed significant differences between respiratory self-gating and non-respiratory gating in sections II, III, and IV. Differences in velocity, WSS, and WSSA were observed in sections II, III, and IV. There were also statistical differences for viscous dissipation, kinetic energy, and energy loss at sections II and IV.

**Table 3.** The Wilcoxon signed-rank test *p*-values for the mean values of the hemodynamic parameters at each section of the vessel, between the self-gated and non-gated data.

\* Existing significative differences.

In Appendix A, Table A1 shows the mean differences and 1.96 times the standard deviation of the Bland–Altman analysis. The 4D flow data with self-gating are used as a reference. The Bland–Altman plots for the velocity and WSS in section I (Figure 2) showed a greater dispersion than sections II, III, and IV. They also show the underestimated values from the data acquired without self-gating. The Bland–Altman plots for velocity, WSS, WSSA, WSSC, vorticity, viscous dissipation, kinetic energy, and energy loss plots are shown in Figure A1. In general, all computations of the hemodynamic parameters acquired without self-gating were underestimated when no respiratory gating was performed. There was a greater difference in velocity, WSS, WSSA, and kinetic energy, and this bias increased along the aortic sections. The standard deviation was greater for most of the parameters in section I compared to other sections.

**Figure 2.** The Bland–Altman plots for velocity and Wall Shear Stress at each section. The difference in Bland–Altman plots was computed as *SG* − *NG*, where *SG* is the data with self-gating and *NG* is the data without self-gating. Appendix A shows the Bland–Altman plots of all hemodynamic parameters considered in this work.

Finally, the mean and range of the DICE coefficient between the segmentation of the self-gated and non-gated data were 0.918 ± 0.056. This result, which considered the whole segmented aorta, showed a high intersection between both segmentations.

A graphical analysis of the local differences [22] between the meshes that were generated from both segmentations shows that the disagreements between them are located at the extremes of the aorta (Figure A2A in Appendix B for one volunteer). The local differences of the segmentations at each section had small variations (Table A2 in Appendix B) and the histogram shows that 88.36% of the differences at each section corresponded to distances lower than or equal to 1 mm (Figure A2B). The local differences were also similarly distributed around the 0% error and thus, there was not a consistent bias between both segmentations.

#### **4. Discussion**

In this work, we compared the hemodynamic parameters computed from thoracic aorta 4D flow MRI data that were obtained with and without respiratory self-gating. The statistical descriptors showed that the data obtained with self-gating had less variability than the data obtained without self-gating (Table 2). The Wilcoxon signed-rank test showed significant differences in most of the parameters calculated between both acquisitions at sections II, III, and IV and minimal differences at section I (ascending aorta). Nevertheless, the computations in section I showed a greater variability.

The Bland–Altman analysis showed that the parameters computed without respiratory self-gating were underestimated in comparison to the parameters computed with self-gating and that the underestimation increased along sections II, III, and IV. Section I did not show bias; however, it presented a greater dispersion of the differences for each hemodynamic parameter (see the Bland–Altman plots in Figure A1 in Appendix A). Therefore, for section I, there was not a clear trend of the under- or overestimation of the parameters that were computed. The differences for each volunteer were greater in section I compared to the other sections, which is an important result that is masked by the averaging of all volunteers.

The only parameter that did not show any statistical difference was the WSSC, which is probably due to the near to 0 value of the WSSC in this group of volunteers. An assessment of the impact of respiratory motion on this parameter in a cohort of volunteers or patients with a greater value of WSSC is necessary. As shown by Sotelo et al., there are significant differences in WSSC between healthy volunteers and BAV patients [19].

We also evaluated the impact of the segmentation process as a potential cause of the statistical differences. However, the DICE coefficient and analyses of the local differences at each section showed high agreement between both segmentations and the differences were minimal and unbiased. Therefore, the segmentation did not influence the differences that were observed in the hemodynamic parameters.

A limitation of this study is that we only included healthy volunteers with laminar flow patterns. The next step is to perform a similar study in a group of patients with non-laminar flow patterns or high Reynolds numbers, such as patients with a bicuspid aortic valve (BAV) or with aortic coarctation. Unfortunately, due to the long acquisition time, acquiring two datasets of 4D flow MRI in patients is challenging. Another limitation of the study is that we only assessed the impact of respiratory gating on the hemodynamic parameters at peak systole. The main reason for this was that the poor signal-to-noise ratio of the other cardiac phases makes the segmentation process very challenging. Although it is expected that respiratory motion artifacts may have a stronger effect at lower velocities, this still needs to be proven. Finally, flow parameters as net or regurgitant flow were not included because these results have been presented before using the same datasets [4]. The purpose of this paper was to assess how other hemodynamic parameters are affected by respiratory motion.

In conclusion, we found that the computations of the hemodynamic parameters from 4D flow data using self-gating acquisition showed statistical differences from the computations from acquisitions without respiratory compensation for the aortic arch and the descending and diaphragmatic aorta in healthy volunteers. The Bland–Altman plots support this fact, showing an underestimation of the hemodynamic parameters acquired without any respiratory compensation. Our results suggest that the use of respiratory motion compensation techniques is mandatory for the accurate evaluation of 4D flow hemodynamic parameters.

**Author Contributions:** Conceptualization, S.U.; methodology, J.S., M.E.A., S.U. and E.D.; software, J.S.; validation, J.S. and M.E.A.; formal analysis, E.D.; resources, S.U.; data curation, E.D.; writing original draft preparation, E.D.; writing—review and editing, S.U., J.S. and C.A.; visualization, J.S. and E.D.; supervision, S.U.; project administration, S.U.; funding acquisition, S.U. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by ANID—Millennium Science Initiative Program—ICN2021\_004, ANID-PIA-ACT192064, ANID—Millennium Science Initiative Program—NCN17\_129 and FONDE-CYT # 1181057. Sotelo J. would like to thank CONICYT-FONDECYT Postdoctorado 2017 #3170737 and FONDECYT de Iniciación en Investigación 2020 #11200481. Andia M. would like to thank FONDECYT # 1180525. Arrieta C. was partially funded by CONICYT FONDECYT Postdoctorado 2019 #3190763, CONICYT PCI REDES 180090, and ANID—Millennium Science Initiative Program— NCN17\_129.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki and was approved by the Ethics Committee at Pontificia Universidad Catolica (protocol code 171127001, 20 June 2019).

**Informed Consent Statement:** Informed consent was not obtained for this study as the data were used retrospectively and obtained in a previous study.

**Data Availability Statement:** Data are available upon request to the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

*Bland–Altman Plots*

**Table A1.** The Bland–Altman plots for the mean and standard deviation (Bland–Altman plots show 1.96 × SD).

**Figure A1.** The Bland–Altman plots for comparing the hemodynamic parameters computed in this research, with and without self-gating: velocity; WSS; WSSA; WSSC; vorticity; viscous dissipation; energy loss; and kinetic energy. The average axis of the Bland–Altman plots corresponds to the average of the data with self-gating (SG) and without self-gating (NG), *Average* = *SG*+*NG* 2 . The difference axis corresponds to the difference between the data with self-gating (SG) and without self-gating (NG), *Di f f erence* = *SG* − *NG*. The difference axis shows the lines plotting the mean ± 1.96 times the standard deviation, and the zero value for comparison purposes.

#### **Appendix B**

*Local Differences*

**Figure A2.** (**A**) The local differences of the respiratory self-gated and non-gated segmentations in a representative volunteer. The differences were calculated by measuring the Euclidean distance between each triangle of the self-gated surface mesh to the nearest triangle of the non-gated mesh along the normal direction of the first one. The region where the aorta was segmented from the left ventricle showed more significant differences between the segmentations, with and without self-gating. (**B**) A histogram of the local differences calculated at each section for all volunteers. The bar plots show the fraction of differences calculated in each section and classified by its distance in millimeters. The local differences obtained between ± 1.0 mm were 87.61% in section I, 88.27% in section II, 88.66% in section III, and 89.11% in section IV.

**Table A2.** The mean values of the local differences between the respiratory self-gated and non-gated segmentations, computed at each section using the signed value differences.

#### **References**

## *Article* **Comprehensive Assessment of Left Intraventricular Hemodynamics Using a Finite Element Method: An Application to Dilated Cardiomyopathy Patients**

**Pamela Franco 1,2,3, Julio Sotelo 1,3,4, Cristian Montalba 1,3,5, Bram Ruijsink 6,7, Eric Kerfoot 6, David Nordsletten 8,9, Joaquín Mura 3,10, Daniel Hurtado 3,11,12 and Sergio Uribe 1,3,5,\***

**Abstract:** In this paper, we applied a method for quantifying several left intraventricular hemodynamic parameters from 4D Flow data and its application in a proof-of-concept study in dilated cardiomyopathy (DCM) patients. In total, 12 healthy volunteers and 13 DCM patients under treatment underwent short-axis cine b-SSFP and 4D Flow MRI. Following 3D segmentation of the left ventricular (LV) cavity and registration of both sequences, several hemodynamic parameters were calculated at peak systole, e-wave, and end-diastole using a finite element approach. Sensitivity, interand intra-observer reproducibility of hemodynamic parameters were evaluated by analyzing LV segmentation. A local analysis was performed by dividing the LV cavity into 16 regions. We found significant differences between volunteers and patients in velocity, vorticity, viscous dissipation, energy loss, and kinetic energy at peak systole and e-wave. Furthermore, although five patients showed a recovered ejection fraction after treatment, their hemodynamic parameters remained low. We obtained several hemodynamic parameters with high inter- and intra-observer reproducibility. The sensitivity study revealed that hemodynamic parameters showed a higher accuracy when the segmentation underestimates the LV volumes. Our approach was able to identify abnormal flow patterns in DCM patients compared to volunteers and can be applied to any other cardiovascular diseases.

**Keywords:** 4D flow MRI; flow quantification; finite elements; left ventricle; dilated cardiomyopathy

#### **1. Introduction**

Dilated cardiomyopathy (DCM) is more common than non-ischemic cardiomyopathy and leads to left ventricular dilation and systolic and diastolic dysfunction [1,2]. The

**Citation:** Franco, P.; Sotelo, J.; Montalba, C.; Ruijsink, B.; Kerfoot, E.; Nordsletten, D.; Mura, J.; Hurtado, D.; Uribe, S. Comprehensive Assessment of Left Intraventricular Hemodynamics Using a Finite Element Method: An Application to Dilated Cardiomyopathy Patients. *Appl. Sci.* **2021**, *11*, 11165. https:// doi.org/10.3390/app112311165

Academic Editor: Julio Garcia Flores

Received: 6 October 2021 Accepted: 8 November 2021 Published: 25 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

process that alters the heart's size, geometry, and function is associated with increased hemodynamic demands, which cause abnormal mechanical stress in the muscle [3,4]. The progression is associated with an incremented risk of heart failure and sudden cardiac death [5]. However, poor survival and high mortality rate reveal that effective treatment of DCM-related heart failure remains challenging. Pharmacological and resynchronization therapies have improved DCM treatment by halting disease progression and leading to reverse remodeling [1].

The preferred imaging technique for assessing the heat in DCM patients is cardiovascular magnetic resonance (CMR). CMR allows the acquisition of anatomical, cine, and velocity images, including 4D Flow MR [5–8].

4D Flow allows a qualitative and quantitative analysis of several hemodynamic parameters. It has been applied extensively in the great vessels, particularly in the aorta [9–11] and in the left ventricle (LV) for assessing intraventricular flow in some cardiovascular diseases [12–19]. Previous studies have demonstrated that lower kinetic energy values in diastole are associated with the deterioration of ventricular filling, induced by morphological alteration commonly found in Fontan patients, mitral regurgitation, and LV dysfunction or remodeling [13–15]. Additionally, turbulent kinetic energy has shown a stronger association with the ventricle's remodeling in patients with Tetralogy of Fallot and higher values in DCM patients compared with normal subjects [16]. Vortex formation has been studied qualitatively (vortex size and location) and quantitatively (Lagrangian Coherent Structures and the curl of velocity) [17–20]. These studies suggest that parameters associated with the 3D intraventricular flow may be critical for LV filling and ejection and could be relevant to the development of dilation, dysfunction, and prognosis in patients with heart diseases. While these measures have a potential role in describing intraventricular flow, the difficulties of implementing them have led to the analysis of only a few combinations of these parameters in a single cohort of patients [13–15,17,19].

Due to the multidirectional velocity data, although impressively comprehensive, it may need to be supplemented by more selective flow imaging at high temporal and spatial resolutions or computational fluid dynamics simulation. Reaching conclusions regarding small-scale methodology, which comprehensively describes the characteristics of intraventricular flow, could improve the use of intraventricular 4D Flow for clinical research and potential translation to clinical settings.

In this work, we adapted a method for quantifying 4D Flow in the aorta [9–11,20]. We modified the methodology applied in the left ventricle to obtain several hemodynamic parameters from a single segmentation from a 4D Flow dataset and cine MRI. To show the applicability of this approach, we performed a proof-of-concept study in which we applied the method in a small cohort of DCM patients to find which parameters were different from volunteers. We obtained three-dimensional hemodynamic parameters, including kinetic energy, vorticity, helicity density, viscous dissipation, and energy loss [9,13,21–25].

#### **2. Materials and Methods**

#### *2.1. Population*

A total of 12 healthy volunteers (HV), mean age 40.8 years (range 27–55 years), and 13 DCM patients, mean age 48.7 years (range 29–62 years), matched according to age and gender, were included in this research. Demographical and clinical data are described in Table 1. At the time of diagnosis, DCM was defined as the presence of symptoms and signs of heart failure with echocardiographic signs of ventricular enlargement and systolic myocardial dysfunction in the absence of hypertension, valve diseases, or significant coronary artery diseases sufficient to cause global systolic impairment, by the definition of the European Society of Cardiology [26]. Our DCM cohort all received treatment with an improved LV ejection fraction (range 51–66%) and LV volume indices at CMR imaging. All patients received standard guideline-directed treatment for DCM following the 2008 heart failure guidelines from the European society of cardiology. The details of treatments were not available, as our center is the referral center for several clinics for cardiac CMR. The

HV had normal electrocardiograms and echocardiographic examinations without valvular or ventricular dysfunction. All subjects participated under informed consent, with data collection approved by the Regional Ethics Committee, South East London, UK (REC, 12/LO/1456).

**Table 1.** Demographical and clinical data for healthy and DCM patients. All quantitative data are expressed as the median (range). HR: Heart Rate, EF: Ejection Fraction, LVSV: Left Ventricle Stroke Volume, CO: Cardiac Output, LVEDV: Left Ventricle End-Diastolic Volume, and Left Ventricle End-Systolic Volume. \* indicates statistically significant differences (*p* < 0.05).

#### *2.2. Data Acquisition*

Multi-slice 2D cine balanced steady-state free precession (b-SSFP) and 4D Flow MRI data were acquired in all subjects using a clinical 1.5 T MT Scanner (Philips Achieva, Philips Medical Systems, Best, The Netherlands). During the MRI examination, multi-slice b-SSFP was used to acquire short-axis morphological images in 40 frames with 8 mm slice thickness, using retrospective cardiac gating. Acquisition parameters were echo time (TE) 1.4 ms, repetition time (TR) of 2.8 ms, flip angle (FA) of 60◦ and acquired and reconstructed pixel sizes were 2.47 × 2.53 mm3 and 1.45 mm2, respectively. 4D Flow MRI data were acquired during free-breathing with MR parameters, as follows: TW of 2.3 ms, TR of 4.7 ms, FA of 6◦, velocity encoding of 130 cm/s, and spatial resolution (acquired and reconstructed) 2.5 × 2.5 × 2.5 mm3. These settings gave a temporal resolution of 58 ms. After the acquisition, the 4D Flow MRI data were reconstructed into 24 cardiac phases on the MRI system.

#### *2.3. Data Analysis*

The 4D Flow MRI datasets were processed using an-house MATLAB library (The MathWorks Inc., Natick, MA, USA), which included the registration of the b-SSFP cine and 4D Flow MR images, interpolation of the b-SSFP images, segmentation of the LV, and generation of the finite element mesh (Figure 1).

The Eidolon software was used to perform the registration between the multi-slice b-SSFP and the 4D Flow MRI (King's College London, London, UK) [27]. To obtain a smooth tetrahedral mesh, we doubled the number of slices in the b-SSFP images by using a cubic interpolation of values at neighboring grid points in each respective dimension, obtaining a final voxel size of 1.43 × 1.43 × 4.04 mm3. LV endocardium was automatically segmented throughout all cardiac phases in the short-axis cine b-SSFP images, using the image analysis software Segment v2.2R6410 (Medviso AB, Lund, Sweden) [28–30]. The segmentation was visually inspected and manually corrected if needed. Segmentations c were then used to generate a binary mask. Afterward, we created a tetrahedral mesh using the iso2mesh MATLAB Toolbox [31]. Once the mesh was constructed, we computed the velocity vector at each mesh node from the 4D Flow datasets using a cubic interpolation. 3D maps of vorticity, helicity density, viscous dissipation, energy loss, and kinetic energy fields were then calculated using a previously published finite element approach [9–11,25]. The description of the equations used to calculate each hemodynamic parameter is presented

in Table S1. The parameters were averaged at peak systole, e-wave, and end-diastole using one timeframe before and after to reduce noise in the data.

**Figure 1.** Schematic description of the quantification process. (**a**) First, we performed registration of the 4D Flow with the b-SSFP images. (**b**) Second, we doubled the number of slices in the b-SSFP images, (**c**) then the LV segmentation and tetrahedral mesh were generated. (**d**) We estimated the cardiac phases under study (**e**) and we transferred the velocity information at each node of the mesh from the 4D Flow MRI datasets using cubic interpolation. (**f**) Then, we calculated hemodynamic parameters under study. (**g**) Finally, the mean values of hemodynamic parameters were included in a bullseye plot to compare volunteers and DCM patients.

#### *2.4. Local Hemodynamics*

A 16-segment model was used to divide the LV. In contrast to a standardized nomenclature, a minor adjustment was made [26]. Due to the generally intricate shape of the apical region of the LV, region 17 was excluded from our analysis. Accordingly, LV mesh was divided into 16 segmentations. The centerline of the LV was calculated automatically by detecting the centroid of the LV contour in each slice and connected to create a line. To determine the three sections of the LV, we divided the centerline into three equal parts perpendicular to the long axis of the heart. An additional point was then manually placed at the junction between the right ventricular free wall and the interventricular septum on the LV. Based on these positions, landmarks were uniformly distributed along the boundaries. Each section was then partitioned into six segments of 60◦ each on basal and mid-cavity sections and four segments of 90◦ each on apical section (Figure 1f). Finally, for visualization purposes, we used the scientific software ParaView version 5.3.0 (Kitware, Clifton Park, NY, USA).

#### *2.5. Statistical Analysses*

Normal distribution in population demographics was evaluated using the Shapiro-Wilk test. Differences between groups for continuous parameters were assessed by a Student *t*-test if they presented a normal distribution, and the Mann-Whitney *U* test otherwise. The χ2 test was applied for categorical variables, which were reported as percentages. A *p*-value < 0.05 was considered statistically significant. The statistical

analyses were performed using GraphPad Prism version 6.0.1 (GraphPad Software Inc., San Diego, CA, USA).

These data were displayed in box-whisker and bullseye plots for global and local analyses, respectively. Additionality, a correlation matrix-based hierarchical clustering method was introduced to extract multiple correlation patterns from hemodynamic parameters. This method can effectively identify highly correlated data. The results are described with a tree structure plot called a dendrogram. The present study used Pearson's correlation method to measure the similarity between hemodynamic parameters [32].

Furthermore, a sensitivity study was performed by looking at changes in the hemodynamic parameters subjected to the LV segmentation changes. We increased and decreased the size of the LV cavity from the first segmentations by moving the segmentation contour in 0.5 to 2 pixels of the b-SSFP image, equivalent to 0.72 to 2.89 mm. We compared the results with the original LV segmentation's respective mean value at each cardiac phase studied. We used the Kruskal-Wallis test to compare the variables across the different LV segmentation, with a *p*-value < 0.05 indicating statistical significance. The significance level was adjusted by using Dunn's test correction.

To assess the inter-observer agreement, data were analyzed by two independent observers, one with three years of experience in MR LV quantification and the other a medical technologist with no previous experience in this field. In addition, re-analyzed images with a 1-month interval to evaluate the intra-observer reproducibility. Inter- and intra-observer reproducibility were analyzed using Bland-Altman plots, and the results are shown in the Supplementary Materials.

#### **3. Results**

#### *3.1. Study Population*

There were no significant differences in age and heart rate (Table 1). However, ejection fraction and stroke volume were lower in DCM patients than volunteers, while enddiastolic and end-systolic volumes were larger. These changes indicated that the LV in DCM patients was enlarged and its cardiac function was reduced, which is consistent with the pathological characteristics of DCM [1,2]. Additionally, eight patients still showed significantly impaired systolic function at CMR's time (non-responders), and five patients showed a complete response to treatment (complete responders). Between DCM groups, complete- vs. non-responders, there were differences in ejection fraction, end-diastolic volume, and end-systolic volume.

#### *3.2. Global Hemodynamics*

Assessment of global hemodynamic parameters is shown in Figure 2 and Table 2. Volunteers showed higher hemodynamics values than patients at peak systole and e-wave, except for helicity density. Remarkably, hemodynamic parameters in complete responder DCM patients remained low compared to volunteers. We found statistical differences between HV and DCM patients: non- and complete responders at peak systole and ewave in velocity, vorticity, viscous dissipation, energy loss, and kinetic energy. In all cases, *p*-values were lower or equal to 0.005. There were no statistical differences in the parameters at end-diastole. Furthermore, we did not find statistical differences between DCM groups. In addition, ROC curves showed that previous parameters discriminated between HV and DCM patients (Figure S1).

The total computational time used to process the data, once the multi-slice b-SSFP was segmented and registered, varied between 30–40 s for one cardiac phase, using a standard computer (3.4 GHz Intel ® Core i7TM, 16 GB RAM).

Hierarchical cluster analysis (Figure 3) provides an alternative method for reliable identification of correlation between ejection fraction and hemodynamics parameters from 4D-flow MRI. According to their similarities, they are classified into two clusters identified at peak-systole, e-wave, and end-diastole. At peak-systole and e-wave, cluster 1 (black): helicity density; and cluster 2 (red): ejection fraction, energy loss, vorticity, viscous dissipation, velocity, and kinetic energy. Finally, at end-diastole, cluster 1 (red): helicity density and ejection fraction, and cluster 2 (red): energy loss, viscous dissipation, vorticity, velocity, and kinetic energy. This means that ejection fraction correlates with all parameters except helicity density at peak systole and e-wave.

**Figure 2.** Box whisker plots for hemodynamic parameters (**a**–**f**) in the entire LV cavity of HV and DCM patients groups at peak systole, e-wave, and end-diastole. On each box, the central mark is the median, the bottom and top edges of the box are the 25th and 75th percentiles, respectively, and the whiskers extend to the most extreme data points not considered outliers. \*,+ Indicates statistically significant differences (*p* < 0.05).

#### *3.3. Sensitivity Study, Intra-, and Inter-Observer Reproducibility*

Figure 4 shows a sensitivity analysis at peak-systole. The relative error in LV cardiac volumes and helicity density did not show significant differences between groups across LV segmentation. There were significant differences for some segmentation for the other parameters, particularly for the velocity magnitude, energy loss, and kinetic energy. The hemodynamic parameters showed a relative error proportional to the dilatation or erosion of the contour in the segmentation. When the segmentation was dilated or eroded 1 pixel or less, the relative error differences, with respect to the original segmentation for volunteers and DCM patients, were: velocity magnitude (9.03%, 6.78%), vorticity magnitude (5.49%, 2.74%), helicity density (12.11%, 13.98%), viscous dissipation (6.31%, 3.34%), energy loss (3.59%, 5.89%), and kinetic energy (7.66%, 6.09%). Similar results were obtained at e-wave and end-diastole, as shown in Figures S2 and S3. Those errors were more significant, particularly when the segmentation was dilated or eroded by more than 1 pixel. Helicity density and energy loss showed greater dependency on the segmentation error.

Regarding reproducibility, there was an excellent agreement of inter- and intraobserver analysis of global hemodynamic parameters. Details are given in the Appendix A.

**Table 2.** Global hemodynamics data for HV, complete, and non-responders DCM patients. All quantitative data are expressed as the mean ± standard deviation. \*,+ Indicates statistically significant differences (*p* < 0.05).

#### *3.4. Local Hemodynamics*

Figure 5 shows the bullseye plots of the hemodynamic parameters for volunteers and DCM patients at peak systole. More areas with statistical differences were observed mainly in velocity magnitude and kinetic energy, particularly in anteroseptal, inferior, inferolateral basal, anterior, inferoseptal, inferior, inferolateral mid-cavity, anterior, septal, and lateral apical segments (all *p*-values < 0.033). Additionally, vorticity magnitude showed statistical differences in anteroseptal basal (*p* = 0.045) and inferoseptal mid-cavity (*p* = 0.046) segments. Energy loss showed statistical differences in inferoseptal (*p* = 0.029) and inferolateral (*p* = 0.046) mid-cavity segments and in anterior (*p* = 0.023), septal (*p* = 0.070), and inferior (*p* = 0.077) apical segments. Helicity density and viscous dissipation did not show statistical differences in any parcellation.

Figure S4 shows the comparison at the e-wave. Similar to peak systole, statistical differences were in velocity magnitude and kinetic energy. Statistical differences in viscous dissipation and energy loss were found in inferolateral basal and septal, anterior, and lateral apical segments (all *p*-values < 0.049). Vorticity magnitude showed statistical differences in anteroseptal basal (*p* = 0.049) and septal (*p* = 0.039) and lateral (*p* = 0.039) apical segments. Helicity density did not show statistical differences in any parcellation of the LV.

We did not find statistical differences in any segment at end-diastole.

The mean values of local hemodynamic parameters for both groups under study at peak systole, e-wave, and end-diastole are available in Tables S2–S4, respectively.

**Figure 3.** Dendrogram and hierarchical clustering results based on average linkage method for ejection fraction and hemodynamic parameters. EF: ejection fraction, V: velocity, KE: kinetic energy, Vo: vorticity, HD: helicity density, VD: viscous dissipation, and EL: energy loss.

**Figure 4.** Relative error values of the volume (**a**) and each hemodynamic parameter (**b**–**g**) obtained comparing the reference segmentation with segmentations given by erosion or dilatation for each group of volunteers and patients at peak systole. \* indicates statistically significant differences (*p* < 0.05).

**Figure 5.** Bullseye plots of mean hemodynamic parameters (**a**–**f**) across 16 segments for volunteers (i) and patients (ii) at the peak systole. \* indicates statistically significant differences (*p* < 0.05).

#### **4. Discussion**

We developed a method to characterize the left intraventricular hemodynamics in the LV from 4D Flow MRI using a finite element method, applied in a cohort of DCM patients. This approach estimates vorticity, helicity density, viscous dissipation, energy loss, and kinetic energy fields from a single segmentation. The hemodynamics results indicated that velocity magnitude, vorticity magnitude, viscous dissipation, energy loss, and kinetic energy revealed statistical differences between volunteers and patients, particularly at peak systole and e-wave.

Some of the parameters reported in this study have been reported before. Nevertheless, those parameters have been obtained from different methodologies in different cohorts of patients. In our case, we calculated several parameters from a single segmentation at once from only one 4D flow dataset, which is difficult to determine with other methods.

Some other methods are based on a finite difference approach, as in Lorenz et al. [33]. However, it is well known that finite difference cannot effectively handle complex geometries, such as those found in the cardiovascular system. Neither can impose boundary conditions on irregular surfaces directly but they are both sensitive to noise. Fouras et al. showed that this approach suffers from a loss of accuracy in estimating hemodynamic parameters due to the omission of out-of-plane velocity information [34]. On the other hand, Sotelo et al. demonstrated the convergence and robustness of the finite element method in cardiovascular flow [9,10]. Further, they also showed that the finite element method is both stable and accurate in the presence of noise.

Although DCM mainly affects the systolic function, we evaluated the hemodynamic parameters at systole and diastole, as several papers have shown that diastolic function is also affected by this disease. For instance, Friedberg et al. and Dragulescu et al. reported that diastolic wall-motion abnormalities are prevalent in pediatric DCM. Their presence is associated with diastolic ventricular dysfunction and adverse outcomes [35,36]. Some papers have assessed diastolic function in DCM patients using 4D Flow data [37,38]. They have described alterations in the flow components related to velocity, vorticity, and kinetic energy in a different cohort of patients, consistent with our results [21–25].

As we showed in this study, intraventricular flow in DCM patients was altered compared to healthy volunteers at diastole. In addition, it is interesting to observe that, while end-diastolic volume was significantly larger in patients than healthy subjects, the maximum hemodynamic values for e-wave and end-diastole were smaller in patients than in volunteers. For instance, we found that in the normal LV, kinetic energy values were high. The highest kinetic energy values were observed during early diastole and regionally distributed near basal LV regions. In contrast, early- and end-diastole kinetic energy was lower than normal in a heterogeneous group of DCM patients and decreased with the LV volume. As we found in our results, this decrease of kinetic energy throughout diastole is associated with viscous dissipation and energy loss. That agrees with a previous study, where comparisons of inflow characteristics in healthy subjects and DCM patients showed more differences at e-wave between the two groups [36–41]. These ventricular diastolic function aspects can be influenced by dynamic load and contractility; these may vary within the spectrum of normal conditions. Furthermore, large ventricles lead to weak suction pumps and have high wall tension, which has been previously suggested to cause energy waste and alter vortex ring dynamics [22,23]. These results indicate that, despite the complex nature of ventricular finding factors, clinically useful information regarding left ventricular diastolic function is associated with distinct mitral flow velocity patterns. Therefore, these alterations of blood flow may be a factor in developing systolic and diastolic dysfunctions.

While eight patients showed significantly impaired systolic function at CMR's time, five showed a complete response to treatment. Despite a nearly normalized LV ejection fraction, it showed similar hemodynamic values to non-responders DCM patients that were markedly different from volunteers. These results suggest a significant increase in the ratio of outflow to inflow during systole in responders DCM patients, but the volumes

were significantly smaller. Therefore, LV ejection fraction cannot reflect subtle ventricular dysfunction, which potentially can be better assessed using flow-based parameters because of the sensitivity to abnormal pumping function [36]. Therefore, the problem of using LV ejection fraction as the pivotal risk marker for DCM patients is that this single parameter does not recapitulate the complexity of the disease.

The location and extent of the changes in intraventricular blood flow, for example, the depth (base to apex) of the vorticity changes or the spread of impaired flow through the ventricular cavity, can potentially be a sensitive marker of the severity of diseases or the progress of the treatment but are hard to quantify because of the 3D nature of the flow. We proposed to use the bullseye plots to depict this data. These plots allow us to display the most important regional differences and extend flow changes in a familiar way to many clinicians. These results could facilitate homogeneity among 4D Flow quantifiable analysis for clinical researchers and clinicians.

The sensitivity study showed a significant relative error, particularly in helicity density, when the differences in the segmentations were greater than 1 pixel in dilation and erosion cases. We performed this sensitivity study even under the pixel resolution of the 4D Flow data. In each pixel, there were four or five elements from the mesh, whose flow values were interpolated from neighborhood pixels. When the segmentation error was lower or equal to 1 pixel, the maximum mean relative error was less than 10% in most hemodynamic parameters studied. Previous research has shown that DCM patients have a lower mean value of velocity magnitude than healthy volunteers [21–25]. These results were considered when the segmentation contour fell inside the LV blood pool. In general, we also observed that lower errors were obtained for almost all hemodynamic parameters when the segmentation underestimated the LV volumes.

Inter-observer and intra-observer assessments showed excellent reproducibility of the results with negligible mean differences and small limits of an agreement at peak systole, e-wave, and end-diastole for all the parameters assessed. It is important to note that the high intra- e inter-observer variability was obtained because we performed an automatic segmentation process using the software Segment. Therefore, the difference in segmentations was minimal, as previously reported by Tufvesson et al. [29]. The automated process corrections were also minimal, which led to a high intra- and interobserver variability. On the other hand, the sensitivity study was performed by modifying a reference segmentation in the entire contour by applying erosion or dilatation. This result implies a more significant volume difference concerning the reference, and, as a result, high sensitivity to the segmentation was obtained.

The limited size of this proof-of-concept study did not allow us to investigate the prognostic impact, but this will be our aim in future research. Nevertheless, in this small cohort of patients, we have shown that the velocity, vorticity, kinetic energy, viscous dissipation, and energy loss revealed statistical differences between volunteers and patients. This finding could be relevant to assess changes in a longitudinal study or to study the response to a particular therapy. Additionally, 4D Flow derived parameters showed that, in responding DCM patients, hemodynamics parameters were low, even though they had a recovered ejection fraction. Nevertheless, hierarchical cluster analysis underlined that a moderate correlation may exist between ejection fraction and 4D Flow-based metrics, which needs to be studied further. A clinical study involving more DCM patients should be performed in order to corroborate a prognostic impact—and hence a clinical relevance—of 4D Flow analysis in monitoring DCM patients.

The segmentation of the data need was performed over the multi-slice b-SSFP. Additionally, multi-slice b-SSFP and 4D Flow images need to be registered before analyzing the 4D Flow data. Ideally, the segmentation would be made directly on the 4D Flow data. However, the contrast between the blood pool and myocardium in our 4D Flow data was insufficient to perform accurate segmentation. New sequence developments will likely improve contrast in 4D Flow acquisitions, potentially allowing direct segmentation from the 4D Flow data.

#### **5. Conclusions**

This study describes a methodology for quantitative evaluation of intraventricular hemodynamics using a single segmentation from a 4D Flow dataset. We demonstrate that velocity, vorticity, viscous dissipation, energy loss, and kinetic energy can characterize changes in intraventricular flow in DCM patients compared to healthy volunteers. Further studies should focus on the impact of different treatments of DCM patients on those parameters. Our evidence shows that, although ejection fraction may be recovered, the hemodynamic parameters remain low.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/app112311165/s1. Table S1: Equations used to calculate each hemodynamic parameter; Table S2: Mean parameter values across 16 segments of the LV, during peak systole, where the bold type is statistically significant between volunteers and patients (*p* < 0.05). υ: Velocity magnitude, ω: vorticity magnitude, Hd: helicity density, VD: viscous dissipation, EL: energy loss, and K: kinetic energy; Table S3: Mean parameter values across 16 segments of the LV during e-wave, where the bold type is statistically significant between volunteers and patients (*p* < 0.05). υ: Velocity magnitude, ω: vorticity magnitude, Hd: helicity density, VD: viscous dissipation, EL: energy loss, and K: kinetic energy; Table S4: Mean parameter values across 16 segments of the LV, during end-diastole, where the bold type is statistically significant between volunteers and patients (*p* < 0.05). υ: Velocity magnitude, ω: vorticity magnitude, Hd: helicity density, VD: viscous dissipation, EL: energy loss, and K: kinetic energy. Figure S1: ROC-curves for hemodynamic parameters (a–f) in the entire LV cavity of the groups of volunteers and patients at peak-systole, e-wave, and end-diastole; Figure S2: Relative error values of volume (a) and each hemodynamic parameter (b–g), obtained comparing the reference segmentation with segmentations given by erosion or dilation, for each group of volunteers and patients at e-wave. \* indicates statistically significant differences (*p* < 0.05); Figure S3: Relative error values of volume (a) and each hemodynamic parameter (b–g), obtained comparing the reference segmentation with segmentations given by erosion or dilation, for each group of volunteers and patients at end-diastole. \* indicates statistically significant differences (*p* < 0.05); Figure S4: Bullseye plots of mean hemodynamic parameters (a–f) across 16 segments for volunteers (i) and patients (ii) at the e-wave. \* indicates statistically significant differences (*p* < 0.05).

**Author Contributions:** All authors were actively involved in reviewing and drafting the manuscript. All authors have approved the final version of this manuscript. P.F. developed the methodology to assess intracardiac flow, process all data, and perform statistical analysis and wrote the manuscript. J.S. developed the finite element method to evaluate the hemodynamic parameters and edited the manuscript. C.M. was the second observer in the inter-observer analysis and edited the manuscript. B.R. participated in the recruitment of patients, acquired CMR images, and edited the manuscript. E.K. developed the Eidolon software to register b-SSFP data and 4D Flow for the cohort of DCM patients and edited the manuscript. D.N. participated in the design of the study to recruit patients and edited the manuscript. J.M. participated helped develop the finite element method to assess hemodynamic parameters and edited the manuscript. D.H. participated in the development of the methodology to assess intracardiac hemodynamics and edited the manuscript. S.U. participated in the study's design, developing the method to assess intracardiac flow, helping with statistical analysis, and wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work has been funded by projects PIA-ACT 192064 and the Millennium Nucleus on Cardiovascular Magnetic Resonance NCN17\_129 of the Millennium Science Initiative, both of the National Agency for Research and Development, ANID. The authors also thanks to Fondecyt project 1181057 also by ANID. Franco P. thanks to ANID-PCHA/Doctorado-Nacional/2018-21180391. Sotelo J. thanks to CONICYT-FONDECYT Postdoctorado 2017 #3170737 and ANID-FONDECYT de Iniciación en Investigación #11200481.

**Institutional Review Board Statement:** All subjects participated under informed consent, with data collection approved by the Regional Ethics Committee, South East London, UK (REC, 12/LO/1456).

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The datasets generated during and/or analyzed during the current study are not publicly available due to data privacy according to the rules of King's College London but could be available from Bram Ruijsink<EMAIL_ADDRESS>on reasonable request.

**Acknowledgments:** This work has been funded by projects PIA-ACT 192064 and the Millennium Nucleus on Cardiovascular Magnetic Resonance NCN17\_129 of the Millennium Science Initiative, both of the National Agency for Research and Development, ANID. The authors also give thanks to Fondecyt project 1181057, also by ANID. Franco P. gives thanks to ANID-PCHA/Doctorado-Nacional/2018- 21180391. Sotelo J. gives thanks to CONICYT-FONDECYT Postdoctorado 2017 #3170737 and ANID-FONDECYT de Iniciación en Investigación #11200481.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Intra- and Inter-Observer Reproducibility**

As shown in Figure A1, excellent intra-observer agreement with minimal mean differences and small limits of agreement were found for peak systole. Mean differences were: velocity magnitude −0.0003 ± 0.0118 m/s, kinetic energy (−0.4580 ± 0.7490) × 10−9 J, vorticity magnitude 0.0029 ± 0.0249 1/s, helicity density (−0.0978 ± 0.5345) × 10−3 m/s2, viscous dissipation −0.0070 ± 0.0578 1/s2, and energy loss (−0.0124 ± 0.2563) × 10−9 W. Similar results were obtained at e-wave and end-diastole, as shown in Figures A2 and A3, respectively. Figure A4 demonstrates excellent inter-observer analysis agreement for peak systole. Mean differences were: velocity magnitude −0.0024 ± 0.0124 m/s, kinetic energy (−0.0349 ± 0.1093) × 10−5 J, vorticity magnitude 0.1122 ± 0.7634 1/s, helicity density −0.0001 ± 0.0139 m/s2, viscous dissipation −2.1652 ± 11.8205 1/s2, and energy loss (0.0012 ± 0.3008) × 10−8 W. Figures A5 and A6 show the results obtained at e-wave and end-diastole, respectively, with comparable results at peak systole.

**Figure A1.** Bland-Altman plots represent the intra-observer reproducibility in the measurements of LV global hemodynamic parameters (**a**–**f**) at peak systole. The thick line represents the mean difference, and the thin lines represent the limits agreement (1.96 SD).

**Figure A2.** Bland-Altman plots represent the intra-observer reproducibility in the measurements of LV global hemodynamic parameters (**a**–**f**) at e-wave. The thick line represents the mean difference, and the thin lines represent the limits agreement (1.96 SD).

**Figure A3.** Bland-Altman plots represent the intra-observer reproducibility in the measurements of LV global hemodynamic parameters (**a**–**f**) at end-diastole. The thick line represents the mean difference, and the thin lines represent the limits agreement (1.96 SD).

**Figure A4.** Bland-Altman plots represent the inter-observer reproducibility in the exams of LV global hemodynamic parameters (**a**–**f**) at peak systole. The thick line represents the mean difference, and the thin lines represent the limits agreement (1.96 SD).

**Figure A5.** Bland-Altman plots represent the inter-observer reproducibility in the exams of LV global hemodynamic parameters (**a**–**f**) at e-wave. The thick line represents the mean difference, and the thin lines represent the limits agreement (1.96 SD).

**Figure A6.** Bland-Altman plots represent the inter-observer reproducibility in the exams of LV global hemodynamic parameters (**a**–**f**) at end-diastole. The thick line represents the mean difference, and the thin lines represent the limits agreement (1.96 SD).

#### **References**

## *Article* **4D Flow MRI in Ascending Aortic Aneurysms: Reproducibility of Hemodynamic Parameters**

**Joe F. Juffermans 1,\*, Hans C. van Assen 1, Bastiaan J. C. te Kiefte 1, Mitch J. F. G. Ramaekers 1,2,3,4, Roel L. F. van der Palen 5, Pieter van den Boogaard 1, Bouke P. Adriaans 2,3,4, Joachim E. Wildberger 3,4, Ilona A. Dekkers 1, Arthur J. H. A. Scholte 6, Simon Schalla 2,3,4, Hildo J. Lamb 1 and Jos J. M. Westenberg 1**

**Abstract:** (1) Background: Aorta hemodynamics have been associated with aortic remodeling, but the reproducibility of its assessment has been evaluated marginally in patients with thoracic aortic aneurysm (TAA). The current study evaluated intra- and interobserver reproducibility of 4D flow MRI-derived hemodynamic parameters (normalized flow displacement, flow jet angle, wall shear stress (WSS) magnitude, axial WSS, circumferential WSS, WSS angle, vorticity, helicity, and local normalized helicity (LNH)) in TAA patients; (2) Methods: The thoracic aorta of 20 patients was semi-automatically segmented on 4D flow MRI data in 5 systolic phases by 3 different observers. Each time-dependent segmentation was manually improved and partitioned into six anatomical segments. The hemodynamic parameters were quantified per phase and segment. The coefficient of variation (COV) and intraclass correlation coefficient (ICC) were calculated; (3) Results: A total of 2400 lumen segments were analyzed. The mean aneurysm diameter was 50.8 ± 2.7 mm. The intra- and interobserver analysis demonstrated a good reproducibility (COV = 16–30% and ICC = 0.84–0.94) for normalized flow displacement and jet angle, a very good-to-excellent reproducibility (COV = 3–26% and ICC = 0.87–1.00) for all WSS components, helicity and LNH, and an excellent reproducibility (COV = 3–10% and ICC = 0.96–1.00) for vorticity; (4) Conclusion: 4D flow MRI-derived hemodynamic parameters are reproducible within the thoracic aorta in TAA patients.

**Keywords:** aorta; aneurysm; hemodynamic; 4D flow MRI; reproducibility

#### **1. Introduction**

Patients with a thoracic aortic aneurysm (TAA) have an increased risk for aortic rupture and dissection [1]. Clinical guidelines use a maximal cross-sectional diameter ≥55 mm as the main criterion for the recommendation of preemptive surgical aortic replacement [2,3]. However, the majority (60–96%) of dissections occur in aortas with diameters below this threshold [4,5]. Therefore, there is a need for additional and more sensitive markers to identify patients at high risk of progressive dilation or adverse aortic events [6].

Evaluation of aortic hemodynamics is considered to be highly promising for the prediction of progressive dilatation and adverse aortic events [7–14]; the blood flow over the

**Citation:** Juffermans, J.F.; van Assen, H.C.; te Kiefte, B.J.C.; Ramaekers, M.J.F.G.; van der Palen, R.L.F.; van den Boogaard, P.; Adriaans, B.P.; Wildberger, J.E.; Dekkers, I.A.; Scholte, A.J.H.A.; et al. 4D Flow MRI in Ascending Aortic Aneurysms: Reproducibility of Hemodynamic Parameters. *Appl. Sci.* **2022**, *12*, 3912. https://doi.org/10.3390/app12083912

Academic Editor: Qi-Huang Zheng

Received: 18 March 2022 Accepted: 11 April 2022 Published: 13 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

cardiac cycle through the aorta can be analyzed using three-dimensional time-resolved phase-contrast magnetic resonance imaging (MRI), also known as four-dimensional (4D) flow MRI. From 4D flow MRI data, several relevant patient-specific hemodynamic parameters can be assessed. It has been demonstrated for patients with a tricuspid aortic valve (TAV) that the amount of flow displacement (a measure of eccentricity in the flow profile) is associated with aortic growth [7,8]. Furthermore, TAV patients showed decreased levels of vortical flow, helical flow, and wall shear stress (WSS) compared to healthy volunteers [9,10].

While such hemodynamic parameters could also be important for future clinical risk stratification of aneurysm patients, the majority of the studies assessing their agreements between observers and examinations only included healthy volunteers, patients with bicuspid aortic valves, or patients with mildly dilated aortas [13–22]. Hence, the reported outcomes are not fully translatable to TAA patients with substantially dilated aortas and a TAV. Therefore, the reproducibility of the hemodynamic parameters in such patients remains unknown [1]. Moreover, several previous studies quantified parameters based on (multiple) two-dimensional plane quantifications or ranked them qualitatively [15–18,21]. Thus, part of the spatial information is lost. In addition, the reproducibility of the WSS magnitude (WSSmag) has been assessed well, but little is known about the reproducibility of the axial or circumferential WSS component (WSSax, WSScir, respectively) and the angle between both WSS components (WSSangle).

To allow quantification of hemodynamic parameters in 4D flow MRI data, a cardiac phase-specific 3D lumen segmentation is required [23,24]. The shape of these lumen segmentations can be described by morphologic parameters such as the volume, centerline length, maximal diameter, and curvature radius. Interestingly, a recent study demonstrated that the curvature radius of the ascending aorta is associated with aortic growth [7].

In most commercial and research software tools, these lumen segmentations are constructed (semi-)automatically on the 4D flow MRI data and then manually adjusted when needed. The manual interaction may introduce observer-dependent variability in the quantitation of the hemodynamic parameters [25,26]. This observer-dependent variability of the 3D lumen segmentation can be assessed directly by evaluating the variability in morphologic parameters. The reproducibility of these morphologic parameters from 4D flow MRI data has been primarily reported in healthy volunteers [25].

Consequently, the aim of this study was to evaluate intra- and interobserver reproducibility of hemodynamic parameters, i.e., normalized flow displacement, flow jet angle, wall shear stress (magnitude, axial, circumferential, and angle), vorticity and helicity (absolute and local normalized) quantitatively assessed for six thoracic aortic lumen segments in patients with TAA and normally functioning TAV. In addition, intra- and interobserver reproducibility of morphological parameters, i.e., aortic volume, centerline length, lumen diameter, and curvature radius, were assessed.

#### **2. Materials and Methods**

#### *2.1. Study Population*

The study protocol was conducted in accordance with the Declaration of Helsinki and was approved by the local Medical Ethics Committee of the Leiden University Medical Center (G20.149, 9 October 2020), which waived the patient's informed consent for anonymized clinical data. The clinical database was used to identify TAA patients with a TAV who underwent a 4D flow MRI of the thoracic aorta between October 2018 and August 2021 (*n* = 260). Patients with aortic valve stenosis, regurgitation, or previous aortic surgery were excluded to obtain a selection of TAA patients without other pathologies, which potentially could also affect the reproducibility (*n* = 171). The clinical reports of the patients were evaluated for the maximal aortic diameter, and the 20 patients with the largest reported maximal diameters were included, see Figure 1. The presence of an aortic root or ascending aorta aneurysm was classified according to Della Corte's classification [27].

**Figure 1.** Population selection diagram. Abbreviations: TAA—thoracic aortic aneurysm, TAV tricuspid aortic valve, MRI—magnetic resonance imaging, and 4D—four dimensional.

#### *2.2. MRI Acquisition*

The MRI examinations consisted of a 4D flow scan covering the entire thoracic aorta (for details, see Supplementary Materials Table S1). All subjects were scanned with a 3Tscanner (Ingenia or Elition, Philips Healthcare, Best, The Netherlands) using a FlexCoverage anterior and dStream Torso posterior coil. Concomitant gradient correction was performed using standard available scanner software. All data were visually inspected, and the absence of aliasing in the 4D flow MRI data was confirmed.

#### *2.3. Image Analysis*

The image analysis consisted of three parts: 4D flow MRI lumen segmentation, aorta morphology quantification, and hemodynamic quantification. First, the aortic lumen of the patients was segmented twice by the first observer (J.J.) and once by the second and third observers (M.R. and B.K., respectively). Lumen segmentation was performed using CAAS MR Solutions v5.2.1 (Pie Medical Imaging, Maastricht, The Netherlands). The software generated initial segmentations on the peak systolic phase and two consecutive phases before and after this peak systolic phase (i.e., five systolic phases, see Figure 2 and Supplementary Materials Video S1). Next, the aortic segmentations were manually corrected for all five phases and partitioned into six consecutive segments by manually placing anatomical planes perpendicular to the aortic centerline, respectively: the aortic root (from the aortic valve to the sinotubular junction), proximal ascending aorta (from the sinotubular junction to the mid-ascending aorta, defined as the midpoint on the centerline between the sinotubular junction and brachiocephalic artery), distal ascending aorta (from the mid-ascending aorta to the brachiocephalic artery), aortic arch (from the brachiocephalic artery until past the left subclavian artery), the proximal descending aorta (from the left subclavian artery to the mid-descending thoracic aorta, defined as the midpoint on the centerline between the left subclavian artery and the descending aorta at the level of the aortic valve), and distal descending aorta (from the mid-descending thoracic aorta to the descending aorta at the level of the aortic valve). More details about the aortic segmentation are described by van der Palen et al. [20]. In addition, the stroke volume and cardiac output were assessed at the level of the sinotubular junction by the first observer for characterizing the left ventricular function of the population using CAAS MR Solutions.

**Figure 2.** 4D flow MRI visualization using CAAS MR Solutions v.5.2.1. Example of a patient (male, 45 years old) with aortic root aneurysm with a diameter of 52 mm. (**I**) 4D flow MRI reconstructions. (**II**) Image analysis including the lumen segmentation, the manual segmentation correction, anatomical partitioning visualized on the wall shear stress maps, and stroke volume assessment. (**III**) Streamline visualization of the 4D flow MRI at the 5 systolic phases (60, 90, 120, 150, and 180 ms). Abbrevations: AoR—aortic root, pAA—proximal ascending aorta, dAA—distal ascending aorta, AoA—aortic arch, pDA—proximal descending aorta, and dDA—distal descending aorta. Secondly, the morphology of the thoracic aortas was analyzed using an in-house developed Python software (Python Software Foundation, Wilmington, DE, USA). This tool was used to quantify aortic volume, centerline length, maximal diameter, and curvature radius of each anatomical segment and for each of the five systolic phases. The aortic diameter was determined by first constructing a cross-section perpendicular to the centerline at every millimeter. Next, at each cross-section, the maximal radial spike length was calculated. The curvature radius was derived by fitting a circle through the segments' centerline as previously described [25]. The maximal aortic diameter within each segment and the radius of the fitted circle were used for statistical analysis.

Thirdly, the hemodynamics of the thoracic aorta were quantified using methods described in previous studies [13,15,19]. These methods were implemented in an in-house developed Python software [9,25]. Initially, the mean blood velocity over all subjects was quantified per anatomical segment. The tool was then used to quantify the normalized flow displacement, flow jet angle, vorticity norm, absolute helicity, and absolute LNH for each anatomical segment and each systolic phase. The normalized flow displacement and jet angle were determined by first constructing a cross-section perpendicular to the centerline at every millimeter. Next, at each cross-section, the normalized flow displacement (distance between the centerline point and center of velocity normalized for the vessel diameter times 100%) and flow jet angle (angle between the centerline and mean velocity vector) were calculated as previously described [15]. The maximal normalized flow displacement and flow jet angle within each segment were assessed and used for statistical analysis. For the calculation of vorticity and helicity of the velocity vector fields (→*V* ), a spatial Gaussian differential operator (→∇ ) with a standard deviation equal to the reconstructed in-plane voxel size was used. The vorticity (→*ω* = →∇ × →*V* ) norm and absolute helicity (*H* = |→*V* · →*ω* |) were quantified via methods previously described by Ramaekers et al. [9]. The vorticity norm and absolute helicity were both normalized for the reconstructed voxel size, and the mean values within each segment were used for statistical analysis. Since a previous study demonstrated that a threshold of >0.6 for absolute LNH (*LNH* = →*V* · →*ω* |→*V* ||→*ω* | ) is

the optimal criterion for detecting differences between patients and healthy volunteers [13], this threshold was used in our study to select volumes with increased absolute LNH. The volumes with increased absolute LNH were assessed and used for the statistical analysis. The in-house developed software can be made available for scientific collaboration per reasonable request. To derive 3D WSS maps for each of the five systolic phases, CAAS MR Solutions v5.2 was used. At each point of the 3D WSS map, WSSmag was calculated as well as the WSSax and WSScir components. The WSSangle between the WSScir and WSSax was quantified for each point as previously described [19]. For each anatomical segment, the mean WSSmag, WSSax, WSScir, and WSSangle were used for statistical analysis.

#### *2.4. Observer Training*

All observers were acquainted with CAAS MR Solutions software. The first, second, and third observer had 4, 3, and 1 years' experience with lumen segmentation on 4D flow MRI data in TAA patients, respectively. Since the observers were recruited at two clinical centers and observers had different levels of experience, all observers received feedback on the segmented aorta from an experienced cardiovascular MRI researcher (J.J.M.W. with >15 years of 4D flow MRI experience) after analyzing the training cohort (*n* = 10 TAA patients). Hereafter, the results of the training cohort were not improved. Next, the observers analyzed the validation cohort (*n* = 10 TAA patients).

#### *2.5. Statistical Analysis*

Statistical analysis was performed using the open-source SciPy v1.2.1. software [28]. Parametric and non-parametric data are expressed as mean ± standard deviation (SD) or median (Q1–Q3), respectively. The Shapiro–Wilk test was used to verify the normality of the data. To assess the intra- and interobserver reproducibility, Bland–Altman analysis was performed, and the coefficient of variation (COV), Spearman rank correlation coefficient (r), and intraclass correlation coefficient (ICC) (two-way mixed effects, absolute agreement) were calculated. For the Bland–Altman analysis, the mean difference and limits of agreement (LoA = 1.96 × SD) were computed [29]. In addition, the relative difference in LoA between the training and validation cohort was calculated (LoAdiff = ((validation LoA–training LoA)/training LoA) × 100%). The COV was classified as low (≤10%), intermediate (11–20%), high (21–30%), or very high (>30%). The r and ICC were classified

as poor (<0.50), moderate (0.50–0.69), good (0.70–0.84), very good (0.85–0.94), or excellent (≥0.95) [25]. A *p*-value of <0.05 was considered statistically significant.

#### **3. Results**

The patient characteristics are presented in Table 1. In Table 2, median values of hemodynamic and morphologic parameters are presented. Of the selected TAA patients, 13 presented with an aneurysm in the aortic root and 7 with an ascending aorta aneurysm. The mean TAA diameter and mean blood velocity were 50.8 ± 2.7 mm and 42.5 ± 14.9 cm/s, respectively.

**Table 1.** Population characteristics.

Data notated as the mean ± standard deviation and presented over all subjects. Abbreviation: TAA—thoracic aortic aneurysm.

**Table 2.** *Cont.*

Data notated as the median (Q1–Q3) and presented over all subject and systolic phases. Abbreviations: AoR aortic root, pAA—proximal ascending aorta, dAA—distal ascending aorta, AoR—aortic arch, pDA—proximal descending aorta, and dDA—distal descending aorta.

A total of 2400 aorta segments were analyzed for the intra- and interobserver reproducibility analysis, see Figure 3. The result of these analyses over all anatomical segments and systolic phases for the validation and training cohort are presented in Table 3 and Supplementary Materials Table S2, respectively.

**Table 3.** Intra- and interobserver reproducibility: results of the validation cohort.

**Table 3.** *Cont.*

**Table 3.** *Cont.*

Characteristics are presented per cohort over all subjects, anatomical segments, and systolic phases (*n* = 300). All probability values <0.01. Abbreviations: Mean Diff—mean difference; LoA—limits of agreement (1.96× standard deviation); COV—coefficient of variation; IA-O—intraobserver analysis by the first observer with 4 years' experience; IE-O1—interobserver analysis by the first and second observer with 4 and 3 years' experience, respectively; IE-O2—interobserver analysis by the first and third observer with 4 and 1 years' experience, respectively; IE-O3—interobserver analysis by the second and third observer with 3 and 1 years' experience, respectively; WSS—wall shear stress; and LoAdiff –relative difference in limits of agreement between the training and validation cohort.

Intra- and interobserver analysis of the hemodynamic parameters demonstrated a good reproducibility in the validation cohort for normalized flow displacement and flow jet angle (LoA ≤4.8% and ≤8.7◦ respectively, COV = 16–32%, r = 0.70–0.92, ICC = 0.80–0.94), an excellent reproducibility for vorticity norm (LoA ≤2779 s−1·mL−1, COV = 3–10%, r = 0.96–1.00, ICC = 0.96–1.00), and a very good reproducibility for absolute WSSmag, WSSax, WSScir, WSSangle, absolute helicity and volumes with increased absolute LNH (LoA ≤242 mPa, ≤244 mPa, ≤103 mPa, ≤8.9◦, ≤631 m·s−2·mL−1 and ≤8.4 mL, respectively; COV = 3–26%, r = 0.86–1.00, ICC = 0.86–1.00). See Supplementary Materials Figures S1–S9 for the corresponding Bland–Altman plots. These plots demonstrated no striking dependency of the measured differences relative to the mean for the analyses and hemodynamic parameters.

Intra- and interobserver analysis of the morphologic parameters demonstrated a very good reproducibility in the validation cohort for volume and centerline length (LoA ≤15.7 mL and ≤12.2 mm, respectively; COV = 8–22%, r = 0.89–0.98, ICC = 0.90–0.98), and an excellent reproducibility for diameter (LoA ≤5.0 mm, COV = 3–7%, r = 0.98–0.99, ICC = 0.96–0.99) and a good reproducibility for curvature radius (LoA ≤46.4 mm, COV = 21–35%, r = 0.72–0.87, ICC = 0.90–0.92). See Supplementary Materials Figures S10–S13 for the corresponding Bland–Altman plots. These plots demonstrated no dependency of the measured differences relative to the mean for most analyses and morphologic parameters, except for the maximal diameter, which had a slightly increasing difference relative to the mean for the interobserver analyses.

The intra- and interobserver reproducibility results per anatomical segment are presented in Supplementary Materials Tables S3 and S4, respectively. The intra- and interobserver analysis showed for the validation cohort a similar reproducibility for most morphologic and hemodynamic parameters per segment. However, the reproducibility was slightly decreased within the aortic root for several parameters (i.e., flow displacement, WSSmag, WSSax, absolute helicity, volumes with increased absolute LNH, aortic segment volume, centerline length, and curvature radius). Furthermore, the reproducibility of the maximal diameter was slightly decreased within the proximal part of the aorta (i.e., the aortic root, proximal and distal ascending aorta) compared to the other anatomical segments further distally to the ascending aorta (intraobserver LoA ≤2.7 and ≤1.8 mm; interobserver LoA ≤7.6 and ≤3.3 mm, respectively). Moreover, the reproducibility of WSSangle was slightly decreased within the proximal and distal descending aorta.

**Figure 3.** Aortic segments and hemodynamic parameters. Example of a patient (male, 71 years old) with an ascending aorta aneurysm with a diameter of 48 mm. (**A**) Anatomical segments with centerline (black color) and flow displacement line (blue color). (**B**) Vorticity norm. (**C**) Absolute helicity. Positive and negative magnitudes indicate a clockwise and anti-clockwise rotation, respectively. (**D**) Wall shear stress magnitude. (**E**) Axial wall shear stress. (**F**) Circumferential wall shear stress. Abbreviations: AoR—aortic root, pAA—proximal ascending aorta, dAA—distal ascending aorta, AoR—aortic arch, pDA—proximal descending aorta, dDA—distal descending aorta, and WSS—wall shear stress.

Compared to the training cohort, the validation cohort demonstrated reduced LoA over all analyses for morphologic and hemodynamic parameters (LoAdiff = −12 [−26–3]%). This improvement was more evident for the morphologic parameters compared to the hemodynamic parameters (LoAdiff = −20 [−41–−6]% and −9 [−21–5]%, respectively). Moreover, this improvement was also more evident within the aortic root, proximal and distal ascending aorta compared to the aortic arch, proximal and distal descending aorta (LoAdiff = −10 [−35–11] % and 2 [−18–41]%, respectively).

#### **4. Discussion**

The main findings of the current study are: (1) hemodynamic parameters (i.e., normalized flow displacement, flow jet angle, WSSmag, WSSax, WSScir, WSSangle, vorticity norm, absolute helicity, and volumes with increased absolute LNH) and morphological parameters (volume, centerline length, diameter and curvature radius) can be assessed with very good reproducibility in the thoracic aorta of TAA patients with a TAV, (2) the reproducibility of these parameters is slightly decreased within the aortic root, and (3) the reproducibility is affected by the observer's experience.

The intra- and interobserver analyses generally demonstrated a very good reproducibility for the assessment of hemodynamic parameters in TAA patients with a TAV. Similar results have been described previously for flow displacement, flow jet angle, WSSmag, WSSangle, secondary flow patterns (which are reflected in vorticity and helicity), and volumes with increased absolute LNH in healthy volunteers and in patients with none-to-marginally dilated aortas [11–19,30], but not yet in a clinically relevant patient group, such as patients with TAA. Moreover, the Bland–Altman plots of the current study demonstrated for most hemodynamic parameters and analyses no striking dependency of the measured differences relative to the mean.

Furthermore, the intra- and interobserver analyses demonstrated very good reproducibility of the morphologic parameters in TAA patients with a TAV. Similar reproducibility results have been described previously for volume, centerline length, maximal diameter, and curvature radius in healthy volunteers [25]. However, the current study demonstrated wider LoAs for all morphologic parameters within TAA patients compared to healthy volunteers. These wider LoAs in TAA patients may potentially be associated with the presence of dilation since, for a constant flow rate, the blood velocity is inversely related to the vessel's diameter. Therefore, a larger area of reduced velocity-to-noise ratio is expected for TAA patients compared to healthy volunteers, especially close to the vessel wall. This reduced velocity-to-noise ratio potentially can have a negative influence on the reproducibility of aortic lumen segmentation. However, the velocity-to-noise ratio within these areas may be improved by 4D flow MRI sequences with multiple-velocity-encoding gradients, which consequently may also improve the reproducibility of hemodynamics parameters [31,32].

When analyzing the results per anatomical segment, it was observed that the reproducibility was slightly reduced but still acceptable for most hemodynamic and morphologic parameters within the aortic root. This slight reduction in reproducibility was also seen in the proximal and distal ascending aorta for the maximal diameter. This reduced reproducibility is potentially introduced by the movement of the aortic root and ascending aorta in contrast to the relatively fixed-positioned aortic arch and descending aorta [33–35]. This effect related to the cardiac movement (contraction and relaxation) may further be amplified by the heart rate variation and breathing variation during the acquisition. Alongside, inadequate respiratory compensation or inadequate electrocardiographic gating may also perturbate the 4D flow MRI data. Moreover, within the aortic root, complex secondary flow patterns (recirculating flow) are present [9,36]. Since this recirculating flow has a relatively low velocity compared to the velocity-encoding sensitivity, these recirculating areas have a reduced velocity-to-noise ratio close to the vessel wall and potentially introduce more observer-dependent variability on the aortic lumen segmentation.

The observer-dependent variability of the 3D lumen segmentation was assessed directly by evaluating the variability in morphologic parameters. This analysis also demonstrated, for example, wider LoAs for the maximal diameter within the proximal part of the aorta (i.e., aortic root, proximal and distal ascending aorta) compared to thoracic aorta distal to the ascending part (intraobserver LoA ≤2.7 and ≤1.8 mm; interobserver LoA ≤7.6 and ≤3.3 mm, respectively). In contrast, a recent study assessing the interobserver reproducibility of manually measured diameters imaged with computed tomography demonstrated in the aortic root, ascending and descending aorta an LoA of 2.5, 2.1, and 1.9 mm, respectively [37]. For MRI, it is recommended to measure the maximal diameter on an MRI sequence with an isotropic voxel size ≤1.5 mm to achieve a clinically acceptable accuracy [38]. Since the 4D flow MRI sequence applied in the current study had an isotropic spatial resolution of 2.5–3.0 mm, the morphologic parameters derived from this sequence should not be used to define the maximal aortic diameter of TAA patients. However, for the quantification of hemodynamic parameters, an isotropic resolution of ≤3.0 mm is indicated as sufficient for the 4D flow MRI [23]. Of note, the 4D flow MRI was acquired to assess the patient's aorta hemodynamics and therefore was not meant to describe the patient's aorta morphology (i.e., TAA diameter).

The decreased LoA within the validation cohort compared to the training cohort indicate that the reproducibility likely is affected by the observer's experience, especially for the aortic root, proximal, and distal ascending aorta. Moreover, the observer's experience likely plays a larger role in the assessment of morphologic parameters compared to the hemodynamic parameters, considering the decrease in LoA between the training and validation cohort. This result demonstrates the robustness of the image analysis for the quantification of hemodynamic parameters, despite the observer-dependent variability and experience.

No major limitations were observed for the hemodynamics quantification in the segments of the thoracic aorta. Although no clear reproducibility limitations were observed for TAA patients with a TAV, still some variation is present in the hemodynamic parameters. This variation must be taken into account when these hemodynamic parameters are used for the prediction of progressive dilation and adverse aortic events in TAA patients with a TAV.

A limitation of the present study was the small population size. However, the robustness of the study was improved by including multiple observers, analyzing five systolic phases, and dividing the thoracic aorta into six anatomical segments. This resulted in a total of 2400 lumen segments for which the hemodynamic and morphologic parameters were quantified. Moreover, future research should include patients with other aortic pathologies. Finally, for patients with TAA, no data were available on repeated 4D flow MRI acquisition; thus, interexamination reproducibility could not be assessed.

#### **5. Conclusions**

The current study demonstrated a very good reproducibility for all hemodynamic parameters within the thoracic aorta in TAA patients. This allows precise quantification of the patient's aorta hemodynamics from a 4D flow MRI.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/app12083912/s1, Figure S1: Bland–Altman plots for the normalized flow displacement; Figure S2: Bland–Altman plots for the flow jet angle; Figure S3: Bland–Altman plots for the wall shear stress magnitude; Figure S4: Bland–Altman plots for the axial wall shear stress; Figure S5: Bland–Altman plots for the circumferential wall shear stress; Figure S6: Bland–Altman plots for the wall shear stress angle; Figure S7: Bland–Altman plots for the vorticity norm; Figure S8: Bland–Altman plots for the absolute helicity; Figure S9: Bland–Altman plots for the absolute local normalized helicity volume; Figure S10: Bland–Altman plots for the volume; Figure S11: Bland– Altman plots for the centerline length; Figure S12: Bland–Altman plots for the maximal diameter; Figure S13: Bland–Altman plots for the curvature radius; Table S1: Characteristics of the 4D flow MRI Sequence; Table S2: Intra- and interobserver reproducibility: results of the training cohort; Table S3: Intraobserver reproducibility: results of both cohorts per anatomical segment; Table S4: Interobserver reproducibility: results of both cohorts per anatomical segment; Video S1: Pathline visualization of the 4D flow MRI of a patient with an aortic root aneurysm with a diameter of 52 mm.

**Author Contributions:** All authors were actively involved in reviewing and editing the manuscript. J.F.J. wrote the original draft version of the manuscript; J.F.J., B.J.C.t.K., R.L.F.v.d.P., P.v.d.B., I.A.D., A.J.H.A.S., H.J.L. and J.J.M.W. contributed to the data curation; J.F.J. and H.C.v.A. developed the in-house developed software; J.F.J., M.J.F.G.R. and B.J.C.t.K. performed the image analysis under supervision of J.J.M.W.; H.J.L., J.J.M.W., A.J.H.A.S., S.S., J.E.W. and B.P.A. contributed to the project administration and funding acquisition. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the Dutch Heart Foundation (CVON2017-08-RADAR).

**Institutional Review Board Statement:** The study was conducted in accordance with the Declaration of Helsinki and approved by the Institutional Ethics Committee of Leiden University Medical Center (project code G20.149, 9 October 2020).

**Informed Consent Statement:** Patient consent was waived due to the retrospective inclusion of anonymized clinical data. Patients were not subjected to any research procedures and had the option to withdraw their participation.

**Data Availability Statement:** Data can be made available upon reasonable request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Review* **Positron Emission Tomography in Coronary Heart Disease**

**José de Almeida 1,2,\*, Sofia Martinho 1, Lino Gonçalves 1,2 and Maria Ferreira 1,2,3,\***

**Abstract:** With advances in scanner technology, postprocessing techniques, and the development of novel positron emission tomography (PET) tracers, the applications of PET for the study of coronary heart disease have been gaining momentum in the last few years. Depending on the tracer and acquisition protocol, cardiac PET can be used to evaluate the atherosclerotic lesion (plaque imaging) and to assess its potential consequences—ischemic versus nonischemic (perfusion imaging) and viable versus scarred (viability imaging) myocardium. The scope of this review is to summarize the role of PET in coronary heart disease.

**Keywords:** positron emission tomography; coronary heart disease; plaque imaging; perfusion imaging; viability imaging; 18F-sodium fluoride; 18F-fluorodeoxyglucose; Rubidium-82; 13N-ammonia; 15O-water

#### **1. Introduction**

PET has been increasingly used in the medical field since its discovery, in the sixth decade of the last century, due to both technological advances and the expansion of indications. Fifty years after its development, the number of detectors increased from 64 to 19,000 and the spatial resolution of the exam from 14 mm to 4 mm [1]. These detectors measure the total energy deposited by the two annihilation photons moving in opposite directions produced after positron emission from a radionuclide-tagged tracer molecule [2]. The tracer molecule is part of the biochemical process that we are interested in quantifying and as such, the radionuclide should ideally only interact with that molecule [3]. The raw data acquired is then corrected for attenuation and finally reconstructed to provide an estimate of the in vivo tracer distribution [2].

The evolution of PET has had a huge impact in the study of cardiovascular disease, including coronary heart disease, heart failure, cardio-inflammatory disease, valvular heart disease, and assessment of cardiac devices and cardiac tumors.

Coronary heart disease is defined as an inadequate blood flow to an area of myocardial tissue due to blockage of the blood vessels that supply it, most often by an atherosclerotic plaque [4].

PET can be used to evaluate the atherosclerotic lesion (plaque imaging) and to assess its potential consequences—ischemic versus nonischemic (perfusion imaging) and viable versus scarred (viability imaging) myocardium [5].

#### **2. Plaque Imaging**

Although coronary atheromatous plaques are often present in older patients, in most cases their development is silent and will not result in a cardiovascular event. In contrast, some atherosclerotic plaques will rupture and cause myocardial infarction. Identifying these "vulnerable" plaques before they rupture and differentiating them from their "stable" counterparts is therefore a key objective for a cardiac imaging technique [6].

**Citation:** de Almeida, J.; Martinho, S.; Gonçalves, L.; Ferreira, M. Positron Emission Tomography in Coronary Heart Disease. *Appl. Sci.* **2022**, *12*, 4704. https://doi.org/10.3390/ app12094704

Academic Editor: Julio Garcia Flores

Received: 20 March 2022 Accepted: 4 May 2022 Published: 7 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

PET scanners have the capacity to assess any biochemical process occurring in the atherosclerotic plaque, if the appropriate surrogate tracer molecule is used.

#### *2.1. Imaging Inflammation with 18F-Fluorodeoxyglucose (FDG)*

FDG is overall the most used radioligand in PET imaging. It is a glucose analog that enters cells via facilitated glucose transporter member (GLUT) 1 and 3 and undergoes phosphorylation to become (18)F-FDG-6-phosphate, which cannot exit the cell before radioactive decay. Therefore, signal intensity correlates with cellular glucose uptake and phosphorylation, which indicates high metabolic activity. In the case of the atherosclerotic plaque, metabolic activity appears to be related to the concentration of proinflammatory macrophages, as demonstrated in animal models [7]. In humans, serum levels of myeloperoxidase were associated with carotid plaque FDG uptake [8].

Aortic FDG uptake was independently associated with cardiovascular risk factors such as increased low-density lipoprotein and total cholesterol, and with the presence of metabolic syndrome [9–11]. On average, patients with myocardial infarction (MI) were shown to have higher aortic FDG uptake compared with stable angina patients [12]. Visceral adipose tissue uptake of FDG was associated with systemic inflammatory status, and with the presence of metabolic syndrome components [13]. Additionally, it was independently associated with the severity of CAD and with the occurrence of AMI [14].

Regarding the assessment of pharmacological treatment, statins were shown to reduce the arterial FDG signal in a dose-dependent manner and pioglitazone was shown to attenuate vascular FDG uptake [15–17]. In contrast, the novel antidyslipidemia drugs Dalcetrapib (cholesteryl ester transfer protein inhibitor) and Rilapladib (lipoprotein-associated phospholipase A2 inhibitor) did not have any effect on vascular FDG activity [18,19].

There is, however, no current evidence to support the implementation of FDG-PET for atherosclerotic risk stratification in clinical practice.

#### *2.2. Imaging Microcalcification with 18F-Sodium Fluoride (NaF)*

NaF has been used to detect bone metastases and it is known to replace the hydroxyl group of hydroxylapatite in areas of calcification [20,21]. Regarding vascular calcification, two processes occur in a continuum in the atherosclerotic plaque as a healing response to the necrotic core inflammatory mediators: microcalcification, followed by macrocalcification. The stage of microcalcification renders the plaque unstable and more predisposed to rupture, making its identification an interesting surrogate marker for detecting "vulnerable" plaques. NaF-PET was shown to be able to identify microcalcification in atherosclerotic plaques in a consistently different pattern of uptake to the macroscopic calcium observed on computer tomography (CT) [22] (Figure 1).