Numerical Simulation of Flow Behavior in Basilar Bifurcation Aneurysms Based on 4-Dimensional Computed Tomography Angiography

Tomoaki Yamazaki¹, Gaku Tanaka¹, Ryuhei Yamaguchi^2*, Yodai Okazaki³, Hitomi Anzai², Fujimaro Ishida⁴, Makoto Ohta²

¹Graduate School of Engineering, Chiba University, Chiba, Japan.
²Institute of Fluid Science, Tohoku University, Sendai, Japan.
³School of Engineering, Tokyo Institute of Technology, Yokohama, Japan.
⁴Department of Neurosurgery, Mie Chuo Medical Centre, Tsu, Japan.

DOI: 10.4236/wjm.2021.114006                

Abstract

Initiation, growth, and rupture of cerebral aneurysms are caused by hemodynamic factors. It is extensively accepted that the cerebral aneurysm wall is assumed to be rigid using computational fluid dynamics (CFD). Furthermore, fluid-structure interactions have been recently applied for simulation of an elastic cerebral aneurysm model. Herein, we examined cerebral aneurysm hemodynamics in a realistic moving boundary deformation model based on 4-dimensional computed tomographic angiography (4D-CTA) obtained by high time-resolution using numerical simulation. The aneurysm of the realistic moving deformation model based on 4D-CTA at each phase was constructed. The effect of small wall deformation on hemodynamic characteristics might be interested. So, four hemodynamic factors (wall shear stress, wall shear stress divergence, oscillatory shear index and residual residence time) were determined from the numerical simulation, and their behaviors were assessed in the basilar bifurcation aneurysm.

Keywords

Basilar Bifurcation Aneurysm, 4-Dimensional Computed Tomographic Angiography, Moving Boundary Method

1. Introduction

The initiation, growth and rupture of cerebral aneurysm are thought to be caused by hemodynamic factors. Hemodynamic factors such as the wallshear stress (WSS), the studies have examined the hemodynamics of cerebral aneurysms using numerical simulation to predict their growth and rupture. It had been extensively accepted to use rigid-wall assumption in computational fluid dynamics [1] [2] [3] [4]. On the other hand, to reproduce the wall movement under cardiac cycle, fluid-structure interaction simulation has been recently applied for the deformation of an elastic cerebral aneurysm model [5] [6] [7]. However, it is difficult to validate these results in another modality. Furthermore, it is still difficult to reproduce the wall deformation of the aneurysm wall in pulsatile flow. Torii et al. reported that there is significant difference of WSS distribution between rigid-wall and elastic-wall simulation. In the previous study [8], the influence of wall motion on wall shear stress is estimated by the similar method to the present study and there is small influence of wall motion on wall shear stress.

The aim of the present study was to exhibit the hemodynamics of basilar bifurcation aneurysms (BBA) using numerical simulation in a realistic moving boundary deformation model obtained by high time-resolution 4-dimensional computed tomographic angiography (4D-CTA) data. First, the wall of the realistic moving deformation model was constructed based on 4D-CTA images. A moving boundary wall was then created based on these images, and used as the wall boundary condition for simulation. Four hemodynamic factors, i.e. WSS, WSSD, OSI and RRT, were assessed by numerical simulation using the moving boundary method. Consequently, the authors have clarified the effect of moving boundary wall in comparison with rigid model.

2. Method

Basilar bifurcation is one of predilection sites of cerebral aneurysm formed a lump-like bulge at cerebral artery, and when it ruptures, it leads to subarachnoid hemorrhage. Construction of the basilar bifurcation aneurysm (BBA) morphology model was based on 4D-CTA data, and the deformation behavior of cerebral aneurysm was taken 30 times per beat. CTA data was measured in Mie Chuo Medical Center using 16-detector multislice acquisition CT scanner (Toshiba, Inc., Tokyo, Japan). The morphology of the aneurysm was constructed based on a patient-specific model is shown in Figure 1 for each phase. Amira ver 5.2.2 (Thermo Fisher Scientific Co., Tokyo, Japan) was used as the morphological creative application. In this study, a moving boundary method is created using an actual shape model and given as a boundary condition for the wall surface. However, since the 30 shape morphologies based on 4D-CTA is not enough to be used as a boundary condition, 10 intermediate shape morphologies interpolated between actual shape morphology models at each phase t and t + Δt were also created. As a comparison, rigid model was selected as the morphology at the 8^th time division based on 4D-CTA data comparable to an average volume from 30 morphologies during one cardiac cycle.

STAR-CCM + Ver14.06 (SIEMENS) was used to create the intermediate shape. This displacement amount during time step Δt was multiplied by n/11 (n is integer). First, the displacement $\Delta x_i=\left(\Delta x_i,\Delta y_i,\Delta z_i\right)$ between shape model 1 (time t) and model 2 (time t + Δt) is obtained, and the vector function $x_{n,i}$ expressed in Equation (1) was created.

$x_{n,i}=\frac{n}{11}\Delta x_i=\left(\frac{n}{11}\Delta x_i,\frac{n}{11}\Delta y_i,\frac{n}{11}\Delta z_i\right)$ (1)

By giving the created vector function of the node at wall surface of shape model 1, the node moves and the shape data of intermediate shape 1 (IS-1) is created. This IS-1 has a shape slightly closer to model 2 than model 1. In this study, we set n = 1, 2, ∙∙∙, 9, 10 and created intermediate shapes between each phase, creating a total of 300 intermediate shapes. In the simulation of moving boundary model, the shape morphology is re-meshed in each time step with changing cell number and created at t/T = 0.00 and 0.20 as shown in Figure 2.

Figure 1. DICOM data and morphology of the basilar bifurcation aneurysm.

Figure 2. Mesh form. (Through one cardiac cycle, average mesh number is approximately 110,000.)

3. Simulation

STAR-CCM + Ver14.06 (SIEMENS) was used for the simulation. The analysis conditions were three-dimensional, implicit, unsteady, laminar flow and non-slip condition on wall, and the moving boundary and rigid walls conditions created based on the actual shape model were applied to the wall surface. The properties of working fluid was assumed to bethe same as the property of blood; Newton fluid, density ρ = 1057 (kg/m³), viscosity coeficientthe same as the property of blood; Newton fluid, density ρ = 1057 (kg/m³), viscosity coefficient μ = 3.7 × 10⁻³ (Pa∙s). The mesh reference length was 0.07 - 0.08 (mm) and the discrete time step was Δt= 1/330 sec. The periodic mean Reynolds number Re_mean and the Womersley number α were defined by the following equations and were 300 and 1.81, respectively.

${\mathrm{Re}}_{mean}=\frac{\rho U_{mean}d}{\mu }$ (2)

$\alpha =\frac{d}{2}\sqrt{\frac{2\pi \rho }{\mu T}}$ (3)

where the mean inflow velocity measured in previous study [9] is U_mean = 0.39 (m/s), the inlet vessel diameter is d = 2.69 (mm), and one period is T = 1.0 (s).

Laminar flow model was used because the flow situation is low Reynolds number. The outlet boundary condition was set to maintain the constant flow division ratio into two outlets. The governing equations are the following continuity equation and the Navier-Stokes equation considering the mesh movement velocity.

$\nabla \cdot u=0$ (4)

$\rho \left[\frac{\partial u}{\partial t}+\left\{\left(u-u_{mesh}\right)\cdot \nabla \right\}u\right]=-\nabla p+\mu {\nabla }^{2}u$ (5)

where u (m/s) is the flow velocity vector, $u_{mesh}$ (m/s) is the mesh moving wall velocity vector, p (Pa) is the pressure and t (s) is the time. The simulation was performed through 3 cardiac cycles, and data were collected from the 2^nd cycle.

The pulsatile flow waveform shown in Figure 3 was used as the inflow boundary condition. This is an adjustment of the pulsatile waveform based on previous studiesin middle cerebral artery [9] [10] [11].

The following quantities were used to evaluate the results obtained by simulation under the above conditions. The first is the magnitude of the resultant wall shear stress WSS (Pa), the second is the divergence of wall shear stress WSSD (Pa/m), the third is the oscillation shear index OSI and the fourth is the residual residence time RRT (Pa⁻¹) as defined following expressions.

$\text{WSS}=\sqrt{{\text{WSS}}_x^2+{\text{WSS}}_y^2+{\text{WSS}}_z^2}$ (6)

$\text{WSSD}=\frac{\partial {\text{WSS}}_x}{\partial x}+\frac{\partial {\text{WSS}}_y}{\partial y}+\frac{\partial {\text{WSS}}_z}{\partial z}$ (7)

Figure 3. Pulsatile flow waveform.

$\text{OSI}=\frac{1}{2}\left(1-\frac{{{\displaystyle \int }}_0^T\text{WSSd}t}{{{\displaystyle \int }}_0^T\text{WSSd}t}\right)$ (8)

$\text{RRT}=\text{1}/\left[\left(\text{1}-\text{2OSI}\right)\text{TAWSS}\right]$ (9)

WSSD reflects the gradient of WSS components in 3 dimensions, and is considered to have a major effect on the pathological diseaseon cerebral aneurysm. OSI indicates the change in the temporal flow direction. RRT denotes the residual residence time at stagnant region.

Firstly, in the process of numerical simulation, the dependency of the mesh reference length (cell number) is examined on the pressure drop between inlet and both outlets. The dependency of the pressure drop between inlet and right-left outlets on mesh reference length is shown in Figure 4 at typical cell number. The difference of pressure drop between cell number 1.1 × 10⁶ and 1.5 × 10⁶ was 0.013% - 0.015%. So, in the present simulation the cell number was approximately 1.1 × 10⁶ in moving boundary model.

4. Results and Discussion

4.1. Flow Pattern and WSS

The typical flow behaviors such as streamline, velocity contour, swirling flow and vortex in moving boundary model are shown in Figure 5 at peak systole. At the front median plane (Figure 5(a)), the main inlet flow collides at the right neck of aneurysm dome. Then, the flow bifurcates to enter inside of the aneurysm as a wall jet along the aneurysm wall, and to flow out to the right outlet vessel. Inside the dome, the global swirling flow along the dome wall surface and the small spiral core vortex flow are induced to be similar to those of the cerebral aneurysm [12]. Also, at the right side plane, large swirling and core vortex flow corresponding to the flow pattern at the front median plane is depicted in Figure 5(b). The comparison of streamline and flow pattern in rigid with moving boundary models is shown in Figure 6 at peak systole. Since the deformation ratio is small, there is little difference between rigid and moving boundary models.

Figure 4. Dependency of mesh reference length (cell number) on pressure drop between inlet and both outlets. (a) Pressure drop from inlet to left outlet; (b) Pressure drop from inlet to right outlet.

Figure 5. Streamline and velocity contour at front median plane and right side plane in the basilar bifurcation aneurysm in moving boundary model at peak systole. (a) Front median plane; (b) Right side plane.

The velocity contour at the entrance of aneurysm dome for rigid and moving boundary models is shown in right and left sides in each phase in Figure 7, respectively. At each velocity contour, the flow with large velocity at right side enters into dome. Conversely, the reverse fluid flows out of left side of dome. The ratio of deformation volume in both models is shown in Figure 8. The maximum volume ratio is 3.0%, i.e. the expansion of depth D of 0.2 mm, is comparable to that [8].

Figure 6. Comparison of streamline and flow pattern at front median plane of rigid with moving boundary models at peak-systole. (a) Rigid model; (b) Moving boundary model.

Figure 7. Velocity at entrance section of aneurysm in rigid and moving boundary models at each phase.

Figure 8. Volume deformation ratio during one cardiac cycle.

The comparison of spatially-averaged WSS in rigid with moving boundary models during one cardiac cycle is shown in Figure 9. WSS in moving boundary is slightly smaller than that in rigid models owing to displacement of aneurysm wall through cardiac cycle. The difference between both models is 0.5 Pa at peak-systole.

Figure 9. Comparison of spatially-averaged WSS in rigid with moving boundary models during one cardiac cycle.

The resultant WSS at four typical phases (early-, mid-, peak-systoles, and mid-diastole) in rigid model is shown in Figure 10. The maximum WSS (denoted in numeral) at each phase is large at the right neck where the inlet flow collides, and reaches 78.6 Pa at peak systole in moving boundary model. In contrast, WSS is globally small at other wall surface, i.e. apart for the aneurysm right neck. WSS in moving boundary is slightly smaller than that in rigid models except for mid-diastole and WSS in rigid model has a similar trend. Comparison of the result [8] with those of present result, the trend is similar to WSS around inlet region larger than that on another surface. Surely, the difference of colliding point depends on morphology.

4.2. WSSD, OSI and RRT

WSSD in both models is shown in Figure 11. At each phase, WSSD at the aneurysm right neck is large compared with that on the global dome surface, and reaches 523 Pa/mm at the flow colliding point in rigid model at peak systole. In moving boundary model, WSSD of 578 Pa/mm at the aneurysm right neck is slightly larger than that in rigid model. WSSD denotes the gradient of WSS, and reflects stretching and compression for the aneurysm wall [3] [13] [14]. WSSD in moving boundary is slightly larger than that in rigid models.

It might be related to the wall displacement, i.e. the aneurysm wall displaces and the flow pattern might be much complicated near wall surface.

The distribution of OSI in both models is shown in Figure 12. In both models, the arrows on surface denote the direction of resultant WSS. The maximum OSI appears at the middle of the front median plane. Three direction vectors in both models are observed around the location of maximum OSI; i.e. the vector at the right side induced by the inlet wall jet, and the vector at the front upper and front lower surfaces induced by the swirling flow and the core vortex flow. Although the maximum OSI of 0.430 in moving boundary is slightly 7% smaller than that (0.460) in rigid models, the spatially- and temporally-averaged OSI in moving boundary is 0.0294, i.e. 6.3 times, larger than that (0.0046) in rigid models. In both models, the vectors of WSS converge around the location of maximum OSI.

The distribution of residual residence time (RRT) is shown in Figure 13. There is a little difference RRT between both models and only there appears the maximum value of 78.9 (Pa⁻¹) and 10.7 (Pa⁻¹) for rigid and moving boundary models, respectively. This small RRT in the moving boundary model implies short residence time of blood flow. When RRT in moving boundary model is 10.7, this aneurysm dome at the middle point on front median plane might not be prone to prevent from rupture. Globally, the specially-averaged RRT of 0.562 (Pa⁻¹) in moving boundary is slightly larger than that of 0.551 (Pa⁻¹) in rigid models. The location of maximum RRT is close to that of maximum OSI. The location with maxi-mum OSI and maximum RRT appear in the middle at front median plane where is different from flow colliding point.

Figure 10. Comparison of WSS in rigid with moving boundary models at four phases. (Numerals denote maximum WSS. Spatially- and temporally averaged WSS in moving boundary and rigid model is 0.326 (−4%) and 0.340, respectively.)

Figure 11. Comparison of WSSD in rigid and moving boundary models at each phases. (Numerals at each panel denote maximum WSSD. For comparison of the rigid with moving boundary models, large and small WSSDs are denoted in red and blue, respectively.) (Spatially- and temporally-averaged WSSD in moving boundary and rigid model is 0.573 [−13%] and 0.657, respectively.)

Figure 12. Distribution of oscillatory shear index OSI. The arrows on surface in both models denote the direction of WSS. (The surface arrows on both models denote the direction of WSS.) (Specially-averaged OSI in moving boundary and rigid models is 0.0294 and 0.0046, respectively.) (a) Rigid model; (b) Moving boundary model.

Figure 13. Comparison of residual residence time (RRT) in rigid with moving boundary models. (Specially-averaged RRT in moving boundary and rigid models is 0.562 [Pa⁻¹] and 0.551 [Pa⁻¹], respectively.)

Our findings indicate that in moving boundary model the WSS might be suppressed and small RRT might be prevented from flow stagnation. WSSD in moving boundary model is larger than that in rigid model and WSSD at right dome neck might be related to stretching and compression for the wall tissue. OSI appears at low WSS region where associates with atherosclerosis of aneurysm dome surface [15] [16]. The spatially-averaged OSI in moving boundary is larger than that in rigid models except for maximum (OSI). On the other hand, RRT in moving boundary is approximately same as that in rigid models.

Although these behaviors cannot be associated with pathological evaluations directly, the location of maximum WSS, maximum WSSD, maximum OSI and maximum RRT might be determined by assessing the hemodynamics of cerebral aneurysms.

5. Summary

In the present study, the effect of moving boundary wall on hemodynamic factors is examined in comparison with rigid wall numerically. Consequently, the temporally- and spatially-averaged WSS around the global dome surface was smaller in comparison with rigid model. Although the maximum residual residence time RRT in moving boundary is smaller than that in rigid models, the spatially-averaged RRT is same as that of rigid models. Both the maximum WSSD at the neck and the global OSI in moving boundary model are larger than those in rigid model.

Acknowledgements

We thank Edanz Group (https://en-author-services.edanzgroup.com/) for editing a draft of this manuscript.

Funding

The present work was mainly supported by the JSPS Science Research Budget #19K04163, and Collaborative Research Projects (J18I108 in 2018, and J19I064 in 2019) of the Institute of Fluid Science, Tohoku University.

Ethical declaration

The present morphology of aneurysm was approved in Mie University, School of Medicine and the patient agreed with application to medical research.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Cebral, J.R., Castro, M.A., Burgess, J., Pergolizzi, R.S., Sheridan, M.J. and Putman, C.M. (2005) Characterization of Cerebral Aneurysms for Assessing Risk of Rupture by Using Patient-Specific Computational Hemodynamics Models. AJNR, 26, 2550- 2559.
[2]	Shojima, M., Oshima, M., Takai, K., Torii, R., Nagata, K., Shirouzu, I., Morita, A. and Kirino, T. (2005) Role of the Bloodstream Impacting Force and the Local Pressure Elevation in the Rupture of Cerebral Aneurysm. Stroke, 36, 1933-1938. https://doi.org/10.1161/01.STR.0000177877.88925.06
[3]	Meng, H., Wang, Z., Hoi, Y., Gao, L., Metaxa, E., Swartz, D.D. and Kolega, J. (2007) Complex Hemodynamics at the Apex of Arterial Bifurcation Induces Vascular Remodeling Resembling Cerebral Aneurysm Initiation. Stroke, 38, 1924-1931. https://doi.org/10.1161/STROKEAHA.106.481234
[4]	Dolan, J.M., Kolega, J. and Meng, H. (2013) High Wall Shear Stress and Spatial Gradients in Vascular Pathology: A Review. Annals of Biomedical Engineering, 41, 1411-1427. https://doi.org/10.1007/s10439-012-0695-0
[5]	Torii, R., Oshima, M., Kobayashi, T., Takagi, K. and Tezduyar, T.E. (2005) Influence of Wall Elasticity in Patient-Specific Hemodynamic Simulations. Computers & Fluids, 36, 160-168. https://doi.org/10.1016/j.compfluid.2005.07.014
[6]	Bazilevs, Y., Hsu, M.C., Zhang, Y., Wang, W., Liang, X., Kvamsdal, T., Brekken, R. and Isaksen, J.G. (2010) A Fully-Coupled Fluid-Structure Interaction Simulation of Cerebral Aneurysms. Computational Mechanics, 46, 3-16. https://doi.org/10.1007/s00466-009-0421-4
[7]	Alexander, A.K., Cherevko, A.A., Chupakhin, A.P., Bobkova, M.S., Krivoshapkin, A.L. and Orlov, Y.-K. (2016) Haemodynamics of Giant Cerebral Aneurysm: A Comparison between the Rigid-Wall, One-Way and Two-Way FSI Models. Journal of Physics: Conference Series, 722, Article ID: 012042. https://doi.org/10.1088/1742-6596/722/1/012042
[8]	Dempere-Marco, L., Oubel, E., Castro, M., Putman, C., Frangi, A. and Cebral, J. (2006) CFD Analysis Incorporating the Influence of Wall Motion: Application to Intracranial Aneurysms. MICCAI 2006 LNCS, Vol. 4191, 438-445. https://doi.org/10.1007/11866763_54
[9]	Harders, A. and Gilsbach, J. (1985) Transcranial Doppler Sonography and Its Application in Extracranial-Intracranial Bypass Surgery. Neurological Research, 7, 129- 140. https://doi.org/10.1080/01616412.1985.11739711
[10]	Valen-Sendstad, K., Mardal, K.A., Mortensen, M., Reif, B.A.P. and Langtangen, H.P. (2011) Direct Numerical Simulation of Transitional Flow in a Patient-Specific Intracranial Aneurysm. Journal of Biomechanics, 44, 2826-2832. https://doi.org/10.1016/j.jbiomech.2011.08.015
[11]	Laleh, Z., Khalid, A., Anders, W., Richard, B., Anders, E. and Jan, M. (2015) Blood Flow Distribution in Cerebral Arteries. JCBFM, 35, 648-654. https://doi.org/10.1038/jcbfm.2014.241
[12]	Yagi, T., Sato, A., Shinke, M., Takahashi, S., Tobe, Y., Takao, H., Murayama, Y. and Umezu, M. (2013) Experimental Insights into Flow Impingement in Cerebral Aneurysm by Stereoscopic Particle Image Velocimetry: Transition from a Laminar Regime. Journal of the Royal Society Interface, 10, Article ID: 20121031. https://doi.org/10.1098/rsif.2012.1031
[13]	Yamaguchi, R., Tanaka, G. and Liu, H. (2016) Effect of Elasticity on Flow Characteristics inside Intracranial Aneurysms. International Journal of Neurology and Neurotherapy, 3, 49. https://doi.org/10.23937/2378-3001/3/3/1049
[14]	Yamaguchi, R., Kotani, T., Tanaka, G., Tupin, S., Osman, K., Shafii, N.S., Khudzari, A.Z.M., Watanabe, K., Anzai, H., Saito, A. and Ohta, M. (2019) Effects of Elasticity on Wall Shear Stress in Patient-Specific Aneurysm of Cerebral Artery. Journal of Flow Control, Measurement & Visualization, 7, 73-86. https://doi.org/10.4236/jfcmv.2019.72006
[15]	Goubergrits, L., Schaller, J., Kertzcher, U., van den Bruck, N., Poethkow, K., Petz, Ch., Hege, H.-Ch. and Spuler, A. (2012) Statistical Wall Shear Stress Maps of Ruptured and Unruptured Middle Cerebral Artery Aneurysms. Journal of the Royal Society Interface, 9, 677-688. https://doi.org/10.1098/rsif.2011.0490
[16]	Malek, A.M., Alper, S.L. and Izumo, S. (1999) Hemodynamic Shear Stress and Its Role in Atherosclerosis. JAMA, 282, 2035-2042. https://doi.org/10.1001/jama.282.21.2035