Quantitative Analysis of Vortical Blood Flow in the Thoracic Aorta Using 4D Phase Contrast MRI

Introduction Phase contrast MRI allows for the examination of complex hemodynamics in the heart and adjacent great vessels. Vortex flow patterns seem to play an important role in certain vascular pathologies. We propose two- and three-dimensional metrics for the objective quantification of aortic vortex blood flow in 4D phase contrast MRI. Materials and Methods For two-dimensional vorticity assessment, a standardized set of 6 regions-of-interest (ROIs) was defined throughout the course of the aorta. For each ROI, a heatmap of time-resolved vorticity values ω→=∇v→ was computed. Evolution of minimum, maximum, and average values as well as opposing rotational flow components were analyzed. For three-dimensional analysis, vortex core detection was implemented combining the predictor-corrector method with λ2 correction. Strength, elongation, and radial expansion of the detected vortex core were recorded over time. All methods were applied to 4D flow MRI datasets of 9 healthy subjects, 2 patients with mildly dilated aorta, and 1 patient with aortic aneurysm. Results Vorticity quantification in the 6 standardized ROIs enabled the description of physiological vortex flow in the healthy aorta. Helical flow developed early in the ascending aorta (absolute vorticity = 166.4±86.4 s-1 at 12% of cardiac cycle) followed by maximum values in mid-systole in the aortic arch (240.1±45.2 s-1 at 16%). Strength, elongation, and radial expansion of 3D vortex cores escalated in early systole, reaching a peak in mid systole (strength = 241.2±30.7 s-1 at 17%, elongation = 65.1±34.6 mm at 18%, expansion = 80.1±48.8 mm2 at 20%), before all three parameters similarly decreased to overall low values in diastole. Flow patterns were considerably altered in patient data: Vortex flow developed late in mid/end-systole close to the aortic bulb and no physiological helix was found in the aortic arch. Conclusions We have introduced objective measures for quantification of vortical flow in 4D phase contrast MRI. Vortex blood flow in the thoracic aorta could be consistently described in all healthy volunteers. In patient data, pathologically altered vortex flow was observed.


Introduction
Phase contrast MRI allows for the examination of complex hemodynamics in the heart and adjacent great vessels. Vortex flow patterns seem to play an important role in certain vascular pathologies. We propose two-and three-dimensional metrics for the objective quantification of aortic vortex blood flow in 4D phase contrast MRI.

Materials and Methods
For two-dimensional vorticity assessment, a standardized set of 6 regions-of-interest (ROIs) was defined throughout the course of the aorta. For each ROI, a heatmap of timeresolved vorticity values o ! ¼ r v ! was computed. Evolution of minimum, maximum, and average values as well as opposing rotational flow components were analyzed. For threedimensional analysis, vortex core detection was implemented combining the predictor-corrector method with λ 2 correction. Strength, elongation, and radial expansion of the detected vortex core were recorded over time. All methods were applied to 4D flow MRI datasets of 9 healthy subjects, 2 patients with mildly dilated aorta, and 1 patient with aortic aneurysm.

Results
Vorticity quantification in the 6 standardized ROIs enabled the description of physiological vortex flow in the healthy aorta. Helical flow developed early in the ascending aorta (absolute vorticity = 166.4±86.4 s -1 at 12% of cardiac cycle) followed by maximum values in mid-systole in the aortic arch (240.1±45.2 s -1 at 16%). Strength, elongation, and radial expansion of 3D vortex cores escalated in early systole, reaching a peak in mid systole (strength = 241.2±30.7 s -1 at 17%, elongation = 65.1±34.6 mm at 18%, expansion = 80.1±48.8 mm 2 at 20%), before all three parameters similarly decreased to overall low values in diastole. Flow patterns were considerably altered in patient data: Vortex flow developed late in mid/end-systole close to the aortic bulb and no physiological helix was found in the aortic arch.

Introduction
Phase contrast (PC) pulse sequences allow for the assessment of moving fluids in magnetic resonance imaging (MRI). Using volumetric data acquisition with multi-directional velocity encoding, complex hemodynamics in the heart and the adjacent great vessels can be assessed [1]. Resulting time-resolved volume datasets provide both the exact quantitative assessment of flow velocities and the intuitive visualization of flow phenomena using three-dimensional stream-and pathlines [2,3]. Numerous 4D flow MRI studies have demonstrated that the presence of vortex flow patterns in the heart chambers [4][5][6], the aorta [7], the carotid sinus [8], and pulmonary circulation [9] is physiological, but also relates to certain pathologies including aortic aneurysms [10][11][12][13][14], pulmonary hypertension [15], bicuspid aortic valves (BAV) [16,17], and congenital heart defects [18].
While physiological helices may help to stabilize flow, preserve energy, and protect vessels from atherosclerosis [7,19], pathologically altered vortical flow may have a number of detrimental effects. First, axial flow velocities are reduced by vortical flow, thus less blood is passing through a constant vessel section per unit time, resulting in a possible decrease of organ perfusion [20] or increase of cardiac work [21] (reduction of blood current). Second, non-physiological vortex patterns may exert strain on the vessel wall, resulting in dilatation or in endothelial disorders [16,[20][21][22] (morphological alterations). Third, the velocity of blood platelets is reduced inside a vortex, possibly promoting blood clots or thrombosis (rheological abnormalities).
First attempts of quantitative evaluation of vortical flow have been made. Hope et al. [11] employed a visual approach. Number of vortices, their start/end/peak frame, angle, and size were analyzed. However, such subjective measures might be inconvenient for clinical practice. The Helical Flow Index (HFI) was proposed as a more objective parameter measuring helical flow based on pathlines [23]. The average of local normalized helicity (LNH, equal to relative helicity H r described below) was calculated along each pathline trajectory. To obtain a measure of global helicity, data of all pathlines were averaged. However, pathlines do not necessarily cover the entire aortic lumen and helicity quantification during diastole might be difficult. Bissell et al. [16] quantified rotational flow by use of circulation Γ [24], which is defined as the line integral of fluid velocity on the boundary of a closed surface (i.e., a vessel's cross-section). Using Stokes's theorem, it can be shown that Γ equals the integral of vorticity over the enclosed area. To reduce noise in the data, 5 time frames around peak systole were averaged for rotational flow measurements. A drawback of this approach is, however, that time-dependent information is lost. Bächler et al. [9] proposed a semi-automated approach by placing twodimensional planes inside a vessel, exporting velocity values, and plotting absolute maximum vorticity values over time using a separate software tool. By this approach, however, important contextual and three-dimensional information is lost. Lorenz et al. [25] used equidistantly distributed 2D planes along the entire thoracic aorta to quantify altered blood flow helicity in bicuspid aortic valve (BAV) patients as compared to healthy volunteers. Helicity density H d (being the dot product of the velocity and the correspondent vorticity vector) and relative helicity H r (equaling the cosine of the angle between both vectors) were averaged over each analysis plane, parts of the aortic lumen, or periods of the heart cycle. Limitations of this approach include the insufficient extraction of vortices without any longitudinal flow component by the helicity parameter [18], information loss on minor counterrotating fluid components due to averaging [11,23], as well as omission of important three-dimensional information in the twodimensional analysis. Recently, Elbaz et al. [6] proposed parameters for the quantification of left ventricular (LV) vortex ring formation during early and late diastole. The location and orientation of vortex rings (defined as λ 2 isosurfaces) were characterized using normalized cylindrical coordinates relative to the shape of the left ventricle (i.e., circumferential, longitudinal, and radial coordinates; the ring's orientation relative to the LV axis; and the circularity index defined as the ratio between the vortex's short and long diameter). However, these parameters are not valid for anatomical regions other than the left ventricle and for vortical structures not equaling a vortex ring; in addition they do not provide insight into the strength or temporal evolution of the vortex. Three-dimensional vortex cores have been identified in previous studies, but have not been quantified [8,22,26,27].
The objective of the present work is to define quantitative metrics for objective characterization of complex vortical flow features (both two-and three-dimensionally), which can be calculated in a fully automated way. The metrics are applied to the analysis of the physiological vortex in the aortic arch during systole [7,10,13,22,25,[28][29][30][31]. MR datasets acquired from 9 healthy subjects and 3 patients with dilated aortas are used to demonstrate the value of the method.

Terminology
The terms vortex and helical flow are not consistently defined in the literature. In medical studies, the term helix often describes fluid rotating around an axis while moving forward. The term vortex denotes rotating motion, where stream-or pathlines tend to curl back on themselves [22]. In the fluid dynamics domain, the terms vortex and vortex core are often used in a broader sense including helical movement [8,[32][33][34][35][36]. Explanatory notes on the difficulty of an exact theoretical definition of the term vortex can be found in [32]. In this paper, the term vortex will be used in the wider definition of the fluid dynamics domain.

Vorticity
To provide a basis for further methods outlined below, calculation of vorticity of the timeresolved three-dimensional vector field derived from 4D flow MRI measurements was implemented. The vorticity ω is defined as the curl of the velocity vector field v ! [37] (Fig 1A): ; @v x @z À @v z @x ; @v y @x À @v x @y : Thus, a vorticity vector is given for every point-the norm of which describes the strength, while the direction represents the orientation of swirl. Positive values represent counter-clockwise (CCW) rotation, while negative values stand for clockwise (CW) fluid motion. The unit of vorticity is given as s -1 .
In addition, the vorticity delivers information about the global form of a vortex (becoming maximum in the vortex's center) due to the tangential velocity profile inside the vortex and its first derivative as shown in Fig 1B.

Lambda 2
To facilitate further vortex characterization, a popular approach commonly referred to as the λ 2 -method [38] was chosen. To define λ 2 , let v ! ðx; y; zÞ :¼ ðv x ; v y ; v z Þ T be the underlying velocity vector field. The velocity gradient or Jacobian of the velocity field is a 3x3 tensor containing information on how the velocity is changing in space [39]: J v ¼ @v x @x @v x @y @v x @z @v y @x @v y @y @v y @z @v z @x @v z @y @v z @z : J v can be decomposed into its symmetric part S and its antisymmetric part O. Physically, S describes the strain-rate tensor (containing both bulk and shear components), while O describes the rotation tensor (containing the rotational part of the flow) [38,39]: Because the 3x3 matrix S 2 +O 2 is real and symmetric, it has 3 real eigenvalues. Let λ 1 , λ 2 , and λ 3 be the eigenvalues with λ 1 λ 2 λ 3. According to Jeong and Hussain [38], a vortex can be defined as a region where two of the eigenvalues are negative, that is λ 2 < 0. According to its definition, λ 2 is a scalar quantity. The main advantages of the method are its reliability in vortex identification [6,33,34,36,38] and its superiority in the context of 4D flow MRI [18]. Since derivatives are examined instead of velocity values, constant deltas added to the flow field cancel each other out, making the method Galilean invariant and applicable for time-varying flow fields [32,33]. Since only adjacent grid cells are involved in the calculation, it is also computationally inexpensive [32,33].

Calculation of Vorticity Fields
While the vorticity of a differentiable vector field is defined infinitesimally, the spatial resolution of experimental MR data is limited. Multiple algorithms have been described for the precise computation of the vorticity field including [40,41]: • first order central difference, • higher order central difference.
Originally defined for two-dimensional vector fields, the algorithms were extended in a three-dimensional fashion and implemented for verification and comparison in our study (Fig 2). Resulting vector fields were rated by visual comparison. Evaluation criteria included noise, consistency of the vector field, and identification of vortex features.

2D: Region-of-Interest Based Analysis
We aimed for an efficient integration of ROI based vortex analysis (Fig 3) yielding two-dimensional quantification of vortices (Figs 4 and 5). Therefore, a standardized set of 6 ROIs (Fig 4A) was defined for each dataset, 0: directly above the aortic bulb, 1: horizontal slice through the ascending aorta, 2: proximal to the brachiocephalic artery, 3: between the left common carotid artery and left subclavian artery, 4: distal to the left subclavian artery, 5: axial slice through the descending aorta at the level of ROI 1. All regions-of-interest were defined on cut planes In the bottom left corner, the calculation scheme of each algorithm is given, illustrating adjacent data points that are taken into account for derivative calculation (schemes similar to [40]). While all three algorithms were adopted to a three-dimensional pattern for this work, 2D patterns are shown for simplicity. perpendicular to the vessel lumen and restricted by the vessel wall. Positions were chosen according to standardized image interpretation guidelines [42,43] and slightly adjusted to fit into the specific context. For each ROI, a heatmap of time-resolved through-plane vorticity values was computed. The evolution of minimum, maximum, and average values throughout the cardiac cycle was plotted and analyzed (Figs 4 and 5). Since these values only contain information about vortex strength, the percentage of CW and CCW values exceeding an arbitrary threshold (set to 10 s -1 to suppress noise) was analyzed to estimate the amount of contrary rotational flow (Fig 5).
Resulting data of all healthy subjects were averaged for analysis; quantitative values are given as mean ± standard deviation. Patient data were analyzed individually to highlight pathologically altered flow.

3D: Vortex Core Analysis
Since a vortex has a spatial extent, it should also be studied in its native three-dimensional domain.
Visual exploration. To ease the understanding of flow characteristics, stream-and pathlines can be limited to a subset of lines already containing all relevant information, while irrelevant lines are abandoned to reduce visual distraction [44,45]. Multiple quality criteria for line extraction have been examined in previous studies with λ 2 giving the best results for vortex identification [18]. In this project, characteristic sets of stream-and pathlines dependent on λ 2 were defined to provide an instant overview of vortical flow features contained in the datasets (Fig 3C).
Vortex core detection. A vortex is intuitively understood as the swirling motion of fluid around a central region [32]. Multiple algorithms have been described to find the core of a vortex to accurately characterize its length, formation in space, and temporal evolution [32][33][34][35][36]38]. In the given scope of application, the vortex core algorithm needs to fulfill certain prerequisites [32]: (1) The core should be represented as a single line to enable length measurement. In addition, measuring the radial expansion of the vortex should be possible. (2) Due to Therefore, an extension of the predictor-corrector method [35] was chosen for vortex core identification, since it conforms to all aforementioned pre-conditions. Assuming that a vortex reveals points of minimum pressure in its core, seed points p 0 are defined as local pressure minima. In every step (Fig 6A), the algorithm elongates the vortex core by following the vorticity vector ω i in the current point p i (prediction step). The new position p i+1 is then corrected to the point of minimum pressure in a plane perpendicular to the vorticity vector ω i+1 (correction step). The core line is extended by repeating both steps until a predefined stopping criterion is reached.
The described method is robust due to its self correcting character (Fig 6B and 6C), but also has certain drawbacks: First, it does not work at low Reynolds numbers [33,38]. Second, pressure is hard to compute from a velocity field [33]. To overcome these issues, we used a method combining the predictor-corrector method with λ 2 , which is easy to compute from the flow field yet robust and effective for vortex identification. Instead of searching for the pressure minimum in the correction step, the vortex core is now corrected to the minimum λ 2 value in the perpendicular plane [33].
Unlike the original suggestion by Banks and Singer [35], we applied a fourth order Runge-Kutta scheme [41] to extend the core line along the vorticity vector in the prediction step. In the correction step, we applied pattern search optimization [46] for finding the λ 2 minimum rather than the originally proposed method of steepest descent [35]. This approach has the advantage of being derivative-free and resulted in more robust and faster convergence.
To quantify the cross section of the vortex, we radially sampled the vorticity field in the perpendicular plane in each correction step. Eight radial lines were tracked until |ω| fell under a given threshold (set to 150 s -1 ). The area πr 2 of a circle with radius r (with r equal to the average line length) was calculated as a measure for the vortex's radial extent.
Vortex core quantification. The following three-dimensional parameters of the vortex core were recorded over time: • Vortex core length. This value contains information about the longitudinal extent of the vortex.
• The minimum, maximum, and average vorticity throughout the course of the vortex core.
These values hold information about the strength of the vortex. Results of all healthy subjects were averaged for further analysis.

Implementation Specific Details
All described methods were implemented in C++ as a separate module inside the commercially available 4D flow software package GTFlow (GyroTools LLC, Zurich, Switzerland). For mathematical computation, the linear algebra template library Eigen (version 3.0.b2, eigen.tuxfamily. org) was used. For 3-dimensional visualization, Coin3D (version 3.1.3, www.coin3d.org) was employed. Coin3D is an OpenGL based graphics library compatible with the Open Inventor 2.1 API.

Sample Data
To evaluate the accuracy of the proposed techniques and to verify their implementation, artificial datasets were created. Results of all visualization and quantification methods analyzing these test data were compared to the mathematically exact solution. The test data contained a vortex superimposed on the laminar flow profile through a simple tube (Fig 1C). Given a linear tube with a radial cross-section of radius R, the parabolic flow profile can be derived [47]:

DP Dx
with v m = maximum velocity, ΔP = pressure difference, Δx = tube length, μ = viscosity, r = radius at point under consideration, and R = total tube radius. The Lamb-Oseen equation describes the temporal evolution of a planar vortex dependent on the fluid's viscosity and was used for vorticity simulation in previous studies [40,48,49]. The tangential velocity equals: vðr; tÞ ¼ G 2pr 1 À exp Àr 2 r 2 c ðtÞ with r and μ as defined above, r c ðtÞ ¼ ffiffiffiffiffiffiffi 4mt p = core radius of vortex, and Γ = circulation contained in the vortex. To define a smooth progression of the vortex, v(r,t) was multiplied by 0.5 (cos(πz / l)+1) for |z| l with l = vortex length.
The first order central difference was used for vorticity calculation of the sample data, since both other methods (i.e., 8-point-circulation and higher order central difference) introduce smoothing [40], which is advantageous in noisy datasets, but adds imprecision to numerically exact data.

MR-Datasets
To demonstrate applicability of the newly defined methods to in-vivo data, a cohort of 9 healthy subjects was scanned. 4D flow MRI datasets covering the left-ventricular outflow tract, ascending, and descending aorta were acquired on a Philips 3T Achieva system (Philips Healthcare, Best, The Netherlands). Scan parameters of the respiratory navigated and retrospectively ECG triggered sequence were: FOV = 320x256x46 mm 3 , matrix = 160x160x20, SENSE reduction factor = 2, TR = 3.05-3.94 ms, flip angle = 5°, Venc x,y = 200 cm/s, Venc z = 100 cm/s, cardiac phases = 25, resulting scan resolution = 2.3x2.3x2.3 mm 3 . Next, 2 patients with mildly dilated ascending aortas were scanned (35.6 mm and 36.8 mm diameter of the proximal ascending aorta, respectively). Last, 1 patient with aortic aneurysm was scanned (49.2 mm diameter of the proximal ascending aorta). In this case, adapted scan parameters included: FOV = 350x280x75 mm 3 , matrix = 176x176x30, TR = 4.10 ms, cardiac phases = 20, Venc x,z = 50 cm/s, Venc y = 60 cm/s. Vessels were segmented by applying a threshold on the velocity component product weighted by the magnitude image value. Cubic temporal interpolation of factor 10 was employed for refinement of limited resolution MR data. Cubic spatial interpolation of dynamic degree was applied at points reached during vortex core calculation.

Ethics Statement
All subjects gave written informed consent prior to MR scanning. The study complied with the Declaration of Helsinki and was approved by the ethics committee of the Canton of Zurich.

Calculation of Vorticity Fields
Vorticity fields in the aortic arch as computed by the three alternative algorithms are shown in Fig 2. Visual comparison indicated sufficiently accurate vorticity detection using the first order central difference (2A) and the 8-point-circulation method (2B). As can be seen in the figure, the latter method showed less vulnerability to noise in the evaluated dataset, resulting in a more consistent vorticity field. The second order central difference did not deliver reasonable results (2C).

Method Validation on the Basis of Sample Data
The exact knowledge of flow characteristics in the artificial sample dataset allowed validation of all described methods (Fig 1C). Characteristic stream-and pathlines clearly depicted the vortex in the linear tube. A ROI intersecting the tube showed the correct vorticity profile as expected from the definition of the Lamb-Oseen equation. Also the vortex core was correctly identified in the center of the tube. Numerical data were in good accordance with anticipated results for higher resolution test data ( Table 1). Analysis of low resolution test data showed higher discrepancy relative to the mathematically exact results.

2D and 3D Vortex Visualization
The combined use of characteristic stream-or pathlines, vorticity or λ 2 heatmaps inside the predefined ROIs, and automatic vortex core identification yielded intuitive three-dimensional visualization of the studied flow patterns (Figs 1C, 3C, 4C and 7A). Vortical flow features hidden in the multiplicity of stream-or pathlines could be easily identified and traced.

2D: Region-of-Interest Based Analysis
By comparing the temporal evolution of vorticity values in the 6 standardized ROIs, it was possible to estimate at which point in space and time a vortex developed, how it propagated, and where and when it faded out.   Table 2) developed during early systole. It could be observed first in ROI 0 and 1, which were located directly above the aortic bulb and in the ascending aorta, respectively (ROI 0: minimal vorticity = -166.4±86.4 s -1 at 12% of cardiac cycle and ROI 1: -183.6±65.3 s -1 at 16%). Thereupon, rotational flow propagated to the more distal ROIs 2, 3, and 4 with minimum values seen in mid systole (ROI 2: -240.1±45.2 s -1 at 16%, ROI 3: -229.2±50.0 s -1 at 20%, and ROI 4: -201.6±92.7 s -1 at 16%). The CW rotational flow diminished with overall low values during diastole. In the descending aorta, only minor rotational flow was present over all time points (ROI 5: -125.0±33.2 s -1 at 20%).
Results were different for CCW rotational flow (positive values in Fig 5 and Table 2). Again, vorticity increased during early systole, followed by maximum values in mid-systole and a final decrease in late systole. However, CCW rotational flow components were weaker during the entire heart cycle. They showed noticeable attenuation in the ascending aorta and proximal aortic arch with the progression of the right-handed helix (ROI 1: 118.8±67.5 s -1 at 16%, ROI 2: 106.0 ±43.1 s -1 at 20%, and ROI 3: 89.9±30.7 s -1 at 20%). Percentage plots of opposing flow components (Fig 5, right graphs) underline this finding (higher CW percentage of 27-29% compared to CCW percentage of 13-17% in ROIs 2, 3, and 4). Even though CCW rotational flow components are present over all time points, they did not lead to a fully developed vortex in the aorta.

3D: Vortex Core Analysis
Three-dimensional vortex core analysis is illustrated in Fig 7. The automatically identified vortex core followed the minimum λ 2 values and was surrounded by the right-handed helical flow (7A). The core line exactly passed regions of minimum λ 2 as can be seen for the intersected ROIs. The serrated nature of the core line was possibly caused by noise in the underlying MR velocity data. However, if gone astray, the core line was re-adjusted in every correction step of the identification algorithm so that single noisy data points did not compromise the rest of the vortex core identification process.
In Fig 7B, results of quantitative vortex core analysis are plotted. Vorticity strength increased in early systole and reached its maximum in mid systole (average vorticity = 241.2±30.7 s -1 at 17% of cardiac cycle). In diastole, the vortex strength decreased finally reaching overall low values (> 50% of cardiac cycle). Slightly delayed in the beginning, vortex elongation paralleled the strength values. The elongation reached its peak in mid systole (65.1±34.6 mm at 18%), before decreasing to low values in diastole. The vortex's radial expansion slightly hurried ahead in early systole and declined faster after peaking in mid systole (80.1±48.8 mm 2 at 20%).

Pathologically Altered Flow
Pathologically altered vortex flow was observed in patient data. Fig 8 shows pathlines and average vorticity plots of the 6 standard ROIs for each case.
A comparison of all three patients indicated that changes in vortical flow got more pronounced as the aortic dilatation progressed (patient 1 < 2 < 3). In patient 1, the physiological helix through central ROIs was still visible (ROI 3: minimum of -115.8 s -1 at 24% of cardiac cycle), even though flow was considerably altered by a strong perpendicular swirl in the ascending aorta. In patient 2, alterations were more pronounced and the physiological helix further vanished (ROI 3: -108.0 s -1 at 20%). In patient 3, swirling flow developed late in mid/ end-systole close to the aortic bulb (ROI 1: -229.1 s -1 at 30%), and no physiological helix was found in the aortic arch (ROI 3: -44.1 s -1 at 10%).
A vortex core line was identified in the ascending aorta of all patients by the vortex core detection algorithm, but results were not fully satisfactory and stable.

Discussion
In our study, we implemented different visualization techniques to allow for intuitive exploration of vortical flow features of time-resolved velocity vector fields. 2D and 3D quantitative metrics for the evaluation of vortex flow were introduced. After validation against artificial sample data, the techniques were successfully applied to quantify physiological vortex flow in the aortic arch of healthy volunteers and to identify pathological flow alterations in 3 cases of aortic dilatation.
Flow patterns in the aortic arch have been examined in a number of previous studies. Kilner et al. [7] were the first to describe a predominant right-handed helical flow in the upper aortic arch in mid to late systole in 9 out of 10 subjects, while flow patterns in the descending aorta varied greatly among subjects. By now, the right-handed helical flow is a well acknowledged Quantitative Vortex Analysis in the Thoracic Aorta Using 4D PC MRI physiological flow feature found in most subjects [10,13,22,25,[28][29][30][31], especially in younger individuals with crook-shaped arch geometry [29]. In accordance with this we found a right- handed vortex arising in the ascending aorta early in systole and propagating to the aortic arch with peak vorticity values reached in mid systole. While also counter-clockwise rotational flow components were measurable, no fully left-handed vortex developed in the aorta. Multiple studies have examined the influence of age [29][30][31], geometry [29], and vascular pathologies [20,30] on aortic flow characteristics. Already minor geometric variations in the aortic size and shape were found to considerably alter blood flow characteristics as compared to healthy subjects [20]. In bicuspid aortic valve disease, increased right-handed helical flow has been described to be the most common flow alteration in the ascending aorta, possibly leading to increased rotational wall shear stress and progressive aortic dilatation [16,25]. Nested helical flow was found at peak systole [17]. In ascending aortic aneurysms, pathological findings included helical and vortical flow [10][11][12]14] as well as prolonged retrograde flow [11]. Comparing aortic aneurysms in more detail, a large variety of pathological vortex flow patterns has been reported depending on the individual cause of dilatation [13]. In our small cohort of patients the severity of hemodynamic aberrations seemed to be gradually related to the degree of aortic dilatation. While in patients with mild aortic ectasia physiological flow features were still discernable, the patient with aortic aneurysm showed pronounced alterations of flow patterns and vorticity values.
All these studies suggest an important role of flow alterations in the pathophysiology of aortic disease. The complexity of flow patterns, however, makes it difficult to reliably identify flow features that determine disease progression and prognosis just by means of visual analysis. In order to advance the knowledge of potentially detrimental effects of certain hemodynamic changes a more systematic analysis is needed, which was the aim of our study.
We therefore compared three algorithms for vorticity calculation, differing in the set of surrounding points that are taken into account for computation. The 8-point-circulation method gave the best subjective results in our specific case, being in good accordance with the literature [40]. Poor results of the second order central difference scheme may have been caused by the low resolution of MRI data; instead of adding valuable information, more distant points may hamper precise derivative computation. Extensive theoretical comparison and uncertainty analysis of the alternative algorithms can be found in [40].
In our study, all methods were validated in a set of artificially defined sample data. Arising deviations from mathematically exact solutions were carefully inspected. We found that the implementation of our methods worked correctly with the given data. However, since values of the tangential velocity profile shift quickly inside a vortex center (Fig 1B), data of low spatial resolution tend to miss relevant velocity values and thereby compromise precise vorticity calculation. These findings indicate that highly resolved data play a decisive role in accurate flow pattern analysis.
With the visualization techniques implemented in our study we were able to consistently visualize vortical flow in the aorta of all scanned individuals. The combination of characteristic stream-or pathlines, standardized ROIs, vorticity/λ 2 heatmaps, and automatically identified vortex core lines allowed identification of regions with vortical flow.
In addition to a more consistent visual exploration, objective quantitative metrics are needed to better characterize flow pattern changes occurring under pathological conditions. Limiting the examination of flow data to certain, carefully chosen regions-of-interest is a wellknown technique widely used in the field of 4D flow MRI [50]. It is used for both the quantification of blood flow velocities and the qualitative description of flow patterns. Even though parts of the vector field might not be fully examined, ROIs help the investigator to abstract from the overwhelming complexity of the flow field and to concentrate on certain areas containing most relevant information. In this study, we transferred this commonly used technique to vorticity quantification: Two-dimensional quantitative analysis in consecutive, standardized ROIs was useful for following vortex paths over time. By these means, we were able to consistently detect and quantify the typical development of right handed helical flow in the aortic arch in all our healthy subjects, fully accordant with findings by Lorenz et al. [25] Three-dimensional quantification was realized by analyzing the elongation, strength, and radial expansion of an automatically detected vortex core. A quantitative description of flow characteristics in the aortic arch observed in 9 healthy subjects was performed with results being in good consistency with physiological flow qualitatively described in the literature [10,13,22,[28][29][30][31]. The three-dimensional algorithm led to less robust results in patient data, which is associated with lower velocity-to-noise ratio (VNR) and low vorticity values. In Fig 1, the Lamb-Oseen equation was introduced as a theoretical approximation of the velocity profile inside a vortex. Accordingly, the vorticity reaches maximum values inside the vortex center due to maximum shift of velocity. If overall velocity is low or the vortex broadens, vorticity values decrease, resulting in a less stable output of the identification algorithm.

Study Limitations
While the focus of our study was the definition and validation of techniques for visualizing and quantifying swirling blood in 4D flow MRI data, the number of healthy subjects and patients enrolled in the study remained small. Even though the measures of vorticity proved stable in all healthy volunteers and were in good agreement with findings in the literature, a larger cohort is needed to establish reference values for vorticity in the healthy aorta. Moreover, further studies with larger numbers of patients are needed to test whether alterations of vortical flow in the aorta are related to disease severity or may even help to predict disease progression.
In our study the methods were applied to a single anatomical region only. Future work should include in-depth analysis of different vascular territories in order to elucidate the role that vortical flow may play in other anatomical regions (e.g., heart chambers, pulmonary circulation, or carotid arteries) and pathologies (e.g., pulmonary hypertension or congenital heart defects).

Conclusions
In the present work, two-and three-dimensional metrics for the quantification of vortical blood flow in the thoracic aorta have been proposed and validated. While this study can serve as a proof of concept, future work should assess whether and how these techniques can help to diagnose and grade flow related pathologies and their prognosis.