Non-Kolmogorov turbulence in carotid artery stenosis and the impact of carotid stenting on near-wall turbulence

Blood flow regime and dynamics play important roles in the development of numerous vascular diseases, such as atherosclerosis,^1,2 stenosis,³ and aneurysm.⁴ Endothelial cells (ECs) respond to hemodynamics through numerous pathways that alter the pathophysiology of the arterial wall at many levels.^5,6 Rapid transient variations in blood hemodynamic pattern, which could be associated with quasi-periodic, transitional, or turbulent flow features, trigger different mechanobiological mechanisms, many of which are yet to be identified and characterized.⁷ Therefore, it is crucial for the progress of vascular medicine to explore and comprehend the physics of such flow regimes as they manifest different layers of complexity in different vascular diseases.

Since the emergence of Computational Fluid Dynamics (CFD) in hemodynamics research in the early 1990s, considerable effort has been made to investigate transitional flow (trans-flow) in carotid artery stenosis.^8,9 Most of the CFD models were based on certain physical assumptions that could be inappropriate to unfold the complexity of the trans-flow problem as it exists. For example, a number of critical studies^10–13 inveighed against the use of Newtonian viscosity assumption that has been promoted in the majority of published CFD studies on carotid artery stenosis hemodynamics.^14–19 In two recent studies, our group have evidently shown in vivo that the Newtonian assumption is inappropriate to represent the Wall Shear Stress (WSS) of the internal carotid artery (ICA).^20,21 Wall Shear Stress (WSS), which is the main hemodynamic parameter investigated in mainstream CFD research,^22,23 is strongly dependent on the closure model of the viscous term in the Navier–Stokes equation, and it does not possess intrinsic vector properties.²⁴ Moreover, the use of ideal two-dimensional axisymmetric models to investigate the stenotic flow patterns²⁵ was shown to be inadequate.²⁶ This casts nontrivial uncertainty on the essential elements of the current paradigm of carotid artery stenosis hemodynamics.²⁷ 

Turbulence remains one of the unsolved problems in physics, and even transition to turbulence in a straight pipe poses unanswered questions until today.²⁸ Thus, studying turbulence in patient-specific geometries, such as carotid artery stenosis, induces further complications into the problems brought about by the complex geometry and the complex multiharmonic flow that is not fully characterized nor understood.^29,30 Nevertheless, there are some interesting characteristics of turbulence in carotid artery stenosis that have been explored by high-fidelity CFD techniques that can provide quasi-exact numerical solutions of the Navier–Stokes equation. Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES) are robust CFD techniques that can provide physically accurate simulations of the carotid artery stenosis flow.³¹ 

Literature records report only a handful of studies using these techniques. Lee et al.⁸ conducted DNS using the Spectral Element Method (SEM) to explore the trans-flow features in a stenosed carotid bifurcation. They demonstrated interesting flow features occurring due to the combination of geometric and inflow waveform effects. The work by Lee et al. highlighted the existence of strong Dean vortices, shear layer destabilization, and transition to weak turbulence during systolic acceleration. Tan et al.³¹ conducted a thorough numerical investigation to demonstrate the merit of LES over RANS (Reynolds-averaged Navier Stokes) turbulence models in modeling stenotic flows. By simulating a steady Newtonian flow past an ideal stenosis, they showed how LES can capture the mean-flow phenomena, such as post-stenotic relaminarization, jet separation, and reattachment. One peculiar finding of their study is the argument that claims that LES is better than DNS in replicating experimental results.

Jabir and Lal⁹ carried out a LES study of multiharmonic non-Newtonian stenotic flow in an ideal geometry to investigate the transition to turbulence. They used DNS and experimental velocity measurements from the literature to establish evidence of their model validity. Such evidence was based on the mean flow velocity, and no information was reported concerning the resolved turbulence kinetic energy with respect to the total energy budget of flow. Their study focused on the investigation of mean flow phenomena, such as Reynolds stresses, coherent structures, and the rms fluctuations of the velocity components. They also showed the turbulence energy cascade as a function of the Strouhal number. Such cascades exhibited a Kolmogorov inertial plateau in addition to other plateaus that seem to be of non-Kolmogorov type, both consuming almost equal frequency bandwidths, unlike the case in fully developed turbulent flow (i.e., Kolmogorov turbulence). Their interesting findings were later extended and better observed by Lancellotti et al.³² in a remarkable LES study. Although they used the Newtonian fluid model, they managed to capture the momentum interactions that lead to intermittency events in physiologic flow conditions. Their study is perhaps the first study to openly argue that the Kolmogorov theory is inapplicable in physiologic turbulence.

The findings of Lancellotti et al. were later supported by Mancini et al.³³ as they openly stated that their well-resolved DNS and LES models did not predict Kolmogorov turbulence, and the theory assumptions could not be met. Their study solved grids with up to 50 × 10⁶ cells in resolution for a multiharmonic flow with a mean Reynolds number of 980. Similar results of non-Kolmogorov turbulence in various stenotic flows were recently depicted in a study by Ozden et al.;³⁴ however, they were not properly discussed. Ozden et al. used ideal geometries with large eddy simulation to reproduce turbulence-induced sound patterns that are close to those obtained from Doppler ultrasound measurements. The depicted power spectrum density of wall pressure fluctuations which demonstrated non-Kolmogorov plateaus and the inverse energy cascade in some geometrical conditions. A summary of the important parameters of such studies is given in Table I.

A recent study by Saqr et al.³⁶ proposed a body of evidence suggesting that physiologic blood flow is inherently turbulent. They showed that blood flow exhibits a sensitive dependence on initial conditions, global instability, and significant levels of the kinetic energy cascade. The study argued that the seven-decade old physiologic blood flow paradigm is inadequate to treat multiharmonic pulsatile flow, let alone pathologic mechanisms and underlying mechanobiological pathways. The study also demonstrated that physiologic blood flow has kinetic energy cascades of non-Kolmogorov type in doctrina and in vivo. Therefore, and in principle, the terminology used with fully developed turbulent flow should be redefined when used to describe blood flow and hemodynamic studies. A previous meta-analysis and review by Saqr et al.³⁶ demonstrated that the old paradigm hindered research on intracranial aneurysm hemodynamics for the past 30 years. The meta-analysis critically established such critical argument based on 1733 published studies covering the elemental research topics comprising the aneurysm disease model of hemodynamics and mechanobiology. Moreover, last year, Rashad et al.³⁷ discovered five mechano-miRNAs with potential pathologic roles under physiologic flow conditions that could not have been reproduced in vitro using the old paradigm of vascular hemodynamics.

The purpose of the present work is to investigate the effects of stenting on the transitional hemodynamics of carotid stents by high-resolution large eddy simulation (HRLES). The CFD model is established to avoid all unnecessary physical assumptions that could produce physical errors in the simulation. Blood is treated here as a non-Newtonian shear-thinning fluid subjected to multiharmonic boundary conditions that represent actual waveform in the carotid artery.

The HRLES approach for resolving the turbulence field is established using Pope’s³⁸ resolution quality measures. Such criteria are used to calculate the ratio between resolved and modeled turbulent kinetic energies to evaluate the validity of a LES model. The HRLES model is based on a patient-specific geometry using a case in which pre- and post-stenting geometries were available for comparison. To the best of our knowledge, this is the first study to explore the therapeutic effects of carotid stenting in light of non-Newtonian transitional flow associated with stenosis and its subsequent effects on near-wall hemodynamics.

This study was approved by the ethics committee of the Asahi University Hospital (Gifu, Japan; Approval No. 2020-01-09), and informed consent was obtained from the patient. All patients were handled in accordance with the Declaration of Helsinki.

An 80-year-old male patient was selected for this study. He had a history of myocardial infarction, atrial fibrillation, diabetes mellitus, and dyslipidemia and was referred to our service for detailed evaluation of left carotid artery stenosis. He was receiving medications in the form of warfarin (1.75 mg/day), aspirin (100 mg/day), ezetimibe (10 mg/day), glimepiride (0.5 mg/day), and teneligliptin hydrobromide hydrate (20 mg/day). Neurological examination demonstrated no abnormal findings. Magnetic resonance angiography (MRA) revealed a mild stenosis in the left internal carotid artery (ICA). The patient was followed up using MRA, but the stenosis of the left ICA developed over time [Fig. 1(a)]. Four years after the first examination, computed tomography angiography (CTA) demonstrated 80% stenosis in the left ICA as shown in Figs. 1(b) and 1(c).

The patient was operated on using endovascular stenting with no complications [PRECISE PRO RX (Cardinal Health, Inc. USA)]; stent characteristics: outer diameter, 8 mm; initial length, 30 mm; stent design, peak to valley; configuration, straight; cell type, open-cell; manufacturing, laser-cut; strut thickness, 0.23 mm; stent width, 0.10 mm; open-cell size, 5.7 mm²; and material, nitinol, as indicated in Figs. 1(d) and 1(e). The postoperative course went uneventful, and he was discharged 4 days after the surgery. No in-stent restenosis existed in the CTA 1 year after the stenting, as demonstrated in Fig. 1(f).

The 3D DICOM images of the MRA and CTA were then imported into Materialise Mimics software (Materialise, NV, Belgium) for visualization, stereolithography (STL) model creation, and initial fine-tuning. The initial models were exported to 3-Matic software (Materialise, NV, Belgium) for further correction and surface tuning. The rate of volume change was suppressed to ≤5% during the smoothing process.³⁹ The final models were exported as STL format files.

Two computational models were created for the patient, one for the pre-stenting condition and one for the post-stenting condition. Models were analyzed using our published high-resolution large eddy simulation methodology (HRLES).⁴⁰ The frequency-domain results of both models were analyzed for two cardiac cycles after the initial flow effects were discarded. The initial flow effects were identified at 10% of the first cycle time. The Courant–Friedrichs–Lewy (CFL) condition was kept below unity in all simulations to ensure solution stability and spatial convergence in each time step. The flow variables are compared spatially between the upstream (proximal) and downstream (distal) regions of the stenosis in each case. The primary parametric comparison is conducted between the flow variables before (pre-stenting) and one year after the endovascular procedure (post-stenting).

The STL files generated from patient angiography imaging were discretized to a computational grid using ICEM-CFD by ANSYS, Inc. ICEM produces computational domains optimized for the finite volume method. Figure 2 shows the snapshots of the computational grid domain, and Table II shows its details. The flow exhibits sub-critical turbulence of non-Kolmogorov type, similar to that measured and deduced by Saqr et al.³⁶ Therefore, it was possible neither to evaluate the inertial and dissipative length scales nor to apply Taylor’s frozen turbulence hypothesis. However, to determine if the mesh resolution was sufficient to resolve the flow in space and time, Pope's critera have been used to determine the LES mesh quality index. This is detailed in the next section.

The present HRLES model was established and solved using the FLUENT code by ANSYS, Inc.⁴⁰ The solver uses a central differencing scheme for the spatial discretization of the governing equations and second order implicit time marching scheme. The Pressure-Implicit with Splitting of Operators (PISO) scheme was used to solve the momentum and continuity equations in a segregated solver. In each time step, the iterative solver was set to converge at residuals of the order of 10⁻⁶. These settings are well established for HRLES models in numerous transitional flow applications.31,41 

To ensure that the HRLES solution produces physically accurate results, the grid was established such that it achieves Pope’s criteria for LES,38,42–44 which ensures that the resolved turbulence kinetic energy in LES is not less than 90% of the total turbulence budget of the flow. This was verified by calculating Pope’s resolution measure in the space–time domain, as shown in Fig. 3. Hence, the validity of the HRLES model was established and verified. Pope’s LES quality measure presented in Fig. 3 is time-averaged over two cardiac cycles, and it demonstrates that in both cases, the ratio between resolved and total kinetic energies above is 0.95 in 99% of the grid cells. It has been demonstrated by Tan et al.³¹ that LES produces results that are sufficiently close to DNS at such high-resolution.

The Newtonian assumption used in cerebral hemodynamic research was recently shown to be inadequate.24,45,46 Therefore, the present work adopted the Carreau–Yasuda non-Newtonian viscosity model to capture the shear-thinning properties of blood. The model choice was motivated by its commonly wide use in cerebral hemodynamic CFD models.47,48 The model expresses blood viscosity as follows:⁴⁹ 

where μ_∞ = 0.0022 Pa s, μ₀ = 0.022 Pa s, λ = 0.11 s, a = 0.644, and n = 0.392.

In order to capture the shear-thinning effects, we propose the local viscosity index that is the ratio between local instantaneous viscosity and local time-averaged viscosity (i.e., local Newtonian viscosity). The local viscosity index, ψ, is expressed as follows:

The physical meaning of ψ represents the local level of blood shear-thinning effects over time.

Blood flow was represented by a multiharmonic inflow boundary condition in the form of a sum of harmonics expressed as follows:

The Womersley coefficients are adopted based on averaged values for carotid blood flow. Given the diameter of the artery at the inlet, the boundary inflow waveform and the corresponding Reynolds number are shown in Fig. 4. Since blood is considered as a Carreau–Yasuda non-Newtonian fluid, the Reynolds number was calculated using the generalized Metzner–Reed50,51 form: $ReMR=DnV2−nρλ8n−13n+14nn$⁠.

The Womersley number, used to describe the time-dependent ratio between inertia and viscosity, is difficult to be used to characterize multiharmonic pulsatile flow. This is because the latter has a bandwidth of frequencies rather than a single frequency value, such as in the former. The authors could not find a definition of the Womersley number for multiharmonic flow in the literature and are now working to establish one. In the present work, the Reynolds number is represented as an average value over one cycle and as a field function of time as shown in Fig. 4. The minimum and maximum Reynolds numbers were 64 and 147, respectively. The time-average Reynolds number was 115 ± 26 over one cycle although the application of Reynolds decomposition in multiharmonic pulsatile flow is inaccurate according to Mancini et al.³³ Therefore, the time-dependent Reynolds number was calculated at the inlet plane of both computational models. The variation between Reynolds number calculations in the two models was sufficiently small to neglect.

Figure 5 shows the time-averaged velocity streamlines of the two cases. The bifurcation region is characterized by helical flow that becomes straight in the stenosis regions and recovers downstream the stenosis region producing secondary flow.

Figure 6 shows the contours of normalized wall shear stress (WSS) at different instances of a cardiac cycle in pre-stent and post-stent conditions. Prior to stenting, the maximum normalized WSS occurs downstream the stenosis region during mid-acceleration, peak systole, and mid deceleration instances. Post-stenting, WSS appears to be more homogeneous at all instants of the pulse cycle. In addition, the resistance to flow is identified as the ratio between pressure drop across the stenosis and its length, expressed as $Rf=ΔpstenosisLstenosis$ Pa m⁻¹. The flow resistance in pre-stenting was 119.2 kPa m⁻¹, which was reduced approximately ten-fold in post-stenting to 14 kPa m⁻¹.

To monitor local turbulence in the domain, six points were identified on the artery centerline in both cases, as shown in Fig. 7. Then, as in Fig. 8, the velocity time-series are compared in pre-stenting and post-stenting cases. Downstream the flow inlet plane, at point 1, the flow is periodic inheriting the base harmonics of the boundary conditions in both cases, as shown in Fig. 8. Downstream the bifurcation apex, the flow conserves its periodic nature without developing any secondary harmonics, as shown in point 2.

In point 3 at the middle of the stenosis region, the flow exhibits weak turbulence features characterized by small-amplitude oscillations in the pre-stenting case. In the post-stenting case, the latter oscillations are mostly diminished, and the periodic state is restored. Downstream the stenosis, the latter phenomena are significantly manifested as depicted in the velocity trace at point 4. The high-amplitude oscillations observable prior to stenting disappeared post-stenting, and the chaotic behavior changes to periodic oscillations characterized by the base harmonics. In the non-stenosed branch, a similar behavior is qualitatively depicted. These oscillations are apparently similar to cycle-invariant turbulent-like flows reported by Khan et al.⁵² 

It is essential to investigate the cascade of kinetic energy to characterize hemodynamic turbulence and its variation in both cases. Figure 9 depicts such a cascade at six points. The kinetic energy is calculated using the Fourier transform of the velocity time-series.⁵³ Clearly, the turbulence cascades shown in Fig. 9 is of non-Kolmogorov type. It is also obvious that such turbulence exhibits Kraichnan’s inverse energy cascade plateau⁵⁴ at higher frequencies. A characteristic frequency (f = 50 Hz) is observed to associate with a sudden rise in energy at some locations, as marked by a black square symbol in Fig. 9. After stenting, most of the locations are characterized by the inverse energy cascade compared with pre-stent conditions at f > 50 Hz.

The local viscosity index is traced in both cases at six centerline points, as shown in Fig. 10. The momentum–viscous interactions associated with local turbulence can easily be observed downstream the stenosis (point 4). In pre-stenting, blood shows more periodic variations in ψ compared to the post-stenting. It can also be observed that in the former condition, there are small-amplitude inertial-viscous oscillations at the peak-systole instance compared with the pre-stenting condition.

Near-wall hemodynamics is crucial to find possible relevance between the findings of CFD models and pathobiology of endothelial cells. Transition to turbulence was investigated at two points located at a distance of 10 µm from the wall. Such a distance was chosen based on physical reasoning to represent the influential flow scale with respect to the thickness of the endothelial layer, which varies from 0.1 to 10 µm. One point is located at the bifurcation apex, and the second point is in the middle of the stenosis throat. Figure 11 shows the velocity profiles and oscillations (top row). The post-stenting recovery of velocity signals in the near-wall region was found to be better than that in the mean-flow region. Figure 11 also shows the trace of kinetic energy cascades. It is easily observed that the characteristic energy spike at f = 50 Hz is maintained near the wall at both locations in both cases. It is also evident that the energy magnitude levels are orders of magnitude less in the post-stent case near the stenosis wall. It is also observed that the inverse energy cascade appears in the pre-stent case at both locations, the bifurcation apex and stenosis throat.

In the present work, large-amplitude velocity oscillations were observed as shown in Fig. 8 (points 3 and 4), which represent flow transition and destabilization. The momentum–viscous interactions due to the shear-thinning properties of blood are vitally important in characterizing this peculiar transitional flow. The resemblance between local velocity fluctuations downstream the stenosis (Fig. 8, point 4) and local viscosity fluctuations at the same location (Fig. 10, point 4) demonstrates such importance. It is a fact that the effects of shear-thinning properties on flow stability and transition are not fully understood in pulsatile flow, whether it is monoharmonic or multiharmonic, with the latter being more complex. However, one of the established observations, as shown in the study by Jabir and Lal,⁹ shows that shear-thinning properties create a highly non-linear energy cascade and inertial sub-range with multiple cascade regimes. This is evident in the present work where energy cascades in pre- and post-stenting conditions exhibit a non-linear inertial subrange, as demonstrated in Fig. 9.

The inverse energy cascade phenomenon is very new in the medical literature, and only a few observations relevant to the present work exist.56,57 Therefore, the impact of different turbulence regimes and energy cascades on vascular biology is a largely uncharted territory. We examined the near-wall turbulence cascades. We observed a sudden rise in kinetic energy (energy spike) at f = 50 Hz, which could indicate the hydrodynamic resonance of the harmonic flow waveforms resulting from vortex formation and breakup, relatively similar to the phenomena reported by Cotton.⁵⁵ The origin of kinetic energy resonant interaction in transitional flow was first characterized by Kachanov and Levchenko.⁵⁶ They showed that transitional flows inherently produce synchronized vortex breakdown frequencies that build up resonant interactions. In pulsatile flow, where the flow is driven by the multiharmonic waveform, transitional resonant interactions become more significant.

It is yet unclear how different turbulence regimes would impact endothelial cells as it is an understudied subject in vascular biology. While indeed turbulent flow (dubbed as disturbed flow in the medicine and biology literature) adversely impacts endothelial cells when compared to laminar flow,⁵⁷ no study in the biology literature has compared different turbulence regimes, which offers an interesting opportunity for novel findings to be elucidated. Therefore, the question remains whether the inverse cascade is a favorable flow regime or a predictor for future pathology and whether endothelial cells could sense different turbulent energy cascades in a differential manner. In Fig. 11, the change in turbulence kinetic energy associated with the inverse energy cascade frequency range (10⁻¹¹– 10⁻⁶ m²/s²) is proportional to forces of pico-Newton scale (10⁻¹³– 10⁻¹⁰ N) at such a distance. Mechanobiological functions of endothelial cells are regulated by forces of such a biologic scale.58,59 

In the present work, we opted not to simulate the stent itself. While this might appear a simplistic approach, we reasoned that following stent deployment, within several weeks to a few months, endothelial regrowth and recovery render the geometry of the stent insignificant to the flow structure and characteristics.60,61 Moreover, vascular remodeling after stent deployment occurs at scales of the same order as the characteristic length scale of the flow.⁶² There was almost one year between the pre-stenting and post-stenting cases, making our argument about the insignificant effect of the stent on near-wall hemodynamics valid for the scope of this work. However, it is important to note that within the near-wall region we selected, the blood is mostly cell free plasma. This is a limitation in our work as our simulation does not account for the variation in the blood composition at varying distances from the wall. While simplification is indeed welcomed to understand the complex physical phenomena we wanted to elucidate, it should be kept in mind and hopefully addressed in the future with more complex simulations.

The objective of this work was to explore the long-term hemodynamic changes in carotid artery stenosis following the endovascular stenting procedure. Now, a growing body of knowledge postulates that physiologic and pathologic blood flow exhibits non-Kolmogorov turbulence. The non-Newtonian HRLES approach used in the present work captured over 95% of the turbulent kinetic energy budget in the flow, revealing peculiar observations related to momentum–viscous interactions, turbulence energy cascade, and near-wall turbulence. Due to the limitations of the existing theory of turbulence and the lack of analytical tools to describe and explain these findings, the observations were presented in representative locations in time and frequency domains. These findings can be concluded as follows:

• Endovascular stenting recovers physiologic blood flow in the stenosed carotid as presented by the instantaneous velocity and viscosity fields.

• The cascade of turbulent kinetic energy exhibited inverse plateaus similar in quality to Kraichnan’s theory. A particular interest of this study is the emergence of an energy spike at f = 50 Hz in all upstream locations in the post-stenting case.

• The inverse energy cascade is conjectured to be related to morphologically promoted hydrodynamic resonance; however, more evidence and physical reasoning are required to clarify such a conjecture.

• Near-wall turbulence shows energy levels that correspond to forces at pico-Newton scales, which are known to regulate epigenetic and metabolic functions of endothelial cells.

This work was supported by a grant-in-aid for Young Scientists (A) under Grant No. 17H04745 (to Niizuma) from the Japan Society for the Promotion of Science.

The authors have no conflicts to disclose.

Khalid M. Saqr, Sherif Rashad, Kuniyasu Niizuma, Toru Iwama, and Teiji Tominaga conceived and designed the study. Kiyomitsu Kano and Yasuhiko Kaku collected the patient data. Sherif Rashad created STL models. Khalid M. Saqr performed CFD simulations. Khalid M. Saqr and Sherif Rashad performed post-processing. Khalid M. Saqr, Sherif Rashad, and Kiyomitsu Kano analyzed the data. Khalid M. Saqr, Kiyomitsu Kano, Sherif Rashad, and Kuniyasu Niizuma wrote the manuscript. Kuniyasu Niizuma acquired the funding and administered the project. Yasuhiko Kaku, Toru Iwama, and Teiji Tominaga critically revised the manuscript. All authors approved the final version of the manuscript. Khalid M. Saqr and Kiyomitsu Kano contributed equally to this work.

The data that support the findings of this study are available within the article.

CTA

computed tomographic angiography

DNS

direct numerical simulation

HRLES

high-resolution large eddy simulation

MRA

magnetic resonance angiography

Trans-Flow

transitional flow

WSS

wall shear stress