Abstract
BACKGROUND AND PURPOSE: Computational fluid dynamics simulations of neurovascular diseases are impacted by various modeling assumptions and uncertainties, including outlet boundary conditions. Many studies of intracranial aneurysms, for example, assume zero pressure at all outlets, often the default (“do-nothing”) strategy, with no physiological basis. Others divide outflow according to the outlet diameters cubed, nominally based on the more physiological Murray's law but still susceptible to subjective choices about the segmented model extent. Here we demonstrate the limitations and impact of these outflow strategies, against a novel “splitting” method introduced here.
MATERIALS AND METHODS: With our method, the segmented lumen is split into its constituent bifurcations, where flow divisions are estimated locally using a power law. Together these provide the global outflow rate boundary conditions. The impact of outflow strategy on flow rates was tested for 70 cases of MCA aneurysm with 0D simulations. The impact on hemodynamic indices used for rupture status assessment was tested for 10 cases with 3D simulations.
RESULTS: Differences in flow rates among the various strategies were up to 70%, with a non-negligible impact on average and oscillatory wall shear stresses in some cases. Murray-law and splitting methods gave flow rates closest to physiological values reported in the literature; however, only the splitting method was insensitive to arbitrary truncation of the model extent.
CONCLUSIONS: Cerebrovascular simulations can depend strongly on the outflow strategy. The default zero-pressure method should be avoided in favor of Murray-law or splitting methods, the latter being released as an open-source tool to encourage the standardization of outflow strategies.
ABBREVIATIONS:
- ACA
- anterior cerebral artery
- CFD
- computational fluid dynamics
- OA
- ophthalmic artery
- OSI
- oscillatory shear index
- PcomA
- posterior communicating artery
- TAWSS
- time-averaged wall shear stress
Simulating blood flows by image-based computational fluid dynamics (CFD) has become a prevalent technique for studying the natural history and treatment options for cerebrovascular disorders, notably intracranial aneurysms.1 The accuracy of these nominally patient-specific simulations is subject to not only the CFD methodology2 but also numerous approximations and assumptions. For example, while patient-specific lumen geometries are readily available from routine clinical angiography and used to create the 3D model, properties of blood and the vessel wall are rarely acquired for each patient. Most of the time, rigid walls are assumed and population-averaged constants are used for blood viscosity. When patient-specific flow rates are available, it is typically only at the model inlet.3 More often, inlet flow rates are assumed from literature values, sometimes coupled to a scaling law to account for interindividual differences in flow.4 On the other hand, an assumption that remains particularly latent is that of the outlet flow rates.
Inspection of articles published in the American Journal of Neuroradiology and other clinical and biomedical engineering research journals reveals that the outflow strategy is rarely or only superficially reported. When information is provided, most studies focusing on multi-outlet CFD models report using “zero pressure” or “traction-free” outlet boundary conditions to determine the division of flow among the outlets, a trend also noted by Marzo et al.5 This outflow strategy is often the default (“do nothing”) setting for CFD solvers, and implicit is the assumption that the necessarily truncated cerebrovascular model has all of its outlets connected to the same artery distally, which is rarely the case. (That assumption should not be confused with the more physiological assumption that outlets feed microvascular beds having the same resistance; see the On-line Appendix). As also detailed in the On-line Appendix, subjective—and rarely documented—choices about the CFD models, such as truncation of downstream branches or addition of flow extensions, may therefore play a role in determining the division of outflow rates, which, incidentally, is typically known only after the CFD simulation is completed.6
Another popular outflow strategy is to explicitly apportion the flow rates among the outlets according to their diameters cubed.7 This “Murray law” strategy derives from the eponymous principle that pumping versus metabolic power is optimized for the circulatory system8 and reflects conduit vessels adapting their calibers to the flow rate demanded by the downstream microvascular beds. Murray's law hinges on the simplifying assumption of Poiseuille flow, which does not necessarily hold for larger, conduit arteries9; thus, like the zero-pressure strategy, the predicted flow rates, and hence the CFD simulations, may also be sensitive to the number and extent of outlets retained in the model. (A third, less common approach is to impose microvascular resistances directly at the outlets; however, absent patient-specific flow rates, population-average resistances must be imposed, resulting in population-average flow divisions.)
Therefore, the aim of the present study was 2-fold: 1) to present a novel and more robust alternative, hereafter referred to as the “splitting” method, which uses vessel diameters to estimate patient-specific outflow divisions that are more grounded in physiology and less sensitive to the subjective extent of the CFD model; and 2) to quantify the uncertainty of cerebrovascular CFD models to outflow strategy.
Materials and Methods
Study Cohort
We evaluated outflow strategies for 70 patients with MCA aneurysms (55 ± 11 years of age, 73% women) from a broader cohort of 244 consecutive patients with 358 aneurysms, included from March 2011 to March 2014. These cases were segmented from 3D rotational angiography, from the cervical ICA to at least the M2 branches, using either a threshold segmentation method (Aneufuse; @neurIST, http://www.aneurist.org/) or a gradient-based watershed technique (Matlab; MathWorks, Natick, Massachusetts). Arteries below a threshold diameter of 0.6 mm were not segmented; thus, the ophthalmic artery (OA) was excluded for roughly 25% of the models. The posterior communicating artery (PcomA) was also typically excluded, except for cases with evident fetal posterior cerebral arteries. Detailed morphologic characteristics of this cohort can be found elsewhere.6
These 3D segmentations were used to construct 0D models, on which the various outflow strategies were tested for their impact on branch flow rates. 3D CFD simulations were performed on a subset of 10 representative cases to determine the impact of outflow conditions on derived hemodynamic indices.
Splitting Method
Murray's law is a specific case of a more general principle that the flow rate in a vessel is proportional to the diameter of that vessel raised to some power (ie, Q∝Dn, with n = 3 for Murray's law). As commonly practiced in CFD, the Murray-law strategy is to simply prescribe the flow rate for outlet i as
where Qi is the flow rate of the outlet, Di is the diameter of the outlet, and, customarily, n = 3. The summations are over all outlets, and by conservation of mass, the sum of outflows must equal the inflows. Note that all outlets are treated the same according to their diameters, irrespective of the branch generation to which they belong—that is, no attention is paid to branching within the model.
Our proposed outflow splitting method is based on the same power law relationship between diameter and flow rate, but in the spirit of Murray's original derivation, it is used locally for each bifurcation within the model (eg, see Fig 1, case A). In other words, starting from the inlets, at each bifurcation encountered, we assume the flow divides according to a power law:
where D1 and D2 are the diameters of the 2 daughter branches at the given bifurcation, and Q1/Q2 is the bifurcation flow division. These outflows then serve as the inflows for subsequent bifurcations, until the CFD model outlets are reached, at which point the proportion of the original inlet flow to each outlet is known. A more detailed derivation, with examples, is given in the On-line Appendix. It is important to note that our splitting method does not impose these individual internal flow divisions on the CFD model. Rather, they are used, as detailed below, to compute outlet flow rates that, imposed on the CFD model, will give rise to the expected internal flow divisions.
The 10 representative cases of MCA aneurysms used for 3D CFD. A subset of 3 cases (A–C) was truncated differently, as depicted by A′, B′, C′. These are color-coded to identify their data in subsequent figures. Cases A and A′ highlight their individual bifurcations, as used by the splitting method.
0D Modeling for the Impact of Outflow Strategy on Flow Rate Distribution
Recently, we developed a novel lumped parameter approach for fast and accurate prediction of flow divisions in 3D CFD models when outflow pressures are prescribed.6 Briefly, given a 3D lumen geometry model, a 0D (electrical circuit) model is derived in 2 steps: 1) The centerlines of the model are computed with the Vascular Modeling Toolkit (VMTK; www.vmtk.org); and 2) these centerlines are decomposed in a network of straight, rigid segments parameterized by their mean radius, length, and bifurcation angles. In these reduced-order networks, we assume fully developed, steady, Poiseuille flow so that for each segment, the flow rate is proportional to the pressure gradient along the vessel and the inverse of its resistance. Equations for all segments are assembled and solved iteratively by applying an analog of the Kirchhoff current law to satisfy segment-to-segment mass conservation.10 The method was previously validated for this MCA aneurysm cohort via pulsatile 3D CFD simulations that had prescribed only ICA inflow but used different pressures for each outlet.6 By setting all outlet pressures to zero, this 0D approach was used to calculate the outflow rates for the zero-pressure method.
To determine the outflow rates for the splitting method, we took advantage of the fact that the above-described 3D to 0D process includes the automated identification of bifurcations and measurement of inflow/outflow branch diameters. It was therefore straightforward to adapt this to our proposed splitting method using, for each bifurcation, the power law equation (Equation 2) relating the outflow division to daughter branch diameters and solving simultaneously to obtain the desired outlet flow rates. A power law exponent n = 2 was used at all bifurcation levels. Finally, for the Murray-law method, we simply applied Equation 1, based on the outlet diameters and with a power law exponent n = 3, reflecting common practice.
The 0D models were solved for each case assuming steady flow, with a cross-sectional mean velocity of 0.27 m/s applied to the ICA inlet, under the assumption that inflow rate scales with cross-sectional area4 (ie, a power law with n = 2). The resulting inflow rates are plotted in On-line Fig 1. Outlet flow rates were then calculated according to the 3 different methods (zero-pressure, Murray-law, splitting) and used to calculate the MCA, anterior cerebral artery (ACA), OA, and PcomA flow rates. It is important to remember that, for these MCA aneurysm cases, the parent artery (ie, the MCA) flow rate is effectively determined by the outflow method in addition to the prescribed ICA inflow.
3D CFD Modeling for Impact of Outflow Strategy on Hemodynamic Indices
As shown in Fig 1, a subset of 10 representative cases was selected for 3D CFD to determine the impact of the outflow strategy on nominal hemodynamic predictors of rupture status. Of these, 3 were chosen to exemplify the impact of model extent. Specifically, as also shown in Fig 1, the alterations were either drastic, with branches clipped immediately distal to the aneurysm (case A) or mild, with shortened M2 segments (case B) or a change of the number of outlets and length of the M2 segments (case C).
For these 10 cases and their 3 variants, we performed CFD simulations using a minimally dissipative solver developed and validated within the Open-Source Finite-Element Method Library, FEniCS (https://fenicsproject.org/).11 The segmented lumens were meshed using VMTK, with 4 layers of boundary elements and the mesh density chosen to be highest in the vicinity of the aneurysm sac, where the tetrahedron side length was 0.12 mm on average.12 The number of tetrahedral elements per case was 3.8 million on average, ranging from 2.3 to 4.7 million, reflecting the variability of the aneurysm sizes and extents of the CFD model domains.
We assumed rigid walls, blood kinematic viscosity of 0.0035 m2/s, and blood density of 1060 kg/m3. Like the 0D models, a cycle-averaged cross-sectional mean velocity of 0.27 m/s was applied to the ICA inlet in all cases and was used to scale a representative older adult ICA flow waveform shape13 to impose fully developed Womersley velocity profiles at the ICA inlet. To account for possible dampening along the ICA and MCA,14 we reduced the harmonics of this inflow waveform in amplitude by 10%. The number of time-steps per cardiac cycle was set to 10,000,12 and all simulations were run for at least 3 cardiac cycles to wash away any initial transients. The analysis was based on the output from 1000 uniformly spaced time-steps from the last cycle.
For the zero-pressure CFD simulations, the pressure was set to zero at all CFD model outlets via traction-free boundary conditions. For the splitting method, outlet flow rates were imposed as mass flux conditions to the outlets of the CFD model, which amounts to traction-free boundary conditions (see “Discussion”). CFD simulations for the Murray-law method were not performed.
For each CFD model, the time-averaged wall shear stress (TAWSS) magnitude and oscillatory shear index (OSI) were determined from the time-varying velocity data. The aneurysm sac was digitally isolated,15 and for each permutation of case and outflow method, we computed the following hemodynamic indices often used for rupture status assessment12: sac-averaged TAWSS and sac-maximum TAWSS, both normalized to the parent artery TAWSS; and sac-averaged OSI.
Results
Impact of Outflow Strategy on Flow Rate Distribution
Table 1 demonstrates the overall impact of outflow strategy on flow rates in the major arteries and their accuracy relative to in vivo measurements averaged from the literature. Although all 3 methods predicted MCA flow rates to within <10% on average, the zero-pressure method overestimated ACA flow rates by 23%, causing an imbalance of flow at the ICA terminus. On the other hand, the Murray-law and splitting methods also predicted ACA flow rates to within <10% on average and thus had ACA/MCA flow divisions nearly identical to the in vivo average. For those cases with a PcomA, the zero-pressure method had the largest overestimation relative to in vivo values, while the small OA flow rates were overpredicted by all outflow strategies. SDs (ie, interindividual variances) were largest for the zero-pressure method, followed by the Murray-law method, while the splitting method was closest to, albeit still higher than, the in vivo variances.
For the 70 patients with MCA aneurysms, mean ± SD flow rates (mL/min) in the major arteries and mean flow divisions at the ICA terminusa
Figure 2 shows the marked impact of outflow strategy on flow rates at the MCA, which is the parent artery for these aneurysm cases. Notably, Fig 2A reveals that relative to the more physiological splitting method, the zero-pressure method underestimates MCA flow rates by 14% on average, but >50% for some cases. Figure 2B shows that the Murray-law and splitting methods give closer results on average, but the limits of agreement still have a wide range and thus large differences for individual cases, from +50% to −70%.
Bland-Altman plots reveal non-negligible limits of agreement (upper and lower dotted lines) for cycle-averaged MCA flow rates computed via the different outflow strategies. Colored circles represent the 10 cases subsequently selected for 3D CFD. Note that the modified extents (cases A′, B′, C′) show non-negligible changes in flow rates for both zero-pressure and Murray-law methods. The y-axis is defined as the considered outflow method minus the splitting method flow rate, divided by the average of flow rates obtained with the 2 methods.
Impact of Model Extent on Flow Rate Distribution
Figure 2 also shows that the modified extents for cases A–C resulted in non-negligible changes in MCA flow rates. This impact of modified extents is further highlighted in Table 2, where removal of the M3 branches for case A resulted in 40% and 138% increases in MCA flow for the zero-pressure and Murray-law methods, respectively. Even the mild modifications to cases B and C resulted in non-negligible changes in MCA flow rates: a 22% decrease with the Murray-law method for case B versus B′ and a 73% increase with the zero-pressure method for case C versus C′. With the splitting method, model extent had no impact on flow rates as designed.
Impact of model extent on predicted MCA flow rates (mL/min) for the 3 cases with modified outflow extents
Such impact of model extent can also be inferred from cases that had an ACA aneurysm in addition to the MCA aneurysm (squares in Fig 2). These models extended at least to the A2 branches, whereas those with only MCA aneurysms were truncated at the A1 branch. Per Fig 2A, such imbalance in the number and size of MCA-versus-ACA outlets—more flow resistance is offered by the M1+M2 branches versus the single A1 stub—led to the marked underestimation of MCA flow rates by the zero-pressure method, whereas the ACA aneurysm cases (having more balanced ACA and MCA flow resistances) had flow rates closer to the splitting method. Per Fig 2B, the Murray-law method was less susceptible to such frank imbalances.
Impact on Hemodynamic Indices
Quantitative differences between zero-pressure and splitting methods for nominal metrics of rupture status are shown in Fig 3, to appreciate the practical impact of outflow strategy. (Qualitative differences in the surface distributions from which these sac-integrated values were derived can be seen in On-line Figs 2 and 3.) Normalized TAWSS was affected by the outflow method, showing absolute differences of 21% ± 18%, reflecting differences in the CFD-computed wall shear stress patterns. Without the parent artery normalization, differences in aneurysm TAWSS were 40% ± 25%, reflecting the compounding effect of differences in parent artery (MCA) inflow rates on wall shear stress magnitudes.
Comparison of 3 normalized hemodynamic indices (see “Materials and Methods”) determined from 3D CFD simulations, computed using zero-pressure-versus-splitting methods. The lower right of each panel shows the coefficients of determination (R2) as measures of correlation and agreement. The dotted line is perfect agreement. Note the triangles and different colors used to highlight the different truncations for cases A, B, and C (see symbol legend in Fig 2). These data points are translated only horizontally because the splitting method is, by design, not sensitive to model extent.
For the zero-pressure method, changes in the extent of models A, B, and C resulted in TAWSS differences of −12%, +22%, and −22%, respectively. It is important to note that these differences are observed for the same cases with the same outflow strategy, with the only difference being an alternative but equally plausible choice of 3D model extent. Sac-maximum TAWSS showed absolute differences of 17% ± 13% between the outflow methods and the comparable effect of model extent. Larger differences were observed for the sac-averaged OSI, showing more pronounced changes in the rank ordering of cases from low-to-high OSI. The absolute difference for OSI was 24% ± 17%, and modified extent had a greater impact than for sac-average and sac-maximum TAWSS, with differences from 24% to 63% relative to the OSI from the splitting method.
Discussion
Relationship to Previous Studies: Impact of Outflow or Not?
Previous studies reporting the impact of the outflow rates on aneurysm CFD have shown mixed conclusions but typically relied on single cases and hemodynamic indices not commonly used in large-scale studies. Ramalho et al16 concluded that the choice of outflow strategy “highly influence[d] the hemodynamics inside the aneurysm [and] should be chosen with special caution.” Venugopal et al17 varied the outflow distribution by factors 2–12 times greater than what we reported, yet seeing its impact on TAWSS maps, they concluded that it was “difficult to draw any conclusions regarding the role of shear stress in the aneurysm's growth.” Grinberg and Karniadakis18 showed that using zero-pressure versus lumped parameter outflow conditions led to a “300% difference in flow rates through the outlets” for their patient-specific cerebrovasculature model. On the other hand, Cebral et al19 reported that outflow rate variations of ±25% “did not alter the main characteristics of the intra-aneurysmal flow patterns.” These disparate conclusions have given mixed messages to the CFD and clinical communities and, as noted by Ramalho et al,16 called for “stud[ies] with larger number of geometries.”
The present study, based on 70 cases (10 with CFD), demonstrates that the choice of outflow strategy can alter the physiological plausibility and results of aneurysm CFD, with non-negligible impact on both flow rate distributions and nominal hemodynamic predictors of rupture status. The 2 most popular outflow strategies—zero-pressure and Murray-law—were also sensitive to subjective choices about CFD model extent, which are rarely documented and can vary widely among the CFD community. For example, in the most recent aneurysm CFD Challenge, which focused on MCA aneurysm cases, the presence and extent of proximal ICA, side branches, or distal bifurcations varied dramatically.20 The splitting method we propose is minimally sensitive to such uncertainties and demonstrably offers robust and realistic flow distributions.
The differences we report in hemodynamic indices due to outflow strategy, as shown in Fig 3, are comparable with those reported due to CFD solution strategy,2 which were, in turn, far greater than those due to non-Newtonian viscosity.21 Similarly, the divergent conclusions of Ramalho et al16 and Venugopal et al17 were drawn from TAWSS patterns that appear, qualitatively at least, to be comparable with those presented in studies of geometric22⇓–24 or inflow3 uncertainties. Together these suggest that the outflow strategy should be considered at least as important as other assumptions or uncertainties in aneurysm CFD.
Implications for Clinical CFD: Does Outflow Really Matter?
It could be argued that our findings are exaggerated by our use of extensive models that include the cervical ICA with side (OA, PcomA) and terminal ACA branches, meaning that the aneurysm (ie, MCA) inflow is also effectively determined by the outflow method. Indeed, CFD studies of MCA aneurysms often comprise just the single (M1/M2) bifurcation.2,25,26 Setting aside the thorny question of whether flow in the MCA can be considered fully developed in such studies considering complex flow patterns almost certainly present upstream of the ICA siphon and terminus,24,27 for these cases, the Murray-law and splitting methods would, assuming the same power law exponent, indeed give the same results. Consider, however, our case J, for which all 3 outflow methods predicted the M1 inflow to within <5%. The difference in M2 outflows predicted by the zero-pressure and splitting methods was a modest 15%, yet this resulted in a 50% difference in the normalized sac-averaged TAWSS (0.17 versus 0.36, respectively). This would seem to suggest that our findings about the zero-pressure strategy would hold even for less extensive aneurysm CFD models.
To appreciate the potential implications of this, consider that Takao et al25 and Miura et al26 used a zero-pressure strategy for their large (50 and 106 cases, respectively) studies of (truncated) MCA aneurysm CFD models and reported differences in sac-averaged TAWSS of −13% and −25%, respectively, between ruptured and unruptured cases. Note that these differences are comparable with differences we observed for the zero-pressure method just by changing the 3D model extent. Such differences would only be exacerbated for smaller studies, and emphasize the need for consistency, both between and within studies, regarding the extent of the 3D models.
Recommendations for Clinical CFD: Stop Doing Nothing
The zero-pressure outflow strategy remains popular because it places fewer constraints on the outlet velocity profile. A traction-free boundary condition is used by default, which requires only that the velocity vectors are approximately parallel to the vessel axis. On the other hand, prescribing outflow divisions as in the Murray-law or splitting methods often requires the velocity profile itself to be prescribed, which, in turn, usually necessitates the addition of long flow extensions to avoid errors or instabilities engendered by this more explicit approach.
Some solvers, like ours, offer outlet flux boundary conditions, which allow the user to specify outflow rates without having to prescribe a velocity profile. In our case, outlet flux is implemented by iterating the outlet pressures (which appear in the traction-free outlet equations) until the desired flow split is achieved.28 For solvers without such capability, we note that our splitting method, when coupled with a 0D solver as described in the “Materials and Methods,” can provide an accurate estimate of the outflow pressures required to achieve a given outflow division.6 That therefore makes it possible for any CFD solver to achieve a desired outflow division using more relaxed traction-free (versus restrictive velocity profile) boundary conditions, simply by prescribing those estimated outlet pressures.
Notes of Caution
In this study, we compared the splitting method using n = 2 with the Murray-law method using its customary n = 3. The impact of the choice of the power law exponent on outflow strategy is summarized in On-line Table 1. Average flow rates from the Murray-law method were relatively insensitive to n; however, intraindividual variations were larger with n = 3, as expected.4 The splitting method was more sensitive to the power law exponent, with n = 3 underestimating flow to smaller branches (OA, PcomA) and biasing flow to the larger MCA. This dependence on n, however, underscores the accuracy of the flexibility of the splitting method: Flow division at the crucial ICA terminus is demonstrably closer to a square law29,30; n can be increased toward a cube law for more distal branches, if desired.
The way diameters are measured from the 3D models could also introduce some variability in the calculation of outflow rates. Owing to our use of a 0D model as the basis for the splitting method, our diameters were averaged over each vessel segment. For the Murray-law method, we used the mean diameter of each outlet face, as is customary, which can create a dependency on where the artery segment is truncated.27 For cases B and B′, for example, the number of outlets was the same, just the branch length was reduced, yet there was a 22% change in the MCA flow rates. This difference would undoubtedly have been less had we based the Murray-law outflows on vessel-average diameter rather than the (customary) single-point outflow face. On the other hand, our use of a 0D model makes it possible to account, a priori, for the hemodynamic impact of vascular abnormalities, such as a stenosis, on the outflow rates, in a way that would be more difficult using standard zero-pressure or Murray-law methods.
Conclusions
Outflow strategy is at least as important as other assumptions made for cerebrovascular CFD. The prevalent zero-pressure or do-nothing strategy should be avoided because it is both unphysiological and operator-dependent. The Murray-law strategy is better than nothing, but here we present our more physiological and robust splitting method as an open source tool (www.github.com/ChrisChnafa/aneuTools) to encourage standardization of cerebrovascular CFD. Irrespective of their outflow strategy, CFD studies should always report, or at least demonstrate, the physiological plausibility of the actual outflow (and inflow) rates experienced by their models.
Footnotes
This work was supported by grants from the Heart & Stroke Foundation (G-13-0001830, to D.A.S.) and the University Health Network (AMO AFP 410004349, to V.M.P.). D.A.S. also acknowledges the support of a Mid-Career Investigator Award (MC7455) from the Heart & Stroke Foundation of Canada.
References
- Received July 24, 2017.
- Accepted after revision October 13, 2017.
- © 2018 by American Journal of Neuroradiology
