Introduction

The worldwide prevalence of intracranial aneurysm disease is estimated to be 3.2% [95% (CI), 1.9–5.2%] [1]. Intracranial aneurysms are characterized by high morbidity and mortality rates following rupture [2, 3]. It is generally considered that hemodynamics plays a critical role in aneurysm development and rupture [4,5,6]. While computational fluid dynamics (CFD) studies have the potential to provide important insights into the mechanisms of aneurysm development and rupture, there remains concern as to the poor consistency of CFD results [7]. Several assumptions are commonly made when conducting CFD simulations of flow in intracranial aneurysms, specifically that vascular compliance is negligible, that flow is laminar, and that blood viscosity is Newtonian. These assumptions have been shown to be reasonable [8, 9], and with those assumptions, the accuracy of CFD rests only on an accurate definition of the geometry of the lumenal surface, and of the inlet flow waveform.

The error in geometry segmentation and its effect on hemodynamics has been explored elsewhere [10, 11]. However, the dependence of hemodynamic descriptors on the inlet flow waveform boundary condition is less well established. Many investigators use generalized or idealized flow waveforms [12], and when patient-specific flow rates are used, little attention is paid to measurement variability. Blood flow in the supplying vessels of the brain can be measured using phase-contrast MR (PCMR). PCMR enables the acquisition of geometry and flow information simultaneously [13]. This technique is also noninvasive, and does not require exposure to radiation or contrast material [14].

There are little data on serial studies of intracranial flow in the same patients as patient-specific flow conditions are rarely assessed as part of the routine clinical examination. de Verdier et al. performed a test–retest measurement on a relatively small sample size [15]. However, that study only included healthy volunteers and therefore did not account for the variability that might be expected in patients with vascular disease that include: variations from making measurements at different locations along the feeding artery; and increased physiological variability. These physiological factors are confounded by reader errors in post-processing and analysis with delineation of the lumenal contour being a factor of major concern.

As noted above, it is important to understand the reproducibility of measurement of the inlet flow waveform in assessing the value of CFD in predicting vascular disease progression. This study was limited to the evaluation of patients with known but untreated intracranial aneurysms. To measure the variation of PCMR that is generally encountered in patients with vascular disease, we carried out a longitudinal study in unruptured intracranial aneurysm patients with multiple follow-up scans. We further investigated the effect that likely uncertainties in flow measurement would have on computed hemodynamic metrics.

Materials and methods

Patient selection

A cohort of patients with intracranial aneurysms were recruited and scanned from April 2005 to January 2013. Inclusion criteria included: diagnosis of at least one unruptured intracranial aneurysm, and the ability to undergo an MR study. Exclusion criteria included: metal implants, claustrophobia, or a known allergic reaction to the MR contrast agent. This study was conducted under IRB approval. All subjects involved gave informed written consent for study participation. A total of 51 patients (37 females and 14 males) with studies performed at three time points were included in this study. For patients with multiple aneurysms, we only analyzed flow into the largest one.

MR protocol

A 1.5 T MR (Achieva, Philips Medical System, Best, The Netherlands) with a six-channel head coil was used in this study. First, a three-dimensional contrast-enhanced MRA (CE-MRA) sequence was acquired following institutional clinical practice. A weight-adjusted single dose of GdDTPA diluted with saline to a 20 mL volume was injected and was followed by 15 mL of saline all delivered at 2 mL/s. The CE-MRA sequence used elliptic-centric phase reordering, and data were acquired using parallel imaging with an acceleration factor of 2. Imaging parameters included the following: TR/TE/flip angle × 5/2/30°. Images were acquired from a 54-mm paracoronal slab, with an FOV of 240 mm and an acquisition matrix of 400 × 380 × 45 zero-filled to 512 × 512 × 90. The resultant images had a resolution of 0.6 × 0.63 × 1.2 mm3 and was interpolated to 0.47 × 0.47 × 0.6 mm3. Total acquisition time was of the order of 35 s.

MIP images reconstructed from the CE-MRA data were used to localize the aneurysm and its feeding arteries. The 2D through-plane velocity-encoded sequence was then oriented in a plane transverse to the feeding artery using two orthogonal views (Fig. 1a) and velocity encoding was prescribed for through-plane flow with the encoding gradient perpendicular to the image plane. To facilitate measurement, the 2D through-plane was placed at least five times of the inlet diameter from the proximal of the aneurysm at a relative straight segment without stenosis. Data were acquired using a finger pulse monitor with retrospective gating. Acquisition parameters were as follows: TR/TE/flip angle = 8/5/15°; FOV 160 × 160 mm with a slice thickness of 5 mm. Velocity encoding (VENC) was set to 100 cm/s. The acquisition matrix was 160 × 160, and the images were reconstructed to a 0.6 × 0.6 mm resolution. The data were retrospectively triggered, and reconstructed to either 15 or 32 time frames, providing the time-varying inlet boundary conditions required for the pulsatile flow simulations.

Fig. 1
figure 1

Placement of 2D PC plane, ROI contouring and flow information extraction for a patient with an aneurysm of the left ICA. a MIP image of CE-MRA; yellow line on anterior–posterior and lateral view shows the position of perpendicular plane at C2 segment; ROI on b magnitude and c phase image after automatic refinement; d net-flow plot through entire cardiac cycle; e velocity values from 32 phases of one cardiac cycle which is consistent with a generally parabolic velocity distribution

Full size image

Post-processing

A research Graphical User Interface (GUI) software Segment (2.2R7056) [16] was used to draw the regions of interest and extract velocity information. First, an elliptic region of interest was roughly drawn around the inlet artery using the magnitude image from the PCMR data; second, an automatic tracking and refinement function was used to refine the regions of interest across all phases (Fig. 1b, c). The final ROIs were independent of the initial contour and further reviewed by two neuroradiologists; discrepancy was resolved by discussion. Third, flow values were measured at each time point through the entire cardiac cycle (Fig. 1d, e). The time-averaged flow rate through the cardiac cycle [mean flow (MF)], flow at peak systole [max flow (MAF)], flow at end diastole [min flow (MIF)], peak systolic velocity [max velocity (MAV)], end-diastolic velocity [min velocity (MIV)], velocity averaged over the whole cycle [mean velocity (MV)], and the area of the region of interest (ROI) were extracted for further analysis.

CFD simulation

The 3D geometry of the intracranial aneurysms was segmented and reconstructed from CE-MRA using Mimics 20.0 (Materialise NV, Leuven, Belgium). The surface was remeshed to an edge size of 0.4 mm and then smoothed while preserving the shape using Geomagic DesignX (Geomagic, Research Triangle Park, NC, USA) [17]. An inlet vessel length of at least four diameters of the parent artery was incorporated to ensure fully developed flow [18]. The models were meshed in Ansys with a mixture of prism and tetrahedral elements. As previously reported [19, 20], grid size verification ranging from 0.4 to 0.1 mm was performed on two randomly selected patients using steady flow. The difference in average wall shear stress (WSS) was calculated and the maximum grid size was selected to provide a difference in WSS smaller than 5%—which was for a grid size of 0.2 mm.

CFD transient simulations were performed three times in ten aneurysms, once with the inlet waveform obtained from PCMR and the other with an inlet waveform either artificially elevated or reduced (where an increase or decrease in flow was enforced at every time point through the cardiac cycle) as the inlet boundary condition. The outlet pressure was set to 0 Pa. Blood was simulated as a Newtonian fluid with a density of 1060 kg/m3 and dynamic viscosity of 0.0035 Pa.s, respectively. The simulation was performed using Fluent’s pressure-based solver with the DOUBLE pressure–velocity coupling scheme. The CFD convergence criteria were set at 0.001 for the continuity and velocity residuals. Three cardiac cycles were performed for each patient, and the data of the last cardiac cycle were saved for visualization and statistical analysis.

CFD post-processing was performed with paraview 5.6 (Kitware, Clifton Park, NY) and python 3.7 (https://www.python.org/). Systolic WSS, diastolic WSS, time-averaged wall shear stress (TAWSS), and low wall shear stress area (LSA) were each calculated. LSA is defined as the percentage of the aneurysm area where the WSS is 1 SD below the mean WSS in the parent artery (within 1 cm of centerline node at the aneurysm neck) [21]. A point-by-point comparison of WSS and LSA was made between the different simulated flow conditions.

Statistical analysis

Statistical analysis was performed using Python 3.7. All measurement data are presented as mean ± standard deviation; enumeration data are presented in percentages. The intraclass correlation coefficient (ICC) and coefficient of variance (CV) of MF, MAF, MIF, MV, MAV, MIV, and ROI area of each artery was calculated. The mean value of 51 patients was used to characterize the variation of the included parameters. Bland–Altman plots were performed to analyze the agreement among three time points. A T test or Wilcoxon test was performed to compare the difference between the two groups. A p value of 0.05 was selected as the threshold for statistical significance.

Results

The average age of the 51 included patients was 62 ± 16 years. 29 aneurysms were located on an internal carotid artery (ICA), 11 were located on the basilar artery (BA), 6 on the vertebral artery (VA), 3 at the anterior communicating artery (ACA), and 2 on the middle cerebral artery (MCA). The average time interval between sequential studies was 0.64 ± 0.21 years.

Normal ranges

The mean value and the standard deviation of MF, MAF, MIF, MV, MAV, MIV, and ROI area are presented in Table 1. The distributions of investigated parameters across different artery segments are displayed in online-Fig. 1. The largest ROI area was in the ICA and the smallest in the ACA.

Table 1 The range of flow, velocity parameters, and ROI at different locations
Full size table

Intra-individual variation

The overall reproducibility for all the parameters included is shown in Table 2. The limits of agreement of three measurements are plotted in Fig. 2. The ICCs are excellent for MF, MAF, MIF, MV, MAV, and ROI. The CV is larger in MIF than in MF and MAF (0.16 vs. 0.10 and 0.11, p = 0.007 and p = 0.01). Similarly, the CV is larger in MIV than MV and MAV (0.18 vs. 0.12 and 0.12, p = 0.02 and p = 0.01). The Histogram of frequency distribution of flow and velocity parameters is shown in online-Fig. 2. The CV distribution of flow and velocity parameters based on artery locations is shown in online-Fig. 3. The CV of ROI area is 11%. The median area size of 0.15 cm2 was selected as a cut-off value for dividing the aneurysms into two groups. A T test showed that there is no significant difference in the CV of MF between the group with a larger ROI area than the group with a smaller ROI area (p = 0.88).

Table 2 Reproducibility and variation of included parameters
Full size table
Fig. 2
figure 2

Bland–Altman plots for flow and velocity parameters for paired comparisons from 3 time points. Horizontally, the average of the two measurements is plotted, and vertically, the difference of these two measurements. Middle dash line: mean difference; upper and lower dash lines: limit of agreement (mean ± 1.96 standard deviation, the range between the lines represents least detectable difference). Blue: measurement agreement between first and second time point; Green: measurement agreement between first and third time point; Red: measurement agreement between second and third time point

Full size image

Influence of flow variance on CFD results

CFD simulations were performed three times in the ten aneurysms from five locations, once with the measured flow waveform and then, to assess the effect of typical measurement inaccuracy, with artificially elevated or reduced flow waveform (of + 10% or − 10% based on these experimentally measured CV of MF) as the inlet boundary condition, to explore their impact on CFD derived flow metrics. The aneurysm information and simulation results are summarized in Table 3 and demonstrated in Fig. 3. A 10% error in flow resulted in a number of hemodynamic changes including: 41.41% of WSS at systolic phase, 39.13% of WSS at diastolic phase, 2.79% of LSA at systolic phase, 2.12% of LSA at diastolic phase, 47.57% of TAWSS, and 0.17% of oscillatory shear index (OSI).

Table 3 The change in computed WSS, LSA, TAWSS, and OSI assuming a ± 10% flow error
Full size table
Fig. 3
figure 3

Example of computed WSS and LSA distribution on diastolic and systolic phases assuming a 10% variation in flow

Full size image

Discussion

In this study, we determined the variability over time of flow and velocity parameters derived from 2D PCMR for patients with unruptured intracranial aneurysms where these patients were examined at 3 time points in follow-up. We found a wide range and generally good reproducibility for the flow and velocity parameters. We then examined the impact of this variability in computed hemodynamic descriptors when the original flow waveforms and artificially elevated flow waveforms (original flow plus variation) were used as inlet boundary conditions in CFD calculations. The results indicate that 10% variation in intra-individual flow could lead to a difference of around 40% in WSS but only around 2.5% in LSA. The result suggests that we should be cautious when we use absolute value of WSS in the evaluation of aneurysm growth or rupture risks, on a patient-specific basis.

Normal flow and velocity distribution have been reported previously [15, 22]. The reported mean flow and velocity were generally of the order of 10–20% greater than the mean flow in our study. A number of factors might contribute to these discrepancies such as magnetic field strength, methods for defining the ROI area [23], and patient age [24]. Given that the population in our study is older (62 ± 16 years) than in these studies, it is likely that age is the major contributor to this discrepancy.

The reproducibility and variation of 2D flow in the intracranial arteries has been reported in the previous studies [22, 25,26,27]. Taviani et al. reported that the CV of MF in the aorta is about 0.107 in nine volunteers [28]. Aart et al. reported that the CV of repeated total CBF measurement is about 11% in 15 volunteers [27]. Our study found a similar result for CV based on a large dataset over a long-term follow-up (mean 2.1 years). The overall reproducibility of MF and velocity parameters was good, even though measurements were repeated with long (> 6 months) intervals between studies. The improved determination of reproducibility in our study could be attributed to three reasons: first, a strict methodology was employed to specify the plane perpendicular to the inlet vessel; second, the ROI was drawn with an automated refinement tool; third, we have a larger sample size and number of follow-up time points. We also compared the CV of flow between arteries with large area size and arteries with small area size. Results showed that CV is independent of artery area size. All these results indicate that PCMR measurements are reliable, and that the flow conditions in these aneurysm patients are fairly stable.

CFD has played an important role in intracranial aneurysm research. Specification of inlet boundary conditions has evolved from generalized and idealized flow values to patient-specific flow conditions [12, 29, 30]. A comparison of simulation results showed that both aneurysm WSS magnitude and WSS distribution were different depending on whether generalized flow or patient-specific flow was specified [12]. In our longitudinal study with a relatively large dataset, we showed that the CV of MF is about 10%. We further verified the effect of flow error on computed WSS. A 10% variation in flow was found to generate of the order of a 20% change in WSS, but only a 4–5% change in LSA and OSI. The determination of the distribution of WSS is likely of high importance in aneurysm evolution. Specifically, we have previously reported aneurysm growth is co-located with areas of LSA [5]. These findings suggest that we should be cautious about the absolute value of WSS, when we interpret the results of a patient-specific CFD calculation. However, the LSA and OSI might be reliable tools for evaluating aneurysm growth or rupture risks when a patient-specific inlet boundary-flow waveform condition is used.

There are additional limitations to this study. First, since this is an intracranial aneurysm patient cohort, more data were acquired for locations in the ICA and BA than other territories, since the measurements were performed to determine inlet flow velocities in the parent vessels of intracranial aneurysms, and these locations were the sites of highest prevalence of aneurysms in this study. Second, the time interval between each repeated study is longer than in typical reproducibility studies. Although this implies that measurements in a narrower time interval might yield more accurate estimates of reproducibility, it also would suggest that measurement reproducibility might be even better than what we report. Third, our study used 2D PCMR which only obtained the through-plane flow information. In recent years, the use of 4D flow imaging which enables full coverage of the intracranial vasculature has gained increased interest [19, 31]. However, those studies require much longer scan times (~ 12 min) than a single 2D PCMR acquisition (~ 30 secs) and are therefore subject to patient motion, and are difficult to accommodate in routine clinical practice. Fourth, our study used a 1.5 T MR system. While the influence of magnetic field strength on flow measurement remains unclear, a higher magnetic field strength of 3 T or 7 T would offer higher SNR, an improved velocity-to-noise ratio, and likely a better ability to delineate the ROI. Fifth, we directly used a post-contrast PC-MRI, and the influence of residual Gd on PC-MRI was not evaluated here. Sixth, in this study, we selected saccular aneurysms for CFD analysis that were representative of the larger set of aneurysms that we have analyzed. However, aneurysm geometry type may lead to different sensitivity to flow variation, and the possible geometric configurations are extremely large. Future studies with larger sample size should be performed to focus on this variability. Seventh, it would be preferable in future simulations, to include a greater length of patient-specific vascular geometry proximal to the inlet location to generate inlet boundary conditions that are even closer to the real physiological conditions.

In conclusion, this study identified that intra-individual MF, as measured by PCMR, varied by an amount of the order of 10%. A 10% variation in flow was found to generate of the order of a 40% change in WSS, but only a 2.5% change in LSA and 0.17% change in OSI. These findings suggest that we should be cautious about the absolute value of WSS, when interpreting the results of a patient-specific CFD calculation. However, the LSA and OSI might be reliable metrics for evaluating aneurysm growth or rupture risks when a patient-specific inlet boundary-flow waveform condition is used.