Interplay of Proximal Flow Confluence and Distal Flow Divergence in Patient-Specific Vertebrobasilar System

Approximately one-quarter of ischemic strokes involve the vertebrobasilar arterial system that includes the upstream flow confluence and downstream flow divergence. A patient-specific hemodynamic analysis is needed to understand the posterior circulation. The objective of this study is to determine the distribution of hemodynamic parameters in the vertebrobasilar system, based on computer tomography angiography images. Here, the interplay of upstream flow confluence and downstream flow divergence was hypothesized to be a determinant factor for the hemodynamic distribution in the vertebrobasilar system. A computational fluid dynamics model was used to compute the flow fields in patient-specific vertebrobasilar models (n = 6). The inlet and outlet boundary conditions were the aortic pressure waveform and flow resistances, respectively. A 50% reduction of total outlet area was found to induce a ten-fold increase in surface area ratio of low time-averaged wall shear stress (i.e., TAWSS ≤ 4 dynes/cm2). This study enhances our understanding of the posterior circulation associated with the incidence of atherosclerotic plaques.


Study design
The purpose of this retrospective study was to investigate the hemodynamic changes in the vertebrobasilar system given a disproportionately high rate of stroke mortality in China [15]. Six patients underwent cerebral CTA at the affiliated hospital of Hebei University, China, to evaluate impairments in vision, body movement, and speaking; unconsciousness; problems with coordination and so on. Table 1 summarizes patient demographics. CTA reconstruction and imaging analysis were performed by researchers at Peking University and Radiologists at the affiliated hospital of Hebei University, which showed no stenoses and aneurysms in patient vertebrobasilar system. This retrospective study was approved by the Institutional Review Board (IRB) for the affiliated hospital of Hebei University. Subjects gave the signed informed consent.

Imaging acquisition
Similar to a previous study [16], all patients underwent CTA scanning from the aortic arch to vertex using the Discovery CT750 HD scanner (HDCT, GE Healthcare, Milwaukee, WI, USA). Briefly, non-enhanced CT brain scan was first performed, which was followed by contrast enhanced CTA. CTA images were acquired when contrast agent (Iopromide 370, Bayer Schering Pharma AG) at the dose of 1.0 ml/kg was injected at a rate of 5 ml/s followed by IV injection of saline chase of 40 ml at a rate of 5 ml/s. A bolus tracking method (Smart Prep) was used to monitor the optimal contrast enhancement. Study parameters included the caudo-cranial scan direction, 120 kVp, rotation time of 0.5 s, 50cm × 50cm DFOV (Display Field Of View), 0.625 mm construction thickness at 0.625 mm intervals, and helical pitch of 0.984:1.

Aortic pressure wave
A patient with coronary artery diseases underwent the examination of invasive coronary angiogram by standard catheterization in accordance with the American College of Cardiology Guidelines for Coronary Angiography [17]. The blood pressure wave was measured by a pressure catheter inserted into the ascending aorta, which was monitored by a pressure control unit (Millar INC, Houston).

3D reconstruction
3D geometry and morphometry of vertebrobasilar arteries were extracted from CTA images using the MIMICS software (Materialise, NV, Belgium), as shown in Fig 1A-1F and Table 2. We focused on the intracranial vertebrobasilar system that is comprised of VAs (LVA and RVA), BA, AICA, SCAs (LSCA and RSCA) and PCAs (LPCA and RPCA). The left and right VAs join at the pontomedullary junction forming the BA which bifurcates into the right and left PCAs as well as right and left SCAs at the pontomesencephalic junction. We excluded some side branches of small diameters that were not observed clearly in the reconstruction (a low CTthreshold of 80 HU). After 3D reconstruction of the vertebrobasilar system, centerlines of the VAs, BA, AICA, SCAs and PCAs were first generated. In the MIMICS software, a centerline was formed by a series of center points which was located in the center on the cross-sectional views of the contour of the 3D vessel. Subsequently, the best fit diameter, D fit , was calculated as twice the average radius between the point on the centerline and the contour forming the 3D vessel.

Mathematical model
Similar to a previous study [18], CFD simulations were performed to analyze the flow patterns in the vertebrobasilar system. The vessel wall was assumed to be rigid and impermeable. The equations of continuity and Navier-Stokes can be written as: where v ! , P, ρ, and μ represent the velocity, pressure, blood mass density, and viscosity, respectively. After the flow computation, Reynolds number is determined as: where v mean and D represent the time-averaged velocity and vessel diameter, respectively.

Method of solution
A finite volume method was applied to solve the governing equations in the ANSYS FLUENT (ANSYS Inc., Canonsburg, USA). Based on these morphometric data, geometrical models were created in the Geomagic Studio software (3D Systems, Rock Hill, USA), which were meshed using the ANSYS ICEM (ANSYS Inc., Canonsburg, USA). Fig 1G-1I show a total of approximately 2.6, 1.8, and 1.3 million hybrid (hexahedral/tetrahedral) shaped volume elements with finer meshes near the vessel wall when the maximal element size was set to 0.35, 0.3, and 0.25 mm, respectively. A mesh dependency was conducted such that the relative error in two consecutive mesh refinements was < 1% for TAWSS and OSI (Fig 2B vs Figs 3A and 4A). A total of approximately 1.8 million hybrid (hexahedral/tetrahedral) shaped volume elements (maximal element size = 0.3 mm) with finer meshes near the vessel wall were necessary to accurately mesh the computational domain.
A measured aortic pressure wave in Fig 2A was applied to the inlet of left and right VAs. The resistance boundary condition was set at each outlet (see the Appendix), based on intraspecific scaling laws of vascular trees that were derived theoretically and validated experimentally in various organs and species [19] including the cerebral circulation [20]. The viscosity and density were set to 4.5×10 −3 PaÁs and 1,060 kg/m 3 , respectively, to mimic the incompressible blood flow with a hematocrit of 45% [21]. Three cardiac cycles were required to achieve the convergence for the transient analysis.

Data analysis
The mean±SD (standard deviation) values were computed for TAWSS (averaged over all nodes on the surface of the vertebrobasilar system). Similar to a previous study [18], hemodynamic parameters, i.e., SAR-TAWSS, SAR-OSI, and SAR-transWSS in a vessel and at a junction (see definitions in the nomenclature), were defined to give integral scalar values instead of traditional local parameters. Similar to previous studies [22,23], a dominant VA was defined as: the diameter of the dominant VA was 0.3 mm larger than another; or the dominant VA straightly connected to the BA if both VAs had similar diameters. ANOVA (SigmaStat 3.5) was used to compare these parameters, where p value < 0.05 represented a statistically significant difference.

Results
3D geometrical models of the vertebrobasilar system were reconstructed from CTA images of six normal patients (i.e., subjects A-F), as shown in Fig 1A-1F. Accordingly, Table 2 lists mean D fit (averaged along the entire vessel length), arch length and chord length in each artery. At the flow convergence with LVA and RVA merging into BA, subjects A and B have a tuning fork geometry with LVA diameter approximately equal to RVA diameter while subject C has a similar geometry, but a dominant LVA. Subjects D-F have a walking geometry with a dominant LVA. Joining LVA to BA forms an S-shaped LVA-BA in subjects D and E and a C-shaped LVA-BA in subject F. These human subjects have similar junctional shapes of flow divergence that is comprised of BA, LSCA, RSCA, LPCA and RPCA despite various geometrical parameters (e.g., length and diameter).   Fig 2C shows the distribution of TAWSS and OSI in subject A when the blood was considered as a non-Newtonian fluid (i.e., Carreau fluid). Fig 2D shows the two hemodynamic parameters when the blood was considered as a Newtonian fluid as well as artificial small sized branches were added to the BA. Fig 3 shows the distribution of TAWSS, which has mean±SD values of 12.7±12.2, 4.9±6.8, 9.9±8.9, 11.5±8. 3, 9.3±8.2, and 7.6±8.1 dynes/cm 2 in subjects A-F, respectively. Fig 4 shows the corresponding distribution of OSI. A comparison between Fig 2C and 2D and Figs 3A and 4A showed negligible effects of non-Newtonian behavior and possible small sized branches in the BA. There were low TAWSS ( 4 dynes/cm 2 ) and high OSI (! 0.15) near the region of high curvature in each vessel. Table 3 lists SAR-TAWSS and SAR-OSI in each vessel of subjects A-F. Subject B had the maximal SAR-TAWSS in LVA, RVA, BA and LPCA, which were significantly higher than others (p value < 0.05), while subject D had the minimal SAR--TAWSS in those arteries. SAR-OSI had a low value (< 0.3%) and SAR-transWSS was zero in each vessel. Moreover, Fig 5A-5F show the complex distribution of streamlines at the time instance with the highest flow velocity at the inlet of VAs (i.e., time equals to 168 ms in Fig 2A).
The upstream flow convergence and downstream flow divergence are frequent sites for atherosclerosis and aneurysm, respectively. Low TAWSS and high OSI occurred near the two sites, as shown in Figs 6 and 7. Table 4 lists SAR-TAWSS and SAR-OSI in the two sites. Subject B had the maximal SAR-TAWSS (p value < 0.05) while subject D had the minimal SAR-TAWSS. Moreover, SAR-OSI had a low value (< 0.3%) and SAR-transWSS was zero in the two sites.

Discussion
Based on 3D geometry and morphometry reconstructed from normal patient CTA images, we carried out a hemodynamic analysis in the vertebrobasilar arterial system that is comprised of  LVA, RVA, BA, AICA, LSCA, RSCA, LPCA and RPCA. The blood was considered as a Newtonian fluid in the analysis given negligible non-Newtonian effects (Fig 2C vs. Figs 3A and 4A), which agreed with previous studies [24,25]. Since Re mean < 400 in Table 3, the blood flow was laminar in vertebrobasilar arteries of these normal patients albeit energetic, rapid fluctuations were observed in CFD models of bifurcation aneurysms [26]. The present study focused on the interplay of upstream flow confluence (LVA and RVA merging into BA) and downstream flow divergence (BA bifurcating into SCAs and PCAs) in the vertebrobasilar system. The corresponding findings were discussed below in turn. Low TAWSS and high OSI (i.e., TAWSS 4 dynes/cm 2 and OSI ! 0.15) are risk factors for the incidence of atherosclerosis and aneurysm [5,6,27,28], which generally occur near curvatures, bifurcations, anastomoses and so on [21,[29][30][31][32]. At the flow confluence of vertebrobasilar system, Ravensbergen et al. showed atherosclerotic plaques at the anastomotic apex and lateral walls of the BA, which were associated with low TAWSS [9,11]. They also found that a blunted apex and a large confluence angle were two geometric risk factors for the occurrence of an atherosclerotic plaque at the BA apex. Smith and Bellon demonstrated at least two flow patterns within the BA due to the non-admixture of vertebral artery flows [13]. Wake-Buck et al. presented the effects of curvatures and their orientations on the flow patterns [14]. In  agreement with the conclusions in Refs. [9,11,13,14], this study showed the distribution of low TAWSS and high OSI at the anastomotic apex and inner side of the curvature. It was, however, difficult to predict the incidence of atherosclerosis by simple geometrical zones (e.g., lateral walls, apex, the site opposite to the flow divider, etc.) given the complex hemodynamics at the upstream flow confluence. Alternatively, SAR-TAWSS, SAR-OSI, and SAR-transWSS were defined to quantify the hemodynamic environment, as shown in Table 4.
A key finding of the study is that SAR-OSI at the flow convergence was zero for subjects A-D (i.e., OSI values at all positions of the flow convergence were smaller than 0.15) and had a low value of 0.14% for subject F while SAR-transWSS was zero (i.e., transWSS values at all positions of the flow convergence were smaller than 6 dynes/cm 2 ) for all patients. Since high OSI resulted from strong reversed flows [29], the flow reversal had small contributions to the incidence of atherosclerotic plaques at the flow convergence. Furthermore, as compared with the tuning fork geometry in subjects A-C, we found stronger spiral flows in the S-shaped or Cshaped walking geometry in subjects D-F (Fig 5A-5C vs. 5D-5F). The spiral flows initiated from the curvature of the dominant VA, which agreed with the computational results in Ref. [14]. The mixing of flows increased the swirling strength and led to the complex distribution of hemodynamic parameters at the flow convergence. The AICA functioned as a shunt path such that it weakened the secondary flows in the BA. Although the secondary flows were weakened at the flow divergence, they continued to propagate in the distal arteries. Mohamied et al. have shown that high transWSS could characterize the multidirectional flows [8]. Because of zero SAR-transWSS for all subjects, the secondary flows could have negligible effects on the incidence of atherosclerosis. It is known that low TAWSS coincided with stagnated, secondary, and reversed flows [29]. Hence, low TAWSS caused mainly by stagnated flows were a risk factor for the incidence of atherosclerosis at the flow convergence. On the other hand, the flow patterns near the flow divergence of vertebrobasilar system, similar to those at a junction in the cardiovascular system of other organs [21,[29][30][31][32], resulted in the non-regular distribution of hemodynamic parameters. Moreover, low TAWSS and high OSI occurred near the aneurysmprone carina due to the oscillatory flows in the direction perpendicular to the BA centerline [31], which led to higher values of SAR-TAWSS and SAR-OSI at the downstream flow divergence than those at the upstream flow convergence.
Another key finding of the study is that the decrease in total outlet area (= ∑ outlet areas of SCAs and PCAs) could increase the vascular resistance to reduce the flow rate that significantly deteriorated the hemodynamic environment (i.e., a significant increase of SAR-TAWSS) in the vertebrobasilar system, as shown in Tables 3 and 4. Here, PCAs only included the P1 segment, which origins at the BA termination to the posterior communicating artery (PCOM), within interpeduncular cistern. Since a grey-scale threshold method with a low CT-threshold of 80 HU was selected to reconstruct the 3D geometry of vertebrobasilar system, this excluded small vessel segments distal to SCAs in the reconstruction. The artificial errors for determination of outlet surface areas were hence negligible. Subject B had the least total outlet area, but the highest SAR-TAWSS in each vessel and flow convergence and divergence of the vertebrobasilar system, as shown in Tables 3 and 4, as well as the lowest mean value of TAWSS. In contrast, subject D with the largest total outlet area had the smallest SAR-TAWSS. This illustrates the significant effects of cerebral microcirculation on the macrocirculation given the scaling relationship between the total outlet area and distal microvasculature [19].

Implications for the posterior circulation
Subject B, the oldest patient, has the highest value of systolic blood pressure, pulse pressure (the difference between the systolic and diastolic pressures), cholesterol and triglycerides.
These cardiovascular risk factors can impair the posterior microvasculature, which, in turn, exacerbates the hemodynamic environment in large arteries of vertebrobasilar system, as shown in Figs 3-7 and Tables 3 and 4. Although a previous study has shown that asymmetric VA flow could be a hemodynamic contributor of BA curvature and peri-vertebrobasilar junctional diseases [22], this study implies that the microcirculation injury should be a significant factor for the incidence and progression of atherosclerosis in large arteries of vertebrobasilar system, which requires further investigations.

Critique of the model
The CFD simulation did not take vessel compliance into account because a previous study has shown negligible effects of vessel compliance on TAWSS and OSI [30]. Moreover, 3D reconstruction of vessels may be affected by imaging parameters and resolution. A study has recently shown that a low CT-threshold of 80 can satisfy the accuracy of 3D reconstruction [33]. The present study did not include human subjects with fPCA, i.e., a PCA arising from the internal carotid artery [34,35], which needs to be considered in the following computational studies. Furthermore, future perspective studies in patients with stenoses and aneurysms should be carried out to validate predictions of the CFD simulation to the incidence and progression of atherosclerosis and aneurysm in posterior circulation.

Conclusions
This retrospective study performed a hemodynamic analysis in patient vertebrobasilar system. The interplay of upstream flow confluence and downstream flow divergence was found to significantly affect the distribution of hemodynamic parameters (e.g., streamlines, TAWSS, OSI and transWSS) in patients of no stenoses and aneurysms. The outlet resistance resulting from the distal microvasculature should be accurately estimated when a CFD simulation is carried out in large arteries of vertebrobasilar system. This study provides insight for understanding of the posterior circulation relevant to the potential incidence of atherosclerosis and aneurysms.
where D outlet refers to the diameter at each outlet of the vertebrobasilar system. Hence, we can get the following equation as: where P outlet refers to the pressure at each outlet of the vertebrobasilar system. Eq [A4] was used as the resistance boundary condition with v outlet and P outlet being the variables.