Determination of Significant 3D Hemodynamic Features for Post-Embolization Recanalization in Cerebral Aneurysms through Explainable Artificial Intelligence

1. Introduction

Cerebral aneurysms are abnormal focal outpouchings of cerebral arteries that are associated with significant morbidity and mortality [1]. They have a prevalence of 1-5% in the adult population, and rupture of an aneurysm can result in subarachnoid hemorrhage, a type of hemorrhagic stroke with a high mortality rate [2,3]. There are three options for treating intracranial aneurysms: observation, craniotomy with clip ligation (clipping), and endovascular occlusion using detachable coils (coiling) [1]. Coil embolization has been increasingly used due to its reduced physical burden on patients compared to clipping surgery [4]. Notwithstanding this advantage, recanalization, or the reopening of an aneurysm after coiling, is a significant challenge [5,6,7]. Therefore, the establishment of predictive models for the post-embolization recanalization in cerebral aneurysm is crucial for physicians in terms of surgery planning, decision-making and postoperative management.

Statistical analysis has shown that various morphological and hemodynamic features are significantly associated with recanalization in cerebral aneurysms. Morphological parameters, such as maximum size, neck width, height, area of aneurysm inlet, neck, posterior communicating (Pcom), as well as dimensionless parameters, such as aspect ratio, bottleneck, size ratio, area ratio and the ratio of the aneurysm volume to coil volume (VER), were identified as risk factors for recanalization [8,9,10,11,12,13,14]. Hemodynamic parameters acquired from computational fluid dynamics (CFD), such as velocity, pressure, wall shear stress (WSS) and shear rate, along with their related dimensionless parameters, including pressure difference (PD), aneurysmal residual flow volume and maximal force with porous media as coil, were also considered important risk factors [5,6,7,15,16,17,18,19,20,21]. In recent developments, machine learning (ML) models emerged as important tools for the prediction of cerebral aneurysm-related diseases. ML models, encompassing random forest, k-nearest neighbor, and support vector machines, were used to predict treatment outcomes with flow diverters, relying on predetermined hemodynamic parameters derived from CFD simulations [22]. However, both statistical and traditional ML models were notably constrained by their reliance on pre-defined input parameters. These parameters may be overlooked in the initial feature selection phase due to lack understanding of the pathology, or they may be difficult to determine accurately due to the complexity of in-vivo fluid dynamics [23,24]. Additionally, these parameters were typically spatiotemporally maximized or averaged zero-dimensional (0D) extracted from 3D or 4D CFD results (time plus 3D coordinates).

Recent advances in deep learning (DL) techniques were applied to the automated extraction of 3D morphological and hemodynamic features from CFD simulation results, albeit with certain limitations. The PointNet architecture was adopted upon its suitability for processing data from CFD simulation expressed as 3D point clouds, where each data point was described by its spatial coordinates (x, y, z) coupled with relevant hemodynamic attributes [25]. One instance developed an integrated model, employing PointNet to extract features for subsequent classification of aneurysm rupture status with traditional ML algorithm. Notably, impingement zones were predetermined in CFD simulation results for feature extraction [26]. An alternative instance used an End-to-End PointNet-based deep neural network (DNN) for the purpose of predicting aneurysm recanalization. This model extracted and determined automatically the key morphological and hemodynamic features, without the need for prior assumption or data processing. A notable drawback of this approach was the "black-box problem," which hindered its acceptance among medical experts and healthcare practitioners, who required transparency in model decision-making processes for making informed clinical inferences [27].

Recently, the technique of explainable artificial intelligence (XAI) garnered substantial attention in the field of medical diagnosis, providing insights into diagnoses and uncovering previously unknown information [28]. XAI techniques were distinguished based on three criteria [29,30]:

• Model-based versus post hoc explanation: Model-based explanation entails the utilization of intelligible yet sufficiently sophisticated models, exemplified by linear regression, which effectively capture and elucidate the relationships between input and output variables. Conversely, post hoc explanation aims at analyzing a trained model to achieve insight into learned relationships. Saliency mapping, also called visual explanation, is the most common post hoc approach in medical image analysis [31].
• Model-specific versus model-agnostic explanation: Model-specific explanation techniques are constrained to specific model categories, whereas model-agnostic explanation methods operate independently of the neural network's architectural choice, focusing solely on the input and output of the neural network.
• Local versus global explanation: Global explanation offers overarching insights by delivering general relationships. In contrast, local explanation focuses on elucidating the rationale behind individual inputs, providing a detailed explanation for the certain dataset.

The objective of this study was twofold: first, to explore previously unknown hemodynamic features associated with aneurysm recanalization, and second, to develop a robust and interpretable predictive model. To achieve this, hemodynamic data were collected for 59 aneurysms through CFD simulations, and subsequently analyzed employing two distinct methodological approaches. Velocity, pressure and WSS simulated in both pre- and post-coiling models were analyzed. This involved the utilization of a PointNet-based DNN trained with the point-cloud-formatted data in 3D. Additionally, the conventional statistical analysis approach was applied, focusing on spatially averaged scalars in 0D. By comparing the results from two approaches, the significant parameters of post-embolization recanalization were determined. The 3D features of these parameters were further analyzed using saliency mapping, the post hoc, model-specific, and localized XAI technique.

2. Materials and Methods

2.1. Patients

Recanalized aneurysms were defined as aneurysms with an increase in the Raymond–Roy grade or enlargement within a grade resulting in retreatment, whereas the stable type was defined as aneurysms that showed no change or enlargement within a grade with no retreatment [32].

Clinical and radiological data of consecutive patients who underwent endovascular treatment for aneurysms between January 2007 and December 2020 were investigated. The inclusion criteria of this study were as follows: (1) saccular aneurysms; (2) coiling performed without stent placement; (3) class I or II initial aneurysm embolization results according to the Raymond–Roy classification; (4) volume embolization ratio (VER), the ratio of the aneurysm volume to the coil volume, over 20%; (5) a follow-up period of at least 1 year after endovascular coiling; and (6) pretreatment 3-dimensional rotational angiography (3D-RA) images available and of good quality for use in computational modeling [33].

A total of 58 intracranial aneurysms from 57 patients were included in the final analysis, including nine recanalized and 49 stable aneurysms. The included patients were regularly followed up with magnetic resonance angiography every six months after treatment. Digital subtraction angiography was performed when recanalization was suspected on magnetic resonance imaging and additional endovascular treatment was administered if needed.

2.2. Ethics approval

The protocol for this retrospective study was approved by the ethics committee of Kanazawa University (No. 1781). Written informed consent was waived due to the retrospective collection of clinical data; however, all patients had the right to opt-out of the study at any time.

2.3. Workflow

The workflow of this study was shown in Figure 1. Velocity, pressure, and wall shear stress (WSS) in both pre- and post-coiling models were derived from computational fluid dynamic (CFD) simulations for 59 aneurysms. Two approaches were utilized to identify the significant parameters for aneurysm recanalization. First, a PointNet-based DNN was used leveraging their inherent 3D features. Second, the conventional statistical analysis was conducted analyzing the spatially averaged parameters in 0D. Further 3D feature analysis was conducted utilizing the visual saliency mapping proposed in [34]. It assigns each point a saliency score to reflect its contribution to the output. A higher (positive) score was expected to indicate a more (positive) contribution. The scores were normalized to the range between 0 and 1 in this study.

2.4. CFD simulation

Digital Imaging for Communication in Medicine format was utilized for 3D-RA (Philips Healthcare, Best, Netherlands). Blood vessels were extracted and converted into standard triangulated language data via Amira (version 5.6, Maxnet Co, Ltd, Tokyo, Japan). Inlet and outlet extensions were performed to generate fully developed flow profiles. Hemodynamics in the pre-coiling and virtual post-coiling models were considered in this study. In the pre-coiling model, we defined the neck plane as the flat plane that divided the aneurysm and branching vessel (posterior communicating artery) or parent artery (internal carotid artery). Based on the pre-coiling model, the virtual post-coiling was created by removing the aneurysm dome, and we defined the coil plane as the surface of the coil mass in the same plane as the neck plane in the pre-coiling model. The inlet plane was defined as the section plane located 1 mm proximal to the aneurysm. For meshing, the 3D images were imported into ANSYS ICEM CFD software (version 16.2, ANSYS Inc., Canonsburg, Pennsylvania). To enhance the analytical precision of the boundary layer, seven prism element layers were created at the wall surface. We performed a steady-state simulation with ANSYS CFX (version 16.2, ANSYS Inc., Canonsburg, Pennsylvania) to acquire the patient-specific hemodynamic data. Blood was modeled as a Newtonian fluid with a density of 1100 kg/m3 and a viscosity of 0.0036 Pa·s. A no-slip boundary condition was applied to the rigid vessel walls. The outflow boundary condition was set as 0 Pa [5,6,7,27]. The inlet boundary condition was a mass flow rate of 0.003465 kg/s, which was the value at diastole end of pulsatile profile from literature [35].

2.5. Parameter deviation for statistical analysis

Morphological parameters were measured using 3D-RA, including the maximum size, neck width, height, area of aneurysm inlet, neck and posterior communicating (Pcom) [5,8]. Aspect ratio was defined as the ratio of the maximum perpendicular height to the neck diameter [9,10]. Bottleneck factor was defined as the ratio of the dome width to the neck diameter [11]. Size ratio was defined as the ratio of the maximum aneurysm height to the parent vessel diameter [12]. Area ratio was defined as the area ratio of the aneurysm neck to the aneurysm inlet in parent artery [13]. VER was the ratio of the aneurysm volume to the coil volume [14].

Dimensionless hemodynamic parameters were utilized for analyzing, building a model independent on patient-specific inflow rate [5,6,7,21]. Hemodynamic parameters in pre-coiling model included velocity in aneurysm dome (volvel), WSS, static pressure (P) and dynamic pressure (Pdyn) at aneurysm neck plane, which were normalized using surface-averaged velocity, WSS and P at aneurysm inlet in parent artery. Inflow rate ratio (FR) was defined as the ratio of the inflow rate at the neck plane to the flow rate at the inlet plane [5]. In the post-coiling model, PP was defined as pressure change between aneurysm neck and inlet, whereas pressure difference (PD) was defined as the ratio of PP to Pdyn at aneurysm inlet. All hemodynamic parameters were spatially averaged.

2.6. Statistical analysis

All continuous parameter values were subsequently expressed as mean ± standard deviation. Continuous variables were analyzed using the Levene test to determine the normality. The Levene test was employed to examine the equality of variances between two groups with continuous variables. Differences in variables between the two groups were tested using Welch’s t-test for data with unequal variance or t-test for data with equal variance. For categorical variables, a chi-square test was used to examine the significant differences between the groups. A P-value less than 0.05 was considered statistically significant. Statistical analyses were performed using Scipy (version 1.9.3) [36].

2.7. Data preparation and augmentation for DNN

To train the DNN, 3D morphological data including x, y, z coordinates and the relevant hemodynamic data were exported from CFD simulation results. The hemodynamic data included velocity, pressure and WSS in pre-coiling (pre_vel, pre_P, pre_WSS) and post-coiling model (post_vel, post_P, post_WSS) as shown in Figure 2. For all the parameters, only data at the aneurysm part was considered. Velocity and pressure data were exported from the internal field of the aneurysm, whereas WSS data were obtained from surface nodes. Voxelization with a voxel size of 0.5mm was employed on velocity and pressure data to ensure uniform distribution in the dataset. Without this process, data points tended to cluster near the wall due to the presence of mesh wall layers in CFD simulations.

The DL algorithms have trouble learning when some classes dominate others [37]. The minority class within the training dataset was initially augmented through oversampling, involving the replication of instances. Subsequently, online augmentation techniques were applied during the training process. The morphological point cloud was augmented by randomly rotating the object along x, y, z axis in the range of [−180°, 180°] and by jittering its position by a Gaussian noise with a mean of 0 and standard deviation of 0.01 [25]. The hemodynamic data was augmented by adding a random type of noise (Gaussian or white noise) with a random level in the range of [0, 0.025] [27].

2.8. DNN model architecture and training

This study utilized a PointNet-based DNN to extract both 3D morphological and hemodynamic features to predict post-embolization recanalization in cerebral aneurysm [25,27]. It consumed the 3D point clouds with two feature channels (morphology and hemodynamics) as inputs and predicted recanalization probabilities as outputs. The DNN consisted of two key modules: feature extraction and classification. In the first module, we approximated the model by a convolutional network; then its output was aggregated in a latent global feature with a vector of size 1024 by the symmetric function, herein, max pooling was utilized. Based on the features from the first module, the probability of aneurysm recanalization was decoded in the second module. In this stage, multi-layer perceptron (MLP) mapped the output to probability of K classes, herein, 2 events (recanalized or stable) for each patient. The detailed description could be found in our previous work [27].

Before training, the input point clouds were normalized to within [0, 1] with each parameter’s maximum and minimum value. A batch size of 10 and batch-norm was used in this study. Dropout layers were utilized for the MLP. Stochastic gradient descent optimization was used to update the weights and biases, herein, adaptive moment estimation (ADAM). During training, an eight-fold cross validation and repetition of 2000 was used. Parameters of the best epoch with lowest loss on validation dataset were recorded. An adaptive learning rate with an exponential decay was adopted with an initial learning rate of 0.0001, decay step of 200000 and decay rate of 0.7. The technique was implemented using the TensorFlow DNN library and executed on a PC cluster with 32 cores of Intel(R) Xeon(R) Gold 6242 CPU @ 2.80GHz and 1.5T RAM in total. The operating system is CentOS 7.4.1708.

2.9. Evaluation of DNN model

The receiver-operating characteristic (ROC) curves were constructed after training and the area under the curves (AUROC) values were calculated to evaluate the predictive performance of DNN model. Furthermore, to address the imbalanced nature of the dataset, in which only 18% of the aneurysms exhibited recanalization, we utilized the area under the precision recall curve (AUPRC) together with AUROC as a metric for performance evaluation. This method is specially designed to detect rare events and are appropriate in these scenarios as they particularly show a classifier having a low performance if it is misclassifying most or all the minority class [38]. Results were further plotted in confusion matrix to analyze its performance for each class. Confusion matrix plots the number of correct and incorrect model-predicted outcomes (horizontal axis) against the actual outcomes (vertical axis). The performance examination was performed with Scikit-learn (version 1.0.2) [39].

2.10. Data and code availability

The datasets generated for this study are protected patient information. The code is available at https://github.com/saradeuJP/XAI_PointNet_CA.git.

3. Results

3.1. Patient characteristics

Characteristics of patients in each group are listed in Table 1. There was no significant difference in age, sex, rupture status and aneurysm locations between two groups. The recanalized group showed higher mean values for all morphological parameters except area ratio in comparison to the stable group. Max size, Height and Area Neck were significantly larger in the recanalized group. Within the hemodynamic features, volvel and PD demonstrated significant differences between the two groups, whereas P, PP, Pdyn, WSS and FR did not show significant differences.

3.2. Model training and performance

Figure 3 and Figure 4 present ROC and PRC analyses on the training and testing datasets, enabling a performance comparison of DNN models trained with 3D velocity, pressure, and WSS in both pre- and post-coiling models. On the training dataset, the model trained with velocity data from the post-coiling model demonstrated the highest performance with AUROC/AUPRC values of 1.000/1.000, a statistically significant improvement over the other five models (P < 0.05). The model trained with velocity data from the pre-coiling model achieved also commendable results with AUROC/AUPRC values of 0.991/0.992. Especially, it displayed its superior performance on the testing dataset, as detailed in Table 2, distinguishing the two recanalized aneurysms, as highlighted by confusion matrix in Fig. 5. In contrast, DNN models trained with pressure and WSS data exhibited less satisfactory performance on both training and testing dataset. Consequently, the subsequent section focused on elucidating the explainability of the DNN model trained with 3D velocity data.

3.3. Model explainability

Figure 6 depicts saliency maps illustrating the top ten salient points within seven recanalized aneurysms from the training dataset in both pre- and post-coiling models. Within the pre-coiling models, the recanalized cases revealed robust injection flow zones with strong impingement impact points at walls. These points were observed within the aneurysms for cases (a), (c), (e), and (g). For cases (d) and (f), these impact points manifested in secondary aneurysms, whereas case (b) exhibited them at the aneurysm inlet. Within the post-coiling models, our observations revealed the presence of three interacting modes in recanalized aneurysms. These modes encompassed the presence of forceful impact points in cases (a), (c), (d), (e), (f), and (g), revealing the salient points on the aneurysm neck plane, as well as the existence of the Pcom within the flow direction in cases (a), (c) and (e), with salient points observed in this artery. Additionally, substantial contact areas in the direction parallel to the flow were also considered as important features, with salient points found at both aneurysm in- and outlet in case (b). In contrast, these influential effects were conspicuously absent in stable cases shown in Figure 7 within both pre- and post-coiling models.

4. Discussion

In this study, we identified the significant hemodynamic parameter using a DNN with PointNet architecture. Its important 3D features that significantly contributed to the DNN's decision-making process were further analyzed through XAI technique. Firstly, significant hemodynamic parameter, that were automatically determined by DNN, were in alignment with those identified through the conventional statistical analysis. DNN models trained with velocity data demonstrated superior performance, as evidenced by ROC and PRC curves, in both training and testing datasets compared to models trained with pressure and WSS data. The statistical analysis had the same conclusions, where velocity within the aneurysm dome in the pre-coiling model (volvel) yielded the lowest P-value, whereas pressure and WSS exhibited no significant differentiation between the two groups. Secondly, the saliency map visualization technique hinted at the existence of hitherto unknown 3D hemodynamic features associated with recanalization. In the pre-coiling models, it identified strong injection flow zones with impingement impact points in recanalized aneurysms as the important 3D features to classify aneurysm recanalization. While this represents a novel finding in the context of aneurysm recanalization, it aligns with prior research in aneurysm development and rupture emphasizing the significance of injection and impingement patterns [18,40]. In post-coiling models, the saliency map illuminated a pronounced interplay between the flow and coiling plane. The XAI technique effectively identified potential factors contributing to this interaction, including impingement points on the coiling plane, the presence of Pcom in the direction of the flow, and substantial contact areas. This discovery, to the best of our knowledge, marks the first report of such phenomena. However, it was noteworthy that the 0D parameters utilized in this study provided an incomplete representation of the complex 3D flow features. Future research endeavors will encompass a broader spectrum of 0D parameters, including flow jet angle and dispersion index, to corroborate the consistency between DNN and statistical analyses, as well as the features observed from XAI visualization [44]. Additionally, DNN models trained with velocity data from both pre- and post-coiling models exhibited no significant differences on the testing dataset, underscoring the necessity for a larger and external patient dataset for future analyses aimed at selecting the most effective model.

The complexity and inherent instability of flow patterns in cerebral aneurysm ranged from simplistic configurations with a singular recirculation region to intricate arrangements featuring multiple vortical structures that could be stable, dynamic, or intermittent throughout the cardiac cycle. These patterns were not solely contingent upon the size and shape of the aneurysm but were profoundly influenced by the manner in which blood flowed from the parent vessel into the aneurysm, a phenomenon intricately linked to the geometry of the parent artery [19]. One widely employed parameter for quantitatively assessing the impact of inflow impingement on the aneurysm wall was WSS. Locally elevated WSS characterized concentrated inflow jets, whereas a more uniform WSS distribution with lower magnitudes delineated areas of stagnation within the aneurysm [41]. However, the complex interplay between morphology and hemodynamics engendered a contentious debate regarding the role of high versus low WSS in cerebral aneurysm pathophysiology [42,43]. To further describe the inflow and impingement pattern, one categorized inject and impingement sizes manually, typically as small or large, the other determined the quantitative size through predefined thresholds such as exceeding 80% of the maximum WSS within the aneurysm sac [15,26]. In contrast, our approach leveraged automatic feature extraction and determination from raw 3D flow pattern data. This approach eliminated the constraints associated with manual preprocessing and prior assumptions or knowledge. Additionally, the flexibility of the DNN model allows for the incorporation of various 3D hemodynamic data inputs, thereby addressing the confounding factors contributing to recanalization. Furthermore, this framework can be readily extended to accommodate 4D input data, incorporating the temporal dimension to capture the temporal instability of flow patterns.

This study did not consider the coiling configuration. The coil surface was represented by a flat and rigid plane; whereas in clinical practice, the coil surface is typically rough and allows for blood penetration into the coil mass. Further investigation is needed to integrate the realistic coil surface obtained through advanced techniques such as silent magnetic resonance angiography into the analysis [45].

Other than saliency map, we also experimented with GradCAM heatmaps but found that the resolution of these is too coarse to capture the fine details that are important to understand the models [46]. For 1024 input points in this work, the GradCAM heatmap would be 170x coarser resolution than the input. In comparison, saliency map have the same resolution as the input image and allow us to examine the predictions at the pointwise level [47]. Other XAI methods are also planned to be used in the future to deepen our understanding of recanalization pathology [29,30].

There are also other limitations to this study. First, the analysis encompassed a total of 58 aneurysms in this study. To further evaluate the model, a larger cohort of cases will be analyzed as a future step. Second, this was a retrospective study. To validate the effectiveness of developed models as a predictive tool, a prospective study with a larger cohort involving all types of aneurysms is required. Third, the fluid-solid interaction between the vessel walls and blood was not considered. Finally, the boundary conditions were uniform for all the patients. To enhance accuracy, future investigations will explore the utilization of patient-specific boundary conditions derived from 4D-Flow MRI or Transcranial Doppler Ultrasonography, if feasible.

5. Conclusions

PointNet-based DNN model trained utilizing 3D velocity data exhibited significant potential in prognosticating post-embolization aneurysm recanalization. The integration of XAI methodologies provided insight into critical 3D features that played a pivotal role in the decision-making framework, holding the potential to unearth hitherto unexplored 3D hemodynamic characteristics in the future.