Vortex Analysis of Intra-Aneurismal Flow in Cerebral Aneurysms

Sunderland, Kevin; Haferman, Christopher; Chintalapani, Gouthami; Jiang, Jingfeng

doi:https://doi.org/10.1155/2016/7406215

Computational and Mathematical Methods in Medicine

On this page

Abstract Introduction Materials and Methods Results Conclusions Abbreviations Disclosure Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2016 | Article ID 7406215 | https://doi.org/10.1155/2016/7406215

Vortex Analysis of Intra-Aneurismal Flow in Cerebral Aneurysms

Kevin Sunderland,¹Christopher Haferman,¹Gouthami Chintalapani,²and Jingfeng Jiang¹

Academic Editor: Hernán G. Morales

Received17 May 2016

Revised27 Aug 2016

Accepted15 Sept 2016

Published07 Nov 2016

Abstract

This study aims to develop an alternative vortex analysis method by measuring structure ofIntracranial aneurysm (IA) flow vortexes across the cardiac cycle, to quantify temporal stability of aneurismal flow. Hemodynamics were modeled in “patient-specific” geometries, using computational fluid dynamics (CFD) simulations. Modified versions of known and -criterion methods identified vortex regions; then regions were segmented out using the classical marching cube algorithm. Temporal stability was measured by the degree of vortex overlap (DVO) at each step of a cardiac cycle against a cycle-averaged vortex and by the change in number of cores over the cycle. No statistical differences exist in DVO or number of vortex cores between 5 terminal IAs and 5 sidewall IAs. No strong correlation exists between vortex core characteristics and geometric or hemodynamic characteristics of IAs. Statistical independence suggests this proposed method may provide novel IA information. However, threshold values used to determine the vortex core regions and resolution of velocity data influenced analysis outcomes and have to be addressed in future studies. In conclusions, preliminary results show that the proposed methodology may help give novel insight toward aneurismal flow characteristic and help in future risk assessment given more developments.

1. Introduction

A weak or thin spot on cerebral arteries, resulting in the bulging or ballooning (dilation) of that vessel, is known as an intracranial aneurysm (IA). Most common forms of IAs are saccular and fusiform IAs. These irregular dilations act like cavities around the blood stream and, therefore, induce flow disturbance and subsequent destructive vascular remodeling [1] in and around IAs. IAs are a potentially life-threatening vascular malformation affecting an estimated 3% of the population [2]. According to the US National Institute of Neurological Disorders and Stroke, the incidence of reported IA rupture is about 0.001%. In other words, there are approximately 30,000 people who could suffer a stroke from ruptured cerebral aneurysms each year in the United States alone. The severe consequences of aneurysm rupture, the cost and risk of aneurysm treatment, the low statistical risk of aneurysm rupture (1% or less per year), and the lack of reliable parameters predictive of the risk of aneurysm rupture, all motivate the continued search for a method to differentiate those aneurysms which are likely to rupture from those at a low risk of rupture.

Currently, geometrical and morphological parameters of IAs (e.g., size and irregular shape) have been used as the basis for clinical predictions of aneurysm rupture, subsequently directing the course of clinical intervention(s), according to recommendations by multiple healthcare professional organizations [3–5]. However, those predictions are not highly accurate. For instance, based on receiver operating characteristic analysis, the size ratio and aneurysm angle had an area under the curve values of 0.83 and 0.85, respectively [6]. Recently, attention has been devoted to aneurysmal hemodynamics, since it is well known that the origin and natural history of IAs are closely associated with disturbed hemodynamics [7, 8]. Consequently, phase-contrast magnetic resonance imaging (PC-MRI) [9, 10] has been used to measure aneurysmal flow and examine factors relevant to the development and progression of IAs in vivo. But disturbed flow in and around IAs makes this task difficult. Recall that individual protons in complex and disturbed aneurysmal flow have incoherent velocities (at the subgrid level) and those components cannot be resolved by one “averaged” velocity measurement from a relatively large resolution cell (at 1 mm scale). Therefore, such measurement errors and potential artifacts can adversely affect the accuracy of PC-MRI results. In parallel to research efforts of PC-MRI flow imaging in and around IAs, blood flow parameters simulated from “patient-specific” computational fluid dynamics (CFD) simulations [11] have also drawn a lot of interest by the clinical and research community [12, 13]. Hemodynamic factors such as wall shear stress [14], oscillatory shear index [15], flow impingement [13], and flow stability [16] have emerged as potential parameters correlating to the risk of IA rupture. Among those efforts, Byrne et al. [16] found that aneurysmal flow, classified by its spatial complexity and temporal stability using vortex core lines, is closely correlated with the risk of aneurysm rupture. Particularly, they concluded that “ruptured aneurysms had more complex and more unstable flow patterns than unruptured aneurysms.” In their work, Byrne et al. have used proper orthogonal decomposition [17] of time-resolved velocities to characterize temporal flow stability. In another study, Kohler et al. [18] presented another approach to visualizing vortex flow using line predicates. However, their technique was primarily for visualization of swirling flow in relatively large arteries such as the aorta and pulmonary artery.

In this study, we intended to develop an alternate technique, also based on the flow vortex core analysis, to characterize the presence, destruction, and motion of spatial swirl flow patterns during the cardiac cycle. More specifically, we implemented an image processing algorithm to first delineate the vortex core region(s) and then quantitatively investigate temporal flow stability by tracking segmented vortex cores over the cardiac cycle. The vortex extraction method was based on velocity data on a rectilinear grid; potentially, this extraction method could be used for PC-MRI measured velocity fields, as well. Toward this end, the primary purpose of this work is to explore and analyze the structural characteristics of flow vortexes as a possible means to complement current assessment methods of temporal flow stability.

The rest of the paper is structured as follows. In Materials and Methods, using original 3D velocity estimates, the proposed image processing method is described in detail. To demonstrate its feasibility, the proposed method is applied to “patient-specific” CFD simulated velocity data using two distinct solvers (a commercial CFD solver using unstructured meshes and a research prototype solver using the Lattice-Boltzmann Method [LBM]). As part of this preliminary study, two CFD solvers were selected in order to investigate whether or not results of the proposed method of vortex core analysis may be influenced by CFD solvers. Outcomes of this feasibility study are reported in Results, followed by Discussions and Conclusions.

2. Materials and Methods

Extraction and tracking of vortex cores are built around time-resolved velocity fields obtained from “patient-specific” CFD simulations. Figure 1 shows two major steps: (1) hemodynamics simulations and (2) vortex core extraction and analysis. The hemodynamic flow was first simulated in vascular models reconstructed from 3D-Digital Subtraction Angiography (DSA) scans and then those time-resolved 3D velocity fields were used to determine vortex core regions within the aneurysmal sac. Of note, the identification of vortex core regions was derived from two known methods: method [19] and the -criterion method [20]. All identified vortex core regions were then segmented out using the classic marching cube algorithm [21] to form 3D triangulated surfaces. Degree of vortex overlap (DVO) among the segmented vortex cores over the cardiac cycle and variations of the number of vortex core regions over the cardiac cycle were used to assess destruction and motion of spatial swirl flow patterns within the aneurismal sac over the cardiac cycle.

Figure 1

Main steps in the generation and analysis of vortex core structures: (a) arterial geometry generated and segmented from a DSA scan; (b) volumetric mesh generation and CDF simulation on the segmented structure; (c) manual extraction of an aneurysm surface from the artery (gray denotes the area of an aneurysm and red denotes the parent artery); (d) extracted aneurysm surface masked over the simulated data to isolate aneurysm flow velocity and resampled at various voxel sizes (0.1–0.8 mm); (e) the and -criterion methods identify vortexes for 21 equally spaced data points over the cardiac cycle; (f) analysis of vortex: characteristics and DOV (stability) between each cardiac step’s vortex core and the cycle-averaged vortex core. In (e), identified vortex core regions are segmented out using the classic marching cube algorithm. In (e), while the white color denotes streamlines of the aneurysmal flow, the black color is used to indicate vortex structures using the two above-mentioned methods— and -criterion methods.

2.1. Modeling of “Patient-Specific” Hemodynamics

2.1.1. Image Segmentation

As shown in Figure 2, 10 cases of IAs were arbitrarily selected from our internal database: in five cases, each has one single terminal aneurysm (TA); in the other five cases, each contains single sidewall aneurysm (SA). Initial 3D-DSA data files were loaded into a commercially available image segmentation package (Mimics Innovation Suite, version 17, Materialise Inc., Leuven, Belgium) where intensity-based image segmentation was first used to isolate the regions of interest and create a volumetric surface of patients’ vascular structures. During this process, a first-order Laplacian smoothing filter in the 3-Matic Software (version 9, Materialise Inc., Leuven, Belgium) was used to reduce irregularities caused by imaging artifacts, while preserving surface shape. Then, manual editing (e.g., removal of kissing vessels) and visual inspection were performed to ensure the integrity of all vessel geometries if necessary. Downstream vessels were shortened in order to reduce computing time. The longest possible upstream section, proximal to an aneurysm, was left intact to help reduce the effects of the inlet(s) and plug-flow alterations [22]. Cylindrical flow extensions (at least 6 times of the vessel diameter) were added to each model (both inlets and outlets) to mitigate influences of boundary conditions using the open-source Vascular Modeling Toolkit (VMTK) software (version 1.2) [23].

2.1.2. Mesh Generation

Each processed vascular surface model (see Figure 2) was converted into an unstructured, 3D, tetrahedralized volumetric mesh using an open-source mesh generator, Tetgen (version 1.4.2) [24]. The mesh generation process was done by an in-house Python script derived from the VMTK program (version 1.2) [23]. Approximately, 1 million computing cells were used per case, with the average mesh size as 0.0022 mm³. The typical edge length of resultant tetrahedrons was roughly 0.12 mm. Mesh sensitivity tests were performed to ensure that results were not sensitive to the mesh size selected.

2.1.3. CFD Simulations

To compute velocity waveforms in and around an aneurysm using the ANSYS-FLUENT solver (version 14.0; ANSYS-FLUENT Inc., Lebanon, NH), we solved the time-dependent incompressible, 3D Navier-Stokes equations [i.e., (1)] for the 3D meshed vessel geometry. The equations for velocity are written aswhere is the three-dimensional velocity vector, is the blood density, is the pressure, and is the viscosity. In the ANSYS-FLUENT CFD solver, the pressure-velocity coupling in (1) is obtained using the SIMPLEC algorithm [25]. The explicit time-marching second-order scheme with a time step of 1 × 10⁻³ seconds (approximately 1000 steps per cardiac cycle) was used for the computations. Although this time step is still relatively coarse, it seems adequate to capture gross flow features. The tetrahedralized volumetric meshes described above were used by the ANSYS-FLUENT solver. More details can be found in our previous publications [26, 27].

The Siemens research prototype CFD solver (version 4.0, Prototype: not for diagnostic use, Siemens Medical Solution, Inc., IL) is based on the Lattice-Boltzmann Method (LBM) with a Multiple Relaxation Time (MRT-LBM) approximation [28]. More details can be found in a previously published paper [28]. In the current implementation, a graphic processing unit (GPU) was used to accelerate the flow simulations. Prior to CFD simulations using the Siemens solver, IA models defined by water-tight 3D surface triangles were automatically discretized with cubical voxels; a smooth level-set function was created to differentiate the vessel and nonvessel regions for each IA model. For each IA, its complex geometry was represented by the zero-crossing of the above-mentioned smooth level-set function, which was used for the accurate imposition of the boundary conditions. All subsequent LBM computations were performed on a Cartesian grid, where the global resolution and time step size were automatically adapted by the Siemens LBM solver for ensuring computational stability. Of note, the exact same vessel geometries (STL files) were used to generate volumetric meshes (for the ANSYS-FLUENT solver) and voxel discretization (for the LBM solver). It is important to note that the Siemens research CFD solver is only an investigative prototype and not commercially available.

In both solvers, vessel walls were assumed rigid with a no-slip boundary condition, and blood flow was considered incompressible and Newtonian. The dynamic viscosity was 0.004 kg/m-s and the mass density of blood was 1050 kg/m³. The inlet boundaries of all models were located in the internal carotid artery or the basilar artery. The zero-pressure condition was used for all outlets. Two pulsatile flow rate waveforms, at a heart rate of 60 bpm, derived from magnetic resonance measurements and taken from Gwilliam et al. [29] were used since patient-specific waveforms were not available. Of note, one waveform was adopted to represent the pulsatile flow at the BA (TA1 and TA2; see Figure 2) and the other waveform was selected to represent the pulsatile flow at the ICA (all other cases; see Figure 2). Each case had one velocity waveform scaled according to the vessel’s inlet cross-sectional area, standardizing the mean volumetric flow rate of 280 mL/min (ICA) or 180 mL/min (BA). Our choice of mean flow rates is consistent with physiological flow rates available in the MR literature [30, 31]. Four cardiac cycles were simulated per case as a means to reduce initial transient conditions, and only the final cycle was saved. Twenty-one equally spaced data points across the final cardiac cycle were saved for subsequent data analysis detailed below.

2.2. Vortex Core Extraction and Analysis

All computational methods for extraction of vortex core regions and analysis of those vortex cores were performed using in-house scripts (C++ and Python) that were derived from the open-source VTK/VMTK (version 1.2) software package [23]. Statistical and correlation analysis was performed using MATLAB scripts.

2.2.1. Aneurysm Extraction and Voxelization of Aneurismal Velocity Data

A published method [32] was used to semiautomatically isolate and extract the aneurysm sac (see white colored surface in Figure 1(c)). The isolated aneurysm sac was sealed at the orifice and converted to a binary mask that is spatially registered with volumetric velocity data. The mask allowed us to analyze intra-aneurysmal velocity data only. To verify the intrarater reliability of proper sectioning of aneurysm masks, two separate users (Kevin Sunderland and Christopher Haferman) sectioned all aneurysms, and Bland-Altman plots [33] were performed on the mask volumes and ostium areas to determine the similarity between chosen masks. Once no significant difference was ensured between sectioned masks, one user was chosen at random and all masks from that user were implemented in the rest of the study.

Simulated ANSYS-FLUENT velocity data are typically located on an unstructured grid. Therefore, resampling of the simulated velocity data was transformed onto a (3D) rectilinear grid, while the velocity data obtained from the Siemens research CFD solver was already on a rectilinear grid. Both data sets were resampled onto a rectilinear grid of 0.2 mm voxel size using the vtkProbeFilter function provided by the open-source VTK package (Kitware Inc., NY, USA), as a means to standardize the resolution of each data set.

In addition, two cases of simulated velocity data using ANSYS-FLUENT were resampled at varying voxel sizes (0.1–0.8 mm) to determine how spatial resolution of velocity data may impact the characteristic structure of generated vortex cores and stability over the cardiac cycle. Determining the changes that voxel size could have on vortex structure characteristics would help to establish proper analysis parameters for this technique.

2.2.2. Vortex Identification

Multiple algorithms exist to locate and extract vortexes in the CFD literature [34]. This study employed two classic methods: method developed by Jeong and Hussain [19] and the -criterion by Hunt et al. [20]. More specifically, given a CFD-simulation velocity field , the velocity gradient can be decomposed as follows [35]:where is the rate of strain tensor and is the vorticity tensor.

Hunt et al. [20] defined a vortex as a spatial region where The vortex core defined by Hunt et al. (i.e., (3)) essentially denotes locations where the Euclidean norm of the vorticity tensor dominates.

Jeong and Hussain suggested that the vortexes are regions where In (4), means the second intermediate eigenvalue of the 3 × 3 tensor . The tensor is symmetric, and therefore all three eigenvalues are real.

In this study, we attempted to use [normalized] and values to condense the distributions of original and values so that the subsequent extraction of vortex cores could be reliably performed. Mathematically, this normalization process can be written as follows: where is the amplitude velocity value.

Once all and values, including [normalized] and values, were determined within the dome of the aneurysm, the classic marching cube algorithm [36] was used to extract the vortex core regions and convert them to surface meshes (see Figure 1(e)). In order to extract a vortex core region, the mean value of interest ( [positive] or [negative]) for each case was used as a threshold value.

To determine the ideal vortex identification methodology and desired range of analysis parameters, each methodology (standard , [normalized] , -criterion, and [normalized] -criterion) was run with 5 different threshold values and their resultant vortex structure was analyzed. The 5 threshold values were selected from the following respective ranges: (1) [0.1, maximum] for -criterion and [normalized] -criterion and (2) [minimum, −0.1] for standard and [normalized] methods in each case. In other words, in an ascending order, selected five (5) threshold values were [mean − STD, mean − (STD/2), mean, mean/2, and mean/4] for the standard and [normalized] methods and [mean/4, mean/2, mean, mean + (STD/2), and mean + STD] for the -criterion and [normalized] -criterion methods.

In order to reduce the appearance of small, isolated areas being mistaken for the dominant vortex structure, only extracted vortex core regions which have a greater volume than 0.5 mm³ were counted.

2.3. Data Analysis

2.3.1. Vortex Stability over the Cardiac Cycle

For each CFD simulation, the vortex structure was generated for all 21 data points along the cardiac cycle, as well as a structure based on the cycle-averaged vorticity. The DVO between each data point’s vortex structure and the cycle-averaged vortex structure were calculated. The mean and standard deviation (STD) of the DVO across all data points for the cardiac cycle were calculated to quantify the average temporal stability of the vortex structure. “Stability” is referring to the tendency of a vortex structure to shift or change over the course of a cardiac cycle. A visual representation of variations in vortex structure (temporal and spatial) over a cardiac cycle can be seen in Figure 3.

Figure 3

Example of variations to vortex core structure across 5 data points along the cardiac cycle: marked waveforms represent a point in the cardiac cycle for each data point. The black structure is the extracted vortex core(s) while gray is an aneurysm. The top row is for SA2 which had a lower DVO (less stable) than the bottom row from TA2 which had a higher DVO (more stable). All cores were extracted using the [normalized] -criterion method. The mean [normalized] threshold values were used to extract the vortex core and only vortex cores with a volume ≥ 0.5 mm³ were saved.

2.3.2. Geometric and Hemodynamic Parameters of IAs

To discern the independence of this method in comparison to established methods for aneurysm analysis, it was explored whether these findings correlated with a number of geometrical [6] and hemodynamic [37] characteristics used to predict IA rupture. Six geometrical characteristics and two hemodynamic characteristics were calculated for each IA using in-house VMTK scripts.

Geometric Characteristics(1)Aneurysm volume (mm³)(2)Ostium (neck) area (mm²): surface area of the opening of an aneurysm(3)Ostium circumference (mm)(4)Aneurysm height (mm): midpoint of ostium to furthest point on aneurysm wall(5)Aspect ratio: aneurysm height/ostium diameter (widest distance across ostium)(6)Aneurysm volume/ostium area (mm) [38]

Hemodynamic Characteristics(1)Spatially Averaged Oscillatory Shear Index (SA-OSI)(2)Spatially and Temporally Averaged Wall Shear Stress (STA-WSS)

The calculated geometric values were verified by performing a secondary measurement on the removed aneurysm surface in the 3-Matic Software package. The geometric parameters were collected for all aneurysm masks that were originally segmented from two users (Kevin Sunderland and Christopher Haferman). Once no significant differences were detected between those two users’ masks, their averaged geometrical values were used for this study.

Wall shear stresses [WSSs] obtained from the ANSYS-FLUENT solver were first collected from each aneurysm wall for each data point along the cardiac cycle. Then, the WSSs of each spatial location on the aneurysm wall, for each data point along the cardiac cycle, were temporally averaged to obtain the time-averaged wall shear stress for the spatial location. Spatial averaging of all TA-WSSs on the aneurysm wall was performed to obtain the STA-WSS for each aneurysm. To measure the WSS directional oscillations over a cardiac cycle, the quantitative metric known as oscillatory shear index (OSI) was calculated using the following equation [37]:where is the number of data points along the cardiac cycle and is the instantaneous wall shear stress vector. Spatial averaging of OSI (i.e., SA-OSI) for each aneurysm was calculated.

Pearson’s linear correlation was performed to compare the mean DVO and change of the number of vortex cores over the entire cardiac cycle to the aforementioned geometrical and hemodynamic parameters. Only the cases sampled at a resolution of 0.2 mm (voxel) were used for this analysis.

3. Results

Results presented in Sections 3.1–3.4 were based on the ANSYS-FLUENT solver. In Section 3.5, the ANSYS-FLUENT results were quantitatively compared to those obtained from the Siemens LBM research CFD solver.

3.1. Comparison of Analysis Methods

As shown in Figure 4, choosing high (or [normalized] ) values or low (or [normalized] ) values significantly reduced the sizes of identified vortex core regions. Across all 10 cases studied, setting the threshold larger than resulted in a breakdown in the overall identified vortex volume, with no areas of vortex being identified in some cases using the normalized methods (see Figure 4(b)). Extending the threshold past mean/2 (increasing the range of values) resulted in a significant increase in identified vortex volume, as regions corresponding to weaker swirling flow were included as part of the vortex core. A visual representation of these changes can be seen in Figure 5.

(a)

(b)

Figure 4

Plots representing variations of identified vortex volumes due to the selection of 5 different threshold values in (a) TA1 and (b) SA1. Threshold values were tested for four vortex extraction methods: standard -criterion, standard method, [normalized] -criterion, and [normalized] method. Geometries of SA2 and TA1 are displayed in Figure 2. Selected threshold values were mean/4, mean/2, mean, mean ± (STD/2), and mean ± STD (positive sign for -criterion and [normalized] -criterion methods and negative sign for and [normalized] methods).

Figure 5

Plots illustrating the impact on extracted vortex cores (i.e., the black surface in each plot) due to changes of isosurface threshold values. From left to right: mean/4, mean/2, mean, mean + (std/2), and mean + std. All images were from TA case 1, [normalized] with a voxel resolution of 0.2 mm. At increased threshold values (≥mean + std/2) a reduction of identified vortex core structures occurred, sometimes not identifying any areas of vortex. All images came from the [normalized] -criterion methodology, case TA1 (see Figure 2), and cardiac cycle-averaged vorticity data.

It is worth noting that, in terms of the absolute volume change, the [normalized] -criterion and [normalized] method provided more consistent results (i.e., smaller changes due to the threshold variation), as shown in Figures 4(a) and 4(b). This is because the standard -criterion and standard methods resulted in a broader area of swirling flow identified as vortex core in relation to their normalized counterparts.

Comparison of the extracted vortexes with velocity streamlines in TA1 and SA2 showed that the normalized methods better identified the central region of swirling flow as part of the vortex without including large portions of the weaker, outer swirling flow, as shown in Figure 6. For subsequent data analysis, the [normalized] -criterion method was chosen using each case’s mean [normalized] value as the threshold value to identify the vortex core while excluding weaker flow patterns from analysis.

Figure 6

Comparison of identified vortex cores from the four mentioned methods: (a) standard , (b) [normalized] , (c) standard -criterion, and (d) [normalized] -criterion. Top row is from SA2 case, and bottom row is from the TA1 case. Velocity streamlines were added to represent simulated flow patterns. Geometries of SA2 and TA1 can be found in Figure 2. Threshold values for each case were the mean (for their representative value), and only extracted vortex cores with a volume > 0.5 mm³ are shown. Figures are from each method’s vorticity averaged data.

3.2. Voxel Size on Vortex Characteristics

Two aneurysm cases studied in the previous subsection (i.e., SA2 and TA1) were resampled at varying voxel sizes and analyzed with the [normalized] -criterion method to determine how spatial resolution of velocity data impacts vortex core analysis. In order to reduce the appearance of small, isolated areas being mistaken for the dominant vortex structure, only connected vortex core regions which have a greater volume than 0.5 mm³ were counted. Through a visual inspection, we found that as voxel size changed, structural characteristics of the vortex core(s) were altered, which could lead to a misinterpretation of the temporal flow stability within an aneurysm. As shown in Figure 7, a high resolution (<0.2 mm) causes the vortex structure to become more fragmented, which may subsequently cause areas of vortex core to be overlooked as their total structural volume may now be smaller than 0.5 mm³. By contrast, in lower resolutions (≥0.6 mm), the overall broad structure of extracted core(s) may dwarf subtle changes that could occur in flow patterns leading to a seemingly more stable vortex structure. Figure 8 shows how changes to voxel sizes impact the mean DVO and the mean number of cores over the cardiac cycle.

(a)

(b)

For the purposes of a concise presentation of the rest of study, only the data with a resolution of 0.2 mm were analyzed in an attempt to analyze the finer details of the hemodynamic flow patterns, all while still identifying the dominant areas of vorticity.

3.3. Vortex Core Parameters in 10 IA Cases

Table 1 shows a summary of the numbers of vortex cores and DVO values among 10 IA cases investigated. Overall, the mean DVO values were similar in TAs (), as compared to those in SAs () with no statistical difference ( using a -test) between aneurysm types. The mean numbers of vortex cores seemed slightly higher in SAs (2.35 ± 1.89) as compared to those in TAs (1.88 ± 1.26) but lacked statistical significance ( using a -test).

3.4. Vortex Correlation with (Aneurysm) Geometrical Parameters

Pearson’s linear correlation was used to determine the connection between vortex characteristics (number of vortex cores and DVO) and aneurysm geometric parameters. Table 2 shows that there are high positive correlations between the number of vortex cores within the aneurysm sac and the volume (, ), height (, ), and volume to ostium ratio (, ) of SA, whereas no statistically significant correlation exists in TAs between the number of cores and any geometrical parameters. In contrast, no statistically significant correlation exists in SAs or TAs between the DVO values and any geometrical parameters.

To account for the possibility (or lack thereof) of correlations between vortex characteristics and geometrical parameters being attributed to the low number of cases for TA and SA, the data was combined and correlations were recalculated. When data from both aneurysm types were combined, correlations were found between the mean number of cores in relation to 3 geometrical parameters: volume (, ), height (, ), and volume/ostium ratio (, ). Slight negative correlations were also found between the mean DVO of the combined data and aneurysm ostium area (, ) and the ostium circumference (, ).

Due to the significant volume difference between individual aneurysms (and number of vortex cores), the coefficient of variation was calculated for each case to normalize their relative change of cores over the cardiac cycle (Figure 9). The large spike of the above-mentioned variation in SA3 and TA5 is due to the fact that no coherent cores (>0.5 mm³) were detected for a number of data points during the cardiac cycle, thereby causing a reduced mean and increased standard deviation, resulting in a coefficient of variation near or above 100. As shown in Table 3 there only exists a correlation between SAs and their aspect ratio (, ) when dealing with the coefficient of variations (vortex cores). All other geometrical parameters held no significant correlation with stability in number of cores across the cardiac cycle. The lack of significant correlations indicates that using the aneurysm geometric parameters (outside of possibly using aspect ratio when dealing with SA) to infer the construction or destruction of vortex cores over the cardiac cycle may not be an ideal methodology.

3.5. Comparison with Hemodynamic Parameters

Pearson’s linear correlation was also performed to compare the vortex core characteristics to the hemodynamic properties of STA-WSS and SA-OSI present in the aneurysms for this study. Upon initial analysis, no statistically significant correlation was seen in either aneurysm type when comparing the mean OSI across a whole aneurysm, to both the coefficient of variation (number of cores) and the mean DVO. When comparing the STA-WSS values, a significant correlation was only seen with the coefficient of variation for SAs (, ). Yet when all data from both IAs were combined, no statistically significant correlative values could be seen for WSS when dealing with coefficient of variations (cores) or the mean DVO across the cardiac cycle. Calculated values for all correlations can be seen in Table 4.

3.6. Comparison of Two Different CFD Solvers

Vortex core results generated from the ANSYS-FLUENT and Siemens CFD solvers are summarized below. Only results from the [normalized] -criterion method at a resolution of 0.2 mm were used. A visual comparison between the vortex cores generated from those two solvers was performed to ensure both their spatial position and quality. It was seen that both solvers generated results that led to the identification of the main vortex structure in a similar location. Two representative examples are provided in Figure 10. Visual inspection of the surface quality of identified vortex cores showed that cores from the Siemens LBM simulation had a smoother, less nebulous structure than that of the ANSYS-FLUENT simulation data.

(a)

(b)

(c)

(d)

Figure 10

Vortex cores identified from two (SA4 and TA2) ANSYS-FLUENT-simulated cases (a and c) and the same cases using the Siemens LBM-simulated velocity fields from the Siemens CFD solver (b and d) under the same boundary conditions. Visual inspection ensured that main vortex core structures occurred in the same generalized location while Siemens CFD solver (b and d) resulted in a smoother core surface. Both cases had cores extracted using the same parameters: [normalized] -criterion; threshold is their mean of [normalized] , and only vortex cores with a volume > 0.5 mm³ were saved. Each vortex core was extracted from the cycle-averaged vorticity data (per respective case).

Analyses of Bland-Altman method [33] were performed as a means to identify the difference between the two CFD platforms (Figure 11). In each plot of Figure 11, solid and dashed horizontal lines represent the limits of agreement (95% confidence interval) and the biases, respectively. Middle (black) lines in Figures 11(a) and 11(b) calculated by the Bland-Altman method indicate there were minimal biases found between the two above-mentioned CFD solvers in terms of the number of vortex cores and the DVO.

(a)

(b)

4. Discussions

Recently, using “image-based” CFD simulations to characterize intra-aneurysmal blood flow, with the primary goal of searching for correlations between local hemodynamics and the risk of aneurysm rupture, has drawn significant interests. In order to impact the clinical management of IAs, aneurysmal hemodynamics have to be processed in a clinical workflow. However, visual assessment of 4D hemodynamic characteristics, through pattern recognition using hundreds of cross-sectional images from multiple cardiac phases, may impose an extensive burden on reviewing physicians. For instance, spatiotemporal characteristics of vorticity patterns, which are the primary objective of this study, are difficult to assess through visual inspection of time-resolved 3D velocity vector fields. Having a computational methodology that can isolate and assess characteristics of vortex cores as a means to understand hemodynamic patterns within an IA may lead to valuable insight toward aneurysm risk assessment and is advantageous as compared to visual assessments of aneurysmal hemodynamics.

For this study, as inspired by an early work by Byrne et al. [16], we proposed an alternate method to identify vortex cores. In the original study by Byrne et al. [16], vortex core lines were extracted from tetrahedral meshes as discrete line segments. They used the length change of vortex core lines during a cardiac cycle as a surrogate for “flow complexity.” However, in the current study, vortex core patterns were derived based on isosurface extraction of the [normalized] values (see (3)). Overall, the initial results (Tables 1 and 2) showed that DVO had no significant correlation with several known geometrical parameters for individual aneurysm types (TA and SA), while the normalized variation (coefficient of variation) of the number of vortex cores only correlated with the aspect ratio of SA. This shows that the two vortex characteristic parameters do not lead to redundant information when compared to geometrical parameters. In addition, limited correlation could be found between the vortex core characteristics and the hemodynamic parameters of STA-WSS and SA-OSI. Therefore, our work introduced two additional parameters, which potentially contribute to a fuller feature space toward assessing swirling aneurismal flow patterns. More work is planned to further test the correlation between the proposed vortex parameters and other important hemodynamic parameters suggested by Chung and Cebral [39].

Our preliminary results (see Table 1) indicated that the mean DVO values were lower (not statistically significant, possibly due to the low number of cases) in SAs as compared to those in TAs. That is consistent with an anecdotal observation that (temporal) swirling flow patterns are more versatile among sidewall aneurysms over a cardiac cycle.

In Section 3.5, we compared results between two different solvers. The primary goal was to demonstrate that the proposed method can work well with velocity data in a rectilinear grid. Our initial visual inspection (Figure 10) showed that the placement and relative shape of vortex cores identified by both solvers were fairly consistent. While a bias toward an increased number of identified cores existed in the Siemens CFD solver over the ANSYS-FLUENT solver, the degree of difference remained minimal (Figure 11(a)). However, comparisons of DVO values showed larger discrepancies (Figure 11(b)) which were not surprising. Recall that vortex core region in this study was derived from the local velocity gradient matrices. Small differences in a 3D velocity field may be amplified during the process of calculations of and [normalized] values (see (3)). Nevertheless, these discrepancies require more careful investigations in future studies.

Although the proposed method could theoretically work with PC-MRI data, which are in rectilinear grids, more work is needed in order to explore whether or not the proposed method is directly applicable. This is because measurement noise in PC-MRI-measured velocity fields could be amplified when the velocity gradient matrices are being evaluated (see (2)). Furthermore, the spatial resolution of PC-MRI is still limited and is typically between 0.5 mm and 1.0 mm voxel size. The relatively large voxel size may compromise our ability to accurately extract clinically relevant vortex core regions (see Figure 7) as this method is susceptible to variations in voxel sizes. Additional studies will be needed to understand the proposed characteristics of the vortex, in which whether a broader view (larger voxel size) of the vortex structure is still able to help identify clinically relevant information toward aneurysm rupture risk.

Our results in Section 3 indicate that both the threshold values used for vortex identification and the resolution of velocity data impact the DVO and number of identified cores in an IA. This limitation is similar to what is known in image segmentation. For this feasibility study, in order to mitigate this shortcoming posed by threshold selection, we adopted the [normalized] -criterion method for two reasons. First, we found that the adaptation of [normalized] -criterion condensed the overall range of the or [normalized] values. This condensing of threshold values helps determine an acceptable range, while mitigating some of the large-scale vortex changes (less change in terms of the absolute volume change as shown in Figures 4(a) and 4(b)) that occur with the standard -criterion. Second, our preliminary results suggested that using the [normalized] -criterion method helped identify the central structures of vortex patterns (see Figure 6), while limiting the identification of outlying weaker vortex patterns unless very low threshold values are used (Figure 5). Inclusion of weak boundaries of the vortex structure may indeed result in large, improper variations in estimated DVO values. In future studies, we intend to develop more automated methods based on machine learning [40]. The integrated machine learning approach will directly connect how the choices of vortex threshold value(s) and the resolution of velocity data impact the extraction of clinically relevant information. Nevertheless, the above-mentioned limitations have to be addressed through additional studies, particularly, through clinical studies, to determine clinical values of the proposed method of vortex core analysis.

Other limitations of our study include the small sample size and the use of standardized flow blood waveforms. However, our study design is still appropriate study design because our primary objective was to evaluate the proposed method of vortex core analysis using velocity data obtained from “patient-specific” CFD simulations. To keep this study concise, we limited the scope of the current study.

5. Conclusions

A vortex core analysis algorithm for time-resolved 4D velocity fields was presented as a means to gain novel insight into aneurysmal hemodynamic characteristics. The proposed method was designed to extract vortex characteristics from velocity data in a rectilinear grid and therefore could be used for both CFD simulated data and, in theory, PC-MRI measured velocity field. Although this technique probably only represents an alternative to other plausible approaches, it does, in our opinion, represent a feasible path to make intra-aneurysmal hemodynamic assessments more quantitative, thereby enabling studies and comparisons of large populations both initially and over time. While this current method gives no additional direct insight toward risk of aneurysm rupture and is susceptible to variations in threshold values and resolution of velocity data, the preliminary results are encouraging in that our findings do not have a high degree of correlation with commonly used geometrical and hemodynamic variables. We believe that the method warrants further studies to explore its full clinical utility and significance.

Abbreviations

3D-DSA:	Three-dimensional-Digital Subtraction Angiography
BA:	Basilar artery
CFD:	Computational fluid dynamics
DVO:	Degree of vortex overlap
IA:	Intracranial aneurysm
ICA:	Internal carotid artery
LBM:	Lattice-Boltzmann Method
MR:	Magnetic resonance
PC-MRI:	Phase-contrast magnetic resonance imaging
OSI:	Oscillatory shear index
SA-OSI:	Spatially Averaged Oscillatory Shear Index
STA-WSS:	Spatially and Temporally Averaged Wall Shear Stress
STD:	Standard deviation
STL:	Stereolithography
TA-WSS:	Temporally Averaged Wall Shear Stress
WSS:	Wall shear stress
VTK:	Visualization Toolkit
VMTK:	Vascular Modeling Toolkit
3D:	Three-dimensional.

Disclosure

The concepts and information presented in this paper are based on research and are not commercially available.

Competing Interests

Researchers from Michigan Technological University (Jingfeng Jiang, Kevin Sunderland, and Christopher Haferman) were partially supported by a research contract from Siemens to Michigan Technological University. Ms. Gouthami Chintalapani is an employee of Siemens Medical Solution USA, Inc.

Acknowledgments

This work is partially supported by a research contract from Siemens Medical Solution USA, Inc., to Michigan Technological University and the Research Excellent Fund from Michigan Technological University. The authors also want to thank their colleague Professor Charlie Strother, M.D. degree holder, from the University of Wisconsin-Madison for supplying some DSA data used for this study and Dr. Thomas Redel from Siemens Healthcare Inc. for technical support.

References

D. M. Sforza, C. M. Putman, and J. R. Cebral, “Hemodynamics of cerebral aneurysms,” Annual Review of Fluid Mechanics, vol. 41, pp. 91–107, 2009.
View at: Publisher Site | Google Scholar | Zentralblatt MATH
M. H. M. Vlak, A. Algra, R. Brandenburg, and G. J. E. Rinkel, “Prevalence of unruptured intracranial aneurysms, with emphasis on sex, age, comorbidity, country, and time period: a systematic review and meta-analysis,” The Lancet Neurology, vol. 10, no. 7, pp. 626–636, 2011.
View at: Publisher Site | Google Scholar
J. B. Bederson, I. A. Awad, D. O. Wiebers et al., “Recommendations for the management of patients with unruptured intracranial aneurysms: a statement for healthcare professionals from the Stroke Council of the American Heart Association,” Stroke, vol. 31, no. 11, pp. 2742–2750, 2000.
View at: Publisher Site | Google Scholar
T. Steiner, S. Juvela, A. Unterberg, C. Jung, M. Forsting, and G. Rinkel, “European Stroke Organization guidelines for the management of intracranial aneurysms and subarachnoid haemorrhage,” Cerebrovascular Diseases, vol. 35, no. 2, pp. 93–112, 2013.
View at: Publisher Site | Google Scholar
B. G. Thompson, R. D. Brown Jr., S. Amin-Hanjani et al., “Guidelines for the management of patients with unruptured intracranial aneurysms: a guideline for healthcare professionals from the American Heart Association/American Stroke Association,” Stroke, vol. 46, pp. 2368–2400, 2015.
View at: Google Scholar
S. Dhar, M. Tremmel, J. Mocco et al., “Morphology parameters for intracranial aneurysm rupture risk assessment,” Neurosurgery, vol. 63, no. 2, pp. 185–197, 2008.
View at: Publisher Site | Google Scholar
B. Weir, “Unruptured intracranial aneurysms: a review,” Journal of Neurosurgery, vol. 96, no. 1, pp. 3–42, 2002.
View at: Publisher Site | Google Scholar
J. L. Brisman, J. K. Song, and D. W. Newell, “Cerebral aneurysms,” The New England Journal of Medicine, vol. 355, no. 9, pp. 928–939, 2006.
View at: Publisher Site | Google Scholar
G. Benndorf, E. Wellnhofer, W. Lanksch, and R. Felix, “Intraaneurysmal flow: evaluation with Doppler guidewires,” American Journal of Neuroradiology, vol. 17, no. 7, pp. 1333–1337, 1996.
View at: Google Scholar
S. Meckel, A. F. Stalder, F. Santini et al., “In vivo visualization and analysis of 3-D hemodynamics in cerebral aneurysms with flow-sensitized 4-D MR imaging at 3 T,” Neuroradiology, vol. 50, no. 6, pp. 473–484, 2008.
View at: Publisher Site | Google Scholar
D. A. Steinman, “Image-based computational fluid dynamics modeling in realistic arterial geometries,” Annals of Biomedical Engineering, vol. 30, no. 4, pp. 483–497, 2002.
View at: Publisher Site | Google Scholar
J. Xiang, S. K. Natarajan, M. Tremmel et al., “Hemodynamic-morphologic discriminants for intracranial aneurysm rupture,” Stroke, vol. 42, no. 1, pp. 144–152, 2011.
View at: Publisher Site | Google Scholar
J. R. Cebral, F. Mut, J. Weir, and C. M. Putman, “Association of hemodynamic characteristics and cerebral aneurysm rupture,” American Journal of Neuroradiology, vol. 32, no. 2, pp. 264–270, 2011.
View at: Publisher Site | Google Scholar
H. Baek, M. V. Jayaraman, P. D. Richardson, and G. E. Karniadakis, “Flow instability and wall shear stress variation in intracranial aneurysms,” Journal of the Royal Society Interface, vol. 7, no. 47, pp. 967–988, 2010.
View at: Publisher Site | Google Scholar
D. A. Steinman, J. S. Milner, C. J. Norley, S. P. Lownie, and D. W. Holdsworth, “Image-based computational simulation of flow dynamics in a giant intracranial aneurysm,” American Journal of Neuroradiology, vol. 24, no. 4, pp. 559–566, 2003.
View at: Google Scholar
G. Byrne, F. Mut, and J. Cebral, “Quantifying the large-scale hemodynamics of intracranial aneurysms,” American Journal of Neuroradiology, vol. 35, no. 2, pp. 333–338, 2014.
View at: Publisher Site | Google Scholar
Y. C. Liang, H. P. Lee, S. P. Lim, W. Z. Lin, K. H. Lee, and C. G. Wu, “Proper orthogonal decomposition and its applications—part I: theory,” Journal of Sound and Vibration, vol. 252, no. 3, pp. 527–544, 2002.
View at: Publisher Site | Google Scholar | MathSciNet
B. Kohler, R. Gasteiger, U. Preim, H. Theisel, M. Gutberlet, and B. Preim, “Semi-automatic vortex extraction in 4D PC-MRI cardiac blood flow data using line predicates,” IEEE Transactions on Visualization and Computer Graphics, vol. 19, no. 12, pp. 2773–2782, 2013.
View at: Publisher Site | Google Scholar
J. Jeong and F. Hussain, “On the identification of a vortex,” Journal of Fluid Mechanics, vol. 285, pp. 69–94, 1995.
View at: Publisher Site | Google Scholar | MathSciNet
J. C. Hunt, A. Wray, and P. Moin, “Eddies, streams, and convergence zones in turbulent flows,” in Proceedings of the Studying Turbulence Using Numerical Simulation Databases, pp. 193–208, December 1988.
View at: Google Scholar
W. E. Lorensen and H. E. Cline, “Marching cubes: a high resolution 3d surface construction algorithm,” Computer Graphics, vol. 21, no. 4, pp. 163–169, 1987.
View at: Publisher Site | Google Scholar
J. Poethke, A. Spuler, C. Petz et al., “Cerebral aneurysm hemodynamics and a length of parent vessel,” in Proceedings of the World Congress on Medical Physics and Biomedical Engineering, O. Dössel and W. C. Schlegel, Eds., vol. 25/4 of Image Processing, Biosignal Processing, Modelling and Simulation, Biomechanics, pp. 1608–1611, Munich, Germany, September 2009.
View at: Google Scholar
M. Piccinelli, A. Veneziani, D. A. Steinman, A. Remuzzi, and L. Antiga, “A framework for geometric analysis of vascular structures: application to cerebral aneurysms,” IEEE Transactions on Medical Imaging, vol. 28, no. 8, pp. 1141–1155, 2009.
View at: Publisher Site | Google Scholar
H. Si and K. Gärtner, “3D boundary recovery by constrained Delaunay tetrahedralization,” International Journal for Numerical Methods in Engineering, vol. 85, no. 11, pp. 1341–1364, 2011.
View at: Publisher Site | Google Scholar | MathSciNet
H. K. Versteeg and W. Malalasekera, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Longman Scientific & Technical, Harlow, UK; John Wiley & Sons, New York, NY, USA, 1995.
J. Jiang, K. Johnson, K. Valen-Sendstad, K.-A. Mardal, O. Wieben, and C. Strother, “Flow characteristics in a canine aneurysm model: a comparison of 4D accelerated phase-contrast MR measurements and computational fluid dynamics simulations,” Medical Physics, vol. 38, no. 11, pp. 6300–6312, 2011.
View at: Publisher Site | Google Scholar
J. Jiang, C. Strother, K. Johnson et al., “Comparison of blood velocity measurements between ultrasound Doppler and accelerated phase-contrast MR angiography in small arteries with disturbed flow,” Physics in Medicine and Biology, vol. 56, no. 6, pp. 1755–1773, 2011.
View at: Publisher Site | Google Scholar
D. D'Humieres, I. Ginzburg, M. Krafczyk, P. Lallemand, and L. S. Luo, “Multiple-relaxation-time lattice Boltzmann models in three dimensions,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 360, no. 1792, pp. 437–451, 2002.
View at: Publisher Site | Google Scholar
M. N. Gwilliam, N. Hoggard, D. Capener et al., “MR derived volumetric flow rate waveforms at locations within the common carotid, internal carotid, and basilar arteries,” Journal of Cerebral Blood Flow and Metabolism, vol. 29, no. 12, pp. 1975–1982, 2009.
View at: Publisher Site | Google Scholar
D. R. Enzmann, M. R. Ross, M. P. Marks, and N. J. Pelc, “Blood flow in major cerebral arteries measured by phase-contrast cine MR,” American Journal of Neuroradiology, vol. 15, no. 1, pp. 123–129, 1994.
View at: Google Scholar
M. Zhao, S. Amin-Hanjani, S. Ruland, A. P. Curcio, L. Ostergren, and F. T. Charbel, “Regional cerebral blood flow using quantitative MR angiography,” American Journal of Neuroradiology, vol. 28, no. 8, pp. 1470–1473, 2007.
View at: Publisher Site | Google Scholar
J. Jiang and C. M. Strother, “Interactive decomposition and mapping of saccular cerebral aneurysms using harmonic functions: its first application with ‘patient-specific’ computational fluid dynamics (CFD) simulations,” IEEE Transactions on Medical Imaging, vol. 32, no. 2, pp. 153–164, 2013.
View at: Publisher Site | Google Scholar
D. G. Altman and J. M. Bland, “Measurement in medicine: the analysis of method comparison studies,” Journal of the Royal Statistical Society. Series D (The Statistician), vol. 32, pp. 307–317, 1983.
View at: Google Scholar
M. Jiang, R. Machiraju, and D. Thompson, “Detection and visualization of vortices,” The Visualization Handbook, pp. 295–309, 2005.
View at: Publisher Site | Google Scholar
G. Haller, “An objective definition of a vortex,” Journal of Fluid Mechanics, vol. 525, pp. 1–26, 2005.
View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
W. E. Lorensen and H. E. Cline, “Marching cubes: a high resolution 3D surface construction algorithm,” ACM Siggraph Computer Graphics, vol. 21, no. 4, pp. 163–169, 1987.
View at: Publisher Site | Google Scholar
X. He and D. N. Ku, “Pulsatile flow in the human left coronary artery bifurcation: average conditions,” Journal of Biomechanical Engineering, vol. 118, no. 1, pp. 74–82, 1996.
View at: Publisher Site | Google Scholar
R. Yasuda, C. M. Strother, W. Taki et al., “Aneurysm volume-to-ostium area ratio: a parameter useful for discriminating the rupture status of intracranial aneurysms,” Neurosurgery, vol. 68, no. 2, pp. 310–317, 2011.
View at: Publisher Site | Google Scholar
B. Chung and J. R. Cebral, “CFD for evaluation and treatment planning of aneurysms: review of proposed clinical uses and their challenges,” Annals of Biomedical Engineering, vol. 43, no. 1, pp. 122–138, 2015.
View at: Publisher Site | Google Scholar
M. Wernick, Y. Yang, J. Brankov, G. Yourganov, and S. Strother, “Machine learning in medical imaging,” IEEE Signal Processing Magazine, vol. 27, no. 4, pp. 25–38, 2010.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2016 Kevin Sunderland et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1568

Downloads

1612

Citations