Patient-Specific Cerebral Blood Flow Simulation Based on Commonly Available Clinical Datasets

Cerebral hemodynamics play an important role in the development of cerebrovascular diseases. In this work, we propose a numerical framework for modeling patient-specific cerebral blood flow, using commonly available clinical datasets. Our hemodynamic model was developed using Simscape Fluids library in Simulink, based on a block diagram language. Medical imaging data obtained from computerized tomography angiography (CTA) in 59 patients with aneurysmal subarachnoid hemorrhage was used to extract arterial geometry parameters. Flow information obtained from transcranial Doppler (TCD) measurement was employed to calibrate input parameters of the hemodynamic model. The results show that the proposed numerical model can reproduce blood flow in the circle of Willis (CoW) per patient per measurement set. The resistance at the distal end of each terminal branch was the predominant parameter for the flow distribution in the CoW. The proposed model may be a promising tool for assessing cerebral hemodynamics in patients with cerebrovascular disease.

1 Introduction

Fifteen percent of the total cardiac output is distributed to the brain (Burton, 1965). The circle of Willis (CoW) is a circulatory anastomosis at the base of the brain and serves as the central hub that distributes blood flow to the brain. The morphology (Kayembe et al., 1984; Chuang et al., 2008; Krasny et al., 2014) and local hemodynamic flow patterns (Li et al., 2009; Can and Du, 2015; Tuenter et al., 2016) within the CoW are associated with cerebrovascular disease, such as carotid artery stenosis and intracranial aneurysms. Numerical modeling of the blood flow through the CoW can help understand the relationship between hemodynamics in the CoW and cerebrovascular disease. However, due to the complexity of the CoW-structure and the large anatomical variation of the CoW in the population, developing a reliable patients-specific numerical model for assessing cerebral hemodynamics is a great challenge.

Various numerical models have been proposed for investigating the hemodynamics in the CoW. The zero-dimensional (0D) model is the simplest model that can provide great insight into the compensation mechanisms of the CoW (Hillen et al., 1988). One-dimensional (1D) modeling has been widely employed for studying hemodynamic effects in the CoWdue to anatomical variations and occlusions (Alastruey et al., 2007; Huang et al., 2018). The more complex three-dimensional (3D) model became popular in recent decades, thanks to the fast development of computational fluid dynamics (CFD) (Alnæs et al., 2007). Detailed flow patterns in the intracranial aneurysms can be obtained with CFD models, which could be used to investigate the rupture risk of IAs (Asgharzadeh et al., 2020). Although the numerical models are based on the rigorous derivation of the fundamental governing equations of fluid dynamics, it was highly sensitive to the model input, particularly the inlet/outlet conditions (Venugopal et al., 2007; Balossino et al., 2009; Boccadifuoco et al., 2018).

In a recent review by Berg et al. (2019) about the hemodynamics modeling in intracranial aneurysms, it was concluded that the accuracy of the pre-and post-processing (e.g., geometry reconstruction, boundary condition setting, flow field analysis) were crucial in numerical simulation. Although many high-end numerical models have been proposed for cerebral hemodynamic studies with prospective potential for clinical practice, to the best of our knowledge, no existing numerical models have been universally accepted by clinicians to directly serve as clinical tool, due to the great uncertainty in patient-specific simulation. For patient-specific modeling, the required input not only includes patient-specific geometry, but more importantly appropriate physiological parameters, as well as proper boundary conditions per patient.

Many attempts have been made to develop patient-specific cerebral blood flow models. Zhang et al. (2016) proposed a 0D-1D coupling model and iteratively adjusted the peripheral cerebral resistance to match the measurement. They demonstrated that combining multiple sources of measurement by individual patients within the hemodynamic model can provide more comprehensive flow information. Lal et al. (2017) presented a more advanced approach for non-invasive estimation of the blood flow in the cerebral arteries, using an ensemble Kalman filter (EnKF) as an optimization tool. They utilized a 3-point clinical measurement of the transient blood flow to tune 21 input parameters. As such, they provided a feasible approach to assess blood flow at non-accessible locations in the cerebral arterial tree. Helthuis et al. (2020) introduced a patient-specific cerebral blood flow model with a structured tree and a simple auto-regulatory model at the distal boundary conditions. They validated their model with 20 healthy subjects and achieved a good agreement. Although these models can provide valuable hemodynamic information for patients, their clinical applicability was limited. The difficulties are two-fold: 1) This interdisciplinary research needs close collaboration between clinicians, biomedical engineers and computer science engineers. The clinicians commonly lack awareness on the role of hemodynamics but prefer traditional medical indices for clinical decision making (Singh et al., 2009). 2) A patient-specific model requires input that may be difficult to be obtained in clinical practice. Thus, many assumptions have to be made, which reduces the reliability of the hemodynamic assessment.

The objective of this study is to investigate cerebral hemodynamics in patients with cerebral vasospasm following aneurysmal subarachnoid hemorrhage (aSAH). A numerical model of patient-specific cerebral flow simulation is proposed, using clinical data obtained from a previously published study that investigating cerebral vasospasm in subarachnoid hemorrhage (TACTICS) study van der Harst et al. (2019). Our hemodynamic model is based on a lumped model in block diagram, which can be intuitively and easily used by clinicians. TCD measurements and common medical records are integrated to correct the input parameters of the model with a three-steps calibration procedure. The proposed model is applied to reproduce 66 cases of patient-specific hemodynamics at the time that vasospasm usually may occur. The validity of the proposed model is demonstrated by comparing the results with and without calibration procedures. The importance of the parameter calibration for the patient-specific simulation is discussed.

2 Methodology

The patients’ specific data were collected from the medical record of the prospective TACTICS study, including basic medical file, CTA images and TCD reports. The CTA and TCD were performed within 24 h of each other. The brachial mean, systolic and diastolic blood pressure (P_m, P_s and P_d), and the heart rate HR were obtained from the basic medical file. The geometry of the arterial network of the CoW, e.g., the diameter of arterial segments ϕ, was extracted from CTA images. The measured mean, systolic, and diastolic blood flow velocity (V_m, V_s and V_d) at specific segments were obtained from TCD reports. Figure 1 illustrates the steps to integrate the collected clinical data to assess patient-specific cerebral hemodynamics.

FIGURE 1

FIGURE 1. Flowchart of the numerical modelling framework.

The workflow consists of three steps:

(1) determine the peripheral resistance R_p at each terminal branch and the periodical blood flow rate at the heart Q_h with a linear steady model of cerebral arterial network;

(2) correct the peripheral compliance of the body artery network and implement the tentative simulation with the Simulink model;

(3) correct the compliance of the cerebral arterial segments by adjusting the Young’s modulus of the cerebral segments and implement the final simulation with the Simulink model.

The first step is to obtain the optimal blood flow distribution at the terminal branches of the cerebral arterial network, which can match the measured blood velocity by TCD. The second step provides the background pulse pressure at the entrance of the cerebral network that is in concordance with the medical record. The third step corrects the pulse waveform from the previous step to match the systolic and diastolic blood velocity by TCD measurement. The three-steps calibration procedure aims to tune the personalized model parameters to improve the fit of model predictions against measured data. The calibration is significantly important, in order to perform patient-specific simulation to reproduce the cerebral hemodynamics at the time of measurement, as will be shown in the result section.

2.1 Physiological Data

From 59 adult patients, diagnosed with aSAH within 4 days after onset, 66 sets of data were obtained. The CTA and TCD measurements have been further detailed in Shen et al. (2020).

In each CTA, the diameters of 17 cerebral segments were measured by semi-automated carotid lumen segmentation. The diameters of begin and end points of each segment were measured, while the simulation was based on an uniform diameter per segment which was calculated based on equivalent circular truncated cone volume. The thickness of arterial wall was assumed to be 25% of its radius. This was not critical, as wall thickness in conjunction with wall elasticity determine the compliance of each segment, while wall elasticity per patient per segment in the model calibration procedure will be adapted.

In each TCD, flow velocities at 11 locations on 9 cerebral segments were measured. Although patients with aSAH are prone to experience cerebral vasospasm, the measured velocity by TCD was thought to reflect the hemodynamics in the CoW at the time of CTA, since the delay between the examinations was short. Thus, patient-specific cerebral hemodynamics under vasospasm condition could be assessed by integrating the TCD and CTA measurements in the hemodynamic model.

For other required but unknown/unavailable physiological properties, such as elastic modulus of the segments, peripheral resistance, and compliance of the terminal branches, we used the same values as that in Alastruey et al. (2007), but would make corrections for each case as described in the following sections.

2.2 Simulink Model

Our hemodynamic model was developed based on Simulink, using Simscape Fluids library, which can provide an intuitive way to model the blood circulation system as a hydraulic system. The model is based on the block diagram environment that allows to easily customize the connectivity of the arterial network, as shown in Figure 2A.In this model, each block consists of a piece of governing equations that represents a simplified hydraulic component. The blocks are linked with connectors that form a complex hydraulic system. The system is interpreted as a set of differential or algebraic equations that are solved through symbolic formula manipulation.

FIGURE 2

FIGURE 2. Simulink model of arterial network with 33 segments.

In the present model, the arterial system originates from the heart (the red block) and ends at the green blocks, which depict the terminal branches. Each segment is represented with an individual blue block with the numbering and name referring to Table 1 in Alastruey et al. (2007), while the segments in the CoW are highlighted in yellow. The polylines represents the interconnection between the segments, which passes the pressure and flow rate data in the network.

TABLE 1

TABLE 1. ICC of various celebrate segments under various simulations.

A simple heart model directly enforces a periodical inlet flow rate using the “Hydraulic Flow Rate Source” block, which describes a cardiac cycle with a half-sinusoidal flow rate wave (systole) and zero in the rest of the period (diastole) that is expressed as:

$Q_h\left(t\right)=\left\{\begin{array}{cc}\frac{\pi }{2{\alpha }_T}{Q}_{h,m}\sin \left(\frac{\pi t}{{\alpha }_T{T}_h}\right),\hfill & \text{if}t<t_s\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}\right.(1)$

where α_T is the ratio of systole duration over the duration of heartbeat, α_T = 0.3 is chosen in the study; T_h is the duration of heartbeat that is obtained from the medical report; Q_h,m is the volumetric flux of the heart, as estimated in Section 2.3.

The arterial segments are described with “Segmented Pipeline” block, which accounts its resistance, fluid inertia, and wall compliance, as shown in Figure 2B. In this block, the deformation of the segmental wall is quantified with a static pressure-diameter coefficient K_p, which establishes the relationship between pressure and internal diameter of the arterial segment at steady-state condition. K_p is determined with the common properties of each arterial segment as:

$K_p=\frac{{\varphi }^2}{2Eh}\left(1-{\sigma }^2\right),(2)$

where ϕ is the internal diameter of the artery at the reference pressure, E is the Young’s modulus of the artery, h is the thickness of the artery; σ is the Poisson’s ratio, σ = 0.5. The pressure loss of the segment due to the friction is computed with the flow regime-dependable friction factor, which can take into account the possible flow separation due to the cerebral vasospasm. The friction factor in turbulent regime is determined with the Haaland approximation. In addition, the sensor blocks are used to monitor the pressure and flow rate at both ends of each segment.

A three-element Windkessel model (WK3) is embedded in the terminal branch to minimize the artificial reflection, see Figure 2C. A WK3 consists of two two resistances and a compliance. The “Linear Hydraulic Resistance” block is used to represent the resistance and The “Constant Volume Hydraulic Chamber” block is used to represent the peripheral compliance. Note that describing the peripheral compliance C in “Constant Volume Hydraulic Chamber” block is not straightforward, but using an equilibrium chamber length, which is proportional to the peripheral compliance, described as:

$l_{eq}=\frac{2C}{\pi \varphi K_p}.(3)$

The peripheral compliance and two resistances in WK3 are customized per terminal branch per patient.

From the description above, we can distinguish the present model from other lumped models. The present model directly uses the hydraulic description (pipeline network) instead of the common hydraulic-electrical analogue (electrical network). It is an object-oriented model, so the interpretation of the model layout is straightforward and intuitive. We can easily build patient-specific artery network by modifying the geometrical parameters and the connectivity of the network that makes this model suitable for investigating the cerebral hemodynamics with considerable CoW variation.

The Simulink model uses an arterial network with 33 segments (Alastruey et al., 2007) as the basis. For each patient-specific case, the diameters of the cerebral segments are altered by the measurement from the CTA images. The resistance and the compliance in each WK3 and the wall compliance of each segment are determined by calibration procedure using TCD measurement, as described in Section 2.3, 2.4.

2.3 Peripheral Resistance Correction

Analytical model by Hillen et al. (1988) has shown that the peripheral resistances dominate the flux distribution in the efferent segments and strongly influence the flux distribution in the afferent vessels. Hence, the peripheral resistance should be adapted to reflect the real flow situation at the time of the measurement for the patient-specific simulation. In the present study, a simple linear steady model was used to estimate the peripheral resistance, as shown in Figure 3.

FIGURE 3

FIGURE 3. Linear steady model of CoW arterial network.

The model only consisted of the CoW-segments, including three afferent arteries (left-right ICA and BA), six efferent arteries (M1L, M1R, A2L, A2R, P2L and P2R) and their interconnection. Several assumptions were made in this model: 1) The pressure drop between inlet and outlet are assumed to be equal as the difference between the mean blood pressure and the pressure at the entrance of the venous system (ΔP = P_m − P_v, with P_v = 5 mmHg); 2) the flow is fully developed (steady and laminar), the resistance of the arterial segments R_a follows the Hagen–Poiseuille law, R_a = 128η/πϕ⁴ with η blood viscosity; 3) the peripheral resistance R_p is considered a variable resistance, which may be adapted to alter the blood flow distribution in the CoW. Mathematically, the model is equivalent to an electric circuit with constant resistance and voltage source, which can be solved with node voltage analysis. Given specific peripheral resistance at the efferent arteries, the flow rate in each segment can be numerically obtained. TCD generally measures blood flow velocity of the target segments, while the exact location is unknown. Because the diameter varies along the vessel, the flow rate of the target segment cannot be directly determined. We measured the diameters at the start point ϕ_s and end point ϕ_e of the segment on CTA, thus the flow rate of the segment Q_m should fall within the interval $\left[0.25\pi {\varphi }_e^2{V}_m,0.25\pi {\varphi }_s^2{V}_m\right]$, where V_m is the measured mean flow velocity. The peripheral resistance can be obtained by searching for the minimum of the norm of the relative error of the flow rates,

$f\left({\mathit{R}}_p\right)=\mathrm{m}\mathrm{i}\mathrm{n}\sqrt{\sum _{i=1}^N{\left(\frac{Q_n^i\left({\mathit{R}}_p\right)-Q_m^i}{Q_m^i}\right)}^2},(4)$

where R_p is a vector of the peripheral resistance of the six outlet arteries, ${\mathit{R}}_p\in {\mathbb{R}}^6$; $Q_n^i$ is the flow rate obtained by the linear model; N is the number of TCD measurement, N = 9. Due to technical limitations of TCD, e.g., inadequate acoustical temporal bone window, it is possible that there is missing data on part of the locations for some cases. In these cases, N equals to the number of locations with available data. Eq. 4 essentially describes an optimization problem about the cerebral blood flow distribution, based on the measured velocity.

Summing up the flow rate of the three afferent arteries, we can obtain the overall blood flow rate of the CoW Q_CoW. Assuming 15% of the cardiac output perfuses the head (with 12% in CoW), 5% supplies each arm and 75% supplies to the rest of the body through the thoracic aorta, the mean inlet flow rate of the heart is determined as:

$Q_{h,m}=\frac{Q_{CoW,m}}{0.12}.(5)$

Consequently, the peripheral resistance of the rest terminal branches can be determined by the ratio of the pressure drop and the flow rate. For example, the peripheral resistance of thoracic aorta is calculated by:

$R_{body,m}=\frac{\mathrm{\Delta }P}{0.75Q_{h,m}},(6)$

To minimize the reflection at the terminal branch, R1 in WK3 is set as the characteristic impedance of the terminal branch,

$R1=\frac{\rho c_0}{A_0}(7)$

with c₀ is the speed of pulse wave propagation; A₀ is the cross-section area of the segment. Then R2 is calculated as:

$R2=R_p-R1.(8)$

By using Eq. 4, we observed that for some cases, R_p may be smaller than R1, this may be caused by the uncertainty of TCD and CTA measurement, as well as the possible unexpected reflection due to the cerebral vasospasm. To avoid a negative resistance, we altered the resistance for those cases with:

$R1=0.8R_p\mathrm{a}\mathrm{n}\mathrm{d}R2=0.2R_p.(9)$

With optimal peripheral resistance and proper flow rate at the heart, the obtained mean flow rate and mean blood pressure matches the measurement very well.

2.4 Compliance Correction

In a cerebral arterial network system, the mean flux distribution is mainly governed by the peripheral resistance, whereas the oscillation pattern is dominated by wall compliance and peripheral compliance. The wall compliance can be adjusted by varying Young’s modulus of the segment. It is more difficult to determine the peripheral compliance because of the complex wave interaction. In the present study, the objective of the compliance correction is to match the systolic and diastolic blood pressures of the brachial artery and the systolic and diastolic cerebral blood flow velocity. Since Alastruey et al. (2007)’s model is employed as the reference model, the pressure pulse calibration can be achieved by scaling the reference model. Given a volumetric flux of the heart in the model, the amplitude of the pressure pulse is inversely proportional to the peripheral compliance. First, we computed the ratio of mean blood pressure and flow rate between patient and reference model with:

${\alpha }_1=\frac{\mathrm{\Delta }P{Q}_{ref,m}}{\mathrm{\Delta }{P}_{ref}{Q}_{h,m}};(10)$

where ΔP is the pressure difference between P_s and P_d, Q_ref,m and ΔP_ref are the reference mean flow rate and pressure difference by the reference model, and Q_ref,m = 92.6 ml s⁻¹ and ΔP_ref = 55 mmHg that is obtained from the Simulink model. The Young’s modulus and the peripheral compliance of the segments can be altered by equation:

$E_1={\alpha }_1{E}_{ref}\mathrm{a}\mathrm{n}\mathrm{d}{C}_{p,1}={\alpha }_1{C}_{p,ref}(11)$

where E_ref and C_p,ref represent the Young’s modulus and the peripheral compliance in the reference model. With the first calibrated peripheral compliance and Young’s modulus, the tentative simulation of the Simulink model is implemented. This simulation can result in a good approximation of the pulse in the aorta arteries, but may not derive a satisfactory waveform in the cerebral arteries, i.e., poor agreement with V_s and V_d in CoW, which requires a second calibration. In intra-aortic hemodynamics simulation, multiple TCD measurement can be approached by iteratively adjusting the peripheral compliance with WK3, as demonstrated by Pant et al. (2014), Bonfanti et al. (2019), while wall compliance is kept constant during the iteration. However, the intracranial arteries are stiffer than extracranial arteries, the waveform in the ICA and BA are therefor not sensitive to the peripheral compliance of the terminal branch of the CoW. Hence, the second calibration is carried out by only adjusting the Young’s modulus of the cerebral arteries by an uniform factor α₂, which is estimated by:

${\alpha }_2=\frac{\mathrm{\Delta }{Q}_{ICAL,m}+\mathrm{\Delta }{Q}_{ICAR,m}+\mathrm{\Delta }{Q}_{BA,m}}{\mathrm{\Delta }{Q}_{ICAL,n}+\mathrm{\Delta }{Q}_{ICAR,n}+\mathrm{\Delta }{Q}_{BA,n}},(12)$

where ΔQ is the flow rate difference between Q_s and Q_d, with subscripts m and n refers to the measurement and numerical results respectively. The Young’s modulus and the peripheral compliance of the cerebral arterial segments can be altered by equation:

$E_2={\alpha }_2{E}_1\mathrm{a}\mathrm{n}\mathrm{d}{C}_{p,2}={\alpha }_2{C}_{p,1}.(13)$

Although individual calibration per arterial segment can be done with Kalman filter Pant et al. (2014), it required many iterations that made the computation inefficient. We applied an uniform factor because the flow velocity results by the tentative simulation in all cerebral segments was observed having similar trend comparing to the TCD measurement. Satisfactory results can be obtained without any complicated iteration.

2.5 Model Validation Metrics

Model validation is quantified with a number of statistical measures. The intraclass correlation coefficient (ICC) was used to assess the agreement between the results obtained by numerical simulation and the TCD measurement. Following the instruction by McGraw and Wong (1996), we selected the type of ICC(A,1) in the study, which assessed the degree of absolute agreement among simulation and measurement. The Taylor diagram (Taylor, 2001) is used to evaluate the validity of the proposed calibrations by comparing the numerical results to that by conventional boundary condition. Taylor diagram is a mathematical diagram to quantify the degree of correspondence between the modelled and observed behavior with three statistics: the Pearson correlation coefficient ρ, the centred-root-mean-square error E′, and the standard deviation σ, which was intensively used in climate study (Kärnä and Baptista, 2016). Based on the definition of the three statistics, they have the following relation:

$E^{\prime 2}={\sigma }_m^2+{\sigma }_n^2-2{\sigma }_m{\sigma }_n\rho .(14)$

where n demotes the variable by the numerical model; m demotes the measured counterpart. Making use of the law of cosines, Eq. 14 can be interpreted as geometrical relationship of a triangle, with E′, σ_m, and σ_n are the length of the sides of the triangle, and ρ is the angle between sides σ_m, and σ_n. In order to compare the data sets at multiple arterial segments in the one diagram, the normalized Taylor diagram was employed by scaling Eq. 14 with σ_m. Thus, the measurement perfectly locates at $\left(R=1,\theta =0\right)$. The numerical results are displayed in polar coordinate with radial coordinate r = σ_n/σ_m and angle θ = arccos ρ.

3 Results

3.1 Validation of Simulink Model

The Simulink model was firstly validated by comparing the results with that of 1D model by Alastruey et al. (2007), as presented in Figure 4.

FIGURE 4

FIGURE 4. Comparison of flow velocity at brachiocephalic bifurcation (top) and in CoW (bottom) by Simulink model (solid lines) and 1D model (dash lines).

The same parameters as that in Alastruey et al. (2007) were applied in the Simulink model. The flow velocity at the brachiocephalic (BrA), right subclavian (SCA) and right common carotid (CCA) showed very good agreement with that by 1D model. Some discrepancy can be found for the flow velocity in CoW. There may be two main causes: 1) wave propagation in the CoW was complex and the waveform was vulnerable to wave interference; 2) the Simulink model is a lumped model, which can be regarded as the first order discretisations of 1D systems (Milišić and Quarteroni, 2004), omitting the (nonlinear) convective acceleration term. Nevertheless, the Simulink model can capture the main waveform similar to that by 1D model, and the waveform can be further calibrated with the TCD measurement. Hence, the Simulink model can be used to investigate the patient-specific cerebral hemodynamics.

3.2 Patient-specific Simulation

In this study, 66 sets of patient-specific simulation were implemented. Figure 5 presents the log-log scatter plot of the optimized peripheral resistance obtained by Eq. 4 against the diameter of afferent segments of CoW, with the black line showing the log-log regression.

FIGURE 5

FIGURE 5. The log-log scatter plot of the optimized peripheral resistance against the diameter of the afferent segments of CoW.

It can be seen that the optimized peripheral resistance showed a trench of inverse-proportional to the diameters (R_p ∝ ϕ^{−0.995 3}), but it was highly scattering. Alastruey et al. (2007) assumed that the peripheral resistance was inverse-proportional to the cross-sectional areas (R_p ∝ ϕ⁻²). This suggested that this simple assumption was not valid for the patient-specific simulation as the peripheral resistance dominated the flow distribution in the CoW. In addition, the resistance of A2 segment was larger than that of P2 with similar diameter, implying that the resistance in the anterior cerebral circulation was larger than that in the posterior cerebral circulation, that was not surprise as the blood velocity in ICA was commonly larger than that in BA. Figure 6 showed the comparison of blood pressure after the first compliance correction, i.e., Eq. 13.

FIGURE 6

FIGURE 6. Comparison of diastolic, systolic and mean blood pressure at brachial artery by simulation and measurement, the cases were sorted in order of low to high measured mean blood pressure.

The excellent agreement was obtained between the numerical and measured P_m, and the very good agreement was observed for P_s and P_d. Since these three parameters were commonly used to sketch the blood pressure waveform, we considered that the obtained pressure waveform by the Simulink model at the entrance of CoW represented the hemodynamic condition at the time of TCD measurement.

Figure 7 showed the comparison of flow rate in the afferent and efferent arterial segments between the measurement and the results by three simulations. The measured flow rate is indirectly obtained by solving the optimization problem of peripheral resistance, i.e., Eq. 4. In the “non-calibrated” simulation, we calculated the peripheral resistance and compliance in CoW of each patient based on the assumption that the blood flux of each part in the body followed a given flow distribution, and the outflow of the CoW was proportional to their initial cross-sectional areas. The agreement between the measurement and the “non-calibrated” simulation was poor, indicating that the simple assumption was not valid for the patient-specific study. The customized parameter must be applied per patient to represent the patient’s hemodynamics at the time of measurement. In the “first calibrated” simulation, the calibrated peripheral resistance and compliance in the body segments was applied in the simulation. The mean flow rate in CoW by simulation agreed well with the measurement, and the results of the efferent segments outperformed that of the afferent segments. The simulation slightly underestimated the diastolic flow rate while overestimated the systolic flow rate, implying that the amplitude of the waveform was lower in simulation. In the “second calibrated” simulation, the compliance of the cerebrate segments were calibrated based on the first calibrated results. The calibration did not alter the mean flow rate in the previous simulation, while the agreement of the diastolic and systolic flow rate was improved. Since TCD-measurement has been used to calibrate the model as well as to validate the model output, it is not surprising that the flow rate by the “calibrated” simulation yields the best agreement. Based on the calibrated model, we can reproduce the cerebral hemodynamic condition for the patient at the time of measurement and provide more comprehensive cerebral hemodynamics for the clinician.

FIGURE 7

FIGURE 7. Comparison of diastolic (1st row), mean (2nd row) and systolic (3rd row) flow rate of various cerebral segments by measurement and various simulations.

The statistic of the results in Figure 7 were summarized in Table 1. The agreement of the second calibrated results was good (0.75 < ICC < 0.90) to excellent (ICC > 0.90), except the BA. The current arterial network did not take into account the minor branches on the CoW that would introduce more error on the BA than that other segments as the BA is a relatively short segment but has the most minor branches. Q_m always presented better agreement comparing to Q_d and Q_s, this was because Q_m was strongly dominated by the peripheral resistance (the peripheral resistance was generally two order larger than the segmental resistance of the CoW). Q_d and Q_s was not only influenced by the peripheral compliance but also the compliance of the segment in CoW. Nevertheless, the agreement was satisfactory, particularly considering that only twice calibration were applied, which was significantly less iteration than using EnKF.

Figure 8 presented the normalized Taylor diagram of various segments that gave an impression of importance of the calibration procedures. The results by the three simulations were plotted in different colour, and three marker symbols were used to distinguish the results of Q_d, Q_m and Q_d. In the Taylor diagram, good results were indicated by having relatively high correlation coefficient and low RMSE (as close to the measured point as possible). The agreement of the “un-calibrated” simulation was poor as it has very low correlation coefficient and very high RMSE, while the first and second calibrated results are much closer to the measurement.

FIGURE 8

FIGURE 8. Normalized Taylor diagram of various celebrate segments by three simulations.

The improvement by the calibration procedures can be seen from Figure 9, the normalized Taylor diagram of flow rate at 9 cerebral locations per patient. The data points cluster of the second calibrated results located at around the $\sim 0.25$ RMSE circle, while the first calibrated results shifted further from the measurement, indicating that the second calibration can describe the waveform pattern better than the first calibration.

FIGURE 9

FIGURE 9. Normalized Taylor diagram of each patient by three simulations.

4 Discussion

When mentioning patient-specific cerebral simulation, many studies refer to the patient-specific geometry of arteries but do not take patient-specific boundary conditions into account. Our results have demonstrated the importance of boundary conditions and wall properties. Comparison by three simulations confirmed the findings of Hillen et al. (1988) that the peripheral resistance plays a predominant role in the flow distribution of the efferent arteries of the CoW. The blood flow waveform in segments in the CoW is not significantly influenced by the peripheral compliance, but would be adapted to any change of the segmental elasticity in the cerebral arterial network. To reproduce the hemodynamic environment in cerebral vasospasm, including the pulse waveform, patient-specific boundary conditions and wall properties must be corrected by integrating flow measurement. On the other hand, the model calibration relies on flow measurement on cerebral segments, which may be missing in retrospective studies, sometimes difficult to be obtained for some patients. The uncalibrated model cannot derive a reliable absolute blood flow rate, as shown in Figure 8 that the correlation coefficient is too low. Therefore, the absolute quantitative comparison should be avoided in the statistical analysis. Alternatively, using the flow rate ratio between the segments may alleviates some of this problem and yield more meaningful results, as it can be seen from Figure 9 that although the universal boundary condition leads to the worst results, there is not any significant difference of the correlation coefficient by three simulations. This idea is similar to the study by Lindegaard et al. (1988) that used the ratio of blood velocity in MCA and ICA by TCD to diagnose vasospasm.

The present model with personalized calibration has potential clinical applications. The model integrates the patient’s morphology information and available flow measurement to provide the comprehensive cerebral hemodynamics, which may help the clinician to establish the correlations between hemodynamic indices and cerebrovascular disease. Since cerebral vasospasm due to the aSAH can result in delayed cerebral ischemia, the model can provide insight in cerebral hemodynamic development during the course of vasospasm. It may also be easily employed to predict the improvement of the blood supply after intracranial bypass operation or treatment of the carotid artery stenosis, thanks to the intuitive Simulink model. In addition, the model can also provide the boundary conditions for a patient-specific CFD model to investigate the detailed local flow in the CoW, such as the flow pattern in intracranial aneurysms. We consider the present work as the first step toward a numerical computational framework serves as a clinical decision-making tool for hemodynamic related cerebrovascular disease. The major advantage of the present model is that the model was developed with a block diagram, which is also user-friendly to the clinician. Also, the required patient’s data are provided by CTA and TCD, which is commonly accessible clinical data. As such, it is easy to build up a patient-specific simulation. The calibration procedure is implemented to correct the model input with a moderate computational burden.

The present study also has limitations. The most significant limiting factor of the model is the use of indirectly obtained flow rate for calibration and validation. TCD measures flow velocity at unknown location of the target segment, the flow rate is numerically determined with a linear model of the CoW network. On one hand, the linear model obeys the law of conservation of mass and integrates the TCD velocity and dimension of segments by CTA, which estimates the flow distribution in the CoW network with minimal mass conservation error. On the other hand, the linear model uses simplified CoW network omitting the less important branch, which introduces uncertainty into the model. Finally, while we have identified the importance of the calibration procedure for patient-specific simulation, additional patient information (arterial blood pressure, TCD records) is required for calibration, which may limit its application. Nevertheless, we assume that the required data is commonly available for patients with cerebrovascular disease.

5 Conclusion

This paper presents a three-steps approach for implementation of a patient-specific cerebral blood flow simulation based on commonly available clinical datasets. The new approach integrates TCD and CTA data to correct the input parameters of an object-oriented hemodynamic model. Simulation with 66 sets of patient data was carried out to evaluate the validity of the model. The results demonstrated the importance of the peripheral resistance in the patient-specific cerebral blood flow simulation, and the satisfactory pulse waveforms can be reconstructed by the personalized model calibration.

The proposed method can reproduce the blood flow in the CoW where the blood flow in afferent and efferent segments was in good agreement with the TCD measurement. The model is potentially a promising tool for developing clinical understanding of cerebrovascular disease. The obtained results will be further investigated by clinical researcher to find the potential relationship between the vasospasm and the configuration of the CoW.

Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Ethics Statement

The studies involving human participants were reviewed and approved by the Institutional Ethical Review Board of the University Medical Center Groningen. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

Author Contributions

YS designed the study, collected and analyzed the data, partially drafted the manuscript. YW designed the study, developed the numerical model and performed the simulation, partially drafted the manuscript. RB assisted the data analysis, critically revised the manuscript. MU provided the patients’ data, technically assisted the collection of the data, and critically revised the manuscript. JD as corresponding researcher, is responsible for the quality of the data, supervised the study, and critically revised the manuscript.

Funding

YS receives financial support from China Scholarship Council (CSC, File No. 201706320024).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Alastruey, J., Parker, K. H., Peiró, J., Byrd, S. M., and Sherwin, S. J. (2007). Modelling the circle of willis to Assess the Effects of Anatomical Variations and Occlusions on Cerebral Flows. J. Biomech. 40, 1794–1805. doi:10.1016/j.jbiomech.2006.07.008

CrossRef Full Text | Google Scholar

Alnæs, M. S., Isaksen, J., Mardal, K.-A., Romner, B., Morgan, M. K., and Ingebrigtsen, T. (2007). Computation of Hemodynamics in the circle of willis. Stroke 38, 2500–2505. doi:10.1161/STROKEAHA.107.482471

PubMed Abstract | CrossRef Full Text | Google Scholar

Asgharzadeh, H., Shahmohammadi, A., Varble, N., Levy, E. I., Meng, H., and Borazjani, I. (2020). A Simple Flow Classification Parameter Can Discriminate Rupture Status in Intracranial Aneurysms. Neurosurg. 87, E557–E564. doi:10.1093/neuros/nyaa189

PubMed Abstract | CrossRef Full Text | Google Scholar

Balossino, R., Pennati, G., Migliavacca, F., Formaggia, L., Veneziani, A., Tuveri, M., et al. (2009). Computational Models to Predict Stenosis Growth in Carotid Arteries: Which Is the Role of Boundary Conditions? Comput. Methods Biomech. Biomed. Eng. 12, 113–123. doi:10.1080/10255840802356691

CrossRef Full Text | Google Scholar

Berg, P., Saalfeld, S., Voß, S., Beuing, O., and Janiga, G. (2019). A Review on the Reliability of Hemodynamic Modeling in Intracranial Aneurysms: Why Computational Fluid Dynamics Alone Cannot Solve the Equation. Neurosurg. focus 47, E15. doi:10.3171/2019.4.focus19181

CrossRef Full Text | Google Scholar

Boccadifuoco, A., Mariotti, A., Celi, S., Martini, N., and Salvetti, M. V. (2018). Impact of Uncertainties in Outflow Boundary Conditions on the Predictions of Hemodynamic Simulations of Ascending Thoracic Aortic Aneurysms. Comput. Fluids 165, 96–115. doi:10.1016/j.compfluid.2018.01.012

CrossRef Full Text | Google Scholar

Bonfanti, M., Franzetti, G., Maritati, G., Homer-Vanniasinkam, S., Balabani, S., and Díaz-Zuccarini, V. (2019). Patient-specific Haemodynamic Simulations of Complex Aortic Dissections Informed by Commonly Available Clinical Datasets. Med. Eng. Phys. 71, 45–55. doi:10.1016/j.medengphy.2019.06.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Burton, A. C. (1965). Physiology and Biophysics of the Circulation. Chicago: Year Book Medical Publishers.

Google Scholar

Can, A., and Du, R. (2015). Association of Hemodynamic Factors with Intracranial Aneurysm Formation and Rupture. Neurosurgery 78, 510–520. doi:10.1227/NEU.0000000000001083

PubMed Abstract | CrossRef Full Text | Google Scholar

Chuang, Y.-M., Liu, C.-Y., Pan, P.-J., and Lin, C.-P. (2008). Posterior Communicating Artery Hypoplasia as a Risk Factor for Acute Ischemic Stroke in the Absence of Carotid Artery Occlusion. J. Clin. Neurosci. 15, 1376–1381. doi:10.1016/j.jocn.2008.02.002

CrossRef Full Text | Google Scholar

Helthuis, J. H. G., van Doormaal, T. P. C., Amin-Hanjani, S., Du, X., Charbel, F. T., Hillen, B., et al. (2020). A Patient-specific Cerebral Blood Flow Model. J. Biomech. 98, 109445. doi:10.1016/j.jbiomech.2019.109445

CrossRef Full Text | Google Scholar

Hillen, B., Drinkenburg, B. A. H., Hoogstraten, H. W., and Post, L. (1988). Analysis of Flow and Vascular Resistance in a Model of the Cricle of willis. J. Biomech. 21, 807–814. doi:10.1016/0021-9290(88)90013-9

CrossRef Full Text | Google Scholar

Huang, G. P., Yu, H., Yang, Z., Schwieterman, R., and Ludwig, B. (2018). 1d Simulation of Blood Flow Characteristics in the circle of willis Using Thinks. Comput. Methods Biomech. Biomed. Eng. 21, 389–397. doi:10.1080/10255842.2018.1468439

PubMed Abstract | CrossRef Full Text | Google Scholar

Kärnä, T., and Baptista, A. M. (2016). Evaluation of a Long-Term Hindcast Simulation for the columbia River Estuary. Ocean Model. 99, 1–14. doi:10.1016/j.ocemod.2015.12.007

CrossRef Full Text | Google Scholar

Kayembe, K. N., Sasahara, M., and Hazama, F. (1984). Cerebral Aneurysms and Variations in the circle of willis. Stroke 15, 846–850. doi:10.1161/01.STR.15.5.846

PubMed Abstract | CrossRef Full Text | Google Scholar

Krasny, A., Nensa, F., Sandalcioglu, I. E., Göricke, S. L., Wanke, I., Gramsch, C., et al. (2014). Association of Aneurysms and Variation of the A1 Segment. J. Neurointervent Surg. 6, 178–183. doi:10.1136/neurintsurg-2013-010669

CrossRef Full Text | Google Scholar

Lal, R., Nicoud, F., Bars, E. L., Deverdun, J., Molino, F., Costalat, V., et al. (2017). Non Invasive Blood Flow Features Estimation in Cerebral Arteries from Uncertain Medical Data. Ann. Biomed. Eng. 45, 2574–2591. doi:10.1007/s10439-017-1904-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z.-Y., Taviani, V., Tang, T., Sadat, U., Young, V., Patterson, A., et al. (2009). The Mechanical Triggers of Plaque Rupture: Shear Stressvspressure Gradient. BJR 82, S39–S45. doi:10.1259/bjr/15036781

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindegaard, K.-F., Nornes, H., Bakke, S. J., Sorteberg, W., and Nakstad, P. (1988). “Cerebral Vasospasm after Subarachnoid Haemorrhage Investigated by Means of Transcranial Doppler Ultrasound,” in Proceedings of the 8th European Congress of Neurosurgery Barcelona, September 6–11, 1987 (Barcelona: Springer), 81–84. doi:10.1007/978-3-7091-8975-7_16

CrossRef Full Text | Google Scholar

McGraw, K. O., and Wong, S. P. (1996). Forming Inferences about Some Intraclass Correlation Coefficients. Psychol. Methods 1, 30–46. doi:10.1037/1082-989x.1.1.30

CrossRef Full Text | Google Scholar

Milišić, V., and Quarteroni, A. (2004). Analysis of Lumped Parameter Models for Blood Flow Simulations and Their Relation with 1d Models. ESAIM Math. Model. Numer. Anal. 38, 613–632. doi:10.1051/m2an:2004036

CrossRef Full Text | Google Scholar

Pant, S., Fabrèges, B., Gerbeau, J.-F., and Vignon-Clementel, I. E. (2014). A Methodological Paradigm for Patient-specific Multi-Scale Cfd Simulations: from Clinical Measurements to Parameter Estimates for Individual Analysis. Int. J. Numer. Meth. Biomed. Engng. 30, 1614–1648. doi:10.1002/cnm.2692

CrossRef Full Text | Google Scholar

Shen, Y., Wei, Y., Bokkers, R. P. H., Uyttenboogaart, M., and van Dijk, J. M. C. (2020). Study Protocol of Validating a Numerical Model to Assess the Blood Flow in the circle of willis. BMJ open 10, e036404. doi:10.1136/bmjopen-2019-036404

CrossRef Full Text | Google Scholar

Singh, P. K., Marzo, A., Coley, S. C., Berti, G., Bijlenga, P., Lawford, P. V., et al. (2009). The Role of Computational Fluid Dynamics in the Management of Unruptured Intracranial Aneurysms: a Clinicians' View. Comput. Intell. Neurosci. 2009, 760364. doi:10.1155/2009/760364

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, K. E. (2001). Summarizing Multiple Aspects of Model Performance in a Single Diagram. J. Geophys. Res. 106, 7183–7192. doi:10.1029/2000JD900719

CrossRef Full Text | Google Scholar

Tuenter, A., Selwaness, M., Arias Lorza, A., Schuurbiers, J. C. H., Speelman, L., Cibis, M., et al. (2016). High Shear Stress Relates to Intraplaque Haemorrhage in Asymptomatic Carotid Plaques. Atherosclerosis 251, 348–354. doi:10.1016/j.atherosclerosis.2016.05.018

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Harst, J. J., Luijckx, G. R., Elting, J. W. J., Bokkers, R. P. H., van den Bergh, W. M., Eshghi, O. S., et al. (2019). Transcranial Doppler versus CT-angiography for Detection of Cerebral Vasospasm in Relation to Delayed Cerebral Ischemia after Aneurysmal Subarachnoid Hemorrhage: A Prospective Single-center Cohort Study: The Transcranial Doppler and CT-angiography for Investigating Cerebral Vasospasm in Subarachnoid Hemorrhage (TACTICS) Study. Crit. Care Explor. 1, e0001. doi:10.1097/CCE.0000000000000001

PubMed Abstract | CrossRef Full Text | Google Scholar

Venugopal, P., Valentino, D., Schmitt, H., Villablanca, J. P., Viñuela, F., and Duckwiler, G. (2007). Sensitivity of Patient-specific Numerical Simulation of Cerebal Aneurysm Hemodynamics to Inflow Boundary Conditions. J. Neurosurg. 106, 1051–1060. doi:10.3171/jns.2007.106.6.1051

CrossRef Full Text | Google Scholar

Zhang, H., Fujiwara, N., Kobayashi, M., Yamada, S., Liang, F., Takagi, S., et al. (2016). Development of a Numerical Method for Patient-specific Cerebral Circulation Using 1D-0D Simulation of the Entire Cardiovascular System with SPECT Data. Ann. Biomed. Eng. 44, 2351–2363. doi:10.1007/s10439-015-1544-8

PubMed Abstract | CrossRef Full Text | Google Scholar