Non-Newtonian Endothelial Shear Stress Simulation: Does It Matter?

Patient-specific coronary endothelial shear stress (ESS) calculations using Newtonian and non-Newtonian rheological models were performed to assess whether the common assumption of Newtonian blood behavior offers similar results to a more realistic but computationally expensive non-Newtonian model. 16 coronary arteries (from 16 patients) were reconstructed from optical coherence tomographic (OCT) imaging. Pulsatile CFD simulations using Newtonian and the Quemada non-Newtonian model were performed. Endothelial shear stress (ESS) and other indices were compared. Exploratory indices including local blood viscosity (LBV) were calculated from non-Newtonian simulation data. Compared to the Newtonian results, the non-Newtonian model estimates significantly higher time-averaged ESS (1.69 (IQR 1.36)Pa versus 1.28 (1.16)Pa, p < 0.001) and ESS gradient (0.90 (1.20)Pa/mm versus 0.74 (1.03)Pa/mm, p < 0.001) throughout the cardiac cycle, under-estimating the low ESS (<1Pa) area (37.20 ± 13.57% versus 50.43 ± 14.16%, 95% CI 11.28–15.18, p < 0.001). Similar results were also found in the idealized artery simulations with non-Newtonian median ESS being higher than the Newtonian median ESS (healthy segments: 0.8238Pa versus 0.6618Pa, p < 0.001 proximal; 0.8179Pa versus 0.6610Pa, p < 0.001 distal; stenotic segments: 0.8196Pa versus 0.6611Pa, p < 0.001 proximal; 0.2546Pa versus 0.2245Pa, p < 0.001 distal) On average, the non-Newtonian model has a LBV of 1.45 times above the Newtonian model with an average peak LBV of 40-fold. Non-Newtonian blood model estimates higher quantitative ESS values than the Newtonian model. Incorporation of non-Newtonian blood behavior may improve the accuracy of ESS measurements. The non-Newtonian model also allows calculation of exploratory viscosity-based hemodynamic indices, such as local blood viscosity, which may offer additional information to detect underlying atherosclerosis.

Introduction

Fundamentally, computational fluid dynamics (CFD) is based on the idea that, given certain assumptions, the mechanics of fluid motion can be accurately described by physical principles and mathematical equations. The computational solution of these equations allows the determination of various hemodynamic indices such as blood velocity and pressure throughout the artery, from which other parameters such as endothelial shear stress (ESS) can be further derived. With continued advances, the underlying computational methods are frequently re-evaluated to optimize the shifting balance between accuracy and complexity. One such computational model concerns the variable viscosity of blood at high and low shear rates.

Earlier studies suggest that for laminar flow in medium to large arteries, blood may be assumed a Newtonian fluid with a constant viscosity independent of shear rate (1–5). However, due partly to its dual solid and liquid phases, blood exhibits non-Newtonian behaviors (6). This includes properties such as shear-thinning, the apparent thinning of blood at high shear rates and thickening at low shear rates. While the Newtonian assumption is generally acceptable in healthy straight segments of larger arteries, it may not be as accurate as non-Newtonian rheological models in the setting of complex flow patterns (7, 8). As artery anatomy changes, high fluctuations in local shear rate are significant enough that the non-Newtonian behavior of blood may emerge (Figure 1). In these near-wall regions, local blood viscosity (LBV) is expected to change from location to location and from instant to instant over the cardiac cycle. However, because the Newtonian model assumes constant viscosity, changes in LBV are not detected. The non-Newtonian model encompasses the variable viscosity of blood and thus provides this information. Our underlying hypothesis is that, if non-Newtonian behavior is negligible in coronary arteries, the two models should present nearly identical results. Divergent results, however, would suggest otherwise.

FIGURE 1

Figure 1. Normal and abnormal viscosity profiles. (A) In a straight unobstructed artery, blood velocity takes on a parabolic profile with low shear rate (high viscosity) in the centre of the artery and low viscosity at the wall. (B) In the presence of stenoses, curvatures, or bifurcations, the blood velocity profile becomes distorted. Localized regions of low shear rate can develop at the wall, resulting in high viscosity near the endothelial surface.

Materials and Methods

Study Design and Patient Selection

This study compares blood flow characteristics generated by CFD analysis in patient-specific coronary arteries using Newtonian and non-Newtonian blood models under pulsatile flow. Patients were retrospectively selected from a previous multicentre randomized clinical trial (NCT01776567). Inclusion criteria for the current study were the presence of an unstented, non-obstructive (diameter stenosis <50%) non-culprit lesion in the culprit vessel. Major exclusion criteria were ST-elevation myocardial infarction within the preceding 48 h, left ventricular ejection fraction <25%, and bifurcation lesions.

Three-Dimensional Reconstruction

Patient-specific 3D arterial models were reconstructed through the fusion of OCT and angiography as previously described (9). The arterial centreline was extracted from two end-diastolic angiographic images with a >25° difference in viewing angles (QAngio XA 3D, Medis Specials Bv, Netherlands). Side branches outside the region of interest were used to co-register OCT images with angiograms. OCT lumen contours were semi-automatically detected and manually corrected as necessary (QCU-CMS, Leiden University Medical Center, Netherlands). OCT contours were placed onto the angiographic centreline (MATLAB R2017b, Mathworks Inc., Natick, MA, United States) and the vessel surface was generated (Solidworks, Dassault Systèmes, Velizy, France). Vessel volume was discretised into tetrahedral elements with an average mesh size of approximately 1 to 2 million depending on the geometric complexity of individual patient-specific arteries (Pointwise v18.2R2). All discretised models included a graduated 10-prism boundary layer to further enhance the resolution of flow phenomena near the arterial wall.

Given the uniqueness and complexity of an individual patient’s coronary arteries, it is difficult to isolate the specific geometric features impacting non-Newtonian blood rheology. Therefore, an idealized 2D model of intermediate stenosis was virtually created to assess the generalizability of the patient-specific results. The 2D model had a 3 mm diameter with diameter stenosis (DS) of 40%. The stenosis anatomy was assumed by the following mathematical equation:

$D=\frac{1}{2}\left[D_{min}+\left(D-D_{min}\right){\left(\sin \left(\mathrm{\pi }x\right)\right)}^2\right]$

where, D is the artery’s diameter, D_min=DS×D is the minimum lumen diameter (MLD) and x is the longitudinal location of the stenotic segment with x = 0 at the MLD. The lesion length was assumed 3mm (i.e., x ranges from −1.5 to 1.5 mm). The idealized geometry was divided into 4 segments (Figure 2). Separating the proximal and distal stenotic geometries allows inspection of common flow phenomena – “favorable” and “adverse pressure gradient” – which are responsible for abnormal flow patterns such as flow separation and reversal.

FIGURE 2

Figure 2. An idealized stenotic model with an artery’s diameter of 3mm and 40% DS. Segment 1: at least 1 diameter from the proximal stenotic segment; Segment 2: 1 diameter from the proximal shoulder of the stenosis to minimal lumen diameter; Segment 3: minimal lumen diameter to 1 diameter distal to the stenosis shoulder; Segment 4: remaining distal segment.

Computational Fluid Dynamics Simulation

Computational fluid dynamics analysis was accomplished through the direct solution of the Navier-Stokes equations describing fluid motion. OpenFOAM, a finite-volume CFD solver, was run on the Magnus supercomputer, consisting of 35,712 Intel Xeon E5-2690V3 “Haswell” processors (Pawsey Supercomputing Centre, Perth, WA, Australia). A time-varying parabolic velocity profile with a mean bulk velocity calculated from the patient-specific TIMI frame count as previously described was applied at the inlet (10). The arterial wall was considered rigid with a no-slip boundary and a non-specific distal vascular resistance was applied at the outlet. A resistance boundary condition assumes linear dependence between the pressure and flow rate at each outlet. It is analogous to a constant pressure boundary condition and is suitable for unsteady simulation with only single outlet (11).

For the Newtonian simulations, blood’s constant dynamic viscosity (μ_Newtonian) was assumed 0.0035Pa⋅s. Blood density was considered 1,060kg/m³ and haematocrit 45%. Simulations were run for 3 cardiac cycles to ensure convergence of results, and all results presented are only from the final cycle. CFD results were post-processed to extract instantaneous and time-averaged hemodynamic indices where appropriate. ESS was calculated as described in Chen et al. (12). ESS gradient (ESSG) was calculated as the spatial gradient of ESS, representing the rate of change in ESS between adjacent spatial points in a local coordinate system (x′, y′, z′):

$\text{ESSG}=\nabla \text{ESSG}=\left[\begin{array}{ccc}\hfill \frac{\partial {\text{ESS}}_{x^{\prime }}}{\partial x^{\prime }}\hfill & \hfill \frac{\partial {\text{ESS}}_{x^{\prime }}}{\partial y^{\prime }}\hfill & \hfill \frac{\partial {\text{ESS}}_{x^{\prime }}}{\partial z}\hfill \\ \hfill \frac{\partial {\text{ESS}}_{y^{\prime }}}{\partial x^{\prime }}\hfill & \hfill \frac{\partial {\text{ESS}}_{y^{\prime }}}{\partial y^{\prime }}\hfill & \hfill \frac{\partial {\text{ESS}}_{y^{\prime }}}{\partial z^{\prime }}\hfill \\ \hfill \frac{\partial {\text{ESS}}_{z^{\prime }}}{\partial x^{\prime }}\hfill & \hfill \frac{\partial {\text{ESS}}_{z^{\prime }}}{\partial y^{\prime }}\hfill & \hfill \frac{\partial {\text{ESS}}_{z^{\prime }}}{\partial z^{\prime }}\hfill \end{array}\right]$

Since ESS represents the tangential force acting on the surface, all normal components (z′) of the tensor are irrelevant to the ESSG calculation. Removing all irrelevant tensor components in the z′ direction, the ESSG was simplified to:

$\text{ESSG}=\left[\begin{array}{cc}\hfill \frac{\partial {\text{ESS}}_{x^{\prime }}}{\partial x^{\prime }}\hfill & \hfill \frac{\partial {\text{ESS}}_{x^{\prime }}}{\partial y^{\prime }}\hfill \\ \hfill \frac{\partial {\text{ESS}}_{y^{\prime }}}{\partial x^{\prime }}\hfill & \hfill \frac{\partial {\text{ESS}}_{y^{\prime }}}{\partial y^{\prime }}\hfill \end{array}\right]$

and the ESSG magnitude is written as the diagonal components of the above matrix. Oscillatory shear index (OSI),

$\text{OSI}=0.5\times \left(1-\frac{\left|{\int }_0^t\overrightarrow{\text{ESS}}\mathit{d}t\right|}{{\int }_0^t\left|\overrightarrow{\text{ESS}}\right|\mathit{d}t}\right),$

indicates the degree of fluctuation in the direction of ESS vectors over the cardiac cycle. It is effectively an index of blood flow recirculation in a pulsating flow environment (13). Hemodynamic variables were calculated at 64 discrete timepoints per cardiac cycle (see Supplementary Table 1 for sensitivity analysis).

Similar procedures were employed in the non-Newtonian simulations except blood was modeled by the Quemada constitutive equation to capture the local variations in blood viscosity (μ_{non–Newtonian}) (14). Local blood viscosity (LBV) was expressed as a ratio of non-Newtonian viscosity to the Newtonian constant viscosity model (μ_{non–Newtonian} /μ_Newtonian) for each tetrahedral element (15). In other words, LBV is a quantitative measure of the impact of non-Newtonian flow relative to that assumed by the Newtonian model. A value of 1 indicates no effect and any values >1 demonstrate the proportional influence of non-Newtonian properties locally in the bloodstream.

Statistical Analysis

Two simulations were carried out for each patient, one using the Newtonian model and the second using a non-Newtonian model. Since the reconstructed arterial models and computational meshes were identical for the Newtonian and non-Newtonian simulations for each patient, a rigorous point-by-point comparison of all vessel nodes was possible. Due to the paired nature of observations, the only variable within each patient was the choice of rheological model. Thus, paired comparisons of the Newtonian and non-Newtonian results were performed, as described below.

Time-averaged ESS and ESSG were calculated for each case. Instantaneous ESS was evaluated at every spatial point on the arterial wall at every step in the cardiac cycle. The maximum ESS at that spatial point was identified. To eliminate the skewing effect of extremely high and low ESS values inherently present in complex geometries, ESS was normalized using the maximum ESS value at every spatial point, yielding a value between 0 and 1. At each time step, a point-by-point comparison of normalized ESS yielded the normalized difference between the Newtonian and non-Newtonian models. The same methods were used to analyze ESSG whereas, by definition, OSI describes the general flow behavior over a cardiac cycle. Similar analyses were also carried out on the four segments of the idealized artery geometry.

Categorical variables are presented as counts and percentages, while continuous variables are presented as a mean ± standard deviation and non-parametric variables in median (interquartile range [IQR]). Because of the paired nature of the simulations within each patient, the only variable was the choice of rheological model, therefore paired t-test or Wilcoxon sign rank test for paired observations (as appropriate) were used. To avoid statistical dependence and to decorrelate data in the point-by-point comparison, bootstrap resampling with replacement was used to randomly select 1.5% of all mesh points and the above-mentioned statistical tests were performed to compare the results of Newtonian and non-Newtonian simulations. This process was repeated 10,000 times for comparison of ESS, ESSG, and OSI. All tests were two-tailed with an α-level of 0.05 to indicate statistical significance. Statistical analysis was performed in R statistical software (R Foundation for Statistical Computing, Vienna, Austria).

Results

Wall-Based Indices: Endothelial Shear Stress, ESS Gradient, and Oscillatory Shear Index

Patient characteristics are shown in Table 1. Qualitative comparison of time-averaged ESS between the Newtonian and non-Newtonian simulations demonstrates broad similarity in its range and distribution for all cases (Figure 3). However, notable differences are found at the stenosis and curved segments (white arrows). By quantitative comparison, time-averaged ESS in the non-Newtonian simulations was significantly higher (1.69 [1.36]Pa versus 1.28 [1.16]Pa, p < 0.001), translating to a mean normalized percent difference of 21.72% over the entire cardiac cycle. Time-averaged ESSG in the non-Newtonian simulations was also significantly higher than the Newtonian model (1.65 ± 0.92Pa/mm versus 1.37 ± 0.78Pa/mm, 95% CI 0.20–0.37.16, p < 0.001). However, OSI was not significantly different between the models (0.0302 ± 0.035 versus 0.0294 ± 0.039, 95% CI 0.0059–0.0075, p = 0.81) (Table 2). The results of the analysis based on bootstrap resampling showed that the 95% confidence intervals of the mean and standard deviation/median and IQR are fully consistent with the analyses conducted on the raw data (Supplementary Table 2).

TABLE 1

Table 1. Patient characteristics.

FIGURE 3

Figure 3. Qualitative differences in ESS from (A) Newtonian and (B) non-Newtonian models. Although the distribution of high and low ESS is similar, the non-Newtonian model predicts higher ESS throughout the artery. This is visually most apparent from the larger areas exposed to high ESS (white arrows). ESS, endothelial shear stress.

TABLE 2

Table 2. ESS, ESSG, and OSI between rheological models.

The absolute and normalized percent difference in instantaneous ESS and ESSG relate inversely with coronary flow rate. At the higher coronary blood flow rates associated with diastole, the difference between Newtonian and non-Newtonian simulations approaches zero, with a minimum difference in absolute ESS of 0.035 ± 0.036Pa (Figure 4). However, the difference increases at low and decelerating flow rates that characterize systole, with a maximum difference in absolute ESS of 1.05 ± 0.42Pa (Figures 4B,C).

FIGURE 4

Figure 4. Quantitative difference in ESS between Newtonian and non-Newtonian models (Single Representative Case). (A) Coronary flow rate, systole is indicated in black, diastole is indicated in green. (B,C) The % difference and absolute difference are calculated as (non-Newtonian ESS – Newtonian ESS), thus positive values indicate non-Newtonian results were higher whereas negative values indicate higher Newtonian results. The non-Newtonian results show consistently higher percent normalized difference and absolute difference in ESS over the cardiac cycle, except during the momentary transition between systole and diastole at approximately 0.35s. (D) Newtonian simulations predict more of the vessel is exposed to atherogenic levels of ESS over the cardiac cycle.

Although the time-averaged results show that the non-Newtonian simulations estimate higher ESS and ESSG over the cardiac cycle, during the momentary transition between end-systole and early diastole, at approximately 0.35s into the cardiac cycle, the Newtonian ESS results are higher than the non-Newtonian results by 0.89 ± 0.52Pa (p < 0.001), or 946.97 ± 898.72% (p < 0.001). However, this is transient, and the non-Newtonian results again become higher immediately thereafter (Figures 4B,C). The non-Newtonian model predicts the highest ESS relative to the Newtonian model during early diastole, as coronary flow is increasing rapidly. However, as the rate of rise slows in mid-diastole (at approximately 0.4s), non-Newtonian ESS is higher by 0.96 ± 0.25Pa (p < 0.001), or 58.13 ± 22.12% (p < 0.001).

The implications of these differences are most apparent by comparing vessel areas predicted to be exposed to low ESS (<1Pa), a generally accepted threshold for stimulating pro-atherogenic processes (16, 17). The Newtonian model predicts significantly greater vessel area exposure to low ESS than the non-Newtonian model (50.43 ± 14.16% versus 37.20 ± 13.57%, 95% CI 11.28–15.18%, p < 0.001) (Figure 4D).

Results from the idealized arterial geometries are consistent with the 3D patient-specific results (Table 3). ESS was higher in the non-Newtonian simulation in both proximal and distal healthy segments (0.8238Pa versus 0.6618Pa, p < 0.001 and 0.8179Pa versus 0.6610Pa, p < 0.001, respectively); in proximal and distal stenotic segments (0.8196Pa versus 0.6611Pa, p < 0.001 and 0.2546Pa versus 0.2245Pa, p < 0.001, respectively). Both non-Newtonian and Newtonian simulations display similar IQR ESS, except for the distal stenotic segment where non-Newtonian simulation demonstrates a significantly higher IQR than Newtonian (0.2199Pa versus 0.08843Pa, p < 0.001). IQR is a measure of the data dispersion (18). In other words, IQR indicates the range of ESS within the region of interest. High IQR indicates wider spread of the ESS values and hence larger variation of ESS from the median, a condition that can be found with increasing chaotic blood flow due to vortices and flow oscillations.

TABLE 3

Table 3. ESS distribution between rheology models in an idealized artery (as depicted in Figure 1).

Local Blood Viscosity

Like ESS, there is high spatial and temporal heterogeneity in blood viscosity over the cardiac cycle. Localized volumetric regions of high blood viscosity are observed in every case, including at the centre and walls of the artery. Across all cases, time-averaged viscosity was 1.45-fold higher than that assumed by the Newtonian model (95% CI 1.43–1.49, p < 0.001) (Figure 5). Some vessel regions are marked by an average 41.5-fold increase in maximum viscosity compared to the Newtonian model (95% CI 30.1–53.0, p < 0.001). In one case, peak viscosity was more than 70 times higher. The peak viscosity invariably occurs during a nadir in coronary flow – either during peak systole or the transition between end systole and early diastole.

FIGURE 5

Figure 5. ESS and local blood viscosity. (A) Multiple areas of low ESS (<1Pa) are identified (white arrows). (B) Non-Newtonian simulation detects regions of high LBV (>1.45) in the centre of the vessel and at the wall (orange arrows). ESS, endothelial shear stress; LBV, local blood viscosity.

Discussion

This study demonstrates: (1) the non-Newtonian model predicts significantly higher ESS and ESSG than the Newtonian model; (2) the Newtonian model shows significantly greater vessel areas exposed to atherogenic low ESS; (3) OSI is not significantly different between the models; and (4) the non-Newtonian model identifies regions of high LBV up to 70-fold higher than that assumed by the Newtonian model.

Local blood flow disturbances are a potent stimulus of endothelial dysfunction and biological processes underlying atherosclerosis (16, 19, 20). ESS disturbances in vivo have been correlated with changes in plaque composition (21, 22), morphology (23), and vessel remodeling (21, 22). Further, low ESS is independently associated with requiring future intervention or even future clinical events (24, 25). In these studies, many vessel areas were predicted to be exposed to low ESS; of lesions that caused future events, most or all were in previously low ESS areas; however, the vast majority of lesions within low ESS areas did not progress to cause events. Previous studies indicate that, while low ESS has high sensitivity to detect future events at a patient level, it has low specificity and positive predictive value (PPV) (24, 25). Furthermore, no studies have associated specific plaques exposed to low ESS as the culprit lesion of later events.

Although blood is a non-Newtonian fluid, all these studies used the Newtonian model of blood behavior. Since coronary blood flow is assumed to have shear rates above 100–200s^–1, such an assumption has long been valid. We investigated the effects of non-Newtonian blood rheology in a series of coronary arteries reconstructed from high-resolution imaging to test whether the Newtonian and non-Newtonian models present results within a small margin of computational error. In other words, if the shear rate was indeed high enough that blood can be assumed a Newtonian fluid, the two models should predict quantitatively similar results. However, the results of this study show that the Newtonian and non-Newtonian models consistently estimate different ESS values throughout the cardiac cycle, suggesting that there are indeed multiple factors, including the instantaneous pulsatile blood flow velocity and local geometric variations, influencing when and where the non-Newtonian behaviors of blood become apparent.

Pulsatile Flow Factor

Unlike Newtonian CFD simulations where the blood viscosity stays constant throughout the cardiac cycle, blood viscosity constantly changes from moment to moment due to the pulsatile nature of blood flow. Indeed, in a Newtonian simulation, the sudden increase in coronary blood flow rate at early diastole is usually accompanied by a very rapid increase in ESS untamed by the constant viscous forces and vice versa. However, the non-Newtonian fluid responds differently to temporal changes in flow rate and hence leads to a very dynamic pattern of the discrepancy between the Newtonian and non-Newtonian results. While non-Newtonian simulations, in general, predict a higher ESS, there is a time point between end-systole and early diastole where the non-Newtonian model momentarily predicts significantly lower ESS than the Newtonian model. Coincidentally, the difference between the two models is most remarkable during this transition phase. While it can be easily speculated that the rapidly changing blood flow conditions play a role, the underlying haemodynamics are far more complicated. One must consider the heart rate, the historical effects of LBV, and its impact on the local flow environment.

Nevertheless, these findings are consistent with many previous fundamental studies showing that non-Newtonian simulations predict higher ESS (26–31). For instance, non-Newtonian flow decreases the area of low ESS in both straight and bent arterial segments, with the largest difference occurring in the straight rather than the bent segment (32). On the other hand, while blood viscosity affects the magnitude of ESS when the flow is disturbed, it does not affect the spatial and temporal distribution of the ESS (33, 34). Our results demonstrate that those results extend to patient-specific coronary arteries. It should be noted that while this difference was also observed for ESSG, there was no significant difference in predicted OSI values.

Reassuringly, ESS calculated by the non-Newtonian model is, on average, 0.44Pa or 21% higher than the Newtonian model in the 3D patient-specific data and 0.63Pa or 27% in the idealized 2D results. While the absolute value of the difference between Newtonian and non-Newtonian models is low, this ultimately means that the Newtonian model estimates a significantly higher percentage of the vessel area exposed to ESS <1Pa (50.43% versus 37.20%, p < 0.001). This could have significant repercussions in the context of earlier clinical CFD studies showing low specificity and PPV of low ESS. Despite the estimated differences in ESS and ESSG, the clinical significance of such fluctuations within a short time is unclear. However, it is hypothesized that the non-Newtonian model, as a more accurate reflection of actual blood behavior, may ultimately offer higher specificity and PPV than Newtonian simulations.

Geometric Factor

Arterial narrowing and widening are major factors in local variations of blood rheology and flow dynamics. Blood flow accelerates as the artery narrows. In haemodynamics, this segment is termed the “favorable pressure gradient” segment. In contrast, “adverse pressure gradient” refers to flow deceleration as the artery widens. The impact of “favorable” and “adverse pressure gradients” can be isolated by carrying out the Newtonian and non-Newtonian simulations in an idealized arterial geometry where the arterial flow will undergo clearly defined favorable (proximal stenotic segment) and adverse (distal stenotic segment) pressure gradients. Our idealized artery results show that the non-Newtonian model continues to display a higher median ESS value compared to the Newtonian model in all segments. However, it is the significantly lower IQR at the distal stenotic segment marked by an adverse pressure gradient in the Newtonian model that is of particular interest. While low IQR reflects lower ESS oscillation and high IQR indicates larger fluctuations in ESS values, it is unclear whether a higher IQR signifies an increase in turbulent activities or vice versa. Nonetheless, regions with adverse pressure gradients have demonstrated an increased likelihood of flow reversal (8). Clinically, flow reversal and abnormal ESS represent a location for plaque development and potentially a nidus for thrombotic events (35–37). This remarkable difference in Newtonian and non-Newtonian models in regions of adverse pressure gradients might lead to different conclusions and hence warrants further objective analyses.

The Need for a Better Reflection of True Blood Rheology

Non-Newtonian models offer other potential advantages. Although ESS, ESSG, OSI, and other wall-based metrics consider the mechanical effect of blood acting on the vessel wall, these indices inherently neglect the physiological response of whole blood. While low ESS and high OSI indicate that blood may be recirculating and stagnating in these areas, these measures do not directly capture or describe flow phenomena within the blood itself.

Compared to the constant viscosity assumed by the Newtonian model, in non-Newtonian rheological models, blood viscosity is treated as a variable dependent on instantaneous local shear rate, allowing determination of viscosity within the entire fluid domain. In this study, the non-Newtonian model identified localized regions of peak LBV, on average, 40-fold higher than that assumed by the Newtonian model. The possibility of detecting localized regions of increased blood viscosity in vivo is intriguing given that blood is a complex fluid with clinically relevant behaviors such as thrombosis. Further, the pathophysiologic mechanisms underlying plaque development may involve the accumulation of cholesterol, pro-inflammatory cells, and humoral mediators in characteristic vessel regions, perhaps exacerbated by increased LBV, recirculation and stagnation in these areas.

In straight unobstructed vessels, high viscosity is expected in the center of a vessel where shear rate is low, and velocity is high (Figure 1). Conversely, blood viscosity at the wall is low since the shear rate at the wall is high. It is hypothesized that, despite high viscosity at the vessel centre, high velocity convects blood axially downstream, preventing significant erythrocyte aggregation or contact with the endothelial surface. However, low blood velocity and recirculation can develop at the distal inner bend of curvatures, the outer walls of bifurcations, and both proximal and distal to stenoses or stent struts. In such regions of disturbed flow, blood velocity and shear rate may decrease at the wall leading to pockets of high LBV near the endothelial surface, potentially facilitating processes leading to both progressive atherosclerosis and thrombosis (38, 39).

Limitations

There are several limitations to the current study. First is a small study population retrospectively selected from a prior randomized clinical trial. However, in the context of patient-specific CFD studies quantitatively comparing rheological models, this is among the largest cohorts. As a result of retrospective patient selection, the cohort also skewed male due to the original study population characteristics (75% male in the original study). Second, we assumed a generic haematocrit of 45% based on standard reference ranges and a desire to limit confounding variables that might have been introduced by incorporating patient-specific values. However, because haematocrit is a determinant of blood viscosity, it is possible that the observed differences between the Newtonian and non-Newtonian models may also be influenced by changes in haematocrit. Future studies should incorporate patient-specific haematocrit to better assess this possibility. Third, this study does not investigate the effect of axial and secondary flow due to the presence of helical inflow, arterial curvature, bifurcation lesions which significantly affect local flow dynamics (40). It is expected that due to the increased complexity of flow in these settings, the differences between the Newtonian and non-Newtonian models would increase further. To facilitate investigating these effects, it is important to correlate LBV with other ESS-based descriptors and helicity indices (41, 42). In terms of bifurcation lesions, future studies should prospectively image-side branches with intravascular techniques if it is feasible and safe to do so. Fourth, as a technical study, changes in plaque composition were not evaluated. Further investigation of ESS derived by Newtonian and non-Newtonian models concerning plaque composition and change over time is necessary. Ultimately, if there is a significant difference in how the models correlate with atherosclerotic plaque, future studies investigating clinical endpoints may be worthwhile.

Conclusion

Although blood is often assumed to be a Newtonian fluid by CFD simulations of the coronary arteries, this study demonstrates that non-Newtonian behaviors of blood are operational, yielding marked differences in calculated flow indices such as ESS and ESSG. Non-Newtonian simulations also allow the calculation of LBV and related indices, potentially presenting novel markers to detect plaques at risk for progression.

Data Availability Statement

The CFD data supporting the conclusions of this article will be made available by the authors. Request for access to patient data (which are subjected to restriction due to the nature of personal healthcare information) can be made to the corresponding author.

Ethics Statement

The studies involving human participants were reviewed and approved by St Vincent’s Ethics Committees (APPOSE HREC/12/SVH/31). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

VT designed the methodology, performed 3D reconstructions and CFD analysis and wrote the original manuscript, and reviewed and edited the manuscript. JD supported intravascular optical coherence tomography analysis. SZ contributed to the CFD modeling. DS, ER, and WG reviewed and edited the manuscript. PS provided expert opinion, review and edited the manuscript. EP co-designed the methodology and performed CFD analysis, co-supervised the project. AO contributed to CFD methods and co-supervised the project. PB collected all medical imaging data, conceptualized the project, supervised, and co-ordinated the project. All authors contributed to the article and approved the submitted version.

Funding

This research was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government.

Conflict of Interest

PS reports personal fees from Sino Medical Sciences Technology, Philips/Volcano, and Xeltis.

The remaining 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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcvm.2022.835270/full#supplementary-material

References

1. Ballyk PD, Steinman DA, Ethier CR. Simulation of non-Newtonian blood flow in an end-to-side anastomosis. Biorheology. (1994) 31:565–86. doi: 10.3233/bir-1994-31505

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Kleinstrauer C, Nazemi M, Archie JP. Hemodynamics analysis of a stenosed carotid bifurcation and its plaque-mitigating design. J Biomech Eng. (1991) 113:330–5. doi: 10.1115/1.2894891

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Morbiducci U, Gallo D, Massai D, Ponzini R, Deriu MA, Antiga L, et al. On the importance of blood rheology for bulk flow in hemodynamic models of the carotid bifurcation. J Biomech. (2011) 44:2427–38. doi: 10.1016/j.jbiomech.2011.06.028

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Rindt CCM, van de Vosse FN, Steenhoven AAV, Janssen JD. A numerical and experimental analysis of the flow field in a two-dimensional model of the human carotid artery bifurcation. J Biomech. (1987) 20:499–509. doi: 10.1016/0021-9290(87)90250-8

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Cho Y, Kensey K. Effects of the non-Newtonian viscosity of blood on flows in a diseased arterial vessel. Part 1. Steady flows. Biorheology. (1991) 28:241–62. doi: 10.3233/bir-1991-283-415

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Merrill EW, Pelletier GA. Viscosity of human blood: transition from Newtonian to non-Newtonian. J Appl Physiol. (1985) 1967:178–82. doi: 10.1152/jappl.1967.23.2.178

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Bodnár T, Sequeira A, Prosi M. On the shear-thinning and viscoelastic effects of blood flow under various flow rates. Appl Math Comput. (2011) 217:5055–67. doi: 10.1007/s10856-010-4210-6

PubMed Abstract | CrossRef Full Text | Google Scholar

8. van Wyk S, Prahl Wittberg L, Bulusu KV, Fuchs L, Plesniak MW. Non-Newtonian perspectives on pulsatile blood-analog flows in a 180^° curved artery model. Phys Fluids. (2015) 27:071901. doi: 10.1063/1.4923311

CrossRef Full Text | Google Scholar

9. Papafaklis MI, Bourantas CV, Yonetsu T, Vergallo R, Kotsia A, Nakatani S, et al. Anatomically correct three-dimensional coronary artery reconstruction using frequency domain optical coherence tomographic and angiographic data: head-to-head comparison with intravascular ultrasound for endothelial shear stress assessment in humans. EuroIntervention. (2015) 11:407–15. doi: 10.4244/EIJY14M06_11

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Gibson CM, Cannon CP, Daley WL, Dodge JT, Alexander B, Marble SJ, et al. TIMI Frame count : a quantitative method of assessing coronary artery flow. Circulation. (1996) 93:879–88. doi: 10.1161/01.cir.93.5.879

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Grinberg L, Karniadakis GE. Outflow boundary conditions for arterial networks with multiple outlets. Ann Biomed Eng. (2008) 36:1496–514. doi: 10.1007/s10439-008-9527-7

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Chen WX, Poon EKW, Hutchins N, Thondapu V, Barlis P, Ooi A. Computational fluid dynamics study of common stent models inside idealised curved coronary arteries. Comput Methods Biomech Biomed Eng. (2017) 20:671–81. doi: 10.1080/10255842.2017.1289374

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Ku DN, Giddens DP, Zarins CK, Glagov S. Pulsatile flow and atherosclerosis in the human carotid bifurcation. Positive correlation between plaque location and low oscillating shear stress. Arteriosclerosis Thrombosis Vasc Biol. (1985) 5:293–302. doi: 10.1161/01.atv.5.3.293

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Quemada D. Rheology of concentrated disperse systems II. A model for non-newtonian shear viscosity in steady flows. Rheol Acta. (1978) 17:632–42. doi: 10.1007/bf01522036

CrossRef Full Text | Google Scholar

15. Poon EKW, Thondapu V, Hayat U, Barlis P, Yap CY, Kuo PH, et al. Elevated blood viscosity and micro-recirculation resulting from coronary stent malapposition. J Biomech Eng. (2018) 140:051006.

Google Scholar

16. Davies PF, Remuzzi A, Gordon EJ, Dewey CFJ, Gimbrone MA. Turbulent fluid shear stress induces vascular endothelial turnover in vitro. Proc Natl Acad Sci U S A. (1986) 83:2114–7. doi: 10.1073/pnas.83.7.2114

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Malek AM, Alper SL, Izumo S. Hemodynamic shear stress and its role in atherosclerosis. JAMA. (1999) 282:2035–42. doi: 10.1001/jama.282.21.2035

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Manikandan S. Measures of dispersion. J Pharmacol Pharmacother. (2011) 2:315–6.

Google Scholar

19. Dolan JM, Meng H, Singh S, Paluch R, Kolega J. High fluid shear stress and spatial shear stress gradients affect endothelial proliferation, survival, and alignment. Ann Biomed Eng. (2011) 39:1620–31. doi: 10.1007/s10439-011-0267-8

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Gimbrone MA, Topper JN, Nagel T, Anderson KR, Garcia-Cardena G. Endothelial dysfunction, hemodynamic forces, and atherogenesis. Ann NY Acad Sci. (2000) 902:230–9. doi: 10.1111/j.1749-6632.2000.tb06318.x

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Corban MT, Eshtehardi P, Suo J, McDaniel MC, Timmins LH, Rassoul-Arzrumly E, et al. Combination of plaque burden, wall shear stress, and plaque phenotype has incremental value for prediction of coronary atherosclerotic plaque progression and vulnerability. Atherosclerosis. (2014) 232:271–6. doi: 10.1016/j.atherosclerosis.2013.11.049

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Samady H, Eshtehardi P, McDaniel MC, Suo J, Dhawan SS, Maynard C, et al. Coronary artery wall shear stress is associated with progression and transformation of atherosclerotic plaque and arterial remodeling in patients with coronary artery disease. Circulation. (2011) 124:779–88. doi: 10.1161/CIRCULATIONAHA.111.021824

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Yamamoto E, Siasos G, Zaromytidou M, Coskun AU, Xing L, Bryniarski K, et al. Low endothelial shear stress predicts evolution to high-risk coronary plaque phenotype in the future: a serial optical coherence tomography and computational fluid dynamics study. Circ Cardiovasc Interv. (2017) 10:e005455. doi: 10.1161/CIRCINTERVENTIONS.117.005455

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Stone PH, Maehara A, Coskun AU, Maynard CC, Zaromytidou M, Siasos G, et al. Role of low endothelial shear stress and plaque characteristics in the prediction of nonculprit major adverse cardiac events: the Prospect study. JACC Cardiovasc Imaging. (2018) 11:462–71. doi: 10.1016/j.jcmg.2017.01.031

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Stone PH, Saito S, Takahashi S, Makita Y, Nakamura S, Kawasaki T, et al. Prediction of progression of coronary artery disease and clinical outcomes using vascular profiling of endothelial shear stress and arterial plaque characteristics: the Prediction study. Circulation. (2012) 126:172–81. doi: 10.1161/circulationaha.112.096438

CrossRef Full Text | Google Scholar

26. Apostolidis AJ, Moyer AP, Beris AN. Non-Newtonian effects in simulations of coronary arterial blood flow. J Non-Newtonian Fluid Mech. (2016) 233:155–65. doi: 10.1016/j.jnnfm.2016.03.008

CrossRef Full Text | Google Scholar

27. Chen J, Lu XY. Numerical investigation of the non-Newtonian pulsatile blood flow in a bifurcation model with a non-planar branch. J Biomech. (2006) 39:818–32. doi: 10.1016/j.jbiomech.2005.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Choi HW, Barakat AI. Numerical study of the impact of non-Newtonian blood behavior on flow over a two-dimensional backward facing step. Biorheology. (2005) 42:493–509.

PubMed Abstract | Google Scholar

29. Hund S, Kameneva M, Antaki JA. Quasi-mechanistic mathematical representation for blood viscosity. Fluids. (2017) 2:10. doi: 10.3390/fluids2010010

CrossRef Full Text | Google Scholar

30. Janela J, Moura A, Sequeira A. A 3D non-Newtonian fluid–structure interaction model for blood flow in arteries. J Comput Appl Math. (2010) 234:2783–91. doi: 10.1016/j.cam.2010.01.032

CrossRef Full Text | Google Scholar

31. Soulis JV, Giannoglou GD, Chatzizisis YS, Seralidou KV, Parcharidis GE, Louridas GE. Non-Newtonian models for molecular viscosity and wall shear stress in a 3D reconstructed human left coronary artery. Med Eng Phys. (2008) 30:9–19. doi: 10.1016/j.medengphy.2007.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Ferrarini A, Finotello A, Salsano G, Auricchio F, Palombo D, Spinella G, et al. Impact of leg bending in the patient-specific computational fluid dynamics of popliteal stenting. Acta Mech Sin. (2021) 37:279–91. doi: 10.1007/s10409-021-01066-2

CrossRef Full Text | Google Scholar

33. Liu B, Tang D. Influence of non-Newtonian properties of blood on the wall shear stress in human atherosclerotic right coronary arteries. Mol Cell Biomech. (2011) 8:73–90. doi: 10.3970/mcb.2011.008.073

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Anastasios S, Panagiotis N, Christos H. Comparison of blood rheological models in patient specific cardiovascular system simulations. J Hydrodyn Ser B. (2017) 29:293–304. doi: 10.1016/s1001-6058(16)60739-4

CrossRef Full Text | Google Scholar

35. Tenekecioglu E, Poon EK, Collet C, Thondapu V, Torii R, Bourantas CV, et al. The nidus for possible thrombus formation: insight from the microenvironment of bioresorbable vascular scaffold. JACC Cardiovasc Intervent. (2016) 9:2167–8. doi: 10.1016/j.jcin.2016.08.019

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Thondapu V, Mamon C, Poon EKW, Kurihara O, Kim HO, Russo M, et al. High spatial endothelial shear stress gradient independently predicts site of acute coronary plaque rupture and erosion. Cardiovasc Res. (2021) 117:1974–85. doi: 10.1093/cvr/cvaa251

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Yamamoto E, Thondapu V, Poon E, Sugiyama T, Fracassi F, Dijkstra J, et al. Endothelial shear stress and plaque erosion: a computational fluid dynamics and optical coherence tomography study. JACC Cardiovasc Imaging. (2019) 12:374–5. doi: 10.1016/j.jcmg.2018.07.024

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Dintenfass L, Julian DG, Miller GE. Viscosity of blood in normal subjects and in patients suffering from coronary occlusion and arterial thrombosis. Am Heart J. (1965) 71:587–600. doi: 10.1016/0002-8703(66)90308-5

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Wells R, Merrill EW. Shear rate dependence of the viscosity of whole blood and plasma. Science. (1961) 133:763–4. doi: 10.1126/science.133.3455.763

CrossRef Full Text | Google Scholar

40. Li Y, Gutierrez-Chico JL, Holm NR, Yang W, Hebsgaard L, Christiansen EH, et al. Impact of side branch modeling on computation of endothelial shear stress in coronary artery disease: coronary tree reconstruction by Fusion of 3D angiography and OCT. J Am Coll Cardiol. (2015) 66:125–35. doi: 10.1016/j.jacc.2015.05.008

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Calo K, De Nisco G, Gallo D, Chiastra C, Hoogendoorn A, Steinman DA, et al. Exploring wall shear stress spatiotemporal heterogeneity in coronary arteries combining correlation-based analysis and complex networks with computational hemodynamics. Proc Inst Mech Eng H. (2020) 234:1209–22. doi: 10.1177/0954411920923253

PubMed Abstract | CrossRef Full Text | Google Scholar

42. De Nisco G, Kok AM, Chiastra C, Gallo D, Hoogendoorn A, Migliavacca F, et al. The atheroprotective nature of helical flow in coronary arteries. Ann Biomed Eng. (2019) 47:425–38. doi: 10.1007/s10439-018-02169-x

PubMed Abstract | CrossRef Full Text | Google Scholar