## Abstract

Hemodynamics associated with the arteries of the circle of Willis (CoW) is analyzed to identify possible cerebral aneurysm initiation locations using computational methods. A numerical fluid–structure interaction model is developed using an idealized geometry with the linear elastic, isotropic arterial wall. Blood flow is assumed to be laminar, incompressible, and modeled using Navier–Stokes equations, non-Newtonian viscosity, and sinusoidal boundary conditions. Available analytical and experimental results are used for the validation of the model. The highest wall shear stress (WSS) and von Mises stress (VMS) are identified for understanding the most vulnerable sites. The WSS distribution in the entire CoW region shows that ACoA junction has the highest value and risk of aneurysm initiation. The flow patterns created due to the geometrical features of the CoW seem to be the significant factor in the distribution of WSS. It is noticed that a decrease in wall elasticity will reduce the magnitude of WSS, both the temporal and spatial averaged value. The wall weakening effects are more pronounced for the posterior communicating artery. The wall weakening creates changes in core velocity and WSS. Changes in Von Mises stress are also noticed due to wall weakening effects. Highly localized VMS is noticed at ACoA and could possess a higher risk for aneurysm initiation and rupture. Despite the simplifications involved in developing the fluid–structure interaction model, this work demonstrates the critical locations at the CoW region regarding aneurysm initiation.

## 1 Introduction

It is demonstrated that intracranial aneurysm (IA) growth [1,2] and rupture mainly depend upon hemodynamic parameters such as wall shear stress (WSS), and the frictional force offered by the blood flow through the vessel lumen. Many studies [3,4] have claimed that the genesis and rupture risk of aneurysms is location-specific with a highly increased frequency of occurrence within the anterior regions of the circle of Willis. The most common locations for IA are the anterior communicating artery (ACoA), the middle cerebral artery (MCA), posterior communicating artery (PCoA), basilar artery (BA) summit, and bifurcation of the internal carotid artery (ICA) near the circle of Willis [3,5]. Computational fluid dynamics (CFD) simulation can provide valuable information on the local hemodynamic parameters, anatomic alterations of the weakened lumen, and the severity of IAs. This technique can predict the pressure distribution along various vessel regions, and the change of wall shear stresses within a cardiac cycle and therefore identify flow mechanisms influencing the initiation and progression of the aneurysm.

A cerebral aneurysm is characterized by local internal elastic lamina loss, extracellular matrix degradation, and tunica media thinning of the intracranial artery [6]. Stehbens [7] noted that aneurysms arise from the crotch and bifurcations of cerebral arteries with large mass flow rates. The degenerative theory on aneurysms suggests that hemodynamically induced changes such as loss of elastic properties of the vessel wall causes the occurrence of aneurysm at cerebral arterial forks [7]. The higher curvature of flexures present in the internal carotid artery may also cause an aneurysm development [8]. Along with WSS, various flow characteristics such as jet impingement [9], recirculation [10], and secondary flow [11] are also considered significant mechanical factors causing the genesis and progression of various vascular diseases. It is observed that locations with hemodynamic stress variation and subsequent internal elastic lamina damage can initiate an aneurysm. Shyy [12] observed that low WSS associated with disturbed flows at bends/bifurcations are atherogenic, whereas high WSS present in laminar flows in the straight vessel is athero-resistant. Cebral and Raschi [13] suggested that cyclic hemodynamic-induced stresses could generate wall degeneration and weakening, which could lead to aneurysm evolution. It was also argued that the higher hemodynamic stress at arterial bifurcations propel an aneurysm initiation as most of them were located in the anterior circulation. Kyriacou and Humphrey [14] showed that in addition to the shape and size, the material properties and the stress loading were important factors causing aneurysm initiation. In the aneurysm wall, the eutrophic changes such as cell proliferation and extracellular matrix production are balanced by destructive changes, cell death, and extracellular matrix degradation. The changes in the local hemodynamic parameters disrupt this balance and create aneurysm growth and rupture. Meng et al. [15] suggested that high WSS and a positive WSS gradient could initiate intracranial aneurysms. The increased stress could lead to activation and production of proteases by mural cells [16] and local elastic lamina deterioration [17]. Once initiated, the aneurysm growth could progress toward rupture through different pathways viz, (i) “inflammatory-cell-mediated destructive remodeling” at high oscillatory shear index and low WSS and (ii) “mural cell-mediated destructive remodeling” at high WSS [18]. After the aneurysmal bulge enlargement, a recirculation zone with vortices and flow instabilities will be formed inside the sac characterized by oscillatory and low WSS [19].

Hassan et al. [20] proposed aneurysm categories based on its location and geometric features of the parent vessel and conducted CFD studies using idealized geometries. The hemodynamic parameters obtained from the simulations supported the standardization proposed using the statistical approach of aneurysm rupture. Local hemodynamic parameters of ruptured aneurysms of MCA were obtained by Omodaka et al. [21] and found that the rupture locations were having the lowest wall shear stress (WSS). Fukazawa et al. [22] had obtained similar results by simulating twelve cases of ruptured MCA aneurysms. The mean WSS at the rupture location was approximately 13% and 4% compared to the dome and parent artery values, respectively. Rupture points were characterized by flow separation, multiple vortices, and low blood velocity. These results have to be interpreted with caution as the postrupture hemodynamics do not truly reflect the prerupture scenario. Xu et al. [23] analyzed the effect of wall elasticity on WSS inside the aneurysm and compared it with the rigid wall using an experimental aneurysm model. The experiment was conducted with a 0.4-mm thick silicon film with Young's modulus of 1.19 MPa and 53% aqueous glycerin as the working fluid. Peak WSS was reduced to 30% for the elastic wall in comparison to the rigid wall. They argued that the low WSS could release forces against the muscle cell, endothelial cell, and intima and might promote thrombosis formation and platelet aggregation. Xiang et al. [24] and Zhang et al. [25] studied the influence of geometry on parameters such as flow streamlines, WSS distribution, and oscillatory shear index (OSI) of the ruptured and unruptured aneurysms. Low WSS and high OSI created flow alterations and contributed around 84% of the rupture cases. The rupture was attributed to decreased endothelial nitric oxide synthase activity at the sites, which increased endothelial permeability. Zhang et al. [25] emphasized the importance of the location-specific analysis of aneurysms for better rupture risk assessment. Jiang et al. [26] pointed out the importance of aneurysm wall thickness in its progression toward rupture. They observed low WSS and high pressure at thin-walled regions, the high pressure being the prime risk factor for local wall thinning. Saqr et al. [27] have conducted a detailed review of CFD IA hemodynamic research, and the important conclusions are, (i) WSS is the most significant parameter affecting aneurysm fate, (ii) only 10% of the studies have included the non-Newtonian nature of the blood, though it results in significant variation in predicted results, (iii) transition from laminar to turbulent flow is unlikely to happen with non-Newtonian shear-thinning blood, and (iv) various fluid dynamics phenomena associated with IA are shear layers, flow impingement, recirculation and separated flows, and the corresponding vortices associated. CFD studies involving complete CoW have been attempted by few researchers. Cebral et al. [28] devised methodologies to construct CoW patient-derived three-dimensional (3D) geometries for use in computational simulations. The effect of pathological conditions such as fetal (reduced size of P1 of posterior cerebral arteries (PCA)) and missing (absence of A1 segments from an anterior cerebral arteries (ACA)) arteries were computationally simulated by Moore et al. [29]. Alnaes et al. [30] studied the effect of vessel radii and bifurcation angles on the distribution of WSS and pressure. Patient-derived CoW geometries with absence or impairment of connective arteries were simulated by Liu et al. [31] The velocity profiles and flow distributions in various arterial branches were studied in detail. Jahed et al. [32] generated a fluid–structure interaction (FSI) model of a patient-derived CoW with an aneurysm. The effect of size of the aneurysm on the flow hemodynamics was estimated. A transient 3D CFD analysis was conducted by Gaidzik et al. [33] on a healthy subject-specific CoW geometry. The study demonstrated the use of MRI velocity measurements and their implementation in numerical simulations. The above studies demonstrated that computational models could be used to analyze the pathophysiology of CoW, including aneurysm initiation risk.

Literature shows numerous CFD studies on ruptured and unruptured aneurysms, though studies on location-specific aneurysm initiation probability are pretty limited. Most of the studies discuss the rupture risk factors and the evaluation of the surgical treatment with devices. Literature data suggest that progressive degradation of the wall and the associated abnormal hemodynamics characterized by low and high WSS is the major reason for aneurysm evolution. The mechanical and geometrical properties of the arterial wall can affect the local WSS and initiate an aneurysm. In this work, an FSI model of the circle of Willis is analyzed to explore the possibility of aneurysm inception at various locations. An idealized geometry and pulsating inflow conditions through sinusoidal velocity profiles are employed. The results may help identify the most vulnerable locations of the aneurysm initiation in the CoW region. The ideal geometry will help generate standard reference data to assess the severity of possible variations such as underdeveloped or completely absent blood vessels. The reduced ability for collateral blood flow via the CoW in such a case can be easily quantified.

## 2 Details of Fluid–Structure Interaction Model

A 3D FSI analysis is carried out on an idealized Circle of Willis geometry with uniform wall thickness. A transient, laminar, incompressible, non-Newtonian blood flow is considered along with the linear-isotropic arterial wall. FSI model is generated using CFD code ansysfluent 19.0. Finite volume method-based solver is used for the flow domain, whereas the finite element method is used for the structural analysis. FSI solver available in ansys is used for two-way sequential coupling. The solution at the fluid–solid interface is obtained by enforcing the displacement compatibility condition, equal and opposite traction, and no-slip at the interface points. The details of geometry, governing equations, boundary conditions, grid, and the solution methodology are described below.

An ideal model of the circle of Willis is considered for the study. The circle of Willis is a major arterial network that ensures adequate blood supply to different lobes of the human brain. The CoW consists of two arterial systems. The posterior system consists of the vertebral arteries, giving rise to the BA. The basilar artery divides into the PCA. The anterior system has the internal carotid artery dividing into the middle cerebral and ACA. The left and right segments of this network are connected by the anterior communicating artery (ACoA). The arteries of two systems are connected by the posterior communicating artery (PCoA) as shown in the Fig. 1. The various geometrical parameters and mechanical properties considered are shown in Table 1. Symmetry is assumed, and the arterial wall is taken to be of a uniform thickness of 0.45 mm. To simulate the wall weakening effects, geometry is segmented at ACoA, MCA, PCoA, BA, and ICA to yield four different sections (as shown in Fig. 16) where Young's modulus is varied from 1.4 MPa to 0.7 MPa.

Geometric parameters | |
---|---|

Artery | Inner diameter (in mm) |

Vertebral arteries | 3 |

BA | 4 |

PCA | 3 |

PCoA | 2 |

ICA | 5 |

MCA | 3 |

ACA | 2 |

ACoA | 2 |

Mechanical properties | |

Density of wall | 1160 kg/m^{3} |

Poisson's ratio | 0.49 |

Young's modulus | 1.4 MPa |

Young's modulus of diseased sections | 1 MPa and 0.7 MPa |

Geometric parameters | |
---|---|

Artery | Inner diameter (in mm) |

Vertebral arteries | 3 |

BA | 4 |

PCA | 3 |

PCoA | 2 |

ICA | 5 |

MCA | 3 |

ACA | 2 |

ACoA | 2 |

Mechanical properties | |

Density of wall | 1160 kg/m^{3} |

Poisson's ratio | 0.49 |

Young's modulus | 1.4 MPa |

Young's modulus of diseased sections | 1 MPa and 0.7 MPa |

*σ*is given by

The boundary conditions used are as follows: For the solid domain, the lateral end surfaces of all six outlets and four inlets are considered to be fixed. A sinusoidal velocity input with a period similar to a cardiac pulse is considered the fluid domain's inlet boundary condition. The outflow condition is set at the outlets. The transient simulation is performed till the peak velocity is reached at the inlets.

Due to the complex geometry, a tetrahedral mesh is used for both the solid and fluid domains. Inflation is used to refine the mesh near the wall in the fluid domain to capture the boundary-layer effects. The multizone meshing technique is used for the fluid and solid domain due to the complex shape of the circle of Willis. The dynamic meshing option in the solver enabled the capturing of the mesh deformation of the fluid domain in response to the displacement in the solid domain. Figure 2 shows the mesh details at various locations of the circle of Willis.

Simulations are conducted using Newtonian and non-Newtonian blood viscosity assumptions. A dynamic viscosity value of 0.00345 Pa·s is assumed for the Newtonian viscosity model. The non-Newtonian viscosity models considered are the Carreau model, the Carreau–Yasuda model, the Herschel–Bulkley model, the Cross model, and the Power-law model. The details are provided in the validation section. The density of blood is taken as 1050 kg m^{–3}. For all analyses, the arteries of the Circle of Willis are modeled having a uniform density of 1160 kg m^{–3}, while the Young's modulus and Poisson's ratio are taken to be 1.4 MPa and 0.49, respectively. Young's modulus for specific regions prone to aneurysms is varied from 1.4 MPa to 0.7 MPa to study the effects of wall weakening.

The semi-implicit method for pressure linked equations scheme is used for the pressure–velocity coupling as a low time-step does not necessitate pressure-implicit with splitting of operators scheme. The scheme provided faster convergence. Spatial and time discretization are carried out using second-order upwind and implicit schemes. Spatial and time discretization are carried out using second-order upwind and implicit schemes. The solid and the fluid domain are coupled using the System coupling module of ANSYS. The coupling solver achieved stability and convergence with an under-relaxation factor of 1 for incremental displacement and 0.3 for force. Convergence is assumed when all residuals are reduced to six orders of magnitude. In addition, mass flow imbalance is ensured to be less than 0.01%. The HP Z4 G4 Z4 Workstation, with a 24 GB RAM, Intel^{®} Xeon^{®} W-2133 CPU at 3.60 GHz processor, is used for the simulations, with each simulation consuming around 20 h on average.

### 2.1 Grid Independence.

A mesh independence study is carried out separately for the fluid and solid domain to ensure the accuracy of the mesh generated. A sinusoidal velocity input at all six inlets similar to that of the original FSI model is applied for the fluid domain. The radial, axial distribution at a critical location is compared for different mesh size models. Constant pressure is used inside the arteries with no blood flow for the solid domain, and the steady-state solution is calculated. The Von mises stress on the inner wall along a certain length of the artery is extracted and compared for different mesh size models. Based on this study, a mesh with an element size equal to 2 × 10^{−4} m is chosen with 25,43,827 elements is selected for the fluid domain. A mesh of element size equal to 2.5× 10^{−4} m with 1,53,506 elements is chosen for the solid domain, based on the convergence observed. Results from the grid independence study are shown in Fig. 3.

### 2.2 Validation of the Computational Model Developed.

*c*is given by

At the inlet, a rectangular pressure pulse is applied as the driving potential for wave propagation. A Fourier approximation is used to simulate the rectangular pulse. For the simulation, a pipe thickness = 5 mm, radius r = 25 mm, length = 500 mm, and Poisson's ratio = 0.3 with a fluid of density *ρ* = 1000 kg m^{−3} are chosen with Young's modulus varying from 0.25 MPa to 5 MPa. Figure 4 shows the comparison between analytical solution and simulated pressure wave speed. The difference is less than 1% at lower Young's modulus. The selection of the finite length of the tube and the fixed constrain at both the ends of the tube could be the limitations in the model causing the variation at higher Young's modulus.

The CFD model is validated using the experimental data of Anastasiou et al. [35]. The experiment involved a pulsatile flow of an analogue fluid generated at the inlet of a parent channel of diameter D = 600 *μ*m. The parent channel bifurcates into two equally sized channels at an approximate downstream distance of two-dimensional, and both have a diameter = 0.5 D. The velocity profiles of the flow at three different stations were experimentally measured. For the validation of the current CFD model, two pulses are simulated. In the first case, the flow rate pulse varies from 140 to 340 ml/h, whereas the second case changes from 20 to 60 ml/h. Several researchers have used various viscosity models to analyze the non-Newtonian behavior of blood flow accurately. Some of the most frequently used models are the Newtonian, the Carreau model, the Carreau–Yasuda model, the Herschel–Bulkley model, the Cross model, and the power-law model, etc. The control parameters for various model considered are shown in Table 2 [36]. Carreau Yasuda model is implemented using a separate program while other models are already available in ansys solver. In this work, simulations are carried out using various viscosity models and the velocity profile of the flow at a particular station for the second pulse at a given time is shown in Fig. 5. The comparison shows very good agreement between the predictions and experimental data. The overall trend is clearly captured by all viscosity models. Similarly, simulation results are extracted at three different time steps and compared with the experimental results. K-S Goodness values are calculated for the various models and the Carreau model is selected for all further simulations. Comparisons are made at various times, and the Carreau model is selected for further simulations as it predicts the experimental results with the highest accuracy.

Model | Equations | Parameters |
---|---|---|

Power–Law model | $\mu p=\lambda \gamma \u02d9n\u22121$ | λ = 0.042 $Pa\u22c5sn$, n = 0.61 |

Hershel–Bulkley model | $\mu hb={\mu 0|\gamma \u02d9|\u2264\gamma 0\u02d9k|\gamma \u02d9|n\u22121+\tau 0|\gamma \u02d9|\u22121|\gamma \u02d9|\u2265\gamma 0\u02d9$ | k = 0.01345 $Pa\u22c5sn$, n = 0.7812 |

τ_{0} = 0.0216 Pa | ||

γ_{0} = 1 s^{−1} | ||

Carreau–Yasuda model | $\mu cy=\mu \u221e+(\mu 0\u2212\mu \u221e)(1+(\lambda \gamma \u02d9)a)n\u22121a$ | λ = 1.902 s, n = 0.22, a = 1.25, |

μ_{0} = 0.056 $Pa\u22c5s,\u2003\mu \u221e$ = 0.00345 Pa·s | ||

Carreau model | $\mu c=\mu \u221e+(\mu 0\u2212\mu \u221e)(1+(\lambda \gamma \u02d9)2)n\u221212$ | λ = 3.313 s, n = 0.3568 |

μ_{0} = 0.056 $Pa\u22c5s,\u2003\mu \u221e$ = 0.00345 Pa·s | ||

Cross model | $\mu cr=\mu \u221e+(\mu 0\u2212\mu \u221e)11+(\lambda \gamma \u02d9)m$ | λ = 1.007 s, m = 1.028 |

μ_{0} = 0.056 $Pa\u22c5s,\u2003\mu \u221e$ = 0.00345 Pa·s |

Model | Equations | Parameters |
---|---|---|

Power–Law model | $\mu p=\lambda \gamma \u02d9n\u22121$ | λ = 0.042 $Pa\u22c5sn$, n = 0.61 |

Hershel–Bulkley model | $\mu hb={\mu 0|\gamma \u02d9|\u2264\gamma 0\u02d9k|\gamma \u02d9|n\u22121+\tau 0|\gamma \u02d9|\u22121|\gamma \u02d9|\u2265\gamma 0\u02d9$ | k = 0.01345 $Pa\u22c5sn$, n = 0.7812 |

τ_{0} = 0.0216 Pa | ||

γ_{0} = 1 s^{−1} | ||

Carreau–Yasuda model | $\mu cy=\mu \u221e+(\mu 0\u2212\mu \u221e)(1+(\lambda \gamma \u02d9)a)n\u22121a$ | λ = 1.902 s, n = 0.22, a = 1.25, |

μ_{0} = 0.056 $Pa\u22c5s,\u2003\mu \u221e$ = 0.00345 Pa·s | ||

Carreau model | $\mu c=\mu \u221e+(\mu 0\u2212\mu \u221e)(1+(\lambda \gamma \u02d9)2)n\u221212$ | λ = 3.313 s, n = 0.3568 |

μ_{0} = 0.056 $Pa\u22c5s,\u2003\mu \u221e$ = 0.00345 Pa·s | ||

Cross model | $\mu cr=\mu \u221e+(\mu 0\u2212\mu \u221e)11+(\lambda \gamma \u02d9)m$ | λ = 1.007 s, m = 1.028 |

μ_{0} = 0.056 $Pa\u22c5s,\u2003\mu \u221e$ = 0.00345 Pa·s |

## 3 Results and Discussions

Variation in the geometry and branch asymmetries may be amplified by the local WSS/flow characteristics in the circle of Willis region and can initiate aneurysms. Therefore, WSS and flow features of various locations are analyzed in detail in Secs. 3.2 and 3.4. In addition, simulations are carried by changing the elastic properties at common locations for IA such as BA, ICA, posterior communicating artery (PCoA) and the MCA. The bifurcation at the tip of the basilar artery, locations near to ACoA and PCoA origin, and a section in the MCA artery away from the MCA-ICA junction are considered in this study. Comparisons are done for the normal and changed wall properties to discriminate the locations based on their hemodynamic properties. The various hemodynamic parameters are compared at the peak of the sinusoidal velocity profile.

### 3.1 Comparison of Newtonian and Non-Newtonian Simulation.

Simulations are done to compare the effects of Newtonian and the Carreau non-Newtonian viscosity model. Near the wall, due to variations in shear, viscosity changes are steeper for the non-Newtonian model and result in boundary layer thickness changes and an increased pressure drop. However, as the vessel diameters are large, the non-Newtonian effects are not significant. This is evident from the mass flow distribution shown in Fig. 6. The mass flow rates in different branches are expressed as a percentage of total mass flow into the geometry. The flow through PCoA is low as most of the blood from ICA flows to MCA directly. No difference is noticed with the paired cerebral arteries, and the highest flow rate is observed for MCA, followed by PCA. It may be noted that the current predicted values may be overestimated as flow through minor arteries branching out from the various locations of CoW are neglected in the present model. The role of communicating arteries seems to be less critical in an ideal geometry, though it might become important with missing or occluded vessels. The direction of flow may be an indication of the CoW side with an occluded or absent artery. Therefore, flow measurement through communicating arteries alone can be used as a marker of such conditions. The difference in blood viscosity is noticed only for PCoA. Figure 7 shows the radial velocity distribution for the PCoA midpoint. The *y*-axis is nondimensionalized using the centerline velocity of the corresponding Newtonian case. The changes in mass flow rates cause changes in velocity distribution. Literature shows that Newtonian models generally underpredict pressure drop, and the effect may be different during various stages of the cardiac cycle [37]. In general, non-Newtonian models generated lower flow rates, weaker but more symmetrical vortices, and large regions of low WSS compared to Newtonian model simulation.

*I*for comparing non-Newtonian viscosity effects. It is defined as

_{L}the ratio of viscosity predicted by the non-Newtonian and Newtonian models. The global non-Newtonian importance factor for the domain, is calculated by considering many nodes for various time steps. It is found that the value is slightly higher than one at the peak fluid flow, though close to 1 as the flow decelerates. Figure 8 shows the local *I _{L}* variation at ACoA, PCoA, and BA terminus during peak flow. It is noted that the non-Newtonian effects are significant at locations where the recirculation zone is expected. The

*I*value near the wall is higher than other nodes due to significant viscosity changes near this region.

_{L}The radial velocity distribution near the bifurcations at PCoA and ACoA are shown in Fig. 9. The Newtonian model predicted higher velocities for both cases. The general flow pattern is not affected much, though the non-Newtonian model predicts a more flat profile. These changes may affect wall shear stress and strength of the recirculation regions present near the bifurcations. The slow flow near the PCoA could cause localized stagnation regions, leading to atherogenesis and other cerebral diseases. A significant change in velocity prediction can be seen in this region from Newtonian and non-Newtonian models. Previous studies [39] have demonstrated that the viscosity model selection is essential to accurately predict the wall shear stress near bifurcations in cerebral arteries. Figure 10 shows the WSS predictions at PCoA and ACoA bifurcations using the Newtonian and Carreau non-Newtonian viscosity model. Higher WSS around the circumference is noticed for the non-Newtonian model. Therefore non-Newtonian blood viscosity models may be more appropriate for modeling the cerebral circulation region. It is also observed that the mass flow imbalance is lesser in non-Newtonian cases. Therefore Carreau non-Newtonian (shear-thinning) model is used for all further simulations.

### 3.2 Wall Shear Stress Analysis.

Wall shear stress, especially at locations such as BA, posterior communicating artery (PCoA), and the anterior communicating artery (ACoA) are discussed in detail. These locations correspond to high bifurcation angles and curvatures, which may influence WSS distribution. Nixon et al. [40] pointed out that high WSS level and WSS gradient can lead to structural aberrations similar to intracranial saccular aneurysm initiation based on in vivo models. WSS distribution at different locations with the typical elastic wall is shown in Fig. 11. The contours are generated from the transient FSI simulation of a sinusoidal velocity input of magnitude equal to 0.2 m/s and a period equal to 1.6 s. The distribution corresponds to *t* = 0.4 s at the first peak of the sinusoid. The locations with the highest and lowest WSS in the figure are also marked as L1 and L2. ACoA has a maximum of 37 Pa followed by 17.2 Pa in ICA bifurcation and 8.2 Pa and 6.4 Pa in BA and PCoA, respectively. In general, high WSS regions mainly appear at reduced cross-sectional areas and sharp turns. Highest WSS is observed at the bifurcation near the ACoA and is around 37 Pa. The cyclic variation with high peak WSS bifurcation walls may affect endothelial cells. Many studies [41,42] have demonstrated that the reorganization of the cellular cytoskeleton may happen due to WSS exposure on endothelial cells. Therefore, the location at ACoA is a potential candidate for aneurysm formation. Low WSS suggests that stagnation of the flow can result in local inflammatory effects leading to wall thickening/thinning. High WSS may have effects on endothelium and can promote degradative processes [43]. Circumferential variation of WSS is seen throughout the region. A comparison of the circumferential variation at various locations of the highest WSS is shown in Fig. 12. The variation is highest at ACoA, followed by MCA and BA. The circumferential variation will create asymmetric WSS distribution and result in deviations in the vessel's shape and high sustained WSS levels at specific sites. These sites could act as an aneurysm initiation and growth. In summary, the high levels and asymmetric variation of WSS are maximum for ACoA followed by MCA. This might be the reason for more occurrences of cerebral aneurysms at ACoA and MCA. These findings are in agreement with the clinical results of Kayembe et al. [44]. The data by Kwak et al. [45] regarding 485 subjects showed that 213 patients had the aneurysm site at ACoA. It was also reported that 80% of the patients had an abnormal circulation in the CoW. In addition, WSS oscillates between a maximum and minimum value during the systolic and diastolic phases of the cardiac cycle with or without direction reversals. This can trigger aneurysm formation in the above-mentioned critical regions. A review by Sadasivan et al. [46] also showed that around 70–75% of intracranial aneurysms occur at three locations, viz., the bifurcation of MCA, ACoA, and ICA junction followed by BA terminus.

Figure 14 shows the RRT variation at the critical locations. RRT shows similar behavior as that of OSI and highly localized. Maximum RRT value is observed in PCoA with 47.93, followed by 26.45 in MCA, 10.57 in ACoA, and 9.42 in BA.

### 3.3 Wall Elasticity Effects.

The effect of changes in arterial wall elasticity is discussed in this section. The Young's modulus of the arterial wall near ACoA, as shown in the Fig. 15, is analyzed. The wall shape changes during the sinusoidal inlet profile's peak velocity are expressed for Young's modulus of $1.4\xd7106$ case. The displacements are magnified by 15 times for clarity in visualization. The percentage denotes the ratio of displacement and radius of the section for 50% of Young's modulus case. It is noted that the cross section decreases with a decrease in Young's modulus. Smooth muscle fibers make up the walls of muscular arteries. The muscles allow these arteries to expand and contract. These changes in size control how much blood moves through the arteries. This change in cross-sectional area is more in a weakened artery when compared to a normal artery. Simulations are conducted using various sinusoidal input profiles. The overall displacement values decrease with a decrease in velocity magnitude; however, the effect of wall elasticity is found similar for each velocity profile. Average WSS at the weakened region decreased with a reduction in wall elasticity. The percentage change in WSS is around 7% and 8% when Young's modulus decreased to 70% and 50%, respectively. This is in agreement with the results obtained by Xu et al. [23], who also observed the decrease of WSS along the aneurysm wall with the reduction of wall elasticity. The results were based on an experimental study conducted with an elastic aneurysm wall.

The young's modulus of the arterial wall at MCA, PCoA, BA, and ICA is varied from 1.4 × 10^{6} Pa to 1 × 10^{6} Pa, and 0.7 × 10^{6} Pa, and its effect on various hemodynamic parameters are analyzed. Figure 16 shows the different regions where the changes in young's modulus are incorporated.

During the cardiac cycle, the radial arterial motion is estimated to be around 3% [49] with respect to the mean radius for the aorta and around 0.5% for carotid arteries. In this study, any additional movement that the changes in wall elasticity might create is examined for the various locations. Table 3 shows the radial wall movement in percentage with variation in Young's modulus of the arterial wall. The changes are not significant for BA and ICA. PCoA showed the highest variation with an average change of around 3% for a 50% change in Young's modulus. MCA showed an average variation of 1%.

Angle with respect to lower-most point, clockwise (deg) | Y = 1.4 × 10^{6} Pa (% change) | Y = 1.0 × 10^{6} Pa (% change) | Y = 0.7 × 10^{6} Pa (% change) | |
---|---|---|---|---|

BA | 0 | −0.25 | −0.23 | −0.22 |

108 | 0.21 | 0.21 | 0.23 | |

216 | 0.32 | 0.33 | 0.36 | |

324 | −0.24 | −0.23 | −0.22 | |

PCoA | 0 | 2.97 | 3.05 | 3.12 |

108 | −2.65 | −2.72 | −2.80 | |

216 | −2.70 | −2.78 | −2.86 | |

324 | 3.015 | 3.09 | 3.17 | |

ICA | 0 | −0.14 | −0.20 | −0.29 |

108 | −0.25 | −0.33 | −0.45 | |

216 | −0.14 | −0.21 | −0.31 | |

324 | −0.23 | −0.31 | −0.42 | |

MCA | 0 | −0.50 | −0.75 | −1.12 |

108 | −0.74 | −0.93 | −1.225 | |

216 | −0.85 | −1.07 | −1.39 | |

324 | −0.31 | −0.54 | −0.87 |

Angle with respect to lower-most point, clockwise (deg) | Y = 1.4 × 10^{6} Pa (% change) | Y = 1.0 × 10^{6} Pa (% change) | Y = 0.7 × 10^{6} Pa (% change) | |
---|---|---|---|---|

BA | 0 | −0.25 | −0.23 | −0.22 |

108 | 0.21 | 0.21 | 0.23 | |

216 | 0.32 | 0.33 | 0.36 | |

324 | −0.24 | −0.23 | −0.22 | |

PCoA | 0 | 2.97 | 3.05 | 3.12 |

108 | −2.65 | −2.72 | −2.80 | |

216 | −2.70 | −2.78 | −2.86 | |

324 | 3.015 | 3.09 | 3.17 | |

ICA | 0 | −0.14 | −0.20 | −0.29 |

108 | −0.25 | −0.33 | −0.45 | |

216 | −0.14 | −0.21 | −0.31 | |

324 | −0.23 | −0.31 | −0.42 | |

MCA | 0 | −0.50 | −0.75 | −1.12 |

108 | −0.74 | −0.93 | −1.225 | |

216 | −0.85 | −1.07 | −1.39 | |

324 | −0.31 | −0.54 | −0.87 |

A representative figure (magnified by a scale of 50:1) showing the cross-sectional variation at location 3 of MCA is shown in Fig. 17 during the peak velocity instant of the inlet sinusoidal wave. The percentage change in the radial direction for the lowest Young's modulus case is also shown. The percentage denotes the ratio of displacement and the initial radius of the section for each young's modulus case. In general, the radial variation and asymetricity increase with a decrease in Young's modulus. As the cross-sectional changes due to elasticity are highest for MCA, coupled with high flow rates compared to other regions, this location may contribute to bigger aneurysms once initiated due to wall weakening. The geometrical variations created by the changes in elasticity can alter the hemodynamic forces around the bifurcations.

Table 4 compares centerline velocity at various locations with different wall elasticity levels. The comparison is done with the rigid wall simulations. It is seen that the core velocity increases as the wall elasticity decreases.

% change in core velocity | |||
---|---|---|---|

Locations | Y = 1.4 × 10^{6} Pa | Y = 1.0 × 10^{6} Pa | Y = 0.7 × 10^{6} Pa |

BA | 0.19 | 0.19 | 0.23 |

ICA | 0.00 | 0.21 | 0.36 |

PCoA | 3.28 | 5.88 | 6.13 |

MCA | 0.56 | 0.63 | 0.74 |

% change in core velocity | |||
---|---|---|---|

Locations | Y = 1.4 × 10^{6} Pa | Y = 1.0 × 10^{6} Pa | Y = 0.7 × 10^{6} Pa |

BA | 0.19 | 0.19 | 0.23 |

ICA | 0.00 | 0.21 | 0.36 |

PCoA | 3.28 | 5.88 | 6.13 |

MCA | 0.56 | 0.63 | 0.74 |

Changes in the arterial wall thickness are also examined for the various elasticity cases. The changes are observed to be less than 1%, though the percentage change doubles as Young's modulus is decreased by 50%. This study demonstrates the importance of considering the deforming walls instead of rigid walls in cerebral artery simulations. Though the gross flow features are not significantly affected, differences in secondary flow structures can be noticed. The wall thickness variation (around 300 *μ*m at parent vessels) may also be considered for more accurate simulation.

### 3.4 Streamlines and Vorticity.

Figure 18 shows the snapshot of streamlines at the peak velocity instant of the sinusoidal profile. Color contours represent the magnitude of flow velocity nondimensionalized by the inlet velocity at the same moment. The highest velocity is seen at MCA. The velocity at both communicating arteries is considerably low.

Figure 19 shows the magnified view of the streamlines and velocity vectors at BA tip, ACoA, and PCoA. BA tip is characterized by high velocities, subsequent rapid deceleration, and local rise in pressure. Sites characterized by high pressure, flow stagnation, and lower shear stress are more vulnerable to aneurysm rupture. Near the impingement zone, the flow bifurcates and then joins the mainstream. This acceleration causes high-velocity gradients and high WSS. A counter-rotating vortex formation can be noticed at ACoA. The shape of the geometry at ACoA is responsible for this vortex formation. This region creates lower shear stress, high oscillation intensity, etc., and is prone to cell destruction [50]. Asymmetry in the A2 segment of the anterior cerebral can also affect the vertical structure and enhance aneurysm formation risk. Many researchers [51,52] have also noticed the vulnerability of bifurcation regions toward aneurysm formation. The reasons for variation in aneurysm neck size might be associated with the resulting WSS generated by these complex flow patterns. The driving force for wide-neck aneurysms could be high WSS associated with wall impingement and degradation, whereas low WSS could initiate narrow-necked aneurysms. The WSS at BA tip increased to approximately 2.9 times when the peak velocity of the sinusoidal input profile is doubled. This indicates that any changes in flow rates can significantly affect the shear rate and the risk of aneurysm initiation.

The non-Newtonian viscosity models showed differences in viscosity compared to the Newtonian model near separation points. Such differences could affect the prediction of the size of the flow circulation regions. Newtonian assumption can also lead to the overestimation of hemodynamic stress. Therefore, it is recommended to use non-Newtonian viscosity models to simulate blood flow in cerebral arteries.

Figure 20 shows the prominent vortex formation locations in the entire circle of Willis. The vorticity contour for the whole of the region is shown. The highest vorticity magnitude is found near the ACoA region, and the contours are nondimensionalized based on the highest vorticity magnitude. The four sites identified are (i) MCA-ICA bifurcation, (ii) ACoA, (iii) PCoA, and (iv) at the bifurcation region of BA. The vorticity changes with lower Young's modulus are examined at these locations to study the effect of wall weakening. It is observed that vorticity increases upstream of the weakened wall, whereas it decreases downstream of the weakened section. The impact of wall weakening on vorticity magnitude is highest at ACoA with a 50% change in Young's modulus resulting in a 3% change in vorticity.

### 3.5 Von Mises Stress.

Previous researchers have noticed that both Von Mises and wall shear stress react with the arterial wall to alterations in wall properties. Von Mises stress (VMS) represents the deformation energy of the arterial wall and can be used as an index to mark the risk regions. VMS is generally used as a criterion for the failure of the arterial wall and represents the global distortion experienced by the wall. The surface distribution of VMS during the peak sinusoidal velocity in the complete circle of Willis region is shown in Fig. 21. L1 represents the highest VMS in the area. ACoA has a maximum of 3505 Pa followed by 1740 Pa in MCA and 1347 Pa and 1192 Pa in PCoA and BA, respectively. Locations shown with maximum von mises stress are close to areas where aneurysms are usually reported. The current trend is similar as ACoA has the highest chances of aneurysm formation followed by MCA, PCoA, and BA.

The highest value of VMS at ACoA and PCoA is highly localized and confined to a small area. Peak values are observed at diametrically opposite points, as evident from the circumferential variation plot (Fig. 22) and suggests that this could be the most likely point of aneurysm initiation. VMS localization results from the wall curvature effects, which change the fluid pressure acting on the wall.

The effect of wall elasticity is analyzed for change in VMS for the ACoA region. Average VMS at the weakened region increased with a decrease in wall elasticity. The percentage change in VMS is around 6% and 7% when Young's modulus decreased to 70% and 50%, respectively. The VMS at MCA increased to approximately 3.4 times when the peak velocity of the sinusoidal input profile is doubled. This shows that any changes in flow rates can significantly affect the rupture risk of aneurysms.

## 4 Limitations of the Study

In this study, numerical simulations of an ideal circle of Willis geometry are analyzed quantitatively in terms of the WSS, VMS, and other parameters. All blood vessels are assumed as ideal tubes with predetermined curves. In reality, the circle of Willis exhibits asymmetry and anatomical variations. Therefore, the accurate analysis may require a patient-derived detail [44,53]. Another drawback regarding the geometry is the neglect of minor arteries branching out from the various locations of CoW. This is done to make the simulation computationally manageable. The mechanical interaction between the blood flow and the arterial wall is taken into account. The wall Young's modulus is varied to accommodate variation in wall elasticity. The variation in wall elasticity that can be present in the different arterial wall layers is not considered. The current analysis has technical limitations, such as its inability to capture wall thickness variations and generate solutions for the thin high elastic walls. Another simplifying assumption used is the inlet blood flow condition. Instead of a pulsatile waveform similar to the actual cardiac cycle, an ideal sinusoidal waveform is applied. This is done to avoid convergence issues associated with the numerical solver used. The simplifying assumptions used in the blood viscosity is as follows. Researchers in hemodynamic simulations have used various rheological models which incorporate viscoelasticity and shear thinning properties of blood. The shear-thinning effect is considered by including the non-Newtonian blood viscosity model, though the viscoelastic models are not considered. Other factors that can be added to the current model are (i) anisotropic properties of the wall, (ii) turbulence, etc.

Therefore, this study can be considered a preliminary exploration of the aneurysm initiation risk at various locations of the circle of Willis. In addition, the results are plotted using nondimensional parameters wherever possible to obtain qualitative information.

## 5 Conclusions

A computational fluid–structure interaction model of the complete CoW is developed using an idealized geometry. The work aims at highlighting the aneurysm initiation risk of various locations of the CoW. Laminar, incompressible, non-Newtonian blood flow is modeled using Navier–Stokes equations. The arterial wall is modeled as a linear elastic, isotropic with different elastic coefficients to capture the wall elasticity effects. Simulations are completed for many sinusoidal velocity pulses to compare the impact of transient pulsatile flow. The study provides a qualitative and quantitative assessment of flow characteristics, arterial wall mechanics, and wall elasticity effects. It is observed that the non-Newtonian model may be more accurate for CoW simulations as there are changes in flow rates, WSS, and velocity profiles with Newtonian and non-Newtonian simulations. The role of communicating arteries is less critical in an ideal geometry. However, changes in blood flow rates and their direction could be an indication of missing or occluded vessels. The highest WSS is noticed at the anterior communicating artery junction followed by the middle cerebral artery and posterior communicating artery. These sites also showed a circumferential variation of WSS, resulting in deviation in the vessel's shape. Therefore, these sites are more prone to aneurysm initiation than any other location of the CoW. Locations near the posterior communicating artery have the highest chances of local thinning and dilatations. As the flow rates through this region are high compared to different locations, the chances of aneurysm growth may be higher once initiated. Another notable feature is the decrease of WSS with the reduction of wall elasticity. Vorticity is found to increase upstream of the weakened wall, whereas it decreases downstream. The highest VMS is noticed at the anterior communicating artery junction, the middle cerebral artery, and the posterior communicating artery. Wall curvature seems to be the major influencing factor for VMS localization. No significant changes are found for the VMS with wall elasticity changes. Results suggest that the risk of an aneurysm is higher due to changes in blood flow rate. The study indicates that the aneurysm initiation risk factor may be higher in the anterior circulation than in the posterior circulation regions. Results suggest that the anterior communicating artery and tip of the basilar artery could be the most common locations in the anterior and posterior regions, respectively. The highly localized VMS at ACoA and PCoA could cause a faster rupture compared to other sites. In summary, the study has highlighted various hemodynamic features of an ideal geometry circle of Willis with normal and weakened arterial walls. The study shows hemodynamic stress variations in the circle of Willis, which plays a role in the initiation and growth of cerebral aneurysms. As the majority of the population differs from a “textbook” circle of Willis, the influence of such deviation may be considerable and may attribute to high aneurysm risk. The use of a regular geometry model of CoW in this study is understandably an idealized scenario. However, the use of this model helps understand the flow characteristics and wall stress responses from an evolutionary viewpoint. The hemodynamic changes that occur secondary to alterations in the flow pattern and blood vessel anatomy will be better understood in comparison to this basic CoW model. This study forms the benchmark for comparing patient-specific anatomy and variations in blood flow patterns.

## Conflict of Interest

There is no conflict of interest.