Study of the Effects of Wall Thickness and Size Variations on the Rupture Risk of Cerebral Aneurysms Using FSI Simulations

1. Introduction

A cerebral aneurysm is a lesion that is presented as a focal bulging in blood vessels, it has a prevalence of about 3% in the population without underlying health conditions [1] and is typically located at a region of the brain known as the circle of Willis [2]. The most severe outcome of this lesion is when the aneurysm ruptures leading to a subarachnoid hemorrhage (SAH) [3,4], which carries a mortality rate exceeding 40% [5].

The serious consequences that a cerebral aneurysm may produce, along with the high-complexity treatments the patient can receive, as neurosurgery, create the need for finding a method to predict the rupture risk of a cerebral aneurysm to make a better decision as to whether the patient needs a complex treatment or not. Nowadays, that method does not exist, but it is known that the interaction between the blood flow and artery walls is the key factor in explaining aneurysm formation, as well as its development and rupture [6].

It is thought that the wall shear stress (WSS) is one of the most important hemodynamic factors to study, where zones with high WSS levels would be more prone to suffer from an aneurysm formation [7,8,9], while zones with low WSS have been related with the development and rupture of an already initiated lesion, finding low-stress levels in those sites where the rupture is produced at the aneurysm dome [10,11], although high WSS levels combined with other factors, could also lead to rupture [12].

Once the aneurysm is initiated, it can have different outcomes, it can rupture young, develop increasing its size to rupture, or be kept in time [13]. The size of an aneurysm does not evolve at a constant rate [14], going through periods where it does not change, and other periods where it can grow, even with different growth rates between ruptured and unruptured aneurysms [15]. It has been shown that aneurysms growth can produce significant morphological changes from at least 23 weeks [16].

Although it may be intuitive to think that an aneurysm will rupture when it is large, it has also been observed that small aneurysms suffer from ruptures [17], which suggests that not only its size is a risk factor. In addition to health conditions such as hypertension or trauma to explain why small aneurysms can rupture, it must be considered that the aneurysm is thinner than the artery [18], and contrary to what could be thought, when the aneurysm grows, it does not necessarily decrease its wall thickness [19] because it is a living structure that can repair itself.

There are several studies modifying hemodynamic parameters to reflect different health conditions that might affect the rupture risk of an aneurysm, and less common are studies where geometrical factors are modified. Sun et al. [20] have shown that when an aneurysm grows and its thickness decreases, the rupture risk increases due to a reduction of the WSS and an increment in the Von Mises stress at the aneurysm, and even keeping a constant wall thickness, the maximum Von Mises stress still increases. Nath et al. [21] studied spherical aneurysm geometries using structural simulations. They found that at the beginning, the aneurysm has a reduction in its maximum Von Mises stress, but then it increases along with the aneurysm size, meaning a higher rupture risk. Other works have used cylindrical/spherical membranes subjected to 1D enlargement to model the development of an aneurysm where the mechanical response is attributed to collagen and elastin [22].

Those previous works have used a uniform distribution of thickness over the aneurysm dome, however, real aneurysms do not have a uniform thickness distribution, and this can be related to stress concentrations [23] or aneurysm expansion for abdominal ones [24]. Nowadays, even with the current medical technology, wall thickness remains a difficult parameter to obtain from patients without performing surgery or removing the aneurysm, especially for cerebral aneurysms. Furthermore, by removing an aneurysm, it has been noted that using a non-uniform thickness distribution in simulations show much higher stress values in the rupture site when compared to a configuration with constant wall thickness [25].

In the present work, the effect of geometrical parameters is studied, specifically, the effects of the aneurysm wall thickness and its size on its rupture risk, separately varying each parameter for better accounting how each parameter affects the rupture risk of cerebral aneurysms. Six patient-specific geometries are chosen and boundary conditions representing healthy conditions for an adult patient are used. Due to the importance of accounting for blood and wall interaction in simulations for an accurate prediction of patient-specific hemodynamics [24,25], a fluid-structure interaction (FSI) approach is used. Besides, as differences in up to 30% in displacement have been observed when comparing two-way, and one-way FSI simulations [28], the former is used.

Additionally, cases where the aneurysm is virtually removed from lateral geometries are included, with the aim of simulating a state before the lesion initiation.

2. Materials and Methods

2.1. Materials and Methods

The obtained geometries belong to a cerebral aneurysm database generated with the methodology exposed by Valencia et al. [29]. That group obtained medical images provided by the Instituto de Neurocirugía Dr. Alfonso Asenjo (INCA) using the tridimensional rotational angiography method with a Philips Integris Allura device. Six geometries were chosen for the study, three of them ruptured and three unruptured. Those geometries can be seen, along with their names in Figure 1, where those with rupture are named RG, and those without rupture are named UG. Table 1 shows each geometry name with its type, location and total geometry volume. Both RG-1 and UG-1 are used to study the aneurysm thickness effect, the rest of the geometries are used for studying the aneurysm size effect. Since the worst outcome for an aneurysm is its rupture, those geometries that presented rupture are modified in such a way that their size is reduced to study a past state, and those with no previous rupture are modified to increase their size to simulate a possible future state. The methodology used to modify the size is shown next.

2.2. Size Variation

The employed algorithm is an adaptation from the one described by Sun in the previously cited work. That algorithm is modified to be used in T-Spline bodies in Autodesk Fusion 360 and to decrease those ruptured aneurysms.

Instead of using offsets as in the original procedure, here, the scaling process uses a scale factor s to expand the aneurysm parallel to its neck plane. s is defined as the ratio between equivalent diameters after (D_f) and before (D_i) scaling:

D_f/D_i = s

(1)

The assumption that the contours in which the aneurysm is sectioned are similar to a circumference is made, with this, the next formula is used:

D_f ≈ D_i + 2Δ

(2)

Where the approximation is taken as equality. Δ is the offset used by Sun et al., i.e., Δ = (α − 1)D_i. Then, the previous expression is written as D_f = D_i + 2(α − 1)D_i and s is written as follows:

s = (2α – 1)

(3)

Expression 2.3 is used for unruptured aneurysms, to make them grow, for ruptured aneurysms, s′ is used instead, and defined as s′ = 1/s.

The main steps of the algorithm are listed below:

• From the original geometry, the aneurysm is extracted to convert its main faces into a T-Spline body. The neck plane is visually identified, and a point is placed at its center, by which a perpendicular axis V passes through.
• A factor α_max = 1.26 is used for all the cases. As in the original algorithm, an α parameter that begins at 0 when V = 0, and linearly varies up to α_max for V > V_max /3, is used. V_max is the distance between the point at the center of the neck and the highest point of the aneurysm dome. The height of each contour in the T-Spline body, approximately parallel to the neck plane, is measured and named V_c.
• Each contour is moved upwards, for unruptured aneurysms, and downwards, for ruptured aneurysms, in a quantity V_c (α − 1). Additionally, each contour is scaled with the factor s or s′ depending on the rupture status of the aneurysm, parallel to the neck plane. The moving and scaling process are repeated with each one of the contour lines of the body.
• When the process is completed, the T-Spline body is converted into surfaces, the top zone of the aneurysm is closed following the curvature of the surrounding faces, and the zone at the bottom of the aneurysm is closed with another surface. Then, the created surface bodies are joined and converted into a solid. Finally, the resultant scaled aneurysm is joined back to the artery. For lateral aneurysms only, a fillet is created in its neck to create an organic transition between the artery and the aneurysm.

An example of the main steps followed can be seen in Figure 2 , where the algorithm is applied to UG-2.

Since it is not known in beforehand at what size the aneurysm will rupture for those unruptured geometries, nor it is known how structural and fluid parameters are going to change before performing the simulations, just two larger variations are created, which can be seen for both original geometries in Figure 3. In a), I size geometries are shown, while in b) and c), II and III sizes are also shown respectively. Here and for ruptured geometries, the smaller the number the smaller the aneurysm size.

For ruptured geometries, its variations can be seen in Figure 4, where a) is the IV size, b), c) and d) are the III, II, and I sizes. The algorithm is applied beginning from geometries a) and b) to get b) and c) respectively, whereas d), which represents an aneurysm recently initiated, has a semi-spherical morphology whose height is approximately the same as two consecutive applications of the algorithm would yield starting from c). Just for size I for ruptured geometries, the aneurysm neck is decreased.

2.3. Thickness Variation

To reflect the aneurysm thickness effect, just in one case that region is set with the same thickness of the artery, 0.35 [mm], whereas in the other two cases, the aneurysm thickness is reduced to 0.2 and 0.10 [mm] since aneurysm walls can be as thin as 0.05 [mm], and even less [30]. Thickness variation is done when creating the mesh, it is set as a property of the structural element, which is mentioned in the following section.

From the original geometry, two surface bodies are created, seen in gray for the aneurysm and green for the artery in Figure 5. The artery is kept at 0.35 [mm] of thickness for all cases, while the aneurysm varies as stated before. To get a compatible mesh between those two zones, the Shared Topology option from SpaceClaim is activated at their boundary.

2.4. Numerical Setup

2.4.1. Fluid Model

Since there is a wide variety of sizes on the selected geometries, it is chosen to use the one that has a total volume similar to the average volume. RG-2 is selected to conduct the mesh tests for both domains. The fluid mesh test is based on the mesh density, the ratio between the total number of elements and the geometry total volume, using computational fluid dynamics simulations in Ansys Fluent to solve the respective equations associated with the fluid using the finite volume method (FVM). Densities go from 250 [elements/mm³] up to 3000 [elements/mm³], using both, linear and quadratic tetrahedral element. The chosen scheme is PISO, with a spatial discretization using least squares cell based on the gradient and second order for pressure and momentum. The transient formulation is first-order implicit, with a time step of 0.0005 [s], 3600 time steps, and 200 maximum iterations per time step.

Blood model is the same as one previous work that uses geometries from the same database [31], where blood is modeled as a laminar and incompressible Casson fluid with a density of 1065 [kg/m³]. Boundary conditions for the entry are set as a Womersley velocity profile, whereas for the exit, a three-element Windkessel model with pressures ranging between 80 and 120 [mmHg] is chosen. The total simulated time is 1.8 [s], slightly more than two complete cardiac cycles, and the considered results are the area-averaged WSS at the diastole (DWSS) and systole (SWSS), from the second cardiac cycle. Simulations were performed using double precision in a machine with an AMD Ryzen 5 5600G processor and 16 GB of RAM.

In Figure 6 the evolution of both parameters with respect to the mesh density can be seen. Since between 500 and 1000 [elements/mm³] the differences are approximately less than 5%, and that two-way FSI simulations need considerable computational resources, it is chosen to use a mesh density of about 700 [elements/mm³], with size refinements at the aneurysm dome and in the zone in contact with the artery wall. Besides, the differences between quadratic and linear elements are negligible, and therefore, the former is chosen.

2.4.2. Structural Model

Mesh tests for the structural domain are based on the element size using the same geometry as the previous tests, but with computational structural dynamics simulations, using Ansys Transient Structural to solve the corresponding associated equations with the finite element method (FEM). Linear shell element is used since the thickness of the artery and the aneurysm are considerably less than other dimensions of the geometry. It is common to assign the element thickness as 10% of the artery diameter, in this case, considering the average diameter of the selected geometries, the thickness is set to 0.35 [mm]. This thickness is used for both the aneurysm and the artery, and it is kept for all the cases for the size effect.

Element size goes from 0.15 up to 0.50 [mm]. An external pressure representing the cerebrospinal fluid action for a healthy adult of 1650 [Pa] [32] is used, and an internal pressure representing the action of the blood flow is set to 15 [kPa]. Inlet and outlets are assumed to have fixed support. A five-parameter Mooney-Rivlin hyperelastic model is used for both mesh tests and FSI simulations. Its elastic energy density function W can be seen in Equation 4:

W = 𝑐₁₀(𝐼₁ − 3) + 𝑐₀₁(𝐼₂ − 3) + 𝑐₁₁(𝐼₁ − 3) (𝐼₂ − 3) +𝑐₂₀(𝐼₁ − 3)² + 𝑐₀₂(𝐼₂ − 3)²

(4)

I₁ and I₂ are the strain invariants, c₁₀, c₀₁, c₁₁, c₂₀ and c₀₂ are model constants obtained from the literature [33]. Although the wall of an aneurysm is anisotropic, in this work it is modelled as an isotropic material, a common practice to do since both approaches yield similar results [34].

Considered results are the total displacement and the Von Mises stress at the aneurysm. Its results and percentage differences with respect to that mesh with the smallest element (D₀₁₅) can be seen in Table 2 where just for the Von Mises there are differences greater than 5%, then, it is chosen to use a maximum size element of 0.20 [mm], since its differences are below that value.

2.4.3. Fluid-Structure Coupling

Simulations performed in this study are two-way FSI using Ansys System Coupling for exchanging information between Fluent and Transient Structural. Boundary conditions to ensure the compatibility of both domains are that the structural and fluid displacement must be compatible with each other (Equation 5), traction at their boundary must be in balance (Equation 6), and the fluid must obey the no-slip condition (Equation 7) [35].

δ_S = δ_F

(5)

$${\stackrel{\mathrm{-}}{\mathsf{\sigma }}}_{\mathrm{S}} · {\widehat{\mathrm{n}}}_{\mathrm{S}}={\stackrel{\mathrm{-}}{\mathsf{\sigma }}}_{\mathrm{F}} · {\widehat{\mathrm{n}}}_{\mathrm{F}}$$

(6)

$$\overrightarrow{\mathrm{u}}={\overrightarrow{\mathrm{u}}}_{\mathrm{g}}$$

(7)

In the previous expressions, δ is the displacement, $\overline{\sigma }$ is the stress tensor, $\widehat{n}$ is the normal vector at the boundary, $\overrightarrow{u}$ is the velocity vector, ${\overrightarrow{u}}_g$ is the mobile velocity coordinate, and subscripts S and F stand for the structural and fluid domain respectively.

Additionally to those settings already mentioned for the fluid domain, dynamic mesh options were activated with smoothing and remeshing methods besides of creating information exchange zones. In the structural domain, the internal pressure is removed, and exchange information zones are also created.

3. Results

A total of 24 FSI simulations were carried out to evaluate the effect of geometrical factors, such as the aneurysm size and wall thickness, and to obtain the structural and hemodynamic environment of the lateral arteries before the lesion initiation. The aneurysm wall thickness effect is first presented, then the size effects, and lastly, the aneurysm removal results.

3.1. Aneurysm Wall Thickness Modification

The aneurysm wall thickness modification produces negligible differences for hemodynamic parameters, where for both, average pressure and WSS, differences are less than 5% and therefore, just the temporal evolution of the average WSS at the aneurysm is presented in Figure 7, for each wall thickness and for each geometry, in the last part of the simulated cycles, where the curves are almost indistinguishable from each other. That figure and following similar ones were created using MATLAB R2023a.

Structural parameters include the displacement and the Von Mises stress. The displacement behavior is opposite between the ruptured and the unruptured aneurysms as can be seen in Figure 8, where the average displacement in the aneurysm is plotted against the wall thickness for each case. In Figure 9 the same plot is presented but for the Von Mises stress, where both geometries increase their average perceived stress as the aneurysm wall becomes thinner. Displacement contours can be seen in Figure 10 and Figure 11.

3.2. Aneurysm Size Modification

The evolution of the aneurysm size makes a large difference for both, hemodynamic and structural parameters, where the average WSS temporal evolution for each geometry size variation is shown in Figure 12, whereas Figure 13 shows the average WSS evolution against the AR of each variation and for each geometry, where except for UG-3, as the AR increases, the WSS reduces.

Figure 14 shows the evolution of the average displacement at the aneurysm against the AR for each case. Once again, UG-3 exhibits a different behavior than the rest of the geometries, which is as their AR increases, the displacement increases as well, except for size I variations of the ruptured geometries, which repeats for the Von Mises stress in Figure 15. Here, the behavior of UG-3 is like the rest of the cases, where an aneurysm of greater AR suffers from greater Von Mises stresses. Examples of the obtained Von Mises contours can be seen in Figure 16 and Figure 17 for both terminal geometries.

3.3. Aneurysm Removal

In order to compare the pre-aneurysm with the aneurysm state, the time-averaged wall shear stress is calculated using the following equation:

$$\mathrm{TAWSS}=\frac{1}{\mathrm{T}}{\int }_0^{\mathrm{T}}\left|{\mathrm{WSS}}_{\mathrm{a}}\right|\mathrm{dt}$$

(8)

Where T is the cardiac cycle time and W SSa is the area averaged WSS at the aneurysm. In Figure 18, TAWSS distribution for the healthy arteries is shown, with an arrow pointing to the future lesion location. In Figure 19 the same can be seen for the pressure distribution.

4. Discussion

4.1. Effect of the Aneurysm Wall Thickness Modification

Although hemodynamic parameters did not significantly perceive the effect of modifying the aneurysm thickness, structural parameters did. As seen in Table 3, if the aneurysm wall becomes thinner, in this case, changing from 0.35 to 0.10 [mm], the average Von Mises stress can increase in more than 200%, which would indicate that the aneurysm has a greater rupture risk due to a structural failure of its wall.

This type of change in the stress magnitude could be expected from Laplace’s law, a principle which relates the pressure, radius, thickness and tension for a sphere, although it could only give an approximate idea of how tension would evolve with the thickness change, whereas a numerical simulation can yield results for the entire dome of an aneurysm with other stresses, such as the Von Mises stress, that could give a better notion of how close the aneurysm is to a structural failure. Besides, it would be not realistic to assume an aneurysm as a sphere, and this has previously shown to give different results [36].

In the other hand, the displacement has a non-intuitive behavior, since it is opposite between both geometries, being the ruptured geometry, RG-1, the one that has the expected trend, where the thinner the wall, the greater the displacement is, furthermore, it is found that this geometry is ten times more sensitive to the wall thickness change than the UG-1 according to Table 4.

The difference in behaviors and sensitivities can be explained due to the parent artery geometry and the blood flow path that it imposes, this can be seen for both chases in Figure 20, where RG-1 has a blood flow that points to the aneurysm, entering that zone more directly, and with more velocity than in UG-1. Once the fluid enters to the RG-1 aneurysm, it travels to its upper zone, where it changes its direction and decreases its velocity magnitude, exchanging momentum with the zone the aneurysm that has the highest displacement values as can be seen in Figure 10 for its contours. In the other case, almost the entire blood flow follows the artery path with very low deviation into the aneurysm, making its displacement very similar to its surrounding zone in the artery, as Figure 11 shows.

4.2. Effect of the Aneurysm Size Modification

As the WSS time evolution shows, there is a clear distinction on its magnitude when the aneurysm evolves in size, except for UG-3, which shows an increase followed by a decrease. These results are similar to those obtained by Sun, in that work, the WSS curves were clearly separated from each other, being the scaled geometry the one with higher levels of WSS. The curves obtained by that research group were closer to each other than in the present case, which can be explained because in the present case of modifying the aneurysm size, the wall thickness remained constant and with a different thickness. Other factors of differences are that the geometries are different and, in this case, a non-Newtonian model for the blood was used.

Table 5 shows that ruptured geometries have an increment in the average WSS for the I size case in between 205 and 276%. Geometry UG-2 shows a reduction in its WSS of around 36% considering both systole and diastole, for its largest variation. UG-3 has an increase of up to 31% at the diastole for its II variation, but a decrease of less than 5% for the largest one, which can be considered as negligible.

The distinctive behavior of geometry UG-3 can be explained by how the path of the blood flow changes as the aneurysm evolves, seen in Figure 21. In the original geometry, the smallest one, there is a recirculation zone for the blood before entering into the aneurysm and then going to the main bifurcation. In the following variation, the size II aneurysm, blood flow enters more directly and with more velocity into the aneurysm, which explains the WSS increment. The blood flow pattern is similar in the largest variation, but just like the rest of the cases, there is a decrease in the WSS levels because now the blood flow reaches the dome with a lower velocity.

In terms of displacement, a larger aneurysm might be expected to have greater displacements because it is farther away from the anchor point to the artery, as it was obtained for most of the cases and the maximum is achieved at the top of the aneurysms, where the blood flow changes its magnitude and direction, transferring its momentum to the aneurysms wall. The distinctive behavior of UG-3 can be explained due to the complexity of its geometry, which includes three bifurcations and whose maximums tend to be located near the neck.

As can be seen in Table 6, the maximum reduction in displacement for both ruptured geometries are achieved for the size II variation, where for RG-2 that reduction is almost 30%, whereas for RG-3, the displacement is reduced both for diastole and systole, in about 54%. UG-2 has the maximum displacement increment in its largest variation, being this of about 60%. For UG3, its size II variation increases its displacement by approximately 10%, while its III variation gets a reduction of about 22% with respect to the original size.

The Von Mises stress evolution with respect to the aneurysm size increment means that larger aneurysms are prone to a higher rupture risk. Percentage differences are presented in Table 7. For both ruptured geometries, their greatest reduction is achieved for II size variations, being approximately 30% for RG-2 and between 23 and 29% for RG-3. For unruptured geometries, in both cases the increase of the Von Mises stress at the aneurysm is between 25 and 29%.

Figure 16 and Figure 17 show examples of the Von Mises stress at the aneurysm evolution for cases RG-3 and UG-3, where it can be noted that the top of the aneurysms shows higher levels of stress than other zones (except for the neck), even in the latter case which has a smooth surface. This is in accordance with the zones where ruptures are typically found. For those ruptured geometries, a reduction of the stresses from size I to size II variation can be seen, to then just increment. This difference of behavior can be explained by the very different morphology that size I variations for ruptured geometries have, which are in both cases, semi-spherical aneurysms, where the change of curvature of their necks can lead to stress concentrations, and since these aneurysms are smaller than the rest of the cases, it has a greater impact on its averaged results.

4.3. Pre-Aneurysm State

The presented contours of the TAWSS and pressure show that the lesion forms in zones where those parameters have medium magnitudes, in no case minimum values. As a comparison, the approximate value of the TAWSS at the zone of the healthy artery where the aneurysm is located afterwards when it occurs (shown in Figure 22), is visually obtained and compared with the same parameter, but averaged at the original aneurysm. Table 8 shows that for all four cases, there is a reduction in the TAWSS when the aneurysm is present, and interestingly, for both ruptured cases, the reduction is approximately 75%, while for the unruptured cases, the reduction is lower, being 45% for UG-1 and 33% for UG-2.

The reduction is explained because when the aneurysm is not present, the blood flow is in contact with the artery with a higher velocity than when the aneurysm occurs, because then it needs to travel a longer path into a cavity that makes it more difficult to the blood to flow due to the change in direction and because there is blood with a lower velocity filling the cavity.

5. Limitations

Real aneurysms exhibit a complex mechanical behavior along with a non-uniform wall thickness distribution and a non-predictable size and shape evolution due to the several factors that can affect its evolution and that can change from patient to patient, nevertheless, assumptions such as constant wall thickness, no pre-stretch, and isotropic wall structure are common practices that have shown to produce realistic structural responses [26,35] to allow investigators to develop models with a complexity that makes it possible to find a solution and predict an structural and hemodynamic evolution to serve as an early clinical aid in a future application of the present methodology.

6. Conclusion

Using computational simulations, this research could study the past evolution of aneurysms and explore potential future developments as well, a task that would otherwise not be possible to perform due to ethical concerns. This study demonstrated that for the studied group, as cerebral aneurysms grow, they trigger significant alterations in hemodynamic and structural parameters which reflect a greater rupture risk environment. For instance, an unruptured geometry can decrease its average WSS by up to 36% when it grows and increase its Von Mises stress by almost 30%. Furthermore, the use of FSI simulations made it possible to assess the effects of aneurysms wall thinning keeping their size unaltered, revealing that even if the aneurysm does not grow, it could be on a high rupture risk state because of a thin wall which is prone to a structural failure since its average Von Mises stress can increase in up to 260%. The use of simulations also enabled the study of the lateral arteries when the aneurysm is removed, finding that for both ruptured aneurysms, the TAWSS was reduced in 75% at the rupture point, whereas for unruptured aneurysms, although they also have a reduction, it was not as high. It is then concluded that even for a healthy patient, either of the geometrical modifications presented in this study implies a significant increment on the rupture risk of the studied cerebral aneurysms, and that there could be an indication of the rupture point of the aneurysms when its TAWSS reduction is obtained, although more samples are needed to give more statistically significant conclusions about this and the geometrical factors effect.