We included all patients operated for a ruptured intracranial aneurysm in whom the surgeon was able to identify the point of rupture on preoperative computed tomography (CT) angiography. Altogether, six aneurysms with a known rupture site were collected within 2 years. All patients s also needed to have a high-quality preoperative CT image available, and these acted as the bases for computational meshes. Their sizes varied from 8.83 to 12.60 mm, and they were located on the middle cerebral artery (5 cases) and right internal carotid artery bifurcation (1 case) as shown in the Table 1. Aneurysm size was calculated as the maximum distance within the aneurysm dome, the diameter of the neck was determined as the maximum distance in the neck plane, and the aspect ratio calculated as the ratio between the height (maximum perpendicular height of the aneurysm) and the diameter of the neck. In aneurysm F, the neck was not defined.

CT angiography was performed at admission using 256-detector multislice CT, with a pixel resolution of 0.6 mm. CT acquisition was synchronized with the contrast bolus throughout a large cranio-caudal, or z-axis, field of view. All CT angiographies were performed after IA rupture and therefore the shape of the aneurysm may have not reflected the pre-rupture conditions. The data were then used for segmentation and 3D visualization of the intracranial blood vessels using ITK-SNAP (11). A 3D model of the aneurysm sac as well as the major inlet and outlet blood vessels was constructed.

Voxel segmentations of patient-specific geometries of cerebral aneurysms obtained from the CT image slice were interactively constructed by the ITK-SNAP program. At least two medical doctors approved this first 3D reconstruction independently. After binary voxel segmentation, the surface meshes were constructed and smoothed using a combination of lowpass and Laplace smoothing filters as implemented in the iso2mesh package (12). The parent vessel was modeled as long as possible, at least six times longer than the inlet diameter of the vessel (13), to achieved fully developed flow. All inlet and outlet branches were extended to obtain cross-sections that were circular and almost perpendicular to the centerline of the vessel. All final surface meshes were approved by medical doctors.

Consequently, three tetrahedral meshes with various average edge lengths were constructed for all six aneurysm cases to test the mesh dependence. A mesh dependency test was performed for steady flow in each case with an input flux of 169 mL/min. Values of velocity magnitude, normal pressure, and wall shear stress were averaged over the aneurysm dome and compared among the three computations.

Navier–Stokes equations with constant dynamic viscosity μ = 3.5 mPa·s and constant blood density ρ = 1,050 kg·m⁻³ were considered to be the governing equations for blood flow. A fast and robust code was used for solving such equations by the finite element method, as described and evaluated by Chabiniok et al. (14). The program is implemented in Python using the FEniCS (15) and PETSc libraries (16). Time derivatives were approximated by the BDF scheme; the finite element method was used for spatial discretization of the velocity and pressure fields, in particular the MINI element was used (17).

The no-slip boundary condition was prescribed on the vessel wall and the zero surface-traction boundary condition on the outlets. The parabolic profile for the inlet velocity was prescribed on the circular inlet surface, with spatially averaged velocity magnitude dependent only on time. The waveform of the velocity magnitude can be characterized by the following three factors: maximum velocity value or peak systolic value (PSV); minimal velocity value or end diastolic value (EDV); and the time-averaged velocity (TAV) over the input plane. Since it was not possible to obtain patient-specific inlet velocity waveforms, the waveform was taken as an averaged curve over a set of patients such that the TAV over the input plane is fixed TAV = 0.3 m/s for all six cases. This resulted in different inlet velocity fluxes as the radii of inlet vessels were different in each case. The velocity and flux waveforms are shown in Figure 1, and the exact values of the flux and radii are collected in Table 2. For all six geometries, the heart rate was taken as 60 bpm.

The following hemodynamic quantities were considered based on Xiang et al. (3) and Axier et al. (18). Time-averaged wall shear stress (TAWSS) is defined as

where T denotes the time period of one cardiac cycle and WSS stands for wall shear stress vector. The oscillatory shear index (OSI) captures the directional changes of wall shear stress vector throughout one cardiac cycle. OSI is a non-dimensional quantity given by the relation

The low shear area (LSA) is the percentage ratio of the aneurysm area with the WSS less than 10% of the spatially averaged WSS and of the entire aneurysm area. Additionally, we focused on the following hemodynamic parameters: the peak wall shear stress (PWSS) and peak LSA—quantities evaluated at the time of the PSV, time-averaged LSA (TALSA)—integrated over one cardiac cycle, normalized TAWSS (NTAWSS), and normalized PWSS (NPWSS)—the percentage ratio between the TAWSS (resp. PWSS) over the dome and TAWSS (resp. PWSS) over the parent artery. All post-processing calculations were performed using Python and visualizations using ParaView (19).

The stereolithographic format (STL) of the aneurysm and connected arteries with patient-realistic dimensions was used to fabricate corresponding silicone models (Figure 2) using a 3D printer (Original Prusa i3MK2, Prusa Research, Czech Republic) and silicone (Sylgard 184, Dow Corning, United States). The lumen of the aneurysm and the corresponding arteries were 3D printed out of water-soluble material [butenediol vinyl alcohol (BVOH)]: verbatim BVOH filament soluble support (Verbatim/Mitsubishi Chemical Holdings Group, Japan). The printout with appropriate tubing was embedded in a silicone cast. The printout was then removed with distilled water. The dimensions of the model lumen were verified using an X-ray computed tomography (horizontal resolution 0.2 mm, vertical resolution 0.5 mm).

Each model was placed in an open-flow loop with steady flow. The inlet flux was fixed at 169 mL/min. All flow parameters are displayed in Table 3.

The inlet of a peristaltic pump (Gilson Miniplus 3, high flow rate head, 8 mm ID tubing for peristaltic loop) was placed in a reservoir. A bubble trap and a pulsation damper were located behind the pump. The damper was made of a “Y” tube connector with a 60 mL syringe connected through a 50 cm long tube followed by the inlet into a silicone model. All model outlets were joined, and the common outlet was placed into the reservoir (Figure 3). The internal diameter of the tubing was kept larger than that of the model inlets/outlets. As a circulating medium served 47% (Vol/Vol) glycerin in water containing 0.1% (Vol/Vol) Tween 20 and supplemented with a suspension of 30 μm polystyrene beads (Cat No. 95531, Merck, Germany). The density of the circulating medium was 1,125 kg m⁻³. To maintain viscosity of 3.75 mPa s the flow loop was placed in a thermostatic chamber (38°C, Okolab, Italy) at Axio Observer Z1 microscope by Zeiss (Germany). The microscope was equipped with a 2.5× Apo Plan objective (NA 0.06) and Orca 4.0 LT digital camera (Hamamatsu, Japan). The focal depth of the setup was 248 μm. The flow velocity was evaluated by means of particle traces length within the exposure time (1–10 ms) using image analysis in Image J (20). The distance of the measurement plane from the base was determined by calibrated vertical movement of the microscope stage.

Kinematic similarity between CFD simulations and PIV measurements is attained by maintaining the Reynolds number. In our case, it was achieved by setting consistent ratios between the fluid density and dynamic viscosity.

The measurement was performed in a plane parallel to the base of the model, the plane was set by the focus of the camera in the PIV experiment. The line segment was set at the beginning of the measurement so that one end of the line lied on the top of the aneurysm. The same plane and line segment in the plane were then found in ParaView using Geometric Shapes and CT scans of the printed models. Computed velocities were projected onto three different line segments with the same slope: one of them lay exactly in the measurement plane and the other two lines were shifted up and down by the distance of depth of field (DOF) of the microscope. Velocity projection was performed in Python using function sample_over_line from PyVista library (21). Moreover, traces of the massless particles were tracked using ParaView. We created randomly distributed particles in the aneurysm dome and tracked their path lines for 10 ms which corresponds to the exposition time of the camera.

Six patients with ruptured IAs were analyzed. All 6 patients were admitted for IA rupture and were operated on in an emergent setting (within 24 h after rupture). The point of rupture was identified during surgical exploration of the aneurysm before clipping and it was marked on the CTA 3D image after surgery. In one case (aneurysm A), a leak of blood from the aneurysm sac was clearly visible on preoperative CT angiography.

The results of the mesh sensitivity test are shown in Table 4. Steady-state quantities, namely the velocity, pressure and WSS, were averaged over the aneurysm dome using three meshes with different average edge lengths. Space-averaged velocity and pressure computed on three different geometries either remain constant or vary on the third significant digit. WSS values depend on mesh refinement in most cases. The finest meshes with an average edge length between 0.173 mm and 0.244 mm were then used for the evaluation of hemodynamics.

The rupture sites were marked by operating neurosurgeons on 3D CT angiography reconstructions taken after rupture. These areas were determined as intersections of computational meshes with a sphere with a radius of 0.5 mm.

The resulting hemodynamic parameters, which can be found in Table 5, were averaged over the aneurysm dome. It can be seen that the PWSS values are nearly two times larger than TAWSS values. Both parameters were additionally normalized with respect to the corresponding parent artery. The PLSA and TALSA in cases A–E in comparison with case F were significantly smaller. However, due to averaging over the aneurysm dome, it is very difficult to distinguish between different hemodynamic cases solely from these values. Thus, we considered the spatial distribution of hemodynamic quantities in the subsequent analysis. We exploit the fact that the rupture site is known and evaluate the OSI and TAWSS averaged both over the dome and over the rupture site and their relative percentage differences, see Table 6. In models A–D there was moderate decrease of TAWSS at the ruptured site compared to the TAWSS over the respective domes, which was even more pronounced in model F. On the contrary, the opposite ratio, significantly higher TAWSS compared to the average TAWSS over the dome, was found in model E.

Figures 4, 5 show the distribution of OSI and TAWSS in each aneurysm, as well as streamlines within the domain. Aneurysms A–D had ruptured near the vortex with high OSI, aneurysm E had ruptured near the flow jet with high WSS, and aneurysm F had ruptured on the bleb. In each figure, the rupture sites are colored white or blue.

Throughout CFD analysis there was some uncertainty about the prescribed boundary conditions and inlet flow profile. Therefore, validation using a laboratory model was carried out. We chose model B and C for this purpose. There was steady flow at a generic value of 169 mL/min with a parabolic profile and zero outlet pressure applied for both models. The same conditions were set for laboratory models (i.e., pulsation damper after peristaltic pump, volumetric control of the flow, open system with tubing diameter exceeding the diameters of the inlet or outlet of the laboratory model). Qualitative characteristics such as flow direction, position of flow jets and central vortexes in the aneurysm dome were very similar between CFD and PIV in the laboratory model shown in Figure 6. In order to provide thorough support to CFD data, quantitative analysis was carried out. Flow profiles were determined by PIV and computed by CFD within the focal depth of the PIV system. Concerning the precise alignment of laboratory models with CFD model meshes (PIV measurement 1 in Figure 7), good conformity was achieved between CFD and PIV data. Repeated positioning of the laboratory model under the microscope stage could result in imprecisions as seen in PIV measurement 2 in Figure 7. Therefore, matching individual PIV measurements could be worse than matching PIV to CFD, if inaccurate positioning of the model occurred. The PIV measurement was charged with a certain variance of data, see Figure 7. This could be attributed to the finite size (30 μm) and density (1,050 kg m⁻³) of particles compared to massless particles used in CFD.

In CFD analysis we utilized the no-slip boundary condition on the vessel wall, zero pressure condition on each outlet part of the boundary, and parabolic profile of inlet flow velocity multiplied by a time-dependent waveform during a cardiac cycle. Concerning the proper choice of boundary conditions, there has been no general agreement on the optimal prescribed parameters. Recent studies found that both inflow and outflow boundary conditions have considerable impact on the hemodynamics (32). Commonly used inflow and outflow conditions are listed, for example (10). In this study, equal outlet pressures were the most frequently used outlet boundary conditions, while on the inlet, most groups prescribed a constant plug profile for either velocity or flow rate. However, in our case we had access to a set of velocity waveforms. Thus, the average time-dependent waveform is likely to model each patient-specific inlet boundary condition more realistically than a constant plug profile or flow rate.

To further address uncertainties about CFD analysis conditions and to evaluate the match of CFD to real flow, a validation using PIV in a laboratory model was carried out. For this purpose, two models (B and C, with the most complex inlet/outlet structures) were chosen. We capitalized on steady inflow with a parabolic profile in order to match the laboratory model with a pump equipped with pulsation damper and connected to tubing. This apparatus in principle develops a parabolic (laminar) flow profile at a flow speed up to several meters per second.

Under these conditions we observed a good qualitative agreement of CFD and PIV as the flow direction, position of flow jets and vortexes matched well. The quantitative agreement of CFD and PIV was well proven by determining flow profiles at line segments from the top of the aneurysm to the center of the aneurysm dome. If the laboratory model was aligned properly with the CFD model mesh there was agreement between CFD and PIV, which was similar or better than reported in literature (33–37). The possible limitation of this approach arose from imprecisions caused by repeated laboratory model positioning for PIV measurement.

Despite providing reasonable data the study has few limitations which should be taken into consideration. Since the primary focus was the relationship of aneurysm rupture and hemodynamics, the study did not consider possible co-morbidities (smoking, hypertension, dyslipidemia, etc.) which can affect the aneurysm development and rupture. Specific limitations of studies in ruptured aneurysms result from the fact that we already assess the condition given by rupture, not the aneurysm at risk, before rupture. Therefore, altered morphology of the aneurysm after rupture can lead to erroneous results or misinterpretations (38). Another limitation in these studies is that we usually do not possess the individual boundary conditions for each aneurysm and need to use literature-based information. In our CFD analysis, we employed a straightforward choice of boundary conditions. We prescribed a generic parabolic velocity profile at the inlet and applied a do-nothing condition on all outlet boundaries. Such an approach may introduce a bias into the CFD analysis if processing aneurysm model based on a patient with anomalous inlet/outlet boundaries.