Source data for mathematical modeling
The original geometric model of the vessel was reconstructed being based on a preoperative CTA examination of the affected left carotid artery of a particular patient. Fig. 1 (a) presents the flow area in the affected part of the vessel (the carotid artery is in the foreground). The curve segment in Fig. 1 (a) shows the supposed vessel wall position under the atherosclerotic plaque. The dotted line marks the plaque inner surface in the projection shown.
The source data for the flow simulation are the patient’s postoperative UDV (Ultrasound Doppler Velocimetry) results. The dependence of the inlet flow rate on time is based on the data from the UDV study of the CCA (Common Carotid Artery). The ratio of flows through ICA (Internal Carotid Artery) and ECA (External Carotid Artery) was calculated using their cross-sectional areas and TAPV (Time-Averaged Peak Velocity). They were also obtained from the UDV data.
Construction of geometric models
The reconstructed three-dimensional model of a healthy vessel was built using SimVascular  and Salome  software. SimVascular was used to build vessel segments (see Fig. 1 (b-d)). Using a Python script of our design, these segments were later imported into Salome to build a geometric model and computational meshes. The reconstructed three-dimensional vessel model is shown in Fig. 1 (e). Further, this model is referred to as the base one and is denoted as m0. The black line on the vessel wall in Fig. 1 (e) indicates the incision line for the subsequent imitation of patch implantation.
The mentioned script was also used to visually simulate the result of a CEA surgery. It visually helps to draw the incision line and patch outlines on the vessel segmentation contours (Fig. 2). The information about the drawn lines is exported by the script to a data file as a list of the patch width values at its intersections with the segmentation contours. The script then uses this data to change the base model geometric shape, simulating the result of the patch implantation.
The modified models are the results of a virtual CEA. They were built by increasing (or decreasing) the perimeters of all segments of the m0 model that intersect with the incision line. This made it possible to simulate any form of patch or incision closure without a patch. Changing the perimeters was achieved by scaling segments relative to their geometric centers following the values recorded in the data file. By default, the scaling coefficient is assumed to be equal to the value of the relative increment of the perimeter stored in the file, divided by 2π. This choice is quite accurate for segments close in shape to ellipses with a small eccentricity. The scaling coefficient can be desirably adjusted manually for irregular segments. As a result of scaling, the upper (in Fig. 2) parts of the ICA proximal contours shifted up and began to intersect with the proximal contours of the ECA (not shown in Fig. 2). To correct this, one performed a small parallel displacement of all the ECA contours along the larger axis of the CCA distal contour. After that, a geometric model of the vessel and computational meshes on it were built. The meshes were then exported to the OpenFoam  for numerical calculations.
For comparative analysis, geometric models m1-m10 were constructed using the method described above. The m1-m9 models simulate the results of CEA surgery on the m0 model with the implantation of p1-p9 patches, respectively (Fig. 3). The m10 model (not shown in Fig. 3) simulates the suture of the incision without patch implantation.
The data on the shapes of the patches are given in Table 1, which contains the widths of the patches in their cross-sections along the incision line. Point 0 of the column "Distance…" corresponds to the proximal end of the patch, point 3.9 does to the distal one. The incision line (shown in Fig. 1 (e)) was the same for all the models. Both its length of 3.9 cm and its location corresponding to the actual incision made during the CEA surgery. The m10 model was constructed by reducing the circumference of the vessel lumen along the incision line. The negative values of the virtual patch p10 width in Table 1 indicate this. Note that Table 1 shows the increments of the vessel perimeter after implantation, while the width of the patch itself before implantation should be slightly larger.
The flow velocity U and pressure p in the constructed geometrical models were described through the three-dimensional unsteady Navier-Stokes equations for the viscous incompressible fluid :
at constant density ρ = 1050 kg/m3 and dynamic viscosity µ = 3.675·10−3 Pa∙s, here τ is the shear stress tensor. The no-slip boundary condition was set for U on the side surface of the flow region, and the parallel flow conditions were set on the inlet and outlets. The flow area boundaries were considered to be rigid. The initial velocity value was set to a constant of 0.15 m/s.
A periodically changing pressure difference was set as the boundary conditions for p at the inlet and outlets. It created a periodic flow with parameters corresponding to the postoperative data of the patient's UDV:
• T = 1.06 s – a period of the cardiac cycle;
• Q = 6.9 ml/s volumetric flow rate through CCA;
• r = 1.72 the ratio of the volumetric flow rate through ICA to the volumetric flow rate through ECA.
The method for constructing boundary conditions was as follows. First, a preliminary CFD simulation was performed where the zero pressure was set at both outlets and a suitable pressure curve was built at the inlet stepwise with a time step of 10 ms. During this simulation, the pressure increments (decrements) at each time step were selected manually. This made it possible to obtain an inlet velocity curve (Fig. 4 (a)) corresponding to the envelope of the UDV spectrum in the patient's CCA. Based on these simulation results, the ratio r was calculated. It turned out to be different from the target value of 1.72. Therefore, the time-dependent pressure was set at the outlets instead of zero to adjust the r-value. Specifically, the same pressure curve was set at the ICA outlet as at the CCA inlet. It was only reduced in magnitude with a coefficient of k = 0.1 and a small phase lag was set. The same pressure curve was set at the ECA outlet as at the ICA one. It was only inverted relative to the abscissa axis. Thus, the ECA flow decreased, and the ICA flow increased. The CCA flow remained practically unchanged. After that, using another series of auxiliary simulations, the value of the parameter k (and the pressure curve shape, if needed) were adjusted so that the Q and r quantities became close to their target values.
Numerical calculations and post-processing
The flow simulation was performed in OpenFoam with the Finite Volume method using the PISO algorithm [34, 35]. Along with commercial software such as Ansys Fluent, OpenFoam is a common tool for performing hydrodynamic calculations and modeling flow in blood vessels (see for example, [4, 18, 36]). In the preliminary numerical calculations described in section 2.3, rough meshes were used to obtain the flow corresponding to the UDV data. After obtaining satisfactory results of preliminary calculations, the final calculations were made on fine meshes. The numerical calculations resulted in the dynamic fields of pressure, velocity, and velocity gradient in the flow domain for several cardiac cycles with a time sampling of 10−2 s. The information about various flow parameters and their derived characteristics was extracted from the calculation results by postprocessing performed in ParaView .
To verify these results, a mesh independence study was conducted. It was found that the results do not significantly change when using meshes with the number of nodes more than 5·105. The mesh cell size on the side surface was set equal to half the cell size inside the computational domain. This was done to improve the accuracy of the calculation of hemodynamic indices which are expressed through the velocity gradient on the vessel wall. The stabilization of pulse oscillations was also studied. As it turns, the process of pulse oscillations can be considered as stabilized one starting with the second cardiac cycle (see also ). In this regard, it is further considered that the time t = 0 corresponds to the beginning of the systolic phase of the second cardiac cycle.
Wall Shear Stress and hemodynamic indices
The WSS index was calculated as the tangential component tw of the shear stress tensor τ on the vessel wall. The TAWSS, OSI, and RRT indices were calculated with the formulas (2) through the averaged value of tw over the one cardiac cycle :
Here T is the duration of the cardiaс cycle and |τw| is the Euclidean norm of the vector τw. For quantitative comparison of the indices in a certain zone σ on the vessel wall their average values were calculated using the formulas:
where S is the area of the zone σ. Dimensionless mean I_RRT and logarithmic maximum M_RRT values of the RRT index for the zone σ were also used, calculated by formulas:
Here TAWSSCCA is the averaged TAWSS value in the cylindrical part of the ССA at a distance of three ССA radii from the vessel bifurcation .
To assess the pathological values of the RRT index, the estimates TAWSS < 0.4 Pa and OSI > 0.3, given in [9, 28], were used. In this case, formula (2) gives the corresponding critical value RRT = 6.25 Pa−1. Further, RRT values will be considered pathological if they are greater than the critical ones.