NON-NEWTONIAN COMPUTATIONAL HEMODYNAMICS IN TWO PATIENT-SPECIFIC CEREBRAL ANEURYSMS WITH DAUGHTER SACCULES*

WANG Sheng-zhang, CHEN Jia-liang, DING Guang-hong

Department of Mechanics and Engineering Science, Fudan University, Shanghai 200433, China,

E-mail: szwang@fudan.edu.cn

LU Gang, ZHANG Xiao-long

Center of Interventional Therapy, Huashan Hospital Affiliated to Fudan University, Shanghai 200040, China

NON-NEWTONIAN COMPUTATIONAL HEMODYNAMICS IN TWO PATIENT-SPECIFIC CEREBRAL ANEURYSMS WITH DAUGHTER SACCULES*

WANG Sheng-zhang, CHEN Jia-liang, DING Guang-hong

Department of Mechanics and Engineering Science, Fudan University, Shanghai 200433, China,

E-mail: szwang@fudan.edu.cn

LU Gang, ZHANG Xiao-long

Center of Interventional Therapy, Huashan Hospital Affiliated to Fudan University, Shanghai 200040, China

(Received September 16, 2009, Revised April 27, 2010)

Hemodynamic factors play important roles in the formation, progression and rupture of cerebral aneurysms, and the Wall Shear Stress (WSS) and Oscillatory Shear Index (OSI) on the aneurysms are considered to be correlated with their growth and rupture. In this article, two computational models based on patient-specific cerebral aneurysms with daughter saccule are constructed from 3D-RA image data, one is lateral aneurysm located in middle cerebral artery (CA1) and the other is terminal aneurysm located in anterior communicating artery (CA2), The corresponding models of the two aneurysms by removing daughter saccule are established in order to investigate the initiation and growth of the daughter saccule. The flow patterns and the distributions of hemodynamic factors in the two aneurysms before and after daughter saccule is removed are obtained by solving the governing equations with the commercial CFD software Ansys CFX11.0 under the non-Newtonian fluid assumption. By analyzing the flow patterns, it is concluded that the aneurysms with daughter saccules have more complex and unstable flow patterns and hence are prone to rupture. By comparing the distribution of OSI, a hypothesis that a high OSI causes the growth of the daughter saccule is presented.

cerebral aneurysm, computational hemodynamics, non-Newtonian, patient-specific

1. Introduction

Cerebral aneurysms are the pathologic dilations of the arteries, generally found in the anterior and posterior regions of the Willis circle. Their most serious consequences are their rupture and intracranial hemorrhage, with an associated high mortality and morbidity rate[1]. Classic treatments of cerebral aneurysms are direct surgical clipping or endovascular coils occlusion. The coils promote blood coagulation inside the aneurysm, thereby avoiding blood flow and excluding the aneurysm from the circulation[2]. Because prognosis of subarachnoid hemorrhage is still poor, preventive surgery is increasingly considered as a therapeutic option. However, every treatment carries a risk, which sometimes matches or exceeds the yearly risk of aneurysm rupture. Therefore, the best patient care would be to treat only those aneurysms that are likely to rupture. Furthermore, the remnant neck of the coiled cerebral aneurysm has the risk of reoccurrence, which is eagerly expected to overcome by the neurosurgeons[3]. Planning elective surgery is essential for the therapies of cerebral aneurysms, but it requires a better understanding of the process of aneurysm formation, progression and rupture so that one can make a sound judge of the risks and benefits of possible therapies. Unfortunately, these processes have not been well understood until now. Hemodynamic factors, such as the blood velocity, Wall Shear Stress (WSS), pressure, particle residencetime and flow impingement, play important roles in the growth and rupture of cerebral aneurysms. Therefore, studying the hemodyanmics of the cerebral aneurysm is significant for assessing the rupture risk and the reoccurrence risk[4]. CFD simulation for blood circulation has many advantages, For example, the model can be modified easily, the cost is low and the computations can be repeated quickly. Therefore, 3-D CFD has been increasingly used to study cerebral blood flow, including flow in cerebral aneurysms.

Steinman et al.[5]reported the image-based computational simulations of the flow dynamics for Newtonian fluid in a giant patient-specific cerebral aneurysm. CFD analysis revealed high-speed flow entering the aneurysm at the proximal and distal ends of the neck, promoting the formation of both persistent and transient vortices within the aneurysm sac. Cebral et al.[6]performed CFD analysis in cerebral aneurysms from CTA and 3DRA image data. They discussed the limitations and difficulties of in vivo image-based CFD. The methodology can be used to study possible correlations between intra-aneurismal flow patterns and the morphology of the aneurysm and eventually the risk of rupture. Hassan et al.[7]reported also a method for the reconstruction of cerebral vessel from images obtained from MRA, CTA, or 3DRA by using iso-surfacing technique, and found that the ruptured areas are correlated with the areas of high fluid-induced WSS. Valencia et al.[8]constructed patient-specific models of cerebral aneurysms of different sizes located in the ophthalmic artery from 3DRA image data. They found a relation between the aneurysm aspect ratio and the mean WSS on the aneurismal sac. Cebral et al.[9]constructed 62 patient-specific models of cerebral aneurysms and performed CFD simulations with these models, and found some association between hemodynamic factors and the ruptured aneurysms.

The influence of non-Newtonian properties of blood has be reported by Perktold et al.[10]and they investigated an idealized arterial bifurcation model with a saccular aneurysm and found there is no essential difference in the results with both fluid models in the regions with relatively low velocities. Valencia et al.[11]reported that the effect of the non-Newtonian properties of blood on the WSS are important only in the arterial regions with high velocity gradients on the aneurismal wall and the predictions with the Newtonian and non-Newtonian blood models are similar. Fisher et al.[12]claimed that the non-Newtonian behavior pf blood has pronounced effects on flow and fluid mechanical forces within the aneurysm.

The effect of arterial compliance has been considered in some recent reported studies of hemodynamics in arteries with aneurysms. Low et al.[13]examined with their models of lateral aneurysms the effects of distensible arterial walls, and found that the increase and decrease of the flow velocity reflect the expansion and contraction of the aneurysm wall where the maximal wall displacement during systolic acceleration is about 6% of the aneurysm diameter. Torii et al.[14]investigated the influence of the arterial-wall deformation on the hemodynamic factors by carrying out computational fluid-structure interaction analysis, and found various patterns of this influence depending on the arterial geometry. Valencia et al.[15]described the flow dynamics and arterial wall interaction with a representative model of a terminal aneurysm of the basilar artery, and compared its wall shear stress, pressure, wall deformation with those of a healthy basilar artery when the arterial wall was assumed to be elastic or hyper-elastic, isotropic, incompressible and homogeneous. Chen et al.[16]studied the variations of hemodynamic factors in one patient-specific cerebral aneurysm after endovascular coiling with the fluid-structure interaction simulation by employing the commercial software systems Ansys and CFX. However, the difficulty is how to reasonably determine the properties of the arterial wall, such as elastic modulus, wall thickness, etc..

According to the clinical observations and statistics, cerebral aneurysms with daughter saccules have more rupture risks than those without daughter saccules, and the rupture always occur on the daughter saccules[16]. In order to investigate the ruptures of the daughter aneurysms, the blood flow dynamics simulations and analysis for cerebral aneurysms with daughter saccule are very necessary. In this article, two patient-specific models of cerebral aneurysms with daughter saccules, located at middle cerebral artery (CA1) and at anterior communicating artery (CA2), are constructed from their 3DRA images data firstly. Table 1 demonstrates the geometrical parameters in the two cerebral aneurysm models, where D denotes the diameter of the parent vessel and AR denote the aspect ratio (depth/neck width) of the aneurysm. Then, computational fluid dynamics is employed to simulate the pulsatile blood flows in the two cerebral aneurysms with the non-Newtonian fluid assumption. Finally, the flow patterns inside the two aneurysms and the distributions of the hemodynamic factors such as wall shear stress, oscillatory shear index, etc., on the wall of the two aneurysms are calculated and analyzed in order to find the relationship between the growth and the rupture ofthe daughter saccule and the hemodynamics.

2. Methods

2.1 Geometrical reconstruction and grid generation

The Phillips Integris Allua system was applied to obtain the cerebral angiography by 180orotational scanning of the cerebral arteries with aneurysms. The 3-D geometry of cerebral arteries with aneurysms was reconstructed from the corresponding 100 projection images at the Phillips 3D-DSA work station. By choosing a suitable Region Of Interest ( ROI ) window, the proper cerebral arteries with aneurysms were obtained and exported in the VRML format to 3DSMax 8.0 software for adjusting the scale, and exported in the STL format. The geometry of arteries with cerebral aneurysm was imported from the STL file to Geomagic Studio 9.0 for smoothing and repairing. On the other hand, in order to compare the hemodynamcis in the cerebral aneurysms with and without daughter saccule, two corresponding models without daughter saccule were created in the following way: (1) Removing the daughter saccule from the cerebral aneurysm, (2) Patching the surface of the cerebral aneurysm in terms of the surface curvature. Then, the geometrical models of the cerebral aneurysms with and without daughter saccule were generated.

The 3-D geometrical model of arteries with aneurysms was imported to the professional mesh generator software ICEM CFD 11.0 in the STL format. An unstructured CFD grid was generated in the fluid domain with tetrahedral cells, and the number of cells was controlled by the maximum element size.

2.2 Computational model

Blood flow in cerebral aneurysms is assumed to be incompressible unsteady flow and satisfies the mass and momentum conservation equations

where ρ is the density of the blood, v the velocity vector, p the pressure andτ the stress tensor equal to

where μ is the dynamic viscocity of the blood, γ˙ijthe strain rate tensor which is defined for incompressible fluid as

Blood behaves as a Newtonian fluid and μ is constant when it flows in tubes that are greater than about 1 mm in diameter and when it flows with rates of shearing strain greater than 100 s-1, however the blood flows very slowly inside the aneurysm sac and the corresponding shear strain rate is much smaller than that inside the arteries, thus the blood inside the cerebral aneurysm behaves is much closer to non-Newtonian fluid[17]. The Carreau-Yasuda blood model predicts decreasing viscosity at higher shear rate and increasing viscosity at lower shear rate and is selected to describe the non-Newtonian behavior of blood in this study and the apparent viscosity varies according to the law

Fig.1 The variation of dynamic viscosity in the Carreau-Yasuda model with the shear rate

where μ0and μ∞r(nóng)epresent the asymptotic values of the viscosity at low and high strain rate and parameters K and n control the transition region. Based on Cho and Kensey’s results, the parameters have typical valuesμ∞=0.00345Pa· s , μ=0.056Pa· s , K=10.976s2and n = ?0 .3216[18].

Figure 1 shows the variation of the blood dynamicviscosity with the shear rate. The density of the blood is assumed to be 1 050 kg/m3. The walls of parent arteries and the aneurysms are assumed rigid mainly due to a lack of information on the material properties of the wall, such as wall thickness and elasticity.

Fig.2 The input velocity waveform

2.3 Boundary conditions

Physiologic boundary conditions are imposed by prescribing fully developed, time-dependent velocity profiles. The pulsating mean velocity profile in the middle cerebral artery was measured by transcranial doppler ultrasound and is prescribed by the Fourier series as

where ω is the fundamental frequency of the flow waveform. Figure 2 displays the input velocity waveform used in this study. The period of the cardiac cycle is 0.75 s. The Womersley solution in a straight pipe is used to prescribe the pulsatile velocity profile at the inlet of the models as

Table 2 Variation of averaged WSS of CA1 for steady blood flow on five grids

The computational workstation in this study is HP xw6400 with Xeon 3.0GHz CPU (dual cores), 8 GB RAM memory, and Windows XP operation system. The computation time for the two cases based on 2 consecutive pulsatile flow cycles were approximately 40 CPU hours (CA1) and 28 CPUhours (CA2), and the results in the second cycle were selected to be analyzed.

2.5 Postprocessing

Besides velocity and pressure obtained by solving the governing equations, other hemodynamic factors are calculated. In view that the blood flow is pulsatile, the Time Averaged Wall Shear Stress (TAWSS) is defined as

According to this definition, the OSI varies between 0 and 0.5, and the higher OSI region indicates flow direction changes over the cardiac cycle.

3. Results

3.1 Results in CA1 aneurysm with and without daughter saccule

Figure 3 depicts the streamlines in the CA1 aneurysms before and after daughter saccule is removed, and the color of the streamlines represents the magnitude of the velocity. The streamlines are spiral in the inflow artery because this artery is curved, and the most streamlines are helixes and flow pattern is complicated inside the sac. Actually, a small part of blood flows to the outflow artery from the inflow artery directly, another small part flows out to the small outflow branch which is located on the aneurysm, and the remaining blood flows into the sac from the distal side of the aneurismal neck, which has a strong impingement on the distal side of the aneurysm and forms impinging region, and flows out from the proximal in side and constructs the helix streamlines. The flow side the daughter saccule has a very weak circulation. The outflow artery bifurcates into two branches and the angle of bifurcation is larger than 120o, therefore the velocity of blood flowing into the right branch is much higher than that into the left one and the maximal velocity is located at the beginning of the right branch.

Fig.3 Streamlines in CA1 aneurysm before (the upper row) and after (the lower row) the daughter saccule is removed at peak systole

Fig.4 Velocity fields and contours in the cross section of the aneurismal neck of CA1 before the daughter saccule is removed. The left figure represents the peak systole and the right one represents the diastole

Fig.5 Velocity fields and contours in the cross section of the aneurismal neck of CA1 after the daughter saccule is removed. The left figure represents the peak systole and the right one represents the diastole

Figure 4 depicts the velocity fields in the cross section of the aneurysmal neck of CA1 at the peak systole and at the diastole (left), and the color represents the magnitude of velocity. There are two vortexes in both figures, the weak one is located on the right-up side of the cross section, close to thesmall outflow branch, and the strong one is located on the left-middle side of the cross section where is close to the impinging region. Although the vortexes at the diastolic instant are much weaker than those at the peak systole, the locations where the two vortexes exist are almost the same. Figure 5 depicts the velocity fields in the cross section of the aneurismal neck of CA1 without daughter saccule at the peak systole (right) and at the diastole (left). Comparing Fig.5 to Fig.4, one can find that they are almost same. Therefore, removing the daughter saccule does not change the flow field at the aneurismal neck.

Figure 6 depicts the distribution of TAWSS on the wall of CA1 aneurysm before (left) and after (right) daughter saccule is removed. TAWSS is larger on the downstream vessel than that on the other surface and the maximum value is about 24 Pa and is close to the bifurcating apex of the downstream vessel. It is lower than 5 Pa on the most surface of the aneurismal sac, and it is higher near the root of the daughter saccule than on the other area. When the daughter saccule is removed, the distribution of TAWSS on CA1 does not change so much.

Fig.6 TAWSS contours on the wall of CA1 before (left) and after the daughter saccule is removed

Figure 7 depicts the distribution of OSI on the wall of CA1 aneurysm before (left) and after (right) daughter saccule is removed. OSI is very small (close to zero) on the most surface of CA1 aneurysm, and the maximum value is located on the top of the daughter saccule and is about 0.5. It is larger on the daughter saccule than on the other surface of the aneurismal sac. When the daughter saccule is removed, the area, where OSI is larger than 0.25, is enlarged on the aneurismal sac.

Fig.7 OSI contours on the wall of CA1 before (left) and after (right) the daughter saccule is removed

3.2 Results in CA2 aneurysm with and without daughter saccule

Figure 8 depicts the streamlines in the In CA2 aneurysm before and after daughter saccule is removed. The streamlines are almost straight in the inflow artery, and the blood flows into the sac from the proximal side of the aneurysm, and the maximal velocity is located at the inflow tract of the aneurysm. Because CA2 is the aneurysm at bifurcated artery and has two outflow arteries, a part of blood flows out to the right outflow arteries directly and the remaining part of blood flows into the sac and flows out from left outflow arteries, which constructs helixes inside the sac. Only very little blood flows into the daughter saccule and it forms a small circulation.

Fig.8 Streamlines in CA2 aneurysm before (the upper row) and after daughter saccule (the lower row) is removed at peak systole

Fig.9 Velocity fields and contours in the cross section of the aneurismal neck of CA2 before the daughter saccule is removed. The left figure represents the peak systole and the right one represents the diastole

Figure 9 depicts the velocity fields in the cross section of the aneurismal neck of CA2 at systolic instant and diastolic instant. There are two vortexes, one is located on the upside of the cross section, close to one outflow artery and the other is located on the downside of the cross section, close to the other outflow artery and the impinging region. The vortex atpeak systole is much stronger than that at diastole. Figure 9 depicts the velocity fields in the cross section of the aneurismal neck of CA2 without daughter saccule at peak systole and at diastole. Comparing Fig.10 to Fig.9, one can see that the two figures are almost the same except that the velocity in the cross section of the aneurismal neck with daughter saccule is a little larger than that without daughter saccule.

Fig.10 Velocity fields and contours in the cross section of the aneurismal neck of CA2 after the daughter saccule is removed. The left figure represents the peak systole and the right one represents the diastole

Fig.11 TAWSS contours on the wall of CA2 before (left) and after (right) the daughter is removed

Figure 11 depicts the distribution of TAWSS on the wall of CA2 aneurysm before (left) and after (right) daughter saccule is removed. TAWSS on the daughter saccule is less than 1 Pa, smaller than that on the aneurysm. The maximum value is about 10 Pa and it is located near the aneurismal neck. When the daughter saccule is removed, the distribution of TAWSS on the aneurysm does not change so much.

Fig.12 OSI contours on the wall of CA2 before (left) and after (right) the daughter saccule is removed

Figure 12 depicts the distribution of OSI on the wall of CA2 aneurysm with (left) and without (right) daughter saccule. OSI on the most surface of CA2 aneurysm with daughter saccule is very small (close to zero) except two small regions, and the maximum OSI is on the top of the daughter saccule and the value is about 0.4. When the daughter saccule is removed, then OSI on the aneurysm increases obviously.

4. Conclusions and discussion

the role of hemodynamics in the initiation and rupture of the cerebral aneurysms.

In the future work, we plan to select more cerebral aneurysms with daughter saccule from the clinic and to perform the computational hemodynamics modeling and simulation. We plan to analyze the numerical results by using the statistical method to draw some quantitative conclusions.

[1] THE CARAT INVESTIGATORS. Rates of delayed rebleeding from intracranial aneurysms are low after surgical and endovascular treatment[J]. Stroke, 2006, 37(6): 1437-1442.

[2] SCHIEVINK W. I. Intracranial aneurysms[J]. The New England Journal of Medicine, 1997, 336(24): 28-40.

[3] SFORZA D. M., PUTMAN C. M. and CEBRAL J. R. Hemodynamics of cerebral aneurysms[J]. Annual Review of Fluid Mechanics, 2009 41: 91-107.

[4] HUMPHREY J. D., TAYLOR C. A. Intracranial and abdominal aortic aneurysms: similarities, differences, and need for a new class of computational models[J]. Annual Review of Biomedical Engineering, 2008, 10: 221-246.

[5] STEINMAN D. A., MILNER J. S. and NORLEY C. J. et al. Image-based computational simulation of flow dynamics in a giant intracranial aneurysm[J]. American Journal of Neuroradiology, 2003, 24(4): 559-566.

[6] CEBRAL J. R., CASTRO M. A. and APPANABOYINA S. et al. Efficient pipeline forimage-based patientspecific analysis of cerebral aneurysm hemodynamics: Technique and sensitivity[J]. IEEE Trans. Med. Imaging, 2005, 24(4): 457-467.

[7] HASSAN T., TIMOFEV E. V. and SAITO T. et al. Computational replicas: Anatomic reconstructions of cerebral vessels as volume numerical grids at three-dimensional angiography[J]. American Journal of Neuroradiology, 2004, 25(8): 1356-1365.

[8] VALENCIA A., BOTTO S. and SORDO J. et al. Comparison of haemodynamics in cerebral aneurysms of different sizes located in the ophthalmic artery[J]. Int. J. Num. Mech. Fluids, 2007 53(5): 793-809.

[9] CEBRAL J. R., CASTRO M. A. and BURGESS J. E. et al. Characterization of cerebral aneurysms for assessing risk of rupture by using patient-specific computational hemodynamics models[J]. American Journal of Neuroradiology, 2005, 26(10): 2550-2559.

[10] PERKTOLD K., PETER R. and RESCH M. Pulsatile non-Newtonian blood flow simulation through a bifurcation with an aneurysm[J]. Biorheology, 1989, 26(6): 1011-1030.

[11] VALENCIA A., ZARATE A. and GALVEZ M. et al. Non-Newtionian blood flow dynamics in a right internal carotid artery with a saccular aneurysm[J]. Int. J. Num. Mech. Fluids, 2006, 50(6): 751-764.

[12] FISHER C., ROSSMANN J. S. Effect of non-Newtonian behavior on the hemodynamics of cerebral aneurysm[J]. Journal of Biomechanical Engineering, 2009, 131(9): 1-9.

[13] LOW M., PERKTOLD K. and RAUNIG R. Hemodynamics in rigid and distensible saccular aneurysms: A numerical study of pulsatile flow characteristics[J]. Biorheology, 1993, 30(3-4): 287-298.

[14] TORII R., OSHIMA M. and KOBAYASHI T. et al. Influence of wall elasticity in patient-specific hemodynamic simulations[J]. Computer and Fluids, 2007, 36(1): 160-168.

[15] VALENCIA A., SOLIS F. Blood flow dynamics and arterial wall interaction in a saccular aneurysm model of the basilar artery[J]. Computers and Structures, 2006, 84(21): 1326-1337.

[16] CHEN Jia-liang, WANG Sheng-zhang and DING Guang-hong et al. Patient-specific blood dynamic simulations in assessing endovascular occlusion of intracranial aneurysms[J]. Journal of Hydrodynamics, 2009, 21(2): 271-276.

[17] OHTA M., WETZEL S. G. and DANTAN P. et al. Rheological changes after stenting of a cerebral aneurysm: A finite element modeling approach[J]. Cardiovasc Intervent Radiol, 2005, 28(6): 768-772.

[18] CHO Y. I., KENSEY K. R. Effects of the non-Newtonian viscosity of blood on flows in a diseased arterial vessel. Part I: Steady flows[J]. Biorheology, 1991, 28(3-4): 241-263.

[19] TAYLOR C. A., HUGHES T. J. R. and ZARINS C. K. Finite element modeling of three-dimensional pulsatile flow in the abdominal aorta: Relevance to atherosclerosis[J]. Ann. Biomed. Eng., 1998, 26(6): 1-14.

[20] HE X., KU D. N. Pulsatile flow in the human left coronary artery bifurcation: Average conditions[J]. J. Biomech. Eng., 1996, 118(1): 74-82.

[21] WANG Qing, WANG Wei-zhe and FEI Zhi-min et al. Simulation of blood flow in intracranial ICA-PCOMA aneurysm via computational fluid dynamics modeling[J]. Journal of Hydrodynamics, 2009, 21(5): 583-590.

[22] CHATZIPRODROMOU I., TRICOLI A. and POULIKAKOS D. et al. Hemodynamics and wall remodeling of a growing cerebral aneurysm: A computational model[J]. Journal of Biomechanics, 2007, 40(2): 412-426.

[23] SHOJIMA M., OSHIMA M. and TAKAGI K. et al. Magnitude and role of wall shear stress on cerebral aneurysm computational fluid dynamics study of 20 middle cerebral artery aneurysms[J]. Stroke, 2004, 35(11): 2500-2505.

[24] MENG H., FENG Y. and WOODWARD S. H. et al. Mathematical model of the rupture mechanism of intracranial saccular aneurysms through daughter aneurysm formation and growth[J]. Neurological Research, 2005, 27(5): 459-465.

[25] PENTIMALLI L., MODESTI A. and VIGNATI A. et al. Role of apoptosis in intracranial aneurysm rupture[J]. Journal of Neurosurgery, 2004, 101(6): 1018-1025.

10.1016/S1001-6058(09)60098-6

* Project supported by the National Natural Science Foundation of China (Grant Nos. 30772234, 30870707), the Shanghai Municipal Natural Science Foundation (Grant No. 08ZR1401000).

Biogaraphy: WANG Sheng-zhang (1976-), Male, Ph. D., Lecturer