Model assumption about circle of Willis

For the CFD simulations, we assumed it sufficient to segment and analyze the anterior part of the CW in order to obtain reliable estimates on the hemispheric MCA pressure laterality. To motivate and exemplify this, we used a circuit analogy (Fig 2). The CW can be represented as a lumped model of an electrical circuit, where the change in pressure across each artery is related to the flow (Q) and the vascular resistance (R) of the artery (Fig 2). If we want to determine the pressure laterality between the ipsi- and contralateral MCA, we can follow the circuit across the anterior pathway, via the proximal anterior cerebral arteries (ACA1) and the anterior communicating artery (ACoA) (i.e., MCA→ACA1→ACoA→ACA1→MCA). Alternatively, across the posterior pathway, via the posterior communicating arteries (PCoA) and the proximal posterior cerebral artery (PCA1) (i.e., MCA→PCoA→PCA1→PCA1→PCoA→MCA). The pressure change across each individual artery is different, but the total sum of pressure changes between the ipsi- and contralateral MCA must be the same for both pathways. Thus, although both the anterior and the posterior circulation contributes to the overall resistance between the left and right hemisphere, analysis of only one pathway is required to determine the pressure laterality. Correspondingly, for estimating the MCA pressure laterality with CFD, the use of only the anterior or posterior pathway should be sufficient. The anterior pathway offers larger vessels and higher flow rates compared to the posterior pathway, which makes the anterior pathway the most feasible choice. In summary, CFD simulations using the anterior pathway for estimating an MCA pressure laterality should be valid as long as the correct flow rates and geometries of the MCA, ACA1 and ACoA arteries are used.

Fig 2

Circuit analogy of the circle of Willis.

The figure illustrates a schematic of the CW, denoting the ipsi (I)- and contralateral (C) side of the ICA stenosis (St). It is modelled as an electrical circuit, where the change in pressure across each artery is related to the flow (Q) and the vascular resistance (R) of the artery. The arrows indicate the normal flow direction in each artery, and the circles indicate the two pressure sites for the ipsi- and contralateral MCA pressure (P_MCA). The two equations for P_MCA,I describe how the pressure changes, starting at P_MCA,C, across the anterior collateral pathway versus the posterior collateral pathway. They show that the pressure laterality between the ipsi- and contralateral MCA must be the same for both pathways. Note that the MCA resistance only refers to the arterial segment from the bifurcation to the location of the pressure measurement.

The correct flow in MCA, ACA1 and the ACoA was accomplished by the boundary conditions of the CFD model (See section about boundary conditions). These conditions required segmentation of the left and right ICA, MCA, ACA1, ACoA and the distal anterior cerebral artery (ACA2). Based on the model assumption, the posterior circulation and the ophthalmic arteries were removed.

Patients

This study was based on a patient cohort previously described, where 38 TIA or ischemic stroke patients with a symptomatic ICA stenosis were included and investigated with 4D flow MRI [3]. For the current study, 10 of these patients were excluded because of no CTA investigation (n = 7), insufficient image quality (n = 2) and a missing ACA1 (n = 1). The final study population thus consisted of 28 patients (73±6 years, range 59–80 years, 21 men) with a symptomatic ICA stenosis ≥50% or occlusion (n = 1) with or without a non-symptomatic contralateral stenosis. The mean degree of stenosis on the ipsi- and contralateral side of the symptomatic stenosis was 76±12% and 38±26%, respectively, graded according to the North American Symptomatic Carotid Endarterectomy Trial (NASCET) [17]. The mean modified Ranking Scale (mRS) score was 0 (range 0–1) and the mean NIHSS score was 1 (range 0–6). The systolic and diastolic blood pressure was 138±17 mmHg and 70±10 mmHg, respectively, with a mean arterial pressure (MAP) of 93±11 mmHg. The heart rate during the MRI was 64±11 bpm. Risk factors were diabetes mellitus (25%), hyperlipidemia (54%), hypertension (75%) and ever smoker (64%).

The ethical review board of Umeå University (Dnr: 2011-440-31M) and the Swedish Ethical Review Authority (Dnr: 2019–05909) approved the study. It was performed in accordance with the guidelines of the Declaration of Helsinki. Oral and written information about the study was given to all participants and written consent was obtained from all participants.

Imaging

For flow rate measurements, MRI scans were performed on a 3T scanner (GE Discovery MR 750, Milwaukee, WI, USA) with a 32-channel head coil. With a scan time of about nine minutes, the 4D flow MRI data was acquired with a balanced 5-point phase contrast vastly undersampled isotropic projection reconstruction (PC-VIPR) sequence with full brain coverage [15, 18]. Scan parameters were: TR/TE 6.5/2.7 ms, Venc 110 cm/s, flip angle 8°, 16000 radial projections, acquisition resolution 300×300×300, imaging volume 220×220×220 mm³, reconstructed resolution 320×320×320, isotropic voxel size 0.7×0.7×0.7 mm³, 20 reconstructed cardiac time-frames. An angiographic complex difference image and velocity images in three directions were reconstructed.

Flow rate quantification was performed according to a previously developed post-processing method [19]. Pulsatile blood flow rates were measured in the left and right ICA, proximal middle cerebral artery (MCA1), ACA1 and the ACA2. There was no measurable flow in three ACA1 and four ICA on the ipsilateral side. Flow rates were obtained by averaging the flow waveforms from 15 consecutive cut planes [19].

Patient-specific vessel geometries were obtained from the clinical CTA images, from different hospitals. The trans-axial image resolution ranged between 0.33 to 0.625 mm, and the reconstructed slice thickness ranged between 0.3 to 0.625 mm.

Geometries and mesh generation

Before segmentation, CTA images were cropped to only include the CW. In axial direction, ICAs were cut at the C3-C4 segment and the distal anterior cerebral arteries at ACA3 level. In left/right direction, the MCAs were cut to include MCA2. We included MCA2, ACA3 and the relatively long ICAs instead of adding flow extensions to the geometry. The segmentations were performed with Synopsys’ Simpleware™ software (ScanIP P-2019.09; Synopsys, Inc., Mountain View, USA).

Images were resampled with linear interpolation to obtain an isotropic voxel size of 0.3 or 0.3125 mm, depending on original resolution. Background noise from surrounding tissue was reduced by a bilateral filter. For vessel wall detection, a mask was generated with a threshold filter, followed by the software specific gradient-based filter ‘Local surface correction’. Before model generation, a volume and topology preserving smoothing filter was applied.

All meshes were generated with the Simpleware FE module and imported as meshes into COMSOL Multiphysics® (COMSOL Multiphysics®, version 5.4, www.comsol.com, COMSOL AB, Stockholm, Sweden). In the refinement region of our simulations, i.e., across the MCA1/ACA1 bifurcations (Fig 3), our target edge length was equal to the interpolated voxel size. Outside the refined region, the target edge length was increased to 1.5 times the voxel size. The meshes included boundary layers with four sublayers and a total target width of approximately 0.14 mm [20]. The total number of elements ranged between 0.9–1.6 million, and the average edge length of all elements was 0.17 mm.

Fig 3

A computational geometry and mesh from a representative patient.

From the posterior view, with left and right denoted in the figure. Note how only the anterior part of the CW was included in the model, according to our model assumption. The black square shows the mesh refinement region of our simulations. Longer segments of the arteries were included as flow extensions.

A coarser and a finer mesh, without any refinement regions, was used to investigate the mesh sensitivity of the stationary solutions. For the coarse mesh, the target edge length was twice the voxel size and for the fine mesh it was equal to the voxel size. The target boundary layer thickness was 0.2 and 0.1 mm, respectively. The total number of volume elements ranged between 0.5–1.1 and 2.7–5.5 million elements for the coarser and the finer mesh, respectively.

Computational fluid dynamics

Navier-Stokes equations were solved with the CFD module in COMSOL Multiphysics®. Blood was assumed to be an incompressible Newtonian fluid with density ρ = 1060 kg/m³ and viscosity μ = 3.45 mPa·s. The flow was assumed laminar and the vessel walls impermeable and rigid. Primary simulations with stationary flows were performed for all of the 28 patients.

Time-resolved simulations with pulsatile flow waveforms were performed in 24 of 28 patients. One patient was excluded because of convergence problems and three because of MR reconstruction problems. Time-dependent simulations were run for two cardiac cycles, with a linear interpolation between each of the 20 cardiac phases, and 4000 time steps per second. The initial flow condition was ramped up from 0 to the flow rate in the first cardiac phase during about 0.5 s, followed by the full cardiac phases. The second cardiac cycle was used in the analysis. A pilot simulation with four cycles indicated two cycles to be sufficient.

Boundary conditions

There were totally six inlet and outlet boundaries in the CFD model. Inlet boundaries were the left and right ICA. Outlet boundaries were the left and right MCA and ACA2. At each inlet or outlet boundary, we applied a laminar parabolic flow rate profile, using COMSOL’s build-in boundary condition for fully developed flow. A no-slip condition was applied to the vessel walls. For mass conservation within the model, and to achieve the measured flow in MCA and ACA1, the flow rate boundary conditions were defined as follows:

The inlet boundary in each ICA was calculated from the summed measured flow rates of MCA1 and ACA1 on the same side.

The outlet boundary in each MCA was set as the measured flow rate.

The right ACA2 outlet was calculated by (ACA2_R / ACA2_tot)*ACA1_tot.

The left ACA2 outlet was calculated by (ACA2_L / ACA2_tot)*ACA1_tot.

The ACA1_tot and ACA2_tot was the sum of the measured flow in the left and right ACA1 and ACA2, respectively. Thus, the proportional difference between the measured flow rates in ACA2 was used as boundary condition.

Since the pressure is only uniquely defined up to a constant when flow rates are applied at all boundaries, i.e., Dirichlet condition, the pressure level was set as a zero-reference pressure point condition at the contralateral ICA inlet boundary.

Pressure and flow laterality

The hemispheric MCA pressure laterality (ΔP_MCA) was defined as the pressure difference between the contra- and ipsilateral MCA1, where ipsilateral is the side of the symptomatic stenosis. The goal was to measure the average pressure at a distance of 5 mm from the center of the MCA1/ACA1 bifurcation. For the time-resolved simulations, the pulsatile pressure waveforms of the contra- and ipsilateral MCA1 were averaged and subtracted to give the ΔP_MCA.

The hemispheric flow rate laterality (ΔQ) was defined as the flow rate difference between the contra-and ipsilateral ICA, MCA1 and ACA1, respectively, measured in the 4D flow MRI data.

Based on findings in a previous modelling study and measurements in canines, it has been shown that a degree of stenosis >50% can generate ICA pressure drops higher than 5 mmHg [21, 22]. This motivated us to set a 5 mmHg MCA pressure laterality as the lower limit of a potentially clinical interesting difference.

Statistical analysis

All mean values were presented as the mean±standard deviation. A one-sample Kolmogorov-Smirnov test was used for normality testing. Wilcoxon signed rank test was used to test for all significant differences. Pearson’s linear correlation coefficient was used for all reported correlations. The significance level was set as p<0.05. All statistical analysis and MRI flow rate measurements were performed with MATLAB (R2019a, MathWorks, Natick, MA).

MCA pressure laterality

In the 28 patients, the hemispheric MCA pressure laterality (ΔP_MCA) between the contra-and ipsilateral side was 6.4±8.3 mmHg (p<0.001) (median 3.9 mmHg, range -1.3 to 31.9 mmHg) (Fig 4A). Degree of ipsilateral ICA stenosis was not correlated to ΔP_MCA (r = 0.23, p = 0.23). When dividing patients into those with or without a contralateral ICA stenosis ≥50%, the ΔP_MCA was 5.5±5.1 mmHg and 6.9±9.6 mmHg, respectively. In Fig 4B, the pressure from the patient in Fig 3 is shown. The pressure result for each patient is found in S1 Fig.

Fig 4

MCA pressure laterality.

(a) Histogram of ΔP_MCA for the patient group. (b) The CFD pressure of a stationary simulation for the case in Fig 3. A zero-reference pressure was set at the contralateral ICA, which was the right ICA inlet in this case. The ΔP_MCA was 3.0 mmHg in this example.

The mean difference in ΔP_MCA between the stationary and the time-resolved simulations was -0.14±0.26 mmHg (p = 0.029), and the correlation was r = 0.99 (p<0.001). The systolic and diastolic ΔP_MCA was 9.9±12.8 (p<0.001) and 3.4±6.2 (p = 0.029), respectively. In the subsequent analysis, only pressures from the stationary simulations are used.

Collateral flow relationship to pressure

There was a correlation between ΔP_MCA and the MCA flow laterality (ΔQ_MCA) (r = 0.47, p = 0.011) (Fig 5). When ΔP_MCA was <5mmHg, there was no MCA flow laterality (6.2±19.8 ml/min, p = 0.14, n = 17). With ΔP_MCA>5mmHg, there was an MCA flow reduction on the ipsilateral side (30.8±27.4 ml/min, p<0.001, n = 11).

Fig 5

MCA flow relationship to pressure.

The ΔQ_MCA versus the ΔP_MCA. The linear correlation was r = 0.47 (p = 0.011).

There was a significant hemispheric flow rate laterality between the contra- and ipsilateral ICA, MCA1 and ACA1, respectively (Table 1). The flow rate laterality between ICA (ΔQ_ICA) and ACA1 (ΔQ_ACA) was highly correlated (r = 0.90, p<0.001). The ΔQ_ICA was correlated to ΔP_MCA (r = 0.85, p<0.001) (Fig 6A). Similarly, there was a strong relationship between the ΔQ_ACA and ΔP_MCA (r = 0.71, p<0.001) (Fig 6B).

Fig 6

ICA and ACA flow relationship to pressure.

(a) The ΔP_MCA versus ΔQ_ICA. The linear correlation was r = 0.85 (p<0.001). (b) The ΔP_MCA versus ΔQ_ACA. The linear correlation was r = 0.71 (p<0.001).

Table 1

Geometry and flow rate features (mean±standard deviation).

Artery	D (mm)	Q (ml/min)	ΔQ (ml/min)	Re	Wo
ICAa	3.8±0.74	193±98	111±137 (p<0.001)	363±169	2.7±0.51
MCA1b	2.6±0.36	130±27	16±26 (p = 0.002)	332±69	1.8±0.26
ACA1a	1.9±0.30	72±68	75±111 (p = 0.001)	291±153	1.3±0.20

D: diameter; Q: flow rate; ΔQ: hemispheric flow rate laterality; Re: Reynolds number; Wo: Womersley number; ICA: internal carotid artery; MCA1: proximal middle cerebral artery; ACA1: proximal anterior cerebral artery. All flow rates were measured from 4D flow MRI data and all diameters from CTA. The significance level was set at p<0.05.

^a Estimated from the whole branches.

^b Estimated from 5 mm long segments.

Mesh sensitivity

The mesh sensitivity analysis revealed a high concordance. The mean difference in ΔP_MCA between the coarser and the intermediate mesh was 0.06±0.36 mmHg (p = 0.68) with a correlation of r = 0.99 (p<0.001). For the finer mesh, the mean difference in ΔP_MCA against the intermediate mesh was 0.16±0.35 mmHg (p = 0.027) with a correlation of r = 0.99 (p<0.001).

Collateral recruitment and perfusion pressure

Without an ICA stenosis, we expect a pressure balance between the left and right MCA. Therefore, ΔP_MCA can be interpreted as the degree of primary collateral system recruitment. Furthermore, this holds with or without an affected ΔQ_MCA, since the ipsilateral flow can be maintained by the autoregulated downstream vasodilatation [23–25]. The ΔP_MCA is visualized as an increased separation between pressures in the MCA Lassen curve (Fig 1). A further increase of the stenosis severity, with an exceeded vasodilatory compensation, could develop a critical blood flow reduction and eventually MCA territory hypoperfusion as well as reduced capacity for washout of emboli [26, 27]. The weak relationship between ΔQ_MCA and ΔP_MCA indicates that this patient group, in general, was on the plateau part of the MCA Lassen curve. However, for the group with ΔP_MCA>5 mmHg there was a significant ΔQ_MCA, indicating that the ipsilateral MCA pressures have decreased towards the flow-reduction region of the curve. For the ΔP_MCA<5 mmHg group, the compensatory capability to maintain MCA flow through collateral recruitment and autoregulatory response was essentially sufficient. A natural MCA flow laterality, due to hemispheric differences in size of vascular territories, challenge the use of MCA flow laterality as an indication of disturbed hemodynamics. A corresponding natural variability in pressure laterality is not expected. Avoiding natural variability and having an interval of a stenosis induced pressure laterality with autoregulatory compensated MCA blood flow, supports that ΔP_MCA has a potential as an early marker for detecting the hemodynamic pressure disturbance.

The ΔP_MCA was highly correlated to the ΔQ_ICA, demonstrating that the stenosis causes the pressure imbalance. This pressure imbalance was expected since a stenosis is causing a recruitment of the collateral system [3, 25]. If ΔP_MCA is an important indicator of disturbed hemodynamics, as suggested in our study, and the CFD analysis to estimate it is unavailable, then the ICA flow laterality is the most reliable indication of this with respect to flow. This is of clinical importance since the degree of stenosis is alone insufficient to characterize the hemodynamic disturbance [23, 28].

Hemodynamic analysis with CFD

A main goal with our CFD model was to incorporate the patient-specific 4D flow MRI measurements and to guarantee a correct flow distribution across the anterior part of the CW. Although patient-specific boundary conditions are important [29, 30], flow outlet boundaries are rare in CFD analysis of CW [8, 16, 31]. A challenge for resistive or pressure outlet boundary condition strategies is to correctly represent the highly resistant capillary bed, which dominates the resistance in the larger arteries. Lumped parameter models or pressures are commonly used instead, but these still require assumptions about territorial resistances [7, 10, 12]. By the model design with flow boundary conditions, we avoided assumptions about territorial vascular resistance and outlet pressures.

Since the posterior circulation of the CW was excluded, the model was limited to patients with an anterior circulation with all arteries visible. Our method should give the correct ΔP_MCA as long as flows and vessel geometries for MCA, ACA1 and ACoA are correct. For patients with bilateral stenoses, the MCA pressure was probably decreased on both sides, potentially reducing the ΔP_MCA, in spite of a distinct bilateral hemodynamic disturbance [32]. For these patients, collaterals might have been recruited from the posterior circulation [28, 33], and a similar CFD analysis between MCA and PCA perfusion pressure could be useful.

The comparison of stationary and time-resolved simulations showed that it was sufficient to use the stationary case to study the mean pressure distribution in the CW. Time-resolved simulations, however, offers the possibility to investigate the pulsatile component of arterial blood flow. Womersley numbers of the major arteries were borderline with respect to using a Womersley profile in the pulsatile analysis [34]. However, we assumed that the effect of these profiles was small and that the inlets, functioning as flow extensions, reduced the boundary effects. A limitation for the time-resolved simulations was the assumption about rigid walls, which excluded dampening effects from the arterial walls [19]. The main advantage of the stationary simulation was that it was performed on a standard PC in less than an hour, while the time-resolved results required up to four days of computational time.

The mesh sensitivity study showed a statistically significant difference in ΔP_MCA of about 2.5% between the intermediate and the finer mesh, but the difference was small and not considered physiologically important in our analysis. Pressure drop is highly correlated to arterial diameters and our results are therefore sensitive to the segmentation. The ACA1 diameters in our study were in agreement with previous studies [32, 35, 36].

Clinical applications

The lower pressure limit for the cerebral autoregulatory function has been shown to have a considerable inter-subject variability [37]. The variability contradicts the use of group defined thresholds in clinical guidance for avoiding hypoperfusion, and rather motivates the use of intra-subject approaches. This is what we have done, when we assumed a normotensive contralateral side and identifies the deviation on the ipsilateral side.

A determined ΔP_MCA can give an insight into the consequences of having a large ICA stenosis or occlusion, where a low cerebral perfusion or hypoperfusion may increase the risk for ischemic stroke [27, 38]. The findings suggest that ΔP_MCA could be an early sign of hemodynamic pressure disturbance in patients with ICA stenosis, observed before MCA flow declines and hypoperfusion symptoms occurs. Specifically, this could be of clinical importance for asymptomatic patients [39], where indications for surgery is not well established, but this has to be confirmed in larger prospective studies.

Another situation for determining perfusion pressure to the brain, is when considering risk for intraoperative hypoperfusion during carotid intervention [40, 41]. Estimation of ΔP_MCA could be important in the planning of surgery, since a reduced perfusion pressure is a risk when the ICA is clamped, to temporarily interrupt the ICA blood flow during the intervention, and ΔP_MCA gives a strong indication of the current hemodynamic pressure disturbance and the collateral capacity.

Limitations

The main limitation was the lack of control group. However, this was a relatively large patient-group for a CFD study with stroke patients and the range of pressure lateralities showed a distribution within the group of patients with ICA stenosis. A control group would have added knowledge about the expected variation without a stenosis. Additionally, the stroke patients recruited for the study were limited to those eligible for ICA surgery or MRI, which excluded asymptomatic cases and patients with the most severe stroke symptoms. Ischemic regions could potentially increase the MCA flow laterality, while the effect on pressure laterality should decrease. From the data in Table 1, we found that the MCA flow laterality was about 12% of the mean MCA flow. For the corresponding average perfusion pressure, we used the approximate MAP of 93 mmHg and assumed an intracranial pressure of 11.6 mmHg [19]. The mean pressure laterality found in our paper was about 8% of this approximated cerebral perfusion pressure. The relative pressure laterality is therefore lower compared to the relative flow laterality. In addition to the effect of autoregulation, the reason for the inconsistency could be that lesioned tissue reduces the needed blood flow and increased the resistance on the symptomatic side. The volumes of the ischemic regions were not measured in the included patients, but the low mRS and NIHSS scales indicated that these volumes should have been small.

The accuracy of the CFD results was primarily limited by the imaging modality [29, 42]. We used two different image modalities, with scans separated in time, which entails a potential risk for geometry discrepancy. An option was to only use the MRI based data. This is, however, sensitive to low or absent flow and would have caused potential vessel-area inaccuracies [43–45]. Therefore, the CTA data was considered a better alternative for our analysis.

The accuracy of the pressure distributions is also limited by the accuracy of the 4D flow MRI measurements and the segmentation method. We can expect some deviations, caused by partial volume effects, small branches along the arteries and random errors between repeated measurements. We tried to reduce these errors by using a flow segmentation method which we have shown is less sensitive to different cross-sectional areas and flow rates [13]. In addition, 4D flow MRI is likely less affected by partial volume effects compared to 2D PC-MRI because of its reduced inflow effects causing more tissue-like magnitude in the blood. The variability along the arteries was managed by averaging multiple of subsequent cut planes, as a representation of the average flow in each artery. We have found this to be important for noise reduction [19]. However, the random variation within repeated measurements is always a remaining source of variability, caused by the true physiological variations of flow in cerebral arteries [13].