Hemodynamic Analysis of Carotid Sinus with OpenFOAM and Machine Learning

Abstract

Computational fluid dynamics (CFD) shows promise in aiding clinical methods in the early detection of atherosclerosis when combined with currently popular machine learning algorithms. In this study, fluid-structure interaction (FSI) analysis of the carotid artery was performed by creating three-dimensional patient-specific pre-operation carotid artery models of four different patients which have vessel stenosis or aneurysms. As a result of numerical simulations, the average flow velocity and average pressure of the patients at 80 specific cross-sections were obtained. The simulation results of three patients’ pre-operation were used for learning in the machine learning algorithm. The training data consists of 80% of the numerical values, while the remaining 20% is used for testing. Then, the algorithm was asked to predict the flow velocity values at different cross-sections of the artery. The values obtained as a result of learning were compared with those obtained from numerical simulation. We found the results promising in terms of guiding the clinical decisions.

1. Introduction

Hemodynamics is the branch of science that explains the blood circulation in the vessels and the physical factors that affect this circulation. Plaque, which disrupts blood flow in the vessels, occurs depending on many parameters such as nutrition, age, gender, high cholesterol, smoking, diabetes, and hypertension. At the point where plaque forms, it causes narrowing of the arteries diameter and disrupts blood flow, which leads to many vital problems, including myocardial infarction, stroke, and acute ischemia of the organs. As a result of chronic inflammation processes of atherosclerosis, the arteries walls thicken, intraluminal spaces become stenotic, and arteries lose their flexibility and become more rigid. The relationship between atherosclerotic plaque formations and bifurcation and the blood flow in the arteries has been theoretically explained by many researchers [1, 2, 3]. Blood flow velocity has a major influence on understanding and analyzing the dynamics of blood circulation in the body. It provides valuable information about the health of the cardiovascular system and can assist in diagnosing conditions such as arterial stenosis and atherosclerosis.

The carotid arteries examined in this study are the vessels on either side of the human neck, and intracranial parts form the circle of Willis polygon that carries blood to the brain. The carotid artery consists of two arteries branching from the common carotid artery (CCA), called the internal carotid artery (ICA) and the external carotid artery (ECA). The bifurcation point of the CCA and the carotid sinus region on the ICA are more prone to plaque formation due to the low velocity and low shear stresses that occur in that region. As a result of plaque build-up, due to contraction blood flow in the vessel is reduced, which can lead to a stroke. The carotid sinus is a small bulb area beginning of the ICA, near where the CCA divides into two branches.

The behavior of blood flow due to plaque formation in the carotid arteries has been investigated by computational fluid dynamics methods in the literature. Because CFD analysis has the potential to be particularly valuable in the diagnostic phase for asymptomatic patients. In the past, analyses with two-dimensional vessel models were replaced by analyses on three-dimensional vessel models. Nowadays, numerical studies are performed on patient-specific vessel models obtained from computed tomography (CT) angiography images. As a result of numerical analysis, parameters such as flow velocity, wall shear stresses (WSS), time-averaged wall shear stress (TWSS), and oscillatory shear index (OSI) can be obtained. These values provide an understanding of the behavior of blood flow in the plaque area. In order to understand the effect of artery stiffness, numerical studies have been conducted considering both rigid and fluid-structure interaction (FSI) models of the carotid artery [4, 5, 6, 7, 8, 9, 10]. FSI takes a more realistic approach to simulate the problem, considering the effect of blood on the artery and deformation in the vessel walls. However, FSI requires modeling of the vessel walls as well as the inner lumen and therefore needs more computational cost. In fact, the rigid models slightly overestimate the WSS compared to FSI values [4]. Those evaluations are trying to give further information that could help guiding the clinical decisions.

Clinically, analysis of hemodynamics and detection of disease is performed using existing radiological imaging modalities such as Doppler ultrasonography, computed tomography angiography (CTA), magnetic resonance angiography (MRA), and digital subtraction angiography (DSA), however, numerical simulations can help to obtain more detailed information about blood flow in patients’ vessels.

Because it has been described in the literature that some of these imaging techniques overestimate and some underestimate the stenosis if the stenosis is >70% [11]. Performing such numerical simulations using a large number of real patient geometries and teaching the hemodynamics of these patients to machine learning algorithms can help to predict the hemodynamics of patients without angiography. Nowadays, the popularity of using machine learning algorithms has increased in the medical area as in every research area. The studies have been conducted in the detection of carotid artery diseases by machine learning techniques in the literature [12, 13, 14, 15]. Machine learning techniques can be utilized to predict blood flow velocity in the carotid artery, considering various factors such as the geometry of the artery, blood viscosity, arterial walls mechanical properties such as elasticity modulus, Poisson ratio, and density. By considering patient-specific geometries and pulsatile wave cycles for blood flow, machine learning algorithms can be trained on a dataset containing information about the blood flow velocity in the carotid artery, as well as associated factors such as age, gender, and presence of cardiovascular diseases. This approach can provide early diagnosis of diseases by making accurate predictions regarding the determination of blood flow velocity in the carotid artery.

This chapter/section gives an idea on how CFD studies and machine learning algorithms can help in the prediction of hemodynamic parameters, considering pre-operation conditions of patients in terms of guiding clinical decisions. CFD analyses were performed with OpenFOAM [16], and machine learning algorithms were created in Python [17].

2. Methodology

2.1 Problem formulation

In this study, blood flow behavior is assumed to be incompressible flow and governed by the continuity and momentum equations. Continuity and Momentum equations are given as;

$$\nabla .\mathbf{u}=0$$E1

$$\mathrm{\rho }\left(\frac{\partial \mathbf{u}}{\mathrm{\partial t}}+\mathbf{u}\mathbf{.}\nabla \mathbf{u}\right)=-\nabla \mathrm{P}+\mathrm{\mu }{\nabla }^2\mathbf{u}\mathbf{+}\mathbf{f}$$E2

The linear momentum equation is the governing equation for artery wall deformation and is expressed as [4],

$${\rho }_w\frac{{\partial }^{2}w}{\partial t^2}-\nabla \sigma ={\rho }_w\overrightarrow{b}$$E3

$$\sigma =2{\mu }_L\epsilon +{\lambda }_L\mathit{tr}\left(\epsilon \right)I$$E4

The first and second Lame parameters are given as ${\mu }_L$ and ${\lambda }_L$, respectively. Lame parameters are related with elasticity modulus and Poisson ratio.

$${\lambda }_L=\frac{\mathit{\upsilon E}}{\left(1+\upsilon \right)\left(2\upsilon -1\right)}$$E5

$${\mu }_L=\frac{E}{2\left(1+\upsilon \right)}$$E6

In general, blood is Non-Newtonian owing to shear rate- and time-dependent viscous properties. Blood viscosity varies as a function of shear rate and blood hematocrit. For Non-Newtonian flow behavior, various models such as Carreau-Yasuda, Power law, Carreau, and Casson can be used. The current study is using the Casson model, which considers shear thinning behavior and the yield stress of the blood. This model is used at low shear rates in small arteries and is suitable for a wide range of shear rates. The shear stress-shear rate relationship for blood can be given by the empirical relation known as the Casson equation [18, 19],

$$\sqrt{\mathrm{\tau }}=\sqrt{{\tau }_y}+\sqrt{{\mu }_0\dot{\gamma }}$$E7

Blood viscosity described as [20],

$$\mu =\frac{\tau }{\dot{\gamma }}$$E8

Wall shear stress is the frictional force per unit area applied by the blood on the surface of the arteries. In this study time averaged wall shear stresses during a pulsatile cycle are calculated as;

$$\mathit{\text{TWSS}}=\frac{1}{T}{\int }_0^T\left|{\overrightarrow{\tau }}_w\right|\mathit{dt}$$E9

OSI is a non-dimensional parameter that shows the deviation of the wall shear stress vector from the flow direction during one pulse cycle. OSI is calculated as,

$$\mathit{OSI}=0.5\times \left(1.0-\frac{\left|{\int }_0^T\overrightarrow{\mathit{WSS}}\mathit{dt}\right|}{{\int }_0^T\left|\overrightarrow{\mathit{WSS}}\right|\mathit{dt}}\right)$$E10

In fluid-structure interaction problems, neither Lagrange nor Euler formulation is fully suitable for solving the problem domain. The Lagrangian formulation, where the traced points are fixed to the mesh domain, is used for solids and the Eulerian formulation, where the points move along a fixed mesh domain, is used for fluids. In fluid-structure interaction problems, Euler formulation for a fluid model cannot accommodate the problem. In this case, a hybrid formulation known as the arbitrary Lagrangian-Eulerian (ALE) method is useful for solving the governing equations of FSI problems. Here the fluid field is allowed to deform arbitrarily according to the deformation of the structural fields. Solving the governing equations of fluid and artery wall deformation with the ALE method is not sufficient. It is also necessary to define the fluid-structure interaction region mathematically. In this study, the fluid-structure interaction region was solved by the Interface Quasi Newton Method using Least Squares (IQN-LS) method placed in the foam-extend-OpenFOAM [21]. Due to the elastic wall assumption on geometry, a free deformation condition was applied thus giving the solid freely deformation capability. The displacement and forces of the fluid is equal to the displacement and forces of the structure, at interaction region.

The accuracy of artificial neural network (ANN) models depends on the optimization of multiple parameters. The parameters consist of determining the optimal hidden layer count, selecting the appropriate number of nodes in the hidden layer, deciding on an effective learning rate, and identifying the ideal momentum factor values.

The number of neurons in the output should match the number of selected parameters used for predicting hemodynamic parameters as flow velocity. In the process of training an ANN, a suitable learning technique is utilized alongside a designated set of input and output data, enabling the network to effectively perform a given task. In order to obtain this desired outcome, it is necessary to adjust the weight coefficients of the neurons. The training process will persist uninterrupted until the network’s output matches the desired output perfectly. Weights and biases are repeatedly adjusted to minimize the difference between the network’s output and the desired output. The training process ends when either the error reaches a predetermined value, or the maximum number of epochs is reached [22, 23, 24].

Deep neural network (DNN) algorithms compute the output of each neuron in the related structure utilizing the following approach [25],

$$y_j=F\left(\sum _{i=1}^m{W}_{\mathit{ij}}{X}_i+B_j\right)$$E11

The input and hidden layers are denoted by subscripts “i” and “j,” respectively. The result of an artificial neural network (ANN) model is commonly denoted as [25],

$$Y=F\left[\sum _{j=1}^n{W}_{\mathit{kj}}\left[F\left(\sum _{i=1}^m{W}_{\mathit{ji}}{X}_i+B_j\right)\right]+B_k\right]$$E12

The rectified linear unit (ReLU) activation function was utilized in the method. The ReLU activation function applies a threshold operation to all input data, setting all values below zero to zero. The formula for this activation function is as follows [26],

$$f\left(x\right)=\mathit{max}\left(0,x\right)=\left\{\begin{array}{c}\begin{array}{c}{x}_i,x_i\ge 0\\ \end{array}\\ 0,x_i<0\end{array}\right.$$E13

$\mathit{\text{ReLU}}\left(x\right)=\mathit{max}\left(x,0\right)$ where x is output of the function. The loss function, referred to as mean squared error (MSE) definition is,

$$\mathit{MSE}=\frac{1}{n}\sum _{i=1}^n{\left(y_i-y_{\mathit{pi}}\right)}^2$$E14

The dataset has been prepared for the purpose of modeling and the inputs for prediction have been selected. Subsequently, the inputs are divided into training and testing datasets in order to normalize all data through scaling. The following step is to define the specification of a deep learning model featuring hidden layers. At this point, the selection of the model’s learning rate, activation function, and loss function is required.

$$\mathit{\text{Normalized data}}\left(\stackrel{\sim }{X}\right)=\frac{\left(\left(X\right)-\underset{\mathit{\text{value}}}{min}X\right)}{\left(\underset{\mathit{\text{value}}}{max}X-\underset{\mathit{\text{value}}}{min}X\right)}$$E15

R² values were evaluated for flow velocity resulting from DNN predictions calculation based on Python.

$$R^2=1-\frac{{\sum }_1^n{\left(y_i-y_{\mathit{pi}}\right)}^2}{{\sum }_1^n{\left(y_i\right)}^2}$$E16

2.2 Modeling

In this study, patient specific geometries of four patients are obtained from CT Angio scans as DICOM files. Then, the solid models of the carotid arteries for pre-operation cases are constructed with Slicer [27] from DICOM files and shown in Figure 1. Patients’ data were anonymized and used after the approval of the local ethics committee approval of the University.

Figure 1.

Solid model of carotid artery for (a) patient 1 (b) patient 2 (c) patient 3 (d) patient 4.

The carotid arteries were modeled as a linear elastic isotropic material. The mechanical properties of the arteries such as density, elasticity modulus, and Poisson ratio are taken as ρ_w = 1160 kg/m³, E = 1.106 × 10⁶ Pa and v = 0.45, respectively. These patients do not all have the same conditions; some have a normal carotid artery, but some have arterial dilation, and some have stenosis.

Pulsatile wave cycle of the patients can be obtained from Doppler ultrasonography. Here, the pulsatile waveform of the velocity is taken from literature [28, 29], then the pressure pulse cycle is generated by using four-element Windkessel model [30]. The pulsatile velocity cycle and pressure cycles are given in Figure 2.

Figure 2.

Pulsatile wave cycles.

In the carotid artery, the blood flow is assumed as laminar and incompressible. The blood density and dynamic viscosity are used as $\rho =1060\mathrm{kg}/{\mathrm{m}}^3$ and $\mu =3.71\times {10}^{-3}\mathrm{kg}/\mathrm{m}\mathrm{s}$, respectively. The Non-Newtonian flow’s Casson model parameters τ_y and μ₀ are taken as 0.0053 Pa and 0.0035 Pa s for normal hematocrit value (45%) in blood, respectively [31, 32, 33].

Mesh independency study was performed for Patient 1 due to its more complex geometry compared to other patients. Figure 3 shows the points where displacement and velocity values were taken for the mesh independency study. Displacement and velocity values were examined at the systolic peak. The mesh sizes for the fluid and structure geometries were calculated in approximately the same way for the other three patients.

Figure 3.

Mesh independency points.

Mesh independency study was performed for 0.18, 0.34, 0.58, 0.65, 0.85, 0.95, and 1.05 million tetrahedral elements, to determine the accuracy of the solution. The mesh constructed with 950,000 elements for fluid and 630,000 elements for structure found to be sufficient according to the results of the mesh independency study, as shown in Figure 4.

Figure 4.

Mesh independency.

The evaluation of the established dataset is performed using a multilayer perceptron (MLP) model. Figure 5 shows the algorithm used. Numerical analysis provided artery data for 80 points, including cross-section area, pressure, and flow velocity for three patients during DNN model construction at pulsatile flow conditions. Analyzes includes 80 data points obtained from the CCA and ICA (sinus and distal) regions of the arterial geometry. The training data consists of 80% of the numerical values, while the remaining 20% is used for testing. Velocity predictions for the fourth patient at various cross-sections were obtained, by deep neural network model. In the deep neural network model, two hidden layers were specified—one with 200 neurons and another with 50 neurons. In this case, the learning rate for the Adam version of stochastic gradient descent is set to 0.001, with momentum values of 0.9. The chosen loss function for fitting the model was mean squared error (MSE) [34]. The ReLU activation function and “he” weight initialization were used to build the model.

Figure 5.

DNN algorithm.

3. Results and discussion

Analyzes were made for all three patients with the OpenFOAM-based foam-extend program. Figure 6 shows the pressure and velocity streamlines from these CFD analysis. At the bifurcation point of the ICA, a flow separation occurs, and circulatory movement begins, which is especially clear in Patient 4.

Figure 6.

Pressure and streamlines at the systolic peak for (a) patient 1 (b) patient 2 (c) patient 3 (d) patient 4.

Reverse flow and circulation movement cause low values of wall shear stress distributions. Pressure distribution is inversely proportional to velocity. ECA of some patients shows higher velocity when compared with ICA. As blood flows through the CCA, it spreads throughout the vessel, but bifurcation, severely impairs the flow. After bifurcation the flow is mostly directed toward the inside wall of the branching arteries. In contrast, lower velocity values are observed in the outer regions. This situation can be clearly observed in Figure 7.

Figure 7.

Close view of streamlines at the systolic peak for (a) patient 1 (b) patient 2 (c) patient 3 (d) patient 4 at bifurcation point.

Figure 8 shows the hemodynamic parameters such as TWSS and OSI for Patient 1. Higher time averaged wall shear stresses were observed at the bifurcation point and at stenosis area. It is known from the literature that TWSS changes inversely with OSI. When TWSS increases, OSI decreases or vice versa. High OSI region could be responsible from atherosclerosis.

Figure 8.

Hemodynamic parameters distribution for patient 1 (a) velocity (b) TWSS (c) OSI.

Figure 9 gives the WSS at the outside wall of the CCA and ICA for pulsatile flow at the systolic peak. After sinus region there is a narrowing part (B) so WSS increases after that there is a dilation part (A) so WSS decreases.

Figure 9.

Wall shear stresses at the outer wall of CCA and ICA for patient 1 at systolic peak.

Figure 10 gives the WSS at the inside wall of ICA for pulsatile flow at the systolic peak. Point B shows the narrowing part at the sinus outlet, so WSS is higher than the carotid sinus region. Normally the wall shear stresses at the outside wall of ICA are lower than the wall shear stresses at the inside wall during the pulsatile cycle due to recirculation and flow separation. However, for Patient 1, the dilation and narrowing parts are located at the outer wall of ICA so because of this the wall shear stresses at the outside wall of ICA higher than inner wall of ICA.

Figure 10.

Wall shear stresses at the inside wall of ICA for patient 1 at systolic peak.

One can see the displacements at systolic peak for three patients in Figure 11. In all patients, maximum deformation occurred at the carotid apex, consistent with the literature, but in two patients, high deformation was also observed in the ICA just after the bifurcation point. The degree of disease and blockage ratio of the artery can affect the deformation. The displacement is zero at the inlet and outlet.

Figure 11.

Displacements at systolic peak for (a) patient 1 (b) patient 2 (c) patient 3.

The DNN model accuracy of the results is achieved as 0.9929 (R²). Table 1 gives the comparison of machine learning results for average velocity with CFD for different cross-sections of Patient 4. Slice locations are numbered from left to right and are shown in Figure 12. The average velocity and average pressure values based on the cross-sectional area of the first three patients were used for DNN learning. Then, Patient 4 was asked to estimate the average velocity values in seven cross-sectional areas. Table 1 shows the comparison of these predicted values with the average velocity values obtained as a result of CFD analysis.

Slice number	Area (A) [m^2]	Average pressure from CFD (P/ρ) [m^2/s^2]	Average velocity from CFD (U) [m/s]	Predicted average velocity values [m/s]
1	2.96479e−05	16.3348	0.389737	0.4183827
2	3.04849e−05	16.2483	0.472065	0.4648768
3	5.01389e−05	16.2795	0.243602	0.2661422
4	2.33626e−05	15.9316	0.560831	0.4901191
5	8.0059e−06	14.0544	1.39763	1.2645323
6	2.83555e−05	15.6028	0.535595	0.604400
7	6.62852e−05	16.2678	0.195996	0.202197

Table 1.

Comparison of machine learning results with CFD for various cross-sections of patient 4.

Figure 12.

Slice locations for patient 4.

The value of the Mean Absolute Percentage Error in Table 1 is 0.08038. The average of the absolute percentage differences between the estimated values and the actual values give the Mean Absolute Percentage Error (MAPE). It evaluates the precision of a forecasting technique. The metric calculates the mean of the absolute percentage errors for each entry in a dataset to determine the accuracy of the anticipated quantities relative to the actual numbers. A Mean Absolute Percentage Error below 5% is deemed to be indicative of a forecast that is sufficiently accurate. A Mean Absolute Percentage Error over 10% but decreasing below 25% denotes a level of accuracy that is low, yet still within an acceptable range. However, a MAPE above 25% implies a far lower level of accuracy.

4. Conclusions

This chapter investigates the blood flow in carotid artery with three-dimensional realistic carotid artery geometry of four patient. The numerical simulations was done with OpenFOAM in order to obtain pressure contours and velocity streamlines. Pulsatile velocity inlet and pulsatile pressure outlet conditions were used at the CCA inlet and the ICA-ECA outlets like as the human body model.

Common carotid artery has two branches; one is external carotid artery, which supplies blood to the superficial structures in the neck and face regions, the other is internal carotid artery, which has no extracranial side branches and supplies blood to the brain. Since ICA supplies blood to the brain, plaque formation and stenosis that effects laminar blood flow on this branch could cause thromboembolic situations and stroke. The low wall shear stress areas are responsible for the formation of atherosclerosis. Due to flow separation, the wall shear stresses at the outside wall of the carotid bifurcation during the pulsatile cycle are lower than the wall shear stresses at the inside wall. Also, due to the recirculation of the flow at the carotid sinus, this region has low wall shear stresses at the outside wall. The results in this study also confirm these inferences in accordance with the literature. The average velocity and average pressure values based on the cross-sectional area of the first three patients were used for learning and testing of machine learning algorithm. After that, the algorithm was asked to estimate the average velocity values in seven cross-sectional areas for Patient 4. The comparison of these predicted values with the average velocity values obtained as a result of CFD analysis gives good accuracy. The accuracy of the results related to machine learning can be increased with the data obtained through numerical analysis using realistic vascular geometries in more patients and different cardiovascular diseases. It can also be used by varying it according to hemodynamic conditions. We found the results promising in terms of guiding the clinical decisions.

Acknowledgments

This research did not receive any funding from funding agencies in the public, commercial, or not for profit sectors.

Conflict of interest

The authors declare no conflict of interest.

Appendices and nomenclature

CCA	common carotid artery
ICA	internal carotid artery
ECA	external carotid artery
WSS	wall shear stress
TWSS	time-averaged wall shear stress
OSI	oscillatory shear index
FSI	fluid-structure interaction
CFD	computational fluid dynamics
ALE	arbitrary Lagrangian-Eulerian
IQN-LS	Interface Quasi Newton Method using Least Squares
ANN	artificial neural network
DNN	deep neural network
CT	computed tomography
CTA	computed tomography angiography
MRA	magnetic resonance angiography
DSA	digital subtraction angiography
ReLU	the rectified linear unit
MLP	multilayer perceptron
MAPE	Mean Absolute Percentage Error
u	velocity vector
ρ	fluid density
$\mathit{\mu }$	dynamic viscosity
P	pressure
f	gravity force
E	elasticity modulus
$\mathit{\upsilon }$	Poisson ratio
${\mathit{\rho }}_{\mathit{w}}$	solid density
σ	Cauchy stress tensor
$\overrightarrow{\mathit{b}}$	body forces
w	wall deformation
${\mathit{\mu }}_{\mathit{L}}$	first lame parameters
${\mathit{\lambda }}_{\mathit{L}}$	second lame parameters
$\mathit{\epsilon }$	strain tensor
tr	trace function
I	identity matrix
τ	shear stress
${\mathit{\tau }}_{\mathit{y}}$	yield stress
${\mathbf{\mu }}_{\mathbf{0}}$	Casson viscosity
$\dot{\mathit{\gamma }}$	shear strain rate
${\overrightarrow{\mathit{\tau }}}_{\mathit{w}}$	wall shear stress vector
T	period of the pulsatile cycle
yj	output from the jth neuron
M	the total number of input variables
Xi	input for the model
Wij	weights
Bj	biases
Y	output of ANN
n	number of hidden layers
k	output layer
F	transfer function used to normalize the neuron’s output
${\mathit{y}}_{\mathit{i}}$	value from numerical analysis
R^2	coefficient of determination
ypi	predicted value