Rupture point is associated with divergent hemodynamics in intracranial aneurysms

Background Understanding the risk factors leading to intracranial aneurysm (IA) rupture have still not been fully clarified. They are vital for proper medical guidance of patients harboring unruptured IAs. Clarifying the hemodynamics associated with the point of rupture could help could provide useful information about some of the risk factors. Thus far, few studies have studied this issue with often diverging conclusions. Methods We identified a point of rupture in patients operated for an IAs during surgery, using a combination of preoperative computed tomography (CT) and computed tomography angiography (CTA). Hemodynamic parameters were calculated both for the aneurysm sac as a whole and the point of rupture. In two cases, the results of CFD were compared with those of the experiment using particle image velocimetry (PIV). Results We were able to identify 6 aneurysms with a well-demarcated point of rupture. In four aneurysms, the rupture point was near the vortex with low wall shear stress (WSS) and high oscillatory shear index (OSI). In one case, the rupture point was in the flow jet with high WSS. In the last case, the rupture point was in the significant bleb and no specific hemodynamic parameters were found. The CFD results were verified in the PIV part of the study. Conclusion Our study shows that different hemodynamic scenarios are associated with the site of IA rupture. The numerical simulations were confirmed by laboratory models. This study further supports the hypothesis that various pathological pathways may lead to aneurysm wall damage resulting in its rupture.


Introduction
In clinical practice there has been a substantial increase in the detection of unruptured intracranial aneurysms (IA) within the population due to growing use of non-invasive neuroimaging diagnostic methods (1).This creates a burden not only for the patients themselves, but for the clinicians as they lack precise indicators for optimal therapy, when consulting these patients.Computational fluid dynamics (CFD) has gained increasing interest in the clinical community, as it could potentially clarify the pathophysiology of IA development, growth, and rupture.Some hemodynamic parameters have already been integrated in newly developed scoring systems (2,3).One of the first grading systems called rupture resemblance score (RRS) system was proposed recently by Meng et al. utilizing a combination of hemodynamic and morphological features of ruptured IAs (4).
Few studies have focused on local hemodynamics with respect to the site of rupture (5,6).Modeling the local hemodynamic environment could reveal key information in understanding its role in aneurysm rupture.On the other hand, a low number of studies is available, which is most likely due to the fact that it is often complicated or even impossible to identify the precise point of rupture in an IAs.Proper identification of the point of rupture can be performed at several points-during surgical exploration of the IA, during clipping or during an endovascular diagnostic procedure, if periprocedural rupture occurs.Thus far, ruptured IAs have been associated with both low and high WSS.The relationship between WSS and aneurysmal rupture is complex and not fully understood (7)(8)(9).The relationship between the point of rupture and hemodynamics is even less clear (5).
Despite a growing body of literature on hemodynamics and IAs, it has little impact on the clinical decision-making regarding patients with unruptured IAs.This may be related to the complexity of the mathematical model as well as the uncertainty surrounding the role of individual hemodynamic parameters (10).Modeling the local hemodynamic environment could reveal key information on its role in aneurysm rupture and possibly bring our knowledge closer to clinical every-day reality (10).

Aneurysm selection/patient population
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.

Image segmentation
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 craniocaudal, 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.

Mesh generation
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.

The model and boundary conditions
Navier-Stokes equations with constant dynamic viscosity μ = 3.5 mPa•s and constant blood density ρ = 1,050 kg•m −3 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.

Hemodynamic quantities
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 postprocessing calculations were performed using Python and visualizations using ParaView (19).

In vitro models
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).

Particle velocimetry
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 −3 .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.

Comparison between CFD and PIV
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   The inlet flux was fixed.
Scheme of the setup for the particle velocimetry.
Hejčl et al. 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.

CFD modeling
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.

CFD validation with particle velocimetry
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 −3 ) of particles compared to massless particles used in CFD.

Discussion
There has been a slow, but gradual increase of interest in hemodynamics of IAs in the field of clinical neurosurgery.The field of hemodynamics has been studied from different perspectives, such as their relationship to IA development, growth or rupture and some hemodynamic parameters have been integrated in newly developed risk factor scoring systems to help physicians estimate the risk of rupture (2,3).In recent years, few studies have focused on local hemodynamic parameters associated with the site of rupture in IAs (5).Due to difficulties in identifying the precise point of rupture, the majority of such studies consisted of case reports or small mini-series.
These mini-series presented divergent hemodynamic scenarios.Interestingly, the site of rupture was characterized by a uniform hemodynamic environment within each study.We believe it is the first time that we described divergent hemodynamic scenarios at the site of rupture within a single study.Most aneurysms (n = 4) in this study, had a rupture site associated with an area of low WSS and increased OSI compared to the aneurysm dome as a whole.In another aneurysm, the site of rupture was associated with a jet flow directed against the aneurysm sac and very low OSI compared to the aneurysm dome as a whole.We have already described this type of hemodynamic environment in our former case report in a ruptured anterior communicating artery aneurysm (22).In the last case, the site of rupture was within a bleb of the aneurysm dome and was associated with low WSS and low OSI.
These typical hemodynamic scenarios have been described in previous studies (5).Our own study consisted of a relatively uniform series of aneurysms.Most of the IA, except for 1, were MCA aneurysms and all of them were large-size aneurysms, ranging from 8 to 13 mm in the largest diameter.Despite the relative homogeneity of aneurysm size and location, hemodynamics at the point of rupture differed significantly.This is of interest, because from a clinical perspective all of these aneurysms, if unruptured, would be viewed similarly when discussed in terms of rupture risk.Specific hemodynamic parameters have been shown to be related to different qualities of the aneurysm wall (23).In that study jet flow and high WSS were associated with thin, translucent regions of the aneurysm wall.Our perioperative photo of the only aneurysm whose point of rupture was defined by high WSS and jet flow demonstrates highly translucent aneurysm wall (Figure 8).We have already described this hemodynamic environment associated with aneurysm sac rupture.
On the other hand, most rupture points in our IAs showed low WSS together with high OSI.Such hemodynamics have been shown to be associated with thick, atherosclerotic regions.Atherosclerosis is often found in large aneurysms, which was the case of all the aneurysms in our study (24).Wall shear stress influences the morphology and orientation of endothelial cell within the vessel lumen.While high WSS may lead to their elongation and alignment, endothelial cells subjected to low WSS are rounded with no preferred alignment (25,26).Moreover, exposure of the arterial wall to a relatively low wall shear stress may increase intercellular permeability and consequently increase the vulnerability of these regions of the vessel to atherosclerosis (5).Gradual enlargement of the LSA over time is associated with IA rupture (27).One of the aneurysms in our series ruptured within the secondary sac or a bleb.Blebs seem to form at or near regions of high flow regions and elevated WSS (28).After bleb is formed it is characterized by low WSS as confirmed in this case.The high WSS associated with bleb development leads to thinning the aneurysm wall making it prone to rupture.In general, IA irregularities including secondary blebs are associated with an increased risk of aneurysm rupture (29-31).
Despite our study is based on a small series of IAs, it confirms the hypothesis that different pathophysiological pathways result in IA rupture (7).More specifically, while low WSS and high OSI may lead to activation of inflammatory cells (atherosclerosis), high WSS may lead to vessel wall cell activation resulting in inadequate growth and rupture (translucent aneurysm wall).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 timedependent 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.

Study limitations
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 Measured velocities along the line segment obtained by two PIV measurements (blue and red circles) together with computed velocities via CFD projected to the same line.The line was projected from the top of the aneurysm dome where rupture occurred through the center of the aneurysm in the measurement plane (black line); upper panel: aneurysm B, lower panel: aneurysm C. In particular, the black line represents computed velocities in the measurement plane, while the grey lines denote velocities at the border of the focal depth (247 μm) of PIV.For clarity, the space between the border lines was highlighted with light gray color.In PIV measurement 1 (PIV1) the laboratory model position was aligned with CFD mesh precisely.Whereas in the 2nd PIV measurement (PIV2) the position of laboratory model was taken from PIV1. 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.

Conclusion
Using a PIV validated CFD model of a ruptured aneurysm, we evaluated hemodynamics of 6 ruptured aneurysms.The most important finding was identification of rupture points associated with various hemodynamic features (low WSS with recirculations contrary to high WSS with a jet flow).Hence our data indicates that divergent biomechanical pathways are involved in the mechanism of aneurysm rupture.

FIGURE 1
FIGURE 1Inflow velocity waveform set for all cases and corresponding inflow flux waveform for radius = 2 mm.

FIGURE 2
FIGURE 2Silicone models: the lumen of the model was 3D printed out of water soluble material which was embedded into silicone together with inlet/outlet tubing and removed with hot water after silicone curing.The lumen of the final models was filled with phenol red solution to enhance visibility.Left column: model B, right column: model C.

FIGURE 4
FIGURE 4 Left column: OSI, middle column: streamlines colored by velocity magnitude, right column: TAWSS.First row: aneurysm A, second row: aneurysm B, third row: aneurysm C.

FIGURE 6
FIGURE 6Left column: PIV, right column: CFD results, first row model A, second row model B. The result of PIV measurement shows sum of 10 consecutive images (focal depth 247 μm) taken at 97 Hz with 10 ms exposure with 30 μm polystyrene particle driven by flow.Thus, the traces were rather nonuniformly distributed and are laid in rows.The CFD shows the 10 ms traces of massless randomly distributed particles.The end of the trace is denoted by a dot.In order to simulate the PIV image with high fidelity the traces of particles within focal depth are displayed and coded with greyscale tones to denote the distance from the center plane of the focal depth.

TABLE 1
Morphology of aneurysms and their location in the brain.

TABLE 2
Inlet information: radius of the inlet artery (mm) and given inlet velocity (m/s).
PSV, peak systolic velocity; EDV, end diastolic velocity; TAV, time-averaged velocity.Inlet flux (mL/min) is computed from the TAV and the inlet area.

TABLE 3
Flow parameters used for the comparison of CFD and PIV.

TABLE 4
Results of the mesh dependency test.
An average edge length h in millimeters, velocity magnitude v in meters per second, normal pressure p and wall shear stress WSS in pascals.Values were obtained as the mean over the aneurysm area after the flow was steady.

TABLE 5
Values of the hemodynamic parameters over the aneurysm dome.Max denotes the maximum over the dome and aver denotes the average value over the dome.PWSS, peak wall shear stress; TAWSS, time-averaged wall shear stress; OSI, oscillatory shear index; NTAWSS, normalized time-averaged wall shear stress; NPWSS, normalized peak wall shear stress; PLSA, peak low shear area; TALSA, time-averaged low shear area.

TABLE 6
The OSI and TAWSS values were averaged over the dome and over the rupture site.Relative percentage differences were obtained with respect to the average value over the dome.