A Novel Patient-Specific Computational Fluid Dynamics Study of the Activation of Primary Collateral Pathways in the Circle of Willis During Vasospasm

The Circle of Willis (CoW) is a redundant network of blood vessels that supplies oxygen to the brain. The primary collateral pathways that connect the major blood vessels within the CoW can be recruited to mitigate the negative effects of vessel occlusion by changing the direction or magnitude of flow in the communicating segment. However, only 40% of the population has a complete and balanced CoW, i.e., all vessels are present and of typical sizes [1].

When blood is introduced into the subarachnoid space during a bleeding event like cerebral aneurysm rupture, it breaks down into reactive oxygen species and oxidation blood products that interact with the outer vessel wall. This interaction is thought to contribute to vasospasm—an involuntary constriction of the blood vessel that can significantly decrease perfusion to the brain [2]. Blood is introduced into the subarachnoid space primarily at the location of aneurysmal rupture. However, vasospasm can appear in any of the blood vessels comprising the CoW. Vasospasm is detected in 50–90% of patients within two weeks of aneurysmal rupture and peaks in incidence between 3 and 7 days [3]. Clinicians monitor all major cerebral vessels daily for vasospasm.

The recruitment of primary collateral pathways during embolic stroke [4,5] and internal carotid artery occlusion [6,7] is well-documented in clinical literature. In contrast, the role of specific collateral pathways during vasospasm has not been systematically reported. While work by Otite et al. suggests that cerebral autoregulation can be impaired following subarachnoid hemorrhage [8], a deeper analysis is required to quantify changes in flow in the CoW during vasospasm. By modeling the entire CoW, hemodynamics in collateral pathways can be investigated in response to different locations and severities of vasospasm.

Earlier computational fluid dynamics (CFD) simulations in the CoW established methods for incorporating clinical measurements into patient-specific models [9,10]. Cebral et al. demonstrated that a high-fidelity CFD model of a healthy patient can be created using a patient-specific segmentation derived from magnetic resonance imaging, with boundary conditions at the inlets and outlets that convert clinically measured flow rates into Womersley profiles [10]. Solving for mass transport of a passive scalar generated a virtual angiogram. Other authors have investigated different pathologies in the complete CoW: aneurysm formation and dynamics [11–13], the transport of emboli to the CoW [14], and changes in collateral flow rates in response to carotid stenosis or occlusion [15,16]. While this previous work has established some of the methods and techniques to study normal and altered blood flow in the CoW, no patient-specific CFD study in the full CoW has been performed for vasospasm patients.

Mukherjee et al. modeled flow in the CoW with the Lagrangian particle tracking to model the transport of emboli [14]. Their work finds that the anterior communicating artery (Acom) acts as a key collateral pathway in the CoW, with the direction and magnitude of flow varying significantly between different patient anatomies. The ratio of flow between each internal carotid artery to the basilar artery increases in patients missing a posterior cerebral artery (PCA) P1 segment.

Previous computational investigations of vasospasm examined the relationship between the degree of stenosis and the magnitude of the velocity on the flow rate in the middle cerebral artery (MCA) [17,18], the effects of hemodilution on flow rate [19], and the sensitivity of physiological parameters to vasospasm dynamics [20,21]. However, the models often used idealized geometries [17–21] without applying patient-specific boundary conditions [17–22]. Despite the availability of patient-specific transcranial Doppler ultrasound (TCD) measurements of velocity in the major blood vessels, only one computational study of a full CoW in a reduced-order model used this type of data to define the boundary conditions [23]. These studies did not account for the uncertainties present in the model inputs, instead creating a single deterministic simulation from the medical imaging.

Ryu et al. investigated the changes in extracranial velocities in the presence of vasospasm in the major arteries in the CoW for different anatomic variants using a reduced-order model of the vasculature between the heart and CoW [24]. The internal carotid arteries (ICA) predominantly supplied the anterior circulation, and the vertebral arteries supplied the posterior circulation. The velocity ratio of the internal carotid arteries to vertebral arteries increased when the vasospasm was centralized in the posterior circulation and decreased when the vasospasm was centralized in the anterior circulation for all anatomic variants. Right-sided vasospasm resulted in an increased velocity from the left carotid artery compared with the right carotid artery. These results suggest that the flow rates in the primary arteries supplying the brain change significantly in patients exhibiting vasospasm.

This study presents a novel technique to quantify the changes in collateral flow rates in vasospasm patients by incorporating patient-specific computed tomographic angiography (CTA) scans and TCD measurements into realistic CFD models. These models were validated by benchmarking diameters, velocities, and flow rates against literature values. Virtual angiograms created using passive scalar transport were visually compared to clinical angiograms. Optimized model parameters, including vessel diameters and velocities, were determined using Bayesian analysis based on metrics of parameter uncertainty and the principle of mass conservation, and then they were applied to the boundary conditions of the CFD model. Uncertainties in the model inputs can propagate into the model outputs, e.g., the collateral flow rates, so a sensitivity analysis of the collateral flow rates was also performed by analyzing the results for a range of input parameters.

Two patients admitted to Harborview Medical Center in Seattle, WA for subarachnoid hemorrhage due to aneurysmal rupture subsequently developed vasospasm. A CTA was performed upon admission, prior to aneurysm treatment. A digital subtraction angiogram (DSA) was also performed between 0 and 2 days after the pre-operative CTA scan. Hourly neurological examinations in the intensive care unit, as well as daily TCD measurements in all major cerebral vessels, were collected postoperatively to detect vasospasm based on standardized guidelines and criteria [25]. Elevated velocities in the TCD exam were a potential indicator for vasospasm. In response to a significantly abnormal TCD exam or a change in neurological condition, a CTA and a DSA were performed to verify the presence of vasospasm prior to intra-arterial treatment.

For each patient, two conditions were modeled using CFD simulations: baseline and vasospasm. For the baseline condition, the pre-operative CTA scan as well as the earliest available TCD exam and angiogram were used. For the vasospasm condition, the same data collected on the date of vasospasm detection were used.

The internal carotid arteries, basilar artery, anterior cerebral arteries (ACA), MCA, PCA, ophthalmic arteries (OA), and superior cerebellar arteries (SCA) were segmented (Fig. 1, left). The primary collateral pathways were defined as the arteries connecting these major blood vessels within the CoW: the Acom, anterior cerebral artery A1 segments (ACA A1), posterior communicating arteries (Pcom), and the posterior cerebral artery P1 segments (PCA P1) (Fig. 1, right). Secondary collateral pathways, such as the pial arteries or more distal collaterals, were not considered in this study.

Patient 1 had a complete CoW with two large Pcoms and a duplicated right SCA (Fig. 1, left), which is a normal anatomic variant. An aneurysm in the left MCA M2 segment was treated with a surgical clip, which was outside of the segmented vasculature. Patient 2 had a fetal-type PCA, i.e., the right PCA P1 segment was missing, and a redundant right PCA was present (Fig. 1, right). The aneurysm was located at the Acom and treated with a surgical clip.

Due to the resolution of the CT scans, the images were interpolated by factor of 4 using a B-spline interpolation method in the open-source tool three-dimensional Slicer². The original pixel size, which was approximately 0.4 mm, was recorded for each CT scan. The resulting images were imported into SimVascular³, an open-source segmentation tool, to trace the centerlines of the vessels. Vessel diameters were measured by applying circular contours closest to where the TCD measurements were recorded. The diameters identified by clinicians as vasospastic on the vasospasm CTA scan, were expected to decrease in size compared with the baseline CT scan.

During the TCD exam, a probe was placed in an acoustic window of the skull and was angulated to align the Doppler beam with different vessels in the CoW [26]. Once the sample volume—the region of the intracranial cavity insonated by the ultrasound beam—was placed inside the vessel of interest, the velocity signal acquired by the TCD sensor was recorded. Measurements for some vessels were recorded at multiple depths, where a lower depth indicated that the sample volume was closer to the transducer than the midline of the brain.

The TCD readings used in the study were collected at the bilateral MCAs, ACAs, PCAs, extracranial ICAs, and basilar artery. The maximum velocity detected from the ultrasound return was determined to most closely correspond to the centerline velocity in the vessel. For each segment, the appropriate waveform was selected based on the depth, the magnitude of the mean velocity, and the quality of the measurement. For the ICAs and PCA P2 segments, there was only one available measurement. For the ACA A2 segment, the deepest measurement was selected as closest to the A2 location. For the basilar artery and MCAs, the highest quality waveforms with the highest mean velocity were typically selected, corresponding to the most favorable angle of insonation.

For vessels with diffuse vasospasm, the velocities were higher than the baseline measurement along the entire segment due to the consistently smaller diameter. In vessels with focal vasospasm, the diameter was significantly smaller more proximal to the CoW than the diameter at the outlet, resulting in lower velocities near the outlet. In this case, the lower velocity measurements were selected for the model boundary conditions since they corresponded more closely with the outlet location. The velocities at the area of focal vasospasm were higher in the CFD due to the constriction and matched the higher TCD velocities at this location.

An in-house automated boundary-detection algorithm in matlab based on pixel intensity thresholding converted the TCD exam tiles into a time-dependent waveform for the centerline velocity (Fig. 2). The resulting data were downsampled to smooth the waveform. Separate cardiac cycles were averaged together, and the standard deviation between the cycles was calculated to represent the intercycle variability. A Fourier transform with 15 modes reconstructed the waveform to ensure periodicity of the cycle. The waveforms were synchronized in time by matching the end diastolic velocities corresponding to the minima between the waveforms. The centerline velocities were used to define Womersley profiles for the boundary conditions in the CFD model as described in Sec. 2.3.2. The Bayesian analysis described in Sec. 2.2 determined the optimal values of the average velocities, and the centerline velocities were scaled accordingly before being applied in the CFD model.

Due to the pixelation of the CT scans and the intercycle variability in the TCD waveforms, a degree of uncertainty was associated with the diameters and average velocities that defined the boundary conditions in the model. In total, there were 18 model parameters: nine diameters measured from the CT scan and nine average velocities associated with TCD measurements of the basilar artery, ICAs, MCAs, ACAs, and PCAs. Using Bayesian analysis, this study combined the knowledge of uncertainties in the clinical measurements with the application of mass conservation to find an optimized set of parameters for the CFD simulation. Performing the Bayesian analysis involved three steps: defining the priors, determining the likelihood, and maximizing the posterior.

Priors accounted for the uncertainty in the original set of parameters. The prior for each parameter was defined with a Gaussian distribution of probability, where the mean value was the initial parameter measurement and the standard deviation was the uncertainty in the parameter. The uncertainties for the diameters and velocities were defined as half of the CT pixel size and the intercycle variability in the TCD waveform, respectively. The priors for the nine diameters and nine average velocities were multiplied together.

The flow rates of one OA and one SCA should be about 10 mL/min [27] and 20 mL/min [28], respectively. The right side of Eq. (1) should evaluate to about 60 mL/min in patient 1 with the complete CoW and about 75 mL/min in patient 2 with a fetal PCA variation with 15 mL/min flowing to the redundant PCA. The likelihood function was prescribed with a Gaussian distribution where the mean value was the target flow rate (60 mL/min or 75 mL/min), and the standard deviation was 10 mL/min.

The posterior was obtained by multiplying the priors and the likelihood together. Figure 3 shows an example of a two-dimensional Bayesian analysis, which is a simplified version of the 18-parameter model used in this study. Maximizing the posterior identifies the parameter set that is both similar to the initial set and that satisfies mass conservation. An in-house script coded in matlab applied an iterative method that maximized the posterior and determined the optimal parameter set, i.e., the scaled diameters and velocities for the three inlets and six outlets. The centerline velocities from the TCD waveforms were scaled to match the optimized average velocities, and the optimized diameters were applied in the final segmentation, summarized in Fig. 4.

In SimVascular, a three-dimensional surface representation was created by lofting circular contours along centerlines. Flow extensions were added to the ICAs (8 mm), basilar artery (4 mm), MCAs (4 mm), ACAs (4 mm), PCAs (4 mm), OAs (1 mm), SCAs (1 mm) and, if present, the duplicate PCA (1 mm). The finalized surface mesh was converted to a volumetric mesh in starccm+ version 12.06 (Siemens, Munich, Germany) consisting of tetrahedral elements (base size = 0.55 mm, surface size = 0.44 mm) and four boundary layer elements. The final meshes consisted of about 4 × 10⁶ elements. A grid resolution study was performed by comparing the average flow rates and maximum velocities in the collateral vessels between the meshes. These metrics for the final mesh exhibited less than a 5% error with the finest mesh comprising about 22 × 10⁶ elements.

The incompressible Navier–Stokes and continuity equations were solved using ansysfluent 18.0 (Ansys, Canonsburg, PA), a finite volume pressure-based solver. Blood was modeled as a Newtonian fluid with a density of 1050 kg/m³ and a viscosity of 0.0035 Pa·s. The vessel walls were modeled as rigid. TCD measurements of the centerline velocity were scaled by the factors determined in the Bayesian analysis, then converted into Womersley profiles [29] applied with user-defined functions for the three inlets (the ICAs and basilar artery) and six outlets (MCAs, ACAs, and PCAs). The OAs were prescribed with a waveform that generated a net flow rate between 10 and 15 mL/min. The remaining boundaries (the SCAs and, in patient 2, the redundant PCA) were prescribed with resistance–capacitance boundary conditions. The flow should be evenly split between the left and right sides in the SCAs (∼20–25 mL/min per side) and to meet a target value of 15 mL/min in the redundant PCA. A CFD simulation was computed with initial estimates of resistance that result in the target flow rates. If the flow rate was higher than the expected value in a vessel, the resistance was increased. The CFD simulations were rerun until all resistances were tuned to generate the target flow rates.

The time-step selected was 0.001 s for the cardiac cycles of 1.59 s and 1.47 s for the patient 1 baseline and vasospasm conditions, respectively, and 1.11 s and 1.05 s for the patient 2 baseline and vasospasm conditions, respectively. These periods were calculated by averaging the periods of the nine TCD waveforms for each patient and clinical condition. The residuals were set to 1 × 10⁻⁴ for mass conservation and 1 × 10⁻⁶ for the three components of velocity. The pressure-implicit with splitting of operators pressure–velocity coupling scheme, with second-order pressure and momentum spatial discretization and second-order implicit time discretization, was employed.

The CFD simulations were computed on four nodes of Hyak—the University of Washington distributed memory massively parallel computing platform—with 28 processor dual Intel Xeon CPUs, or on the Pittsburg Supercomputer Center through the National Science Foundation Extreme Science and Engineering Discovery Environment, for simulation times of about 10 h. The simulations were run for two cardiac cycles to eliminate transient behavior associated with initialization, and the results were analyzed for the third cardiac cycle.

Virtual angiograms were performed by injecting a passive scalar into the basilar artery and the ICAs. An advection–diffusion equation was solved to model the transport of radiopaque contrast. The intensity of the contrast indicated the magnitude of the flow through the vessel, with lighter regions indicating limited flow. Details on the computational implementation are available in the Supplemental Materials on the ASME Digital Collection.

The diameters [30–39], velocities [40–47], and flow rates [44–56] from studies of healthy patients were compared with the simulation values. Literature flow rate measurements in the smaller vessels, such as the ACAs, are likely less accurate than those in the large vessels, such as the ICAs, due to the relatively large voxel size in the phase-contrast magnetic resonance imaging measurements compared with the small vessel diameters [54]. For patient 2 with the fetal PCA, the flow rates corresponding to this variation for the ICA and basilar arteries were plotted [51,53,55].

Some uncertainties in the measurements of the diameters and velocities in the primary inlets and outlets are easily quantified, such as the CT pixel size and intercycle variability in the TCD measurement. However, other uncertainties in the input parameters, such as a nonzero Doppler angle in the collection of the TCD velocities and the determination of the correct measurement location of the vessel diameter during segmentation, cannot be easily quantified. A sensitivity analysis can determine how all types of uncertainty simultaneously propagate into uncertainties in the collateral flow rates.

The collateral flow rates were estimated by solving a linear system of equations representing mass conservation at each bifurcation in the CoW. For the complete CoW anatomy in patient 1, there were seven bifurcations and seven unknown collateral flow rates. However, the solution is not unique. Therefore, a Moore–Penrose pseudo-inverse was used to solve the system of equations. For the fetal PCA anatomy in patient 2, the number of equations exceeded the number of unknowns, creating an over-defined system solved with a least-squares approximation.

One thousand different combinations of initial diameters and mean velocities that could vary up to 10% from their optimal values were introduced into the network. Bayesian analysis was performed on each initial parameter set to generate a final parameter set that satisfies mass conservation. The collateral flow rates associated with the final parameter set were computed by solving the linear system of equations. The resulting changes in the collateral flow rates were represented using violin plots. The collateral flow rates corresponding to the optimized parameter set were associated with the maximum width of the violin. Positive values of collateral flow rate corresponded to the following directions: left to right in the Acom, A1 to A2 in the ACA, anterior to posterior in the Pcom, and P1 to P2 in the PCA.

Figure 4 summarizes each stage of the workflow used to create the computational model: taking initial measurements of vessel diameters and velocities, leveraging Bayesian analysis and benchmarking to determine the final model parameters, applying the final parameters in the simulation, and generating virtual angiograms and violin plots representing the sensitivity analysis for the collateral flows.

For both patients, the locations and severity of vasospasm as well as the direction and magnitude of the collateral flows were determined by reviewing clinical DSAs with experienced neurosurgeons. A communicating vessel was considered to be activated if the direction or magnitude of the flow changed to help supply the cerebral territories that were affected by the vasospastic vessels. A supplying artery to the CoW—an ICA or basilar artery—was considered to be recruited if the flow rate notably increased during vasospasm compared to baseline.

Patient 1 had a complete CoW with two large Pcoms (Fig. 1, left). At baseline, the anterior circulation strongly supplied the posterior territory via the bilateral Pcoms, and there was some reverse flow in the left PCA P1 segment, i.e., flow from left to right (Fig. 5, top left). The locations of vasospasm were centralized in the anterior circulation, as shown by the valve icons (Fig. 5, top right), where the lightest valve represents mild vasospasm and the darkest valve represents severe vasospasm. During vasospasm, the basilar artery was recruited to help supply the anterior circulation through reverse flow in both Pcoms, as shown by the large double arrows (Fig. 5, top right). Additionally, the reverse flow was no longer present in the left P1 segment due to the strengthening of the flow from the basilar artery.

Patient 2 had a fetal-type PCA on the right side (Fig. 1, right). At baseline, most of the collateral pathways exhibited typical directions of flow with the exception of weak reverse flow in the left Pcom (Fig. 5, bottom left). The locations of vasospasm were centralized in the MCA and ACA territories with more severe vasospasm in the right supraclinoid ICA than the left one (Fig. 5, bottom right). During vasospasm, the left ICA was recruited, as shown by the large double arrows, to help supply the right side via bidirectional flow in the right ACA A1 segment.

Figures 6 and 7 plot the diameters, velocities, and flow rates for patients 1 and 2, respectively, as described in Sec. 2.4. The shaded Gaussian distributions represent the literature values, the dashed black line and solid black line represent the Bayesian priors for baseline and vasospasm, respectively, and the black point represents the optimal value used to define the CFD boundary conditions. The standard deviations of the priors were half of the CT pixel size for the diameter and the intercycle variability for the velocity. The valve icons indicate which vessels were in vasospasm, and the bold arrow represents the main inlet artery that was recruited during vasospasm (Fig. 5). The standard deviation for the flow rate reflects how high or low the flow rate could be given an increase or decrease in the diameter by half of the CT pixel size. The flow rates scaled linearly with the velocity and with the square of the diameter (Eq. (2)), making them more sensitive to changes in diameter than velocity. The Gaussian distributions for flow rates were wider than the distributions for the diameter and velocity priors due to this square relationship between changes in diameter and the resulting flow rate. The Gaussian distributions for the vasospasm flow rates (Fig. 6, solid line) were wider than the baseline distributions; changes in flow rate due to the change in diameter were then multiplied by the higher vasospasm velocity, resulting in a wider distribution.

For patient 1 at baseline, most of the diameters, velocities, and flow rates from the simulation (Fig. 6, dashed lines with a black point) fell within typical literature ranges (Fig. 6, shaded Gaussian distributions). While the diameters of the PCAs were larger than the typical literature values, their corresponding velocities were lower, resulting in physiologically realistic flow rates. The velocity in the basilar artery was especially low, resulting in a flow rate significantly below typical literature values. The Bayesian analysis resulted in smaller outlet diameters and larger inlet diameters as well as higher velocities in the ICAs (Fig. 6, black points on dashed lines). The remaining velocities were relatively unchanged due to their low intercycle variabilities.

During vasospasm for patient 1, the diameters decreased in the left MCA and in both ACAs compared to the baseline condition (Fig. 6, distributions with valve icons). The decrease in cross-sectional area of these vessels caused their velocities to increase. This combination of the decreases in diameter and increases in velocity resulted in similar flow rates to baseline. The Bayesian analysis increased the MCA and ACA diameters and decreased the ICA diameters, with other parameters remaining unchanged. In the left ICA, the velocity increased significantly but was accompanied by a decrease in the diameter, resulting in a flow rate similar to baseline. In the basilar artery, both the diameter and velocity increased, resulting in a significantly higher flow rate than baseline (Fig. 6, black arrow at bottom right). The basilar flow rate increased from supplying 8% of the total perfusion before vasospasm to 22% during vasospasm due to the recruitment of the posterior circulation to help supply the anterior circulation in vasospasm via reverse flow in both Pcoms (Fig. 5, top right).

In patient 2 at baseline, most of the diameters (Fig. 7, dashed lines with a black point) were higher than the literature values, which was balanced by lower velocities, resulting in similar flow rates to those reported in the literature. The availability of literature values for the ICAs and basilar arteries for fetal PCA anatomy is limited, but the simulation values agreed with at least one of the referenced studies [51,53,55]. The Bayesian analysis did not noticeably change the parameter values.

Vasospasm in both ACAs and MCAs in patient 2 (Fig. 7, distributions with valve icons) was characterized by smaller diameters and higher velocities. The velocity in the left PCA was higher than literature values, which could be explained by vasospasm near the P1-P2 bifurcation or by the recruitment of the left Pcom as a collateral pathway. Both the diameter and velocity in the left ICA increased, resulting in a significantly higher flow rate (Fig. 7, bottom right). The flow rate in the left ICA increased from supplying 45% of the total perfusion before vasospasm to 55% during vasospasm, helping supply the right side of the CoW through reverse flow in the right ACA A1 segment.

Figure 8 compares the clinical (left column) and virtual (right column) angiograms for patient 1 at baseline (top row) and vasospasm (bottom row). At baseline, the left ICA was injected with contrast with an angiographic view from front to back, i.e., the right side of the CoW appears on the left side of the figure. In both the clinical and virtual angiograms, the presence of contrast in both ACA A2 segments indicated that the flow in the Acom was from left to right (see Fig. 8, top row). In the posterior circulation, the reversal of flow in the left PCA P1 segment in the clinical angiogram was inferred from the presence of contrast in the right PCA and right SCA (Fig. 8, top left). The virtual angiogram exhibited the same result, with reverse flow in the left PCA P1 filling the right PCA (Fig. 8, top right).

For the vasospasm condition in patient 1, the basilar artery was injected with contrast with an angiographic view from the side, i.e., the posterior circulation appears on the left side of the figure (Fig. 8, bottom row). Due to the presence of vasospasm in the anterior circulation, the basilar artery was recruited to help maintain physiological perfusion; blood flowed from the posterior circulation to the anterior circulation via both Pcoms for part of the cardiac cycle. Reverse flow in the right Pcom was stronger than in the left, likely due to vasospasm present in the left Pcom. The strong reverse flow in the right Pcom resulted in a small amount of contrast filling the right MCA, seen in both the clinical and virtual angiograms (Fig. 8, bottom row).

Figure 9 shows the left ICA injection in patient 2 during baseline (top) and vasospasm (bottom), with an angiographic view of front to back, i.e., the right side of the CoW appears on the left side of the figure. At baseline, contrast was present in both ACA A2 segments, corresponding to flow in the Acom from left to right (Fig. 9, top row). No reverse flow was present in the right ACA A1 segment. During vasospasm, blood continued to flow strongly from left to right in the Acom. Unlike at baseline, the flow in the right ACA A1 reversed direction, as evidenced by the presence of contrast in the vessel connected to the Acom (Fig. 9, bottom left). Some flow from the right ICA injection (not presented here) was antegrade, i.e., toward the ACA A2 segments, indicating bidirectional flow in the right ACA A1 segment. The left ICA was recruited to help supply the right side of the CoW due to the severe vasospasm in the right supraclinoid ICA compared to the moderate vasospasm in the left supraclinoid ICA. Due to the fetal PCA, the basilar artery and right-sided anterior circulation were not connected, preventing the basilar artery from being recruited to help supply the vasospastic right MCA (Fig. 5, bottom).

In patient 1 during vasospasm, the velocity in the basilar artery increased while velocities in the ICAs remained similar to baseline, aligning with findings from Ref. [24]. The recruitment of the left ICA in patient 2 in the presence of severe right-sided vasospasm agrees with the trends exhibited by Ryu et al.

Figure 10 shows the results of the sensitivity analysis quantifying the changes in collateral flow rates caused by changes in both the vessel diameters and velocities, as described in Sec. 2.5. The widest part of the violin, seen at the centerline with a white point, represents the median value in the distribution and corresponds to the optimal parameter value used in the CFD simulation. The black bar at the centerline of the violin shows the interquartile range. A violin with a wide range of values indicates that the collateral flow rate is more sensitive to variations in diameter and velocity than a violin with a smaller range.

In patient 1 at baseline (Fig. 10, left), the flow was weakly from left to right in the Acom, strongly from A1 to A2 in both ACAs, strongly from anterior to posterior in both Pcoms, and strongly from P1 to P2 in the right PCA. In the left PCA P1 segment, the net flow rate was close to zero, resulting in bidirectional flow (Fig. 5, top left). During vasospasm, the flow rates in the Acom and ACA A2 segments were almost identical to the baseline values. However, the net flow rates in both Pcoms were close to zero, indicating bidirectional flow. The flows in both PCA P1 segments were strongly from the P1 to P2 segments due to the recruitment of the basilar artery, eliminating the reverse flow seen in the left PCA P1 segment at baseline (Fig. 5, top right). Despite the uncertainties in the input parameters, the same trends for the flow directions hold for different parameter sets and match the clinically observed flow directions (Fig. 5, top).

In patient 2 at both baseline and vasospasm (Fig. 10, right), the flow rates were from left to right in the Acom and antegrade in the left ACA A1, the right Pcom, and the left PCA P1. At baseline, the flow was weakly bidirectional in the left Pcom (Fig. 5, bottom left), which resulted in a low net flow rate in this segment. During vasospasm, the strength of the flow increased in the Acom and left ACA A1 segments, resulting in bidirectional flow in the right ACA A1, and the median flow rate fell to zero. The flow in the left Pcom increased, resulting in antegrade flow. For all collateral pathways, the directions and magnitudes of flow were similar across different parameter sets, agreeing with the clinical results (Fig. 5, bottom).

The present study models only two patients representing two different CoW anatomies. Future work will study a cohort of patients with a broader range of anatomical configurations and with different locations of vasospasm. This expanded study will help develop an understanding of how arteries supply the CoW under both healthy and vasospasm conditions, as well as how the collateral pathways are recruited to mitigate the negative effects of vasospasm.

While this study does not include smaller vessels that branch from the CoW, like the anterior choroidal arteries and anterior inferior cerebellar arteries, these arteries are typically too weakly perfused to detect in the CT scan. In addition, the literature indicates that these arteries are negligible in size compared to the SCAs and OAs [28,57,58] and thus have limited effect on the intracranial hemodynamics [59].

During ICA occlusion or MCA stroke, retrograde flow in the OAs can help supply the CoW [60,61]. The ophthalmic arteries are not insonated during vasospasm, i.e., the direction of flow during vasospasm is unknown. Additionally, flow reversal has not been documented in the literature in patients exhibiting vasospasm. Antegrade flow in the OAs was assumed.

Transcranial Doppler ultrasound measurements are collected by insonating a blood vessel with an ultrasound beam oriented as closely as possible with the centerline of the vessel. For measurements collected at high Doppler angles, the TCD reading could underestimate the velocity by a factor of $(1/cos(θ))$⁠, where $θ$ is the Doppler angle. The literature is not in agreement on whether velocity measurements from TCD and angle-corrected color duplex sonography are significantly different from each other [43,46,62]. In this study, the velocities were not angle-corrected because the uncertainty in the Doppler angle was too large, and the vessel flow rates fell within physiological ranges using the uncorrected TCD velocities, as seen in Figs. 6 and 7.

For the sensitivity study in patient 1 with the complete CoW, the solution to the system of equations is not unique. The linear system of equations does not account for the differences in resistance between the collateral pathways, which could be significant in the presence of vasospasm. The resistance affects through which pathways the blood preferentially flows. Therefore, the collateral flow rates calculated from this method deviated from the CFD flow rates. The percent differences were about 20±10% for the larger flow rates and about 7 mL/min for the smaller flow rates. These percent differences in collateral flow rates fell within the uncertainties associated with the primary inlet and outlet flow rates due to the pixelation of the CT scans, as seen in Figs. 6 and 7. Meaningful trends during baseline and vasospasm could still be detected in the violin plots in Fig. 10. For the sensitivity study in the fetal PCA variation for patient 2, the collateral flow rates found using the least-squares method agreed closely with the flow rates from the CFD.

Additional collateral pathways in the distal circulation such as the pial arteries can be recruited to maintain typical perfusion during ICA occlusion [61,63] and stroke [64]. CT scans and TCD velocity measurements, which could quantify the extent to which these collaterals are recruited during vasospasm, were not available. Therefore, the primary collateral arteries within the CoW were exclusively modeled.

This study assumed that the maximum velocity from the TCD measurement corresponded to the centerline velocity of a symmetric Womersley profile. The velocity profiles in the major arteries could be asymmetric due to the tortuosity of the blood vessels. Therefore, some error could be introduced into the model by assuming a symmetric velocity profile. The sensitivity analysis accounts for all types of uncertainties, including the error of assuming a symmetric velocity profile. Flow extensions were added at the end of each artery where the symmetric Womersley profiles were applied as boundary conditions. These extensions ensured that the symmetric boundary conditions did not inhibit the development of asymmetric velocity profiles in the proximal locations of the CoW.

Clinical data were acquired for the baseline condition upon admission for a subarachnoid hemorrhage and before the development of vasospasm, whereas the data used to benchmark the diameters, velocities, and flow rates were acquired in healthy patients. In the absence of stroke or vasospasm at that time of admission, the vessel calibers and brain perfusion are estimated to represent a baseline for the purpose of comparing against changes during vasospasm, compatible with measurements in a healthy patient.

This study presents a novel technique to incorporate TCD measurements into patient-specific CFD simulations of blood flow in the CoW before and during vasospasm. Bayesian analysis used information about the uncertainties in the vessel diameters and velocities, as well as mass conservation, to identify optimal parameters that were applied as CFD boundary conditions. A sensitivity study demonstrated that the magnitude and direction of collateral flow rates were similar across different parameter sets, indicating that the trends in collateral flow rates seen in the CFD simulations are valid despite uncertainties in the input parameters. In patient 1 who had a complete CoW, the posterior circulation was recruited in vasospasm to help supply blood to the anterior circulation via retrograde flow in both Pcoms. In patient 2 who had a fetal PCA, the left ICA was recruited in vasospasm to help supply the right circulation via the right ACA A1 segment. Virtual angiograms showed good agreement with the clinical ones in both patients. This analysis will be extended to a cohort of patients to relate the activation of collateral pathways with the anatomy of the CoW and the locations of vasospasm.

The authors would like to thank Do Lim for collecting clinical CT scans and TCD measurements. We would like to express our gratitude to Dr. Juan Carlos del Alamo for his advice to explore Bayesian analysis, as well as to solve the system of linear equations for the complete CoW with a pseudo-inverse. We would like to acknowledge Shumaila Ahmad for assisting in the segmentation of the patients' CT scans. This work was facilitated through the use of advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system and funded by the STF at the University of Washington. This work was funded by the generous support of the Catchot family.

• National Science Foundation Graduate Research Fellowship Program (NSF GRFP) AHRQ (Grant No. 1 R18 HS 026690; Funder ID: 10.13039/100000133).

• National Institute of Neurological Disorders and Stroke (Grant No. R25NS079200; Funder ID: 10.13039/100000065).

• National Institutes of Health (Grant No. NIH/NINDS R01NS105692; Funder ID: 10.13039/100000002).

• Student Technology Fund at the University of Washington.

• Extreme Science and Engineering Discovery Environment (XSEDE), National Science Foundation (Grant No. ACI-1548562; Funder ID: 10.13039/100000001).

A =

cross-sectional area of the vessel

d =

vessel diameter

Q =

volumetric flow rate

V_ave =

average velocity in space and time in the vessel

V_centerline =

velocity at the centerline of the vessel averaged over time

$θ$ =

Doppler angle

ACA =

anterior cerebral artery

Acom =

anterior communicating artery

CFD =

computational fluid dynamics

CoW =

Circle of Willis

CTA =

computed tomographic angiography

DSA =

digital subtraction angiography

ICA =

internal carotid artery

MCA =

middle cerebral artery

OA =

ophthalmic artery

PCA =

posterior cerebral artery

Pcom =

posterior communicating artery

SCA =

superior cerebellar artery

TCD =

transcranial Doppler ultrasound