Governing equations.

The governing equations for blood flow in 1D arteries are derived from the principle of conservation of mass and momentum, as follows [5,6]: (1) (2) Herein, t represents the time; x indicates the axial coordinate along the artery; A, Q, and P denote the cross-sectional area of the artery, volumetric flow rate, and internal pressure, respectively; ρ = 1060 kg m⁻³ indicates the blood density; and K_R = 22πμ/ρ represents the resistance parameter [4] with blood viscosity μ = 0.0047 Pa s. The (A, Q) system in Eqs (1) and (2) is closed by the relationship between the pressure and cross-sectional area, derived from Laplace’s law [3, 5] as follows: (3) where A₀ denotes the cross-sectional area at reference pressure P₀ = 85 mmHg, h indicates the arterial wall thickness, E represents Young’s modulus, and σ denotes Poisson’s ratio. In this study, σ was set to 0.5, and Eh in Eq (3) was assigned based on the empirical relationship with the arterial radius [3].

In a stenotic artery, the abrupt changes in the cross-sectional area cause a large pressure loss associated with flow separation and reattachment. As the 1D model alone cannot completely describe such a pressure loss, a stenosis model, which relates pressure loss across the stenosis (ΔP) to geometric parameters, was coupled with the 1D model. We employed the model reported in previous studies [42–44]: (4) (5) where R_v denotes the viscous resistance of the stenosis, evaluated considering the axial diameter change D(x); D_n indicates the maximum diameter distal to the stenosis; SR represents the stenosis ratio defined as the percentage reduction in diameter (1−D_s/D_n, with the minimum stenosis diameter D_s); L_s denotes the stenosis length; and indicates the time derivative of Q. The first, second, and third terms in Eq (4) account for the contribution of viscous friction, flow separation, and pulsatility to the pressure loss, respectively. Although coefficients K_t and K_u rely on the stenosis geometry, they have been empirically set to 1.52 and 1.2, respectively, in the literature [43,45]. In this study, we considered K_t as an uncertain parameter ranging between 1.0 and 2.699 [46], whereas K_u was maintained constant at 1.2, owing to its negligible influence on ΔP.

The 0D closed-loop model comprises the peripheral artery, upper and lower body blocks, and heart (Fig 2). The peripheral arteries distal to the 1D terminal arteries are represented by the three-element Windkessel model (RCR circuit). In each upper or lower body block, the capillaries, venules, and veins are modeled as RLC circuits in series. The heart is modeled based on the time-varying elastance method, which provides the inlet boundary condition to the 1D network, generating a closed-loop system. The governing equations for the 0D model are derived by linearizing and integrating Eqs (1)–(3) along the axial direction, as follows [47,48]: (6) (7) where R, L, and C represent the viscous resistance, inertia of blood, and vascular compliance, respectively; and subscripts 1 and 2 denote the quantities upstream and downstream, respectively.

Numerical methods.

The governing equations for the 1D model were solved using the two-step Lax–Wendroff scheme. Bifurcated 1D segments were coupled by enforcing the conservation of mass and total pressure at the bifurcations. As this yields the coupled nonlinear equations (see [41] for detailed formulas), we used the iterative Newton–Raphson method to solve them [6,41]. Furthermore, simultaneous ordinary differential equations in the 0D model were solved using the fourth-order Runge–Kutta scheme. The 1D, 0D, and stenosis models were coupled using Riemann invariants [41].

Patient-specific modeling.

Initially, all model parameters were assigned as reported by Liang et al. [18,41]. Subsequently, certain parameters associated with cerebral circulation were assigned or adjusted based on the patient’s clinical data. The patient-specific parameters included the diameters and lengths of the carotid and cerebral arteries, stenosis model parameters, and peripheral resistances (PRs), which represent the sums of the two resistances in the three-element Windkessel model at the six outlets of the CoW. Additionally, the stiffness and diameter of the aorta were adjusted based on the patient’s age, and the total PR was adjusted to match the measured pressure.

The diameters and lengths of the carotid and cerebral arteries were extracted from medical images (CT or MRI) and directly assigned to the corresponding parameters in the 1D model. In the case of the stenotic artery, R_v, D_n, SR, and L_s in Eqs (4) and (5) were evaluated based on the acquired geometry. Image processing for arterial lumen segmentation, centerline extraction, 3D reconstruction, and calculation of geometric parameters was conducted using in-house software, namely “V-Modeler” [49]. As routine diagnostic imaging generally involves only the diseased region (head and neck in this case), patient-specific geometries for the remaining 1D segments could not be obtained. Therefore, we used the geometries prescribed in the literature for these segments; however, the stiffness and diameters of the aortic segments were modified to reflect age-related changes [50].

At the six outlets of the CoW, including the left and right anterior, middle, and posterior cerebral arteries, the PRs were adjusted through iterative calculations to match the measured flow rates [20,51]. Outflow rates in these arteries were measured using single photon emission computed tomography (SPECT) combined with phase contrast magnetic resonance imaging (PC-MRI) or ultrasound measurement [20]. Initially, we converted the regional brain perfusion map on SPECT images to the flow rates averaged over a cardiac cycle duration at the six outlets, , using vascular territory templates [52]. Subsequently, these flow rates were corrected as follows based on the measured total inflow rate to the CoW while maintaining constant flow distribution ratios among the outlets: (8) Herein, denotes the summation of the flow rates in the three inlets of the CoW (the left and right ICAs and the basilar artery), which is measured either using PC-MRI or ultrasound. Finally, we used in Eq (8) as the target value for the PR adjustment.

The total PR was adjusted to match the patient’s mean arterial pressure measured at the upper arm [51]. This was implemented by changing the PRs of the terminal arteries, excluding the CoW, with the scaling factor relative to the initial values.

Defining inputs and outputs.

Although the 1D–0D model includes a large number of parameters, only some have a significant impact on cerebral circulation. As described in the “Patient-specific modeling” subsection, we set those parameters in a patient-specific manner based on the patients’ clinical data. The parameters include

• Diameters of 22 carotid and cerebral arteries in the 1D model (22 parameters);
• Lengths of 22 carotid and cerebral arteries in the 1D model (22 parameters);
• R_v, D_n, and SR in Eq (4) for each left and right ICA stenoses (6 parameters);
• PRs at the six outlets of the CoW (6 parameters);
• Scaling factor for the total PR (1 parameter);
• Age (1 parameter);

and an uncertain parameter for the stenosis model, which is

• K_t in Eq (4) for each left and right ICA stenoses (2 parameters).

The aforementioned 60 parameters characterize the patient’s anatomical and physiological conditions and significantly affect the cerebral circulation. We selected all these parameters as inputs to capture their influence on the cerebral circulation. Note that we do not select the stenosis length, L_s, as an input. As shown in Eq (4), L_s has two effects on ΔP: one on the third term and the other on the first term via R_v. The effect of L_s on the third term can be ignored because the third term is negligible compared to the other terms. Furthermore, since we selected R_v as the input representative of the stenosis geometry, encompassing the variations of D(x) and L_s, it is unnecessary to select L_s as a separate input to be varied.

Based on the 1D–0D simulation, A(t, x), Q(t, x), and P(t, x) at each axially aligned grid point of the 1D artery can be obtained as the output. However, according to the definition of CH (percentage increase in time-averaged flow rate), the focus lies on the assessment of cerebral circulation as a “time average” in several clinical situations. Therefore, we aimed to construct a surrogate model that predicts hemodynamic quantities averaged over a cardiac cycle duration and limits the outputs to be predicted, which include

• Cycle-averaged flow rates, , in the middle of the carotid and cerebral arteries (22 quantities);
• Cycle-averaged pressures, , in the middle of the carotid and cerebral arteries (22 quantities);
• Mean arterial pressure, which is the cycle-averaged pressure in the middle of the left subclavian artery (1 quantity).

Here, cycle-averaged flow rate and pressure refer to Q and P averaged over a cardiac cycle duration: (9) (10) where t_s and T_c respectively denote the time to start averaging and cardiac cycle duration (fixed as 1 s). In the axial direction, is constant and decreases almost linearly unless there is a significant axial change in . Therefore, and at the middle grid point of each artery can be regarded as the axially averaged quantities in each artery. The aforementioned 45 outputs are the primary clinically relevant quantities describing the cerebral circulation. Consequently, we constructed a surrogate model that defines a mapping from the inputs x∈ℝ⁶⁰ to the outputs y∈ℝ⁴⁵.

Design of experiments.

The input–output paired learning data can be obtained by randomly sampling x∈ℝ⁶⁰ and performing 1D–0D simulations to obtain the corresponding y∈ℝ⁴⁵. In this step, the sampling ranges for x must be adequately prescribed. If the ranges are extremely narrow, the trained surrogate model would be accurate only in limited input space, restricting the model’s coverage. Particularly, the prediction accuracy of the DNN outside the trained range decreases significantly because of its interpolative nature [53]. Therefore, we inferred physiologically reasonable ranges for by investigating the data of the seven patients (Table 1) and reviewing the literature [54–58]. The basic policy was to calculate the mean and standard deviation (SD) of the data to adopt a range that covers the mean ± 3SD (details in S1 Appendix). Table 2 summarizes the ranges of inputs used to generate the learning data.

Simulations.

Learning data were generated considering four scenarios, namely (i) intact ICAs, (ii) left ICA stenosis, (iii) right ICA stenosis, and (iv) left and right ICA stenoses. In each scenario, x was randomly sampled in the prescribed range (Table 2); however, the stenosis parameters for the intact ICA were set as R_v = 0, D_n = D_ICA (diameter of the ICA), SR = 0, and K_t = 0. We sampled 50 000 sets of x for each scenario and obtained the corresponding y from the simulation. Consequently, learning data were generated with the number of samples N_data = 200 000. We observed that certain inputs resulted in unphysical or non-physiological outputs, such as < 0, or reversed flow in terminal arteries. Such samples were replaced with new samples. All simulations were performed on the Oakforest-PACS supercomputer system provided by the Information Technology Center at The University of Tokyo (Tokyo, Japan). The samples were equally allocated to 31 280 CPU cores (Intel Xeon Phi 7250). The total computation time required was approximately 25 h.

Data splitting.

We split the learning dataset into training, validation, and test data in the ratio of 6:2:2. The training data were used to construct the surrogate model, and the prediction accuracy of the model was evaluated using the validation/test data. The validation data were specifically used to determine the stopping point of model training (see the “Model training” subsection), whereas the test data were used to assess the performance of the trained model.

Deep neural network.

We used a fully connected DNN as a regression model to fit the training data. The DNN comprises a total of N_layer + 2 layers: an input layer, a series of N_layer hidden layers, and an output layer (S1 Fig). The input and output layers include nodes equal to the number of inputs and outputs, respectively. Each hidden layer comprises an equal number of nodes, N_node, and each node is connected to all nodes in the adjacent layers. Initially, the values of serve as input to the nodes in the input layers. Subsequently, each node in the first hidden layer receives the weighted inputs, sums them up, adds a bias, and finally applies the rectified linear unit (ReLU) activation. This process continues for each layer up to the last hidden layer. The nodes in the last hidden layer and the output layer are fully connected without ReLU activation. Consequently, the DNN turns into a recursive function, as follows: (11) where y^l, W^l, and b^l denote the output vector, weight matrix, and bias vector of the l-th layer, respectively.

Model training.

The DNN was trained using the data by adjusting the weights and biases to minimize the loss function, which is defined as the mean squared error of the outputs, as follows: (12) where ‖∙‖ denotes the Euclidean norm (l₂-norm of a vector), N_out = 45 indicates the number of outputs, N_sample represents the total number of samples used for evaluation, denotes the outputs in the training data (outputs from the 1D–0D simulation), and indicates the outputs predicted by the DNN. We used the gradient-based algorithm “Adam” [59] for optimization, with an initial learning rate lr. Furthermore, mini-batch training was employed with a batch size N_batch, and batch normalization [60] was applied between the linear transformation and ReLU activation in individual layers.

During the training, the coefficient of determination (R² score), defined as (13) was evaluated based on the validation data at the end of each epoch to assess the performance improvements. Herein, denotes the mean of . The closer the R² score is to 1, the more precise the prediction; R² = 1 if and are equal for all s. For every 100 epochs, we monitored the R² score averaged over the latest 100 epochs; training was terminated if no improvements were observed in three successive evaluations. The weights and biases with the highest R² scores at the epoch were selected for the trained model. The DNN and its training were implemented using “Chainer” [61], which is a Python-based open-source framework for deep learning.

Training the DNN involves certain hyperparameters, namely N_layer, N_node, N_batch, and lr, which are not trainable through the optimization process as they are chosen arbitrarily. Although the choice of hyperparameters affects the prediction accuracy of the trained model significantly, the optimal values cannot be known in advance as they vary considerably depending on the data. Therefore, we conducted a grid search to identify the best combination of hyperparameters in N_layer∈{5, 7, 10, 13}, N_node∈{50, 100, 200, 400}, N_batch∈{300, 1000, 3000, 10000}, and lr∈{10⁻³, 10^−2.5, 10⁻², 10^−1.5}. After 4⁴ = 256 rounds of training, the R² scores evaluated by the models based on the test data were compared, and the best-performing model was selected as the final surrogate model.

The training data were preprocessed to improve the model performance. We normalized the inputs such that their upper and lower bounds were 1 and −1, respectively, and standardized the outputs to ensure zero mean and unit SD. The inputs of the validation and test data were normalized similar to the training data, whereas the outputs were scaled as y′ = (y−μ_train)/σ_train. Herein, μ_train and σ_train denote the mean and SD of the outputs of the training data, respectively. The normalization and standardization (also known as z-score normalization) applied to the data herein constitute standard preprocessing in supervised learning [62].

Model verification.

The surrogate model was verified by (i) assessing the prediction accuracy using 40 000 samples of test data, and (ii) comparing the surrogate model and 1D–0D simulation in terms of the predicted outputs and adjusted inputs of the seven patients (Table 1). In the second step, the procedure for assigning or adjusting the inputs based on the patient’s clinical data during the prediction performed by the surrogate model was identical to that of the simulation (“Patient-specific modeling” subsection). In both steps, we used the mean absolute error (MAE) (14) in addition to the R² score to assess the prediction accuracy of the surrogate model.

Quantity of interest.

CH is defined as an increase of more than 100% in the flow rate of cerebral arteries due to ICA stenosis surgery [35–38]. Therefore, we focused on predicting the cerebral circulation when the stenosis is dilated, to evaluate the percentage increase in outflows of the CoW as follows: (15) Herein, and denote the cycle-averaged flow rates at the six outlets of the CoW before and after dilating the stenosis, respectively. By definition, > 100% represents CH.

Target patient characteristics.

Three patients (Patients 1–3 in Table 1) were included in the surgical outcome prediction. The imaging data (CT), measurements of inflows (PC-MRI or ultrasound) and outflows (SPECT) of the CoW, and mean arterial pressure collected before the surgery were used for the predictions. The evaluated stenosis ratios (SR) for Patients 1, 2, and 3 were 59%, 83%, and 91%, and the corresponding R_v values (evaluated using Eq (5)) were 0.5 mmHg s mL⁻¹, 11.3 mmHg s mL⁻¹, and 66.6 mmHg s mL⁻¹, respectively, exhibiting the same trend as the SR (see S2 Fig for the stenosis geometry). Patients 1 and 3 each had a complete CoW; however, CT images of Patient 2 suggested hypoplasia (missing) of an anterior communicating artery (ACoA). Patient 2 was identified by the surgeon as being at risk for CH based on the collected data. To minimize the potential risk of CH, Patient 2 underwent staged surgery, where the stenosis was pre-dilated using a balloon, followed by complete dilation with a stent after two weeks.

Uncertainty modeling.

We evaluated the uncertainty in the clinical data that were used to assign or adjust the patient-specific inputs. We focused on uncertainties in the arterial diameters and stenosis parameters, which were used directly as inputs, and those in the CoW inflow and outflow measurements, which were used to obtain the target outputs. The arterial length is more robustly measured than the diameter and has a minor effect on flow resistance; therefore, the uncertainty in length was not considered.

In all three patients, arterial diameters and stenosis parameters were obtained through segmentation of the arterial lumen on CT images. The geometry obtained during the segmentation can vary based on the threshold used to determine the boundary. In the case of CT, the lumen boundary spanned 2–3 pixels, and the diameter changed by ±2 pixels based on the threshold used. Therefore, we assumed uncertainty of ±2 pixels (±0.702–0.936 mm, depending on image resolution) with respect to the arterial diameter obtained from the segmentation. Similarly, uncertainties in the stenosis parameters were estimated by considering a 2-pixel uncertainty in the underlying geometry. However, an exception was made for Patient 2, as the ACoA was not recognized on CT images of this patient, suggesting hypoplasia of the ACoA. Nevertheless, we could not rule out the possibility that the ACoA, hidden between the extremely close presence of the left and right anterior cerebral arteries, might have failed to resolve on the images (Fig 3). Therefore, we assumed uncertainty of 0.1–2.6 mm in the ACoA diameter, thereby including the possibility of its absence as well as presence.

Uncertainties in the measured flow rates were determined based on modality. The uncertainties in the measured values were assumed to be ±16%, ±35%, and ±16% for PC-MRI, ultrasound, and SPECT, respectively, based on the literature [19,63–68] and discussions with surgeons. The uncertainty ranges were intentionally overestimated to maximize the chance of including the “true” (yet unknown) value regardless of the modality used.

Uncertainty propagation.

Uncertain inputs and targets were treated as random variables with uniform distribution on the determined interval. To estimate the statistics of the predicted under uncertainties, we used the most straightforward approach for UQ, namely the Monte Carlo method. In each realization, uncertain inputs and targets were sampled from a specified probability distribution, and was predicted through successive steps of “preoperative adjustment” and “postoperative prediction” (Fig 4). In the first step, the PRs of the CoW and scaling factor for the total PR were adjusted to match the predicted outputs to the targets. The samples were rejected if target convergence was not attained. Subsequently, the stenosis parameters were modified to R_v = 0, D_n = D_ICA, SR = 0, and K_t = 0 to reflect the complete dilation of the stenosis. The modified stenosis parameters and adjusted PRs were used as inputs in the subsequent steps to predict cerebral circulation immediately after the stenosis surgery. Finally, was calculated using the flow rates before and after the surgery. The statistics of under uncertainties were estimated using the collected , where N_MC denotes the number of realizations. N_MC was increased sequentially until the statistics of converged. As a basic policy, we increased N_MC by 10 000 and ensured that the change in mean and variance of was within 0.1%. We also confirmed that there was no significant change in the probability of > 100% when N_MC was increased. A detailed description of the algorithm for uncertainty propagation is provided in S2 Appendix.

In this study, we assumed that the surgery did not alter the arterial geometry (except for stenosis) and PRs. This assumption is justified because we aim to predict the cerebral circulation immediately after the surgery. Additionally, autoregulation and remodeling of the cerebral arteries generally prevent an abrupt change in blood flow. Therefore, our assumption is appropriate for predicting the maximum possible , which is the most dangerous surgical outcome in terms of CH.

Sensitivity analysis.

In addition to UQ, we performed SA to measure the impact of each parameter (uncertain input or target) on . We adopted a variance-based global SA proposed by Soboľ [69] to consider the interaction between the parameters. In this method, the impact of parameter x_n on output y is quantified as the Soboľ sensitivity indices [69,70]: (16) (17) where denotes the conditional expectation of y for a fixed x_n; indicates the variance of y; x_−n represents all parameters except x_n; S_n, the first-order sensitivity index, quantifies the independent contribution of x_n to the measured variability of y; and S_T,n, the total sensitivity index, quantifies the overall contribution of x_n to the variability of y, including indirect contributions through interactions with other parameters. A large S_T,n − S_n indicates that the impact of x_n varies significantly with the values of other parameters.

We used Saltelli’s algorithm with the Monte Carlo method to compute the sensitivity indices [70, 71]. The accuracy of the sensitivity indices in terms of the sampling error was assessed by estimating the 95% confidence interval using the bootstrap method [72] with a sample size of 1000. The SA was implemented using the open-source Python library “SALib” [73].

Effect of hyperparameters.

Based on the grid search for 256 sets of hyperparameters, the highest R² score was achieved when N_layer = 7, N_node = 200, N_batch = 3000, and lr = 10^−2.5 (S3 Fig). Combinations with N_node = 200 yielded an overall higher R² score, indicating that N_node affects the R² score more than the other hyperparameters.

N_layer and N_node determine the total number of trainable parameters (weights and biases) of the DNN, whereas N_batch and lr control the gradient and the rate of parameter update, respectively, during the optimization process. To compare the influences of these effects on prediction accuracy, we plotted the R² score with respect to the number of trainable parameters, as illustrated in Fig 5A, using the following equation: (18) where denotes the number of nodes in the l-th layer. Note that the index of summation starts at 2 instead of 1 since the input layer (l = 1) has no parameters. Fig 5A indicates that the R² score has an inverted U-shaped relationship with N_param. The highest R² score was achieved when the DNN contained 262 445 trainable parameters; increasing or decreasing the number of parameters from this optimal number resulted in lower R² scores. The vertical variations in the R² score indicate that the effects of N_batch and lr were relatively small when the DNN comprised an optimal number of parameters. However, the choice of N_batch and lr significantly affected the prediction accuracy when training the DNN with more than 1 million parameters.

Effect of number of training samples.

To investigate the effect of the number of training samples on the R² score, we trained the DNN with different numbers of training samples while maintaining the hyperparameters constant in the optimal combination. Fig 5B compares the networks’ R² scores evaluated using identical test data. Increasing the number of training samples improved the prediction accuracy significantly, particularly with a smaller number of samples. However, the accuracy reached a plateau with 120 000 samples, indicating that the accuracy cannot be improved further with more samples.

Model performance.

The best-performing DNN, trained with hyperparameters N_layer = 7, N_node = 200, N_batch = 3000, and lr = 10^−2.5, was selected as the final surrogate model and verified. Initially, we assessed the prediction accuracy of the model using 40 000 samples of test data. The overall R² scores for the flow rate and pressure were 0.9959 and 0.9973, respectively. On average, the MAE was 2.617 mL/min for the flow rate and 0.7226 mmHg for the pressure, which correspond to approximately 4% and 0.9% of the flow rate and pressure mean absolute values, respectively. A detailed comparison of the flow rate and pressure in each artery predicted by the surrogate model and 1D–0D simulation are illustrated in S4 and S5 Figs.

Furthermore, to verify the model accuracy using the patients’ clinical data for assigning and adjusting the inputs, we compared the surrogate model and 1D–0D simulation in terms of flow rates, pressures, and adjusted PRs of the CoW for the seven patients (Fig 6). The flow rates at the six outlets of the CoW were excluded from the evaluation, as they matched the measured flow rates. As indicated in Fig 6, the outputs and adjusted PRs from the surrogate model were in agreement with those from the simulation. Even in the case of patient-specific predictions that involved iterative adjustment of inputs, the flow rate and pressure errors were comparable to those evaluated using the test data.

Fig 7 compares the surrogate model and simulation in terms of the time required for a single prediction. On a single CPU core (Intel Core i9-9900K, 3.6 GHz), the surrogate model achieved a prediction time of several milliseconds, reducing the computation time of the simulation by a factor of over 43 000. Furthermore, the surrogate model exhibited excellent parallelization performance, particularly when executed on a GPU (NVIDIA GeForce RTX2080 Ti), and significantly reduced the computation time per prediction. As illustrated in Fig 7, the computation time was only five times longer when the surrogate model performed 10 000 predictions on a GPU than a single prediction on a single CPU core. Parallelization on the GPU was performed using the built-in backend of Chainer [61] for CUDA-based parallel matrix operations. The latest deep learning libraries, including TensorFlow, Keras, PyTorch, and Chainer, support GPU execution using their built-in backends, allowing easy parallelization of matrix operations in training and predictions.

Flow rate increase following the stenosis surgery.

The percentage increase in flow rate () following the ICA stenosis surgery was evaluated for Patients 1–3, considering the uncertainties in arterial diameters, stenosis parameters, and target flow rates. The number of realizations (N_MC) to obtain the statistics of was set to 100 000. The time required for the UQ was a few minutes on a single CPU core, which was shorter than the time required for a single prediction using the 1D–0D simulation. The time was reduced to less than a minute when executed on a GPU.

Although we obtained at each outlet of the CoW, we focused only on the results at the middle cerebral artery (MCA) on the stenosis side, which was subjected to the largest . Fig 8 depicts the probability density of the predicted for Patients 1–3 along with the values from the deterministic 1D–0D simulation (represented as triangles). Additionally, the figure depicts the interval, mean, and mode (the value with the highest frequency) of and the probability of being more than 100%. A negative indicates a decrease in the flow rate following the surgery. This situation may be rare in actual patients with severe stenosis; nonetheless, it is not non-physiological, as observed in certain clinical cases [36].

Overall, uncertainties in the clinical data generated large variations in the predicted . Based on the comparison of patients’ results, we observed that the mode of was close to the predicted by the deterministic simulation and higher when stenosis was more severe (larger SR and R_v). In all patients, the distribution of was skewed to the right, with a higher mean than the mode. The distribution of was spread extensively to large values in Patients 2 and 3, wherein the stenosis was more severe than in Patient 1. The increase in the prediction uncertainty in with higher stenosis severity is attributed to the 2-pixel uncertainty considered for the arterial diameter. With the same variation width of diameter, the uncertainty in R_v (Eq (5)) and SR (= 1−D_s/D_n) increases with a smaller diameter, leading to a larger uncertainty in .

However, the comparison of Patients 2 and 3 indicated that CH ( > 100%) is not caused solely by the severity of stenosis. Patient 2 exhibited a 3.8% chance of CH, whereas the corresponding estimates for Patients 1 and 3 were 0% and 0.001% (only one sample out of 100 000 samples), respectively. In Patient 2, who was assumed to have a possible missing ACoA, the variability of to values above 100% was prominent compared to Patient 3, implying that was significantly affected by this artery.

Patient conditions causing cerebral hyperperfusion.

To clarify the conditions under which CH occurs, we further investigated the characteristics of 3796 samples in Patient 2 and 1 sample in Patient 3 with > 100%. The left column in Fig 9 depicts the relationship between the preoperative PR of the MCA on the stenosis side (PR_MCA) and in each patient. Furthermore, the right column in Fig 9 illustrates the variation in with respect to the diameters of the ACoA and the posterior communicating artery (PCoA) on the stenosis side in each patient.

As indicated in the left column of Fig 9, exhibits an inverse relationship with PR_MCA. This is natural because , where denotes pressure recovery at the MCA on the stenosis side caused by the surgery. Even with the same , a smaller PR_MCA results in a larger . Fig 9C and 9E depict the results of Patients 2 and 3, respectively, wherein the PR_MCA is smaller than 20 mmHg s mL⁻¹ in most samples when exceeds 100%. However, we observed that a small PR_MCA did not always result in > 100%, as varied with respect to some factors. Samples with > 100% were associated not only with a small PR_MCA but also with small diameters of the ACoA and PCoA that form the collateral pathway to the artery on the stenosis side (Fig 9D and 9F). Particularly, an extremely small ACoA diameter (<1 mm) resulted in > 100%, regardless of the PCoA diameter (Fig 9D).

Fig 10 depicts the variation in the preoperative flow rate in the ACoA (left column) and PCoA (right column) with respect to the diameter. Note that the flow rate shown in Fig 10 varies both horizontally and vertically. As indicated by the relationship between ΔP and R_vQ in Eq (4), the flow rate in an artery is proportional to the pressure difference between the two ends and is inversely proportional to the fourth power of the diameter. The horizontal variation in flow rate shown in Fig 10 is attributed to the diameter variation of the communicating artery within the uncertainty range. On the contrary, the vertical variation is caused by variations in the pressure difference between the ends (i.e., the pressure difference between arteries on the normal and stenosis sides) resulting from uncertainties in the diameter of other arteries, stenosis severity, and flow measurements. As seen from the large vertical variations, the flow rate of the communicating artery is strongly influenced not only by the diameter uncertainty of this artery but also by other uncertainties.

As indicated in Fig 10C, the flow rate in the ACoA of Patient 2 is distributed up to 250 mL/min regardless of the diameter when it is ≥1 mm. However, when the diameter < 1 mm, the flow rate decreases rapidly, and samples with > 100% appear frequently near the upper end of the distribution. In other words, exceeds 100% when the collateral flow through the ACoA is limited owing to the small diameter despite the large pressure difference between arteries on the normal and stenosis sides. Conversely, no such condition was observed in Patients 1 and 3. The small amount of collateral flow in Patient 1 indicates that the pressure difference was small, and in Patient 3, the diameters of the ACoA and PCoA were sufficiently large.

Sensitivity of uncertain parameters.

To gain further insight into the factors associated with CH, we quantified the influence of each uncertain parameter on through SA. According to Saltelli’s algorithm, the number of samples required to compute the sensitivity indices was 370 000. Fig 11 depicts the first-order (S_n) and total (S_T,n) sensitivity indices. In most parameters, S_T,n is considerably larger than S_n, indicating a strong interaction between the parameters. In all patients, the diameters of the ACoA and PCoA on the stenosis side influenced the significantly. Additionally, the diameters of the anterior cerebral artery I (ACA I) and posterior cerebral artery I (PCA I) exhibited substantial sensitivity. As depicted in Fig 12, these arteries form collateral pathways to supply blood to the MCA on the stenosis side.

Furthermore, the severity of the stenosis affected the considerably. Although both R_v and SR are measures of stenosis severity, only R_v contributed to the variance of . The impact of R_v on was smaller in patients with highly severe stenosis. In Patient 1, the sensitivity of R_v was higher than that of collateral pathway diameters, whereas the opposite behavior was observed in Patients 2 and 3.