Quantitative analysis of flow vortices: differentiation of unruptured and ruptured medium-sized middle cerebral artery aneurysms

Abstract

Background

Surgical intervention for unruptured intracranial aneurysms (IAs) carries inherent health risks. The analysis of “patient-specific” IA geometric and computational fluid dynamics (CFD) simulated wall shear stress (WSS) data has been investigated to differentiate IAs at high and low risk of rupture to help clinical decision making. Yet, outcomes vary among studies, suggesting that novel analysis could improve rupture characterization. The authors describe a CFD analytic method to assess spatiotemporal characteristics of swirling flow vortices within IAs to improve characterization.

Methods

CFD simulations were performed for 47 subjects harboring one medium-sized (4–10 mm) middle cerebral artery (MCA) aneurysm with available 3D digital subtraction angiography data. Alongside conventional indices, quantified IA flow vortex spatiotemporal characteristics were applied during statistical characterization. Statistical supervised machine learning using a support vector machine (SVM) method was run with cross-validation (100 iterations) to assess flow vortex-based metrics’ strength toward rupture characterization.

Results

Relying solely on vortex indices for statistical characterization underperformed compared with established geometric characteristics (total accuracy of 0.77 vs 0.80) yet showed improvements over wall shear stress models (0.74). However, the application of vortex spatiotemporal characteristics into the combined geometric and wall shear stress parameters augmented model strength for assessing the rupture status of middle cerebral artery aneurysms (0.85).

Conclusions

This preliminary study suggests that the spatiotemporal characteristics of flow vortices within MCA aneurysms are of value to improve the differentiation of ruptured aneurysms from unruptured ones.

Appendix

A computational method based on informational entropy determined the spatially varying direction of the velocity field as a means to identify flow vortices. First, the 3D angular space of the flow velocity field was divided into 360 bins of equal area [25], resulting in cones connecting a unit sphere center and surface patches. Velocity vectors were then designated to a patch when said vector was contained within the connected cone. Assessing the local flow direction (x ϵ (x₁, x₂, x₃,…x_n)) of the velocity field (X), the probability p(x_i) of flow direction can be used to calculate Shannon’s entropy.

$$ {\displaystyle \begin{array}{c}\boldsymbol{H}\left(\boldsymbol{X}\right)=-\sum \limits_{{\boldsymbol{x}}_{\boldsymbol{i}}\in \boldsymbol{X}}\boldsymbol{p}\left({\boldsymbol{x}}_{\boldsymbol{i}}\right){\mathbf{\log}}_{\mathbf{2}}\boldsymbol{p}\left({\boldsymbol{x}}_{\boldsymbol{i}}\right)\\ {}\\ {}\boldsymbol{NE}\left(\boldsymbol{X}\right)=\frac{\boldsymbol{H}\left(\boldsymbol{X}\right)}{{\mathbf{\log}}_{\mathbf{2}}\left(\boldsymbol{N}\right)}\end{array}} $$

(1)

Entropy is seen as minimal as the velocity directions concentrate toward one bin value and tend toward the maximum as the probability of velocity vectors in all directions becomes equally likely. Dividing H(X) by the maximum possible Shannon’s entropy (log2(N)) normalizes the entropy (NE) values between 0 and 1. For this work, a fixed volume of interest (VOI) was selected around each voxel within an IA: N_x × N_y × N_z; N_x = N_y = N_z = 11. Velocity vectors within each VOI, centered at each ith voxel, were used to determine the normalized entropy of each point in the flow field.

Reliance on normalized Shannon’s entropy to identify areas of flow vortices carries with it an inherent limitation: the inability to distinguish rotational flow vortices from random Brownian flow. From the perspective of Shannon’s entropy, both Brownian and rotational flow fields have a high probability that their velocity vectors will orient in any direction. The adaptation of aspects of the λ₂ method of critical point analysis improves the NE’s ability to assess areas of vortices by identifying the central-most region (critical point) of swirling flow, with Brownian flow unlikely to have a centralized critical point [15]. The λ₂ method decomposes the $\mathbf{\nabla}\overrightarrow{\boldsymbol{v}}$ tensor into its strain-rate (symmetric: S) and spin (asymmetric: Ω) tensor:

$$ \boldsymbol{S}=\frac{\mathbf{\nabla}\overrightarrow{\boldsymbol{v}}+\mathbf{\nabla}{\overrightarrow{\boldsymbol{v}}}^{\boldsymbol{T}}}{\mathbf{2}},\boldsymbol{\varOmega} =\frac{\mathbf{\nabla}\overrightarrow{\boldsymbol{v}}-\mathbf{\nabla}{\overrightarrow{\boldsymbol{v}}}^{\boldsymbol{T}}}{\mathbf{2}}. $$

(2)

The critical point of swirling flow is identified where S² + Ω² has two negative eigenvalues and λ₁ > λ₂ > λ₃. When a proper point is identified, the dot product of the normalized velocity vector and eigenvector is calculated, relating the angle between the velocity vector and the degree of change of said vector: 0 as co-directional and 1 as orthogonal vectors.

$$ du\left(\theta \right)=\left\Vert \overrightarrow{v}\right\Vert \cdotp {\uplambda}_{2^{\cdotp }} $$

(3)

The du(θ) value multiplied to NE to help distinguish swirling flow patterns from Brownian patterns: swirling patterns having lower values than Brownian patterns:

$$ CM(X)= NE(X)\times du\left(\theta \right). $$

(4)

This combined method (CM) was used across all voxels within an IA dome for each time step of the final simulated cardiac cycle.

Once the CM was applied to all voxels within the 3D velocity field, the classic marching cube algorithm extracted the regions of vortex flow and mapped them to iso-surfaces [29]. To standardize the identification and extraction of vortex iso-surfaces, a set value > 0.3 of the combined method was used as the marching cube threshold value for all cases. To reduce the appearance of small, isolated areas of disturbed flow being mistaken for vortex flow, iso-surfaces with a volume < 0.5 mm³ were excluded from the analysis.