Published online before print August 2, 2007, doi: 10.1161/STROKEAHA.107.482471
| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
(Stroke. 2007;38:2500.)
© 2007 American Heart Association, Inc.
Original Contributions |
Computation of Hemodynamics in the Circle of Willis
From the Simula Research Laboratory (M.S.A., K-A.M.), Oslo, Norway; the Department of Neurosurgery (J.I., T.I.), University Hospital of North Norway, Tromsø, Norway; the Department of Neurosurgery(B.R., T.I.), Institute for Clinical Medicine, University of Tromsø, Tromsø, Norway; the University Clinic of Neurosurgery (B.R), Neuroscience Centre, Rigshospitalet, Copenhagen, Denmark; and the School of Advanced Medicine (M.K.M.), Macquarie University, Sydney, Australia.
Correspondence to Prof Tor Ingebrigtsen, Department of Neurosurgery, University Hospital of North Norway, N-9038 Tromsø, Norway. E-mail tor.ingebrigtsen{at}unn.no
 |
Abstract |
|---|
 Top Abstract IntroductionMethodsResultsDiscussionReferences |
|---|
Background and Purpose— Wall shear stress (WSS) and pressure are important factors in the development of cerebral aneurysms. We aimed to develop a computational fluid dynamics simulator for flow in the complete circle of Willis to study the impact of variations in vessel radii and bifurcation angles on WSS and pressure on vessel walls.
Methods— Blood flow was modeled with Navier-Stokes equations as an incompressible newtonian fluid within rigid vessel walls. A model of the circle of Willis geometry was approximated as a network of tubes around cubic curves. Pulsatile inlet flow rates and constant outlet pressure were used as boundary conditions.
Results— The simulations confirmed that differences in vessel radii and asymmetric branch angles influence WSS magnitude and spatial distribution. High WSS occurred at locations where aneurysms are frequent and in anatomic variants known to be associated with an increased risk for aneurysm development.
Conclusions— Computational fluid dynamics analysis can be applied to the complete circle of Willis and should be used to study the pathophysiology of this complex vascular structure, including risk factors for aneurysm development. Further development of the method should include simulations with flexible vessel walls.
Key Words: aneurysm • computational fluid dynamics • circle of Willis • hemodynamics • wall shear stress
|
Introduction |
|---|
|
TopAbstractIntroductionMethodsResultsDiscussionReferences |
|---|
Disruption of the internal elastic lamina is required for the creation of saccular aneurysms. Hemodynamic factors play an important role in this process. Saccular aneurysms usually arise at the distal carina of bifurcations, where vessels are exposed to the maximum impact of wall shear stress (WSS).1 The amount of WSS depends on the geometry of the bifurcation.2–5 WSS is minimized when the relation between vessel radii and bifurcation angles follows optimality principles of minimum work.6–8
In the circle of Willis, there is a confluence of flow from 3 vessels: both internal carotid arteries and the basilar artery (BA). Therefore, the hemodynamics in the circle of Willis is anatomically significantly different from the hemodynamics in normal branching situations addressed by the optimality principle. Accordingly, the normal physiology of flow and the likely impact of deviation from normality in the circle of Willis are not fully understood.
In a previous study, we analyzed 3-dimensional digital subtraction angiography images of cerebral vessels with respect to vessel radii and bifurcation angles and concluded that bifurcations beyond the circle of Willis approximated optimality principles, whereas those within the circle of Willis did not.9 In addition, we observed an increased prevalence of aneurysms at bifurcations with large branch angles. Furthermore, studies of this complex vascular structure in patients, animal models, or experimental in vitro models are difficult. Therefore, simulations with computational fluid dynamics (CFD) may contribute to the understanding of this problem. In the present study, we aimed to develop a CFD simulator for flow in the complete circle of Willis to study the impact of variations in vessel radii and bifurcation angles on WSS and pressure on vessel walls.
|
Methods |
|---|
|
TopAbstractIntroductionMethodsResultsDiscussionReferences |
|---|
Geometric Model
Patient-specific images of the circle of Willis were obtained by computed tomography angiography, magnetic resonance angiography, and digital subtraction angiography by standard protocols. A representative model of the geometry was constructed by manually analyzing natural variations in 10 patients and incorporating the known anatomy from textbooks of neurosurgical anatomy.
We have written customized software that generates computational grids from a parameterization of the geometry. Each vessel was approximated as a smooth tube around a sequence of cubic curves, with the potential for radius variation. With a network of such tubes connected by bifurcations, more complicated vascular regions can be described, and geometries with variations in radii and bifurcation angles can be created. Figure 1 shows examples of computational grids.
|
Mechanical Model and Boundary Conditions
Blood was modeled with the Navier-Stokes equations for an incompressible newtonian fluid. The mass density of blood was set to 1 g/cm3 and viscosity was set to 0.0035 Pa · s.10,11 To complete the Navier-Stokes equations, we imposed additional physical conditions on blood flow at the boundary of the computational domain. All vascular walls were considered rigid with a no-slip condition. The inflow boundary was defined as the most proximal cross sections of the 2 internal carotid arteries and vertebral arteries (VAs) and the outflow boundary as the distal cross sections of the 2 first segments of the posterior cerebral artery (P1s), middle cerebral arteries, and anterior cerebral arteries. At the inflow boundary, the velocity was set to a Womersley profile fitted to normalized temporal flow-rate curves.12,13 The average total volumetric flow rate through the circle of Willis was set to 11 mL/s, with 
through the internal carotid arteries (4.125 mL/s per artery) and 
through the VAs (1.375 mL/s per artery).14 On outflow boundaries, the average pressure was set to a constant reference value of 5000 Pa according to natural "do-nothing" boundary conditions.15,16
Numerical Methods
To solve the Navier-Stokes equations, we used the software package FeatFlow (version 1.3, release candidate 3; University of Dortmund, Dortmund, Germany). FeatFlow uses the mixed finite-element method with the Rannacher-Turek element for spatial discretization.17 For time discretization, an adaptive second-order, fractional-step scheme was applied. An algebraic projection method was used, which accurately computes pressure and velocity on the boundary. To resolve the nonlinearities, a fixed-point iteration was applied, which results in a series of linearized algebraic equations. These equations were solved with an efficient geometric multigrid algorithm.18
From the computed velocity field and Newtonian stress description, shear stress on the walls was calculated directly on a finite-element basis. WSS and pressure on the boundary of the geometric model are visualized with color-coded magnitudes according to ParaView software (ParaView software, versions 2.4; Kitware Inc, Los Alamos National Laboratory, Los Alamos, NM).
|
Results |
|---|
|
TopAbstractIntroductionMethodsResultsDiscussionReferences |
|---|
Figure 2 shows the pressure at peak systole from a simulation on the complete 3-dimensional model of the circle of Willis. Blood flow direction is from high to low pressure. Blood flow velocities (not shown) were within normal values.
|
Figure 3 shows results from simulations on 3 different geometric variants of the posterior part of the circle of Willis. The first simulation was based on the complete 3-dimensional model of the circle of Willis (Figure 2), with the model cut at the anterior side of the posterior communicating artery. In the second simulation, the radius of the left P1 was increased from 1.30 to 1.60 mm and the radius of the right P1 was decreased from 1.30 to 1.05 mm. This caused a redistribution of the high-pressure area at the bifurcation apex and increased WSS at the origin of the smaller P1 branch. In the third simulation, the branch angle of the left P1 was reduced from 69° to 26°. This further accentuated the changes observed in the second simulation.
|
Figures 4 and 5
show results from 2 simulations on a simplified model of the posterior part of the circle of Willis. In this 3-dimensional model, the vessels were restricted to 1 plane. In the first simulation, the pressures on the 2 outflow boundaries were equal. In the second simulation, the outflow pressures were set with a difference of 125 Pa. This pressure difference caused a change in flow diversion (Figure 4), resulting in increased WSS on the side of highest flow (Figure 5).
|
|
|
Discussion |
|---|
|
TopAbstractIntroductionMethodsResultsDiscussionReferences |
|---|
The present study shows that simulations with the use of CFD can be applied to studies of hemodynamics in the complete circle of Willis. Specifically, we used this method to investigate the impact of variations in vessel radii and bifurcation angles on pressure and WSS. The simulations indicate that deviations from normal anatomy result in a redistribution of wall pressures and increased WSS at branch points.
The bifurcation of the BA is normally symmetric, both with respect to radii in the 2 P1 segments and to branch angles. In this situation, the simulation showed a symmetric distribution of pressures at the bifurcation angle and low WSS. When we increased the radius of the P1 on 1 side and decreased it on the other side, the simulations showed a redistribution of pressure on the vessel wall and increased WSS at the origin of the smaller branch. The estimated peak values for WSS were >30 Pa, which is considered sufficient to cause disruption of the endothelium and eventual development of an aneurysm.19,20 This was further accentuated in simulations of asymmetric bifurcations. The increased WSS occurred at the origin of the smaller branch and at the origin of the branch with the larger branching angle. These observations are in accordance with our previous findings in a clinical study, wherein a higher prevalence of aneurysms was observed at branch points deviating from optimal bifurcation geometry.9
We also experimented with asymmetric outflow pressures in the 2 P1 segments, thus simulating the situation in patients with arteriovenous malformations, in whom the outflow pressure is lower and the flow higher on the side of the malformation. These simulations indicated increased WSS on the side with the higher flow. This corresponds to the observed increased occurrence of aneurysms on arteries feeding arteriovenous malformations.21,22
Application of CFD methods to studies of the cerebral circulation is in an early phase, and few have attempted to simulate hemodynamics in the complete circle of Willis. Moore et al23 constructed a geometric model similar to ours. They studied blood flow patterns with a focus on the impact of anomalies like a fetal P1 and a missing A1 and concluded that the model might be developed into a clinically useful method. Cebral et al24 developed image-based methods based on magnetic resonance angiography for obtaining accurate patient-specific geometries and applied simplified tree models for the vascular bed supplied by the circle of Willis to improve the outflow boundary conditions. Those studies did not attempt to estimate the impact on WSS and pressures on vessel walls in bifurcations.
Furthermore, CFD methods have been used to simulate intra-aneurysmal hemodynamics in patient-specific geometries based on computed tomography angiography or magnetic resonance angiography.25–28 The results are conflicting concerning the location and magnitude of WSS. Some report simulations indicating a high WSS near the neck of the aneurysm,26,27 whereas others conclude that WSS might be highest near the aneurysm dome.28 A study of patient-specific data obtained by digital subtraction angiography showed that simulations incorporating parent-vessel geometry before the region of interest are preferable.28
In our simulations, flow velocity was defined on each point of the inflow boundaries. The applied Womersley velocity profile is correct only in a long, straight tube. To quantify the impact of the chosen velocity profile, we compared several simulations with a flat ninth-degree profile (results not reported). In a straight tube, flow through 
5 vessel diameters was sufficient to reduce the visible flow difference. In the complete geometric model, there was no flow difference where the 2 VAs meet to form the BA. This suggests that average flow rates are sufficient for accurate simulations, so long as enough geometry is included in the model.13 In the circle of Willis, inclusion of 1 bifurcation before and after the area of interest probably is sufficient. Furthermore, there are numerous small vessels leaving the circle of Willis that were not included in our geometric models. This may affect the boundary layer flow, but the flow through these vessels is low, and we have therefore assumed it would not influence the flow significantly.
Exploratory simulations (not reported) showed that changing the length of the outflow vessels had the same effect as changing outflow pressure. Navier-Stokes equations are sensitive to pressure gradients only and not to absolute values. Differences in pressure between boundaries led to changes in flow division between branches. Thus, the absolute values of the prescribed pressures led to no error, but pressure gradients can do that.
The vascular walls were assumed to be rigid in our simulations. This is a good approximation for smaller vessels (capillaries) but not for larger arteries. For the medium-size vessels in the circle of Willis, this assumption could overestimate WSS but probably not influence its spatial distribution.10,29 Assuming that blood behaves like a newtonian fluid is a good approximation in larger vessels but is incorrect in the smallest vessels, where the size of blood cells approximates the vessel diameter. In medium-size arteries, blood starts to show deviations from this stress description. Cebral et al24 performed simulations with a Casson stress model and obtained the same spatial distribution of WSS with slightly higher values than the newtonian model. During systole, blood velocity is high enough for turbulence to occur (Reynolds’ number >2000). However, turbulence takes some time to develop, and the strong pulsatility of flow decreases the effects of eventual turbulent motion. Turbulence is generally rare in healthy arteries, and for this reason we have not included turbulence models in our simulations. However, turbulence is often associated with atherosclerosis30 and may be of interest when studying particular pathologic conditions.
In conclusion, CFD analysis can be applied to the complex flow-related pathophysiology in the circle of Willis. Our simulations confirm that differences in vessel radii and asymmetric branch angles influence WSS magnitude and spatial distribution. The simulations showed high WSS at locations where aneurysms are frequent and in anatomic variants known to be associated with an increased risk for aneurysm development. To improve the accuracy of the method, further refinement of the model is needed. Of particular importance are more accurate geometry, elastic vessel walls, and nonnewtonian fluid properties, and this work is in progress.
|
Acknowledgments |
|---|
We thank Trygve Karper for participating in writing the first version of the grid software.
Source of Funding
This study was supported by a grant from the Regional Health Authority of North Norway.
Disclosures
None.
Received January 14, 2007; revision received February 27, 2007; accepted February 28, 2007.
|
References |
|---|
|
TopAbstractIntroductionMethodsResultsDiscussionReferences |
|---|
- Stehbens WE. Familial intracranial aneurysms: an autopsy study. Neurosurgery. 1998; 43: 1958–1959.
- Ferguson GG. Physical factors in the initiation, growth, and rupture of human intracranial saccular aneurysms. J Neurosurg. 1972; 37: 666–677.[Medline] [Order article via Infotrieve]
- Forbus WD. On the origin of miliary aneurysms of the superficial cerebral arteries. Bull Johns Hopkins Hosp. 1930; 47: 239–284.
- Roach MR, Scott S, Ferguson GG. The hemodynamic importance of the geometry of bifurcations in the circle of Willis (glass model studies). Stroke. 1972; 3: 255–267.
[Abstract/Free Full Text] - Ujiie H, Liepsch DW, Goetz M, Yamaguchi R, Yonetani H, Takakura K. Hemodynamic study of the anterior communicating artery. Stroke. 1996; 27: 2086–2093.
[Abstract/Free Full Text] - Rosen R. Optimality Principles in Biology. London: Butterworths; 1967: 41–59.
- Rositti S, Löfgren J. Optimality principles and flow orderliness at the branching points of cerebral aneurysms. Stroke. 1993; 24: 1029–1032.
[Abstract/Free Full Text] - Zamir M, Bigelow DC. Cost of departure from optimality in arterial branching. J Theor Biol. 1984; 109: 401–409.[Medline] [Order article via Infotrieve]
- Ingebrigtsen T, Morgan MK, Faulder K, Ingebrigtsen L, Sparr T, Schirmer H. Bifurcation geometry and the presence of cerebral artery aneurysms. J Neurosurg. 2004; 101: 108–113.[Medline] [Order article via Infotrieve]
- Quarteroni A, Formaggia L. Mathematical modeling and numerical simulation of the cardiovascular system. In: Ciarlet PG, Lions JL, eds. Modeling of Living Systems. Amsterdam: Elsevier; 2004.
- Travagli V, Zanardi I, Boschi L, Turchetti V, Forconi S. Evaluation of a torsional-vibrating technique for the hemorheological characterization. Clin Hemorheol Microcirc. 2006; 35: 283–289.[Medline] [Order article via Infotrieve]
- Womersley JR. Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known. J Physiol. 1955; 127: 553–563.
[Free Full Text] - Ford MD, Alperin N, Lee SH, Holdsworth DW, Steinman DA. Characterization of volumetric flow rate waveforms in the normal internal carotid and vertebral arteries. Physiol Meas. 2005; 26: 477–488.[CrossRef][Medline] [Order article via Infotrieve]
- Scheel P, Ruge C, Petruch UR, Schoning M. Color duplex measurement of cerebral blood flow volume in healthy adults. Stroke. 2000; 31: 147–150.
[Abstract/Free Full Text] - Fogarty-Mack P, Pile-Spellman J, Hacein-Bey L, Osipov A, DeMeritt J, Jackson EC, Young WL. The effect of arteriovenous malformations on the distribution of intracerebral arterial pressures. Am J Neuroradiol. 1996; 17: 1443–1449.[Abstract]
- Heywood J, Rannacher R, Turek S. Artificial boundaries and flux and pressure conditions for the incompressible Navier-Stokes equations. Int J Num Methods Fluids. 1996; 22: 325–352.[CrossRef]
- Rannacher R, Turek S. A simple nonconforming quadrilateral Stokes element. Num Methods Part Diff Equations. 1992; 8: 97–111.[CrossRef]
- Turek S. Efficient solvers for incompressible flow problems. In: Lecture Notes in Computational Science and Engineering, vol 6. Berlin-Heidelberg-New York: Springer; 1999.
- Resnick R, Yahav H, Shay-Salit A, Shushy M, Schubert S, Chen L, Zilberman M, Wofovitz E. Fluid shear stress and the vascular endothelium. Prog Biophys Mol Biol. 2003; 81: 177–199.[CrossRef][Medline] [Order article via Infotrieve]
- Fry DL. Arterial intimal-medial permeability and coevolving structural responses to defined shear-stress exposures. Am J Physiol Heart Circ Physiol. 2002; 283: 2341–2355.
- Brown RD Jr, Iebers DO, Forbes GS. Unruptured intracranial aneurysms and arteriovenous malformations: frequency of intracranial hemorrhage and relationship of lesions. J Neurosurg. 1990; 73: 856–863.
- Redekop G, TerBrugge K, Montanera W, Willinsky R. Arterial aneurysms associated with cerebral arteriovenous malformations: classification, incidence, and risk of hemorrhage. J Neurosurg. 1998; 89: 539–546.[Medline] [Order article via Infotrieve]
- Moore S, David T, Chase JG, Arnold J, Fink J. 3D models of blood flow in the cerebral vasculature. J Biomech. 2006; 39: 1454–1463.[CrossRef][Medline] [Order article via Infotrieve]
- Cebral JR, Castro MA, Soto O, Löhner R, Alperin N. Blood-flow models of the circle of Willis from magnetic resonance data. J Eng Math. 2003; 47: 369–386.[CrossRef]
- Oshima M, Torii R, Kobayashi T, Taniguchi N, Takagi K. Finite element simulation of blood flow in the cerebral artery. Comp Methods Appl Mech Eng. 2001; 191: 661–671.[CrossRef]
- Oshima M. A new approach to cerebral hemodynamics: patient-specific modeling and numerical simulation of blood flow and arterial wall interaction. IACM Expressions. 2004; 6: 4–9.
- Shojima M, Oshima M, Takagi K, Torii M, Hayakawa M, Katada K, Morita A, Kirino T. Magnitude and role of wall shear stress on cerebral aneurysm: computational fluid dynamic study of 20 middle cerebral artery aneurysms. Stroke. 2004; 35: 2500–2505.
[Abstract/Free Full Text] - Castro MA, Putman CM, Cebral JR. Computational fluid dynamics modeling of intracranial aneurysms: effects of parent artery segmentation on intra-aneurysmal hemodynamics. Am J Neuroradiol. 2006; 27: 1703–1709.
[Abstract/Free Full Text] - Loots E. Fluid-structure interaction in hemodynamics. PhD thesis, Riksuniversiteit Groningen, Germany. 2003.
- Fung YC. Biomechanics: Circulation. New York: Springer; 1997.
| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||