Influence of right coronary artery motion, flow pulsatility and non-Newtonian rheology on wall shear stress metrics

Introduction

The highly non-uniform distribution of atherosclerosis within the arterial tree (Mitchell and Schwartz, 1965; Grøttum et al., 1983; Friedman et al., 1993) implies the existence of powerful local risk factors. Excessive uptake of plasma macromolecules by the arterial wall is one such factor; it shows excellent spatial correlation with rabbit and human lesions (Sinzinger et al., 1980; Sebkhi and Weinberg, 1994; Barnes and Weinberg, 1998; Weinberg, 2004). Another local factor – haemodynamically induced mechanical stress on endothelial cells – is likely to trigger the elevated uptake but the nature of the forces involved has been a matter of debate for many years.

The low, oscillatory wall shear stress (WSS) theory postulates that disease is triggered where haemodynamic WSS is low on average (Caro et al., 1971) and the Oscillatory Shear Index (OSI) is high (Ku et al., 1985). The OSI captures near-wall flow that deviates from the mean flow direction during the cardiac cycle (Ku et al., 1985). A more recent metric of multidirectional flow is the transverse WSS (transWSS), which captures components of instantaneous flow vectors that are perpendicular to the mean flow vector (Peiffer et al., 2013). Unlike the OSI, it does not include simple reverse flow. Its potential importance is predicated on the assumptions that endothelial cells elongate and align close to the mean flow direction and are adversely affected by flow across their long axis. TransWSS correlates better than low WSS or the OSI with patterns of permeability and lesions around side branches of the aorta (Mohamied et al., 2015). Furthermore, in vitro studies have shown that transverse flow leads to endothelial dysfunction (Wang et al., 2013; Ghim et al., 2018) and enhanced wall permeability (Ghim et al., 2017; Ghim et al., 2021), and that endothelial cells have evolved mechanisms which minimise it (Arshad et al., 2021).

Coronary artery disease has greater clinical importance than disease of the aorta. Conventional WSS metrics have been computed in static coronary models for over 20 years (e.g., Gibson et al., 1993; He and Ku, 1996; Perktold et al., 1998). TransWSS has also been investigated in such models: recent prospective studies have not demonstrated a convincing difference in plaque progression between regions experiencing low, intermediate or high transWSS (Kok et al., 2019; Hoogendoorn et al., 2020). However, the coronaries have a more dynamic mechanical environment than the aorta: there is greater bending, twisting and translation during the cardiac cycle. Only a few studies of coronary haemodynamics have included the effect of arterial motion.

Using a simple curved-tube model with a cyclically changing radius of curvature, Moore et al. (2001) found that WSS magnitudes could fluctuate substantially at individual locations in the tube but that flow patterns were broadly similar to those observed in stationary tubes and that the mean value was also predicted well by the stationary model. Zeng et al. (2003) imposed physiologically realistic motion on a constant-diameter model with an anatomically correct centre-line geometry and reached a similar conclusion: motion caused local fluctuations in WSS but the average value was similar to the static case. The fluctuations were in any case smaller than those produced by a pulsatile inflow boundary condition. Ramaswamy et al. (2004) incorporated cyclical diameter changes as well as motion of the vessel and their anatomically accurate geometry contained a stenosis. A notable finding was a substantial effect of motion on the pattern of OSI and on its average magnitude, which increased. Consistent with this, Torii et al. (2010) found that motion had a greater effect on OSI than on TAWSS in a dynamically and anatomically accurate model of a healthy right coronary artery (RCA).

The effect of arterial motion on transWSS has not previously been examined. The Dean-like secondary flow features that largely determine transWSS (Mohamied et al., 2017) are particularly sensitive to the changes in local curvature produced by coronary bending (Myers et al., 2001). The present study examines the effects of cyclical motion on WSS metrics including the transWSS in an anatomically realistic model of the normal human RCA. It also examines influences of blood rheology. Three previous studies of flow in coronary arteries found little difference in WSS patterns when non-Newtonian rheology was incorporated into the model (Johnston et al., 2004; Johnston et al., 2006; Soulis et al., 2006), but Apostolidis et al. (2016) showed differences as high as 50% at the left coronary bifurcation. A study by Cherry and Eaton showed that errors arising from the assumption of Newtonian rheology were particularly high in the presence of secondary flows (Cherry and Eaton, 2013). Indeed, non-Newtonian rheology can lead to the appearance of additional vortices in models of curved arteries (van Wyk et al., 2015). Hence shear-thinning properties might be particularly critical for transWSS. The present study examines the effect of non-Newtonian rheology on patterns of transWSS under coronary artery motion.

Methods

Data acquisition using magnetic resonance imaging

The dynamic geometry used in this analysis was derived from the 4-D MRI data of Torii et al. (2010). Fluid-structure interaction (FSI) techniques, used for other vessels, are not appropriate because most of the deformation in the coronary arteries is caused by cardiac motion and not by the fluid within the artery itself. An additional advantage of using MRI is that flow waveforms can be obtained at the same time, with the same imaging modality.

Briefly, data for the right coronary artery (RCA) of a healthy 31-year-old male were acquired using an ECG- and navigator-gated interleaved spiral sequence on a Siemens Sonata 1.5 T scanner. Written consent was obtained from the volunteer, and the protocol was approved by the local ethics committee; the study complied with the Declaration of Helsinki. A total of 20 interleaves were required to fill the k-space, with each lasting 10 ms. Imaging was performed during free-breathing and a water excitation pulse was incorporated to eliminate signal from fat. The spatial resolution of the acquired images was 0.9 × 0.9 mm; the images were later reconstructed to a finer resolution of 0.45 × 0.45 mm. Approximately 15–21 cross sections were obtained along the length of the RCA at 14 different time points. Segmented lumen contours from the cross sections were interpolated using cubic splines for 3-D surface reconstruction. The cross-sectional shape of the reconstructed geometry was assumed to be constant throughout the entire cardiac cycle as effects of radial changes are small (Torii et al., 2010).

An ECG- and navigator-gated interleaved spiral phase velocity mapping sequence (Keegan et al., 2004) was used to obtain velocity data in the proximal region of the RCA in the same volunteer. Figure 1A shows how the measured velocity varied throughout the cardiac cycle. The phases of geometry and velocity data were aligned using the timing of the R-wave (i.e. the beginning of the cycle), time from the R-wave and the RR interval as references.

FIGURE 1

FIGURE 1. (A) Velocity waveform for the RCA. Data were obtained from 0.012 s onwards. Earlier values were obtained by extrapolation, using a spline. (B) Flowchart showing how the moving wall boundary condition was generated for flow simulations.

Computational framework

Numerical methods were used to model blood flow in the coronary artery. A finite volume approach implemented in STAR-CCM+ (Siemens Industry Software Ltd., UK) approximated the solution of the continuity and Navier-Stokes equations. More specifically, a segregated flow solver was used to solve the equations in an uncoupled manner using a predictor-corrector scheme. A second order upwind scheme was employed for spatial discretisation and a first order Euler implicit scheme for temporal discretisation. Details of how the transient, convective, diffusive and pressure terms in the Navier-Stokes equation are discretised can be found in the STAR-CCM user guide. Polyhedral meshing was used because of its shorter computational time and greater wall shear stress accuracy (Spiegel et al., 2011); eight prism layers were added to capture the effect of the boundary layer near the wall. The flow was assumed to be laminar, and blood was assumed to be incompressible (with density 1,044 kg/m³) and homogeneous.

The flow waveform shown in Figure 1A and a flat velocity profile (i.e. uniform through-plane velocity) were prescribed at the inlet, and a zero relative pressure condition at the outlet. When the inlet and outlet were extended by 10 vessel diameters in trial computations, the results did not change substantially: the Spearman rank correlation coefficient for transverse wall shear stress between the case with and without extensions was 0.98. However, the morphing method (described below) did not work well near the proximal root and required a slight inlet extension to ensure smooth deformation. Hence, all simulations were performed with inlet extensions having a length equal to the vessel diameter.

The boundary condition for the moving wall was incorporated using a combination of the Vascular Modelling Toolkit (VMTK), Meshmixer, Blender and the mesh morphing solver in STAR-CCM+. Figure 1B summarises the method used to generate wall displacements. First, the centrelines were computed for the RCA geometries acquired at 14 timepoints, using VMTK. Interpolation using cubic splines gave centrelines at intermediate time intervals. Surface deformation in Blender was implemented by using centrelines (50 control points) of geometries as source and target points, and involved three steps: (i) rigging, (ii) binding, and (iii) skinning. During the rigging phase, the armature or bone structure of the 3-D RCA geometry was created using its centreline. Bones were then bound to the mesh using the bone heat algorithm (Baran and Popović, 2007) in which each mesh vertex is assigned a weight based on its distance from a particular bone. The weight controls how much influence the movement of the bone has on the vertex. The RCA was deformed during the skinning phase using linear blend skinning, also known as skeleton subspace deformation (Lewis et al., 2000). Finally, incremental wall displacements were computed for each time step; when combined with the morphing algorithm in STAR-CCM+, this allowed non-rigid deformations of the mesh by redistributing mesh vertices. Deformation of an exemplar geometry throughout the cardiac cycle is shown in Supplementary Figure S1, and values of curvature and torsion along the RCA are given in Supplementary Figure S2. A comparison of the original images and derived geometries in ParaView showed only negligible differences (data not shown).

To validate the methods, flow was computed in a 180-degree curved tube with sinusoidally varying curvature and steady inlet velocity, as modelled by Santamarina et al. (1998). Sensitivity tests were carried out to determine the optimal time step, mesh resolution and number of cardiac cycles.

Wall shear stress metrics

Wall shear stress vectors from the simulation were used to compute the four metrics defined in Equations (1)–(4). Time average WSS (TAWSS) is a magnitude-based metric and not sensitive to flow direction. OSI (the equation in (He and Ku, 1996) was used) and transWSS (Peiffer et al., 2013) have been described above. The Cross Flow Index (CFI) is similar to transWSS in that it represents the directionality of WSS, but it does not take magnitude into account. Mathematically it is the average over the cardiac cycle of the sine of the angle between the temporal mean WSS vector and instantaneous WSS vector (Mohamied et al., 2017).

$TimeAveragedWallShearStress\left(TAWSS\right)=\frac{1}{t}\underset{0}{\overset{T}{{\displaystyle \int }}}\left|\overrightarrow{{\tau }_w}\right|dt(1)$

$OscillatoryShearIndex\left(OSI\right)=\frac{1}{2}\left(1-\frac{\left|{{\displaystyle \int }}_0^T\overrightarrow{{\tau }_w}dt\right|}{{{\displaystyle \int }}_0^T\left|\overrightarrow{{\tau }_w}\right|dt}\right)(2)$

$CrossFlowIndex\left(CFI\right)=\frac{1}{T}\underset{0}{\overset{T}{{\displaystyle \int }}}\frac{\overrightarrow{{\tau }_w}}{\left|\overrightarrow{{\tau }_w}\right|}\cdot \left(\overrightarrow{n}\times \frac{{{\displaystyle \int }}_0^T\overrightarrow{{\tau }_w}dt}{\left|{{\displaystyle \int }}_0^T\overrightarrow{{\tau }_w}dt\right|}\right)dt(3)$

$TransverseWallShearStress\left(transWSS\right)=\frac{1}{T}\underset{0}{\overset{T}{{\displaystyle \int }}}\left|\overrightarrow{{\tau }_w}\cdot \left(\overrightarrow{n}\times \frac{{{\displaystyle \int }}_0^T\overrightarrow{{\tau }_w}dt}{\left|{{\displaystyle \int }}_0^T\overrightarrow{{\tau }_w}dt\right|}\right)\right|dt(4)$

where $\overrightarrow{{\tau }_w}$ is the instantaneous WSS vector, $\overrightarrow{n}$ is the surface normal and T is the period of the cardiac cycle.

In order to investigate the effect of vessel motion on these metrics the following three simulations were undertaken:

i) Static Pulsatile – pulsatile velocity with static geometry, to isolate the unsteady effect of flow pulsatility. (The geometry was obtained at t = 0, defined as 50 ms after the R-wave of the ECG.)

ii) Dynamic Nonpulsatile – constant velocity at the inlet with wall movement, to isolate the unsteady effect of RCA motion. (The velocity was the mean over the cardiac cycle.)

iii) Dynamic Pulsatile – pulsatile velocity at the inlet with wall motion, to give the most physiologically realistic solution.

For (ii) and (iii), it is important that the moving reference frame is taken into account when computing OSI, CFI and transWSS. Since the vessel wall is moving and the mesh elements are constantly reorienting, the conventional postprocessing method of mapping WSS vectors at all time steps on to the t = 0 wall geometry results in out-of-plane vectors and can lead to errors. To avoid this, coordinate transformation was computed between t = 0 and the time under consideration for each surface element, and the vector was transformed appropriately.

Distributions of the shear metrics computed in the RCA geometry were converted into 2-D rectangular maps using the t = 0 anatomy as the reference geometry. Parametrisation implemented in VMTK used the curvilinear abscissa and the angular position of each point on the wall in a plane perpendicular to the local centreline. These longitudinal and circumferential coordinates of the rectangular parametric space were used to make the cut along the inner wall of the vessel (the inner curvature at t = 0) and subsequently open the surface using Python, as shown in Supplementary Figure S3. Spearman rank correlation coefficients were computed to quantify similarity between maps. To facilitate this, maps were divided into individual elements using a square grid of 19 × 19 µm (which is smaller than the mesh resolution); pixels were ranked according to their value of the metric being mapped.

Non-Newtonian models

Non-Newtonian blood rheology was investigated using the generalised Carreau-Yasuda model. Eq. 5 defines the dependence of dynamic viscosity (µ) on shear rate ($\dot{\gamma }$) in this model:

$\mu ={\mu }_{\infty }+\left({\mu }_0-{\mu }_{\infty }\right){\left[1+{\left(\lambda \dot{\gamma }\right)}^a\right]}^{\frac{n-1}{a}}(5)$

where ${\mu }_0$ is viscosity at zero shear rate, ${\mu }_{\infty }$ is viscosity at infinite shear rate and $\lambda$, a and n are Carreau-Yasuda rheological parameters (Mahmood et al., 2020).

The parameters are determined empirically and vary significantly with haematocrit, gender, cholesterol level and fitness index. Four different parameter sets (from Cho and Kensey (Cho and Kensey, 1991), Gijsen et al. (Gijsen et al., 1998), van de Vosse et al. (Van de Vosse et al., 2003) and Chen et al. (Chen et al., 2006)) were used (Supplementary Table S1).

To better understand the importance of non-Newtonian effects, we calculated local and averaged global non-Newtonian importance factors, as defined by Johnston et al. (Johnston et al., 2004) and Ballyk et al. (Ballyk et al., 1994), respectively. The local importance factor is simply the ratio of non-Newtonian viscosity to Newtonian viscosity (Eq. (6)) for each mesh point. Values significantly different from one would indicate a non-Newtonian regime. The global importance factor (Eq. (7)) is a spatially averaged measure of deviation from the Newtonian approximation. Values greater than 0.25 indicate non-Newtonian flow behaviour (Soulis et al., 2008; Soulis et al., 2011).

$I_L=\frac{\mu }{{\mu }_{\infty }}(6)$

$I_G=\frac{1}{N}\frac{{\left[{{\displaystyle \sum }}_N{\left(\mu -{\mu }_{\infty }\right)}^2\right]}^{\frac{1}{2}}}{{\mu }_{\infty }}\times 100(7)$

where N is the total number of mesh points on the surface.

Results

Validation and sensitivity analysis

The computational method, including mesh morphing, used to model the flow in the moving RCA was validated by using it to simulate flow in a tube with dynamic curvature, for which results obtained using experimentally validated numerical methods are already available (Santamarina et al., 1998). Figure 2A compares axial velocity profiles 13 tube diameters from the inlet at t = 0.25 s of a 1 Hz cycle: the average difference between velocity values was 3.5%. Figure 2B compares newly computed and previously reported wall shear rates along the inner and outer wall of the tube at the same time. The average differences were 1.2 and 1.6%, respectively.

FIGURE 2

FIGURE 2. Validation by comparison of simulations for a tube of time-varying curvature with previously published values for the same dynamic geometry. (A) Velocity profiles thirteen tube diameters from the inlet and (B) wall shear rate along the inner and outer wall, both at t = 0.25 s in a period of 1 s.

To assess the parameter values required to obtain convergence, simulations were carried out for different mesh resolutions, time steps and number of cardiac cycles until the average of the absolute difference between successive iterations across all mesh points and for all four metrics was approximately 2.5%. For example, the maximum base size was repeatedly halved (element count increased) until the criterion was satisfied. The histograms in Figure 3A show the difference in transWSS estimates between different mesh resolutions. Differences reached 30% at some mesh points when comparing coarse meshes (90 k elements vs. 278 k elements) but were below 10% for almost all the points when comparing the finer meshes (626 k elements vs. 775 k elements). Average absolute percentage differences for TAWSS, CFI, OSI and transWSS between the meshes with 626 k and 775 k elements were 1.3, 2.1, 2.4 and 2.6%, respectively. Therefore, 775 k elements was considered sufficiently fine for mesh independence.

FIGURE 3

FIGURE 3. Sensitivity analysis. (A) Histogram comparing absolute percentage differences in transWSS at each mesh point between different mesh resolutions: i) 90 k elements vs. 278 k elements, ii) 278 k elements vs. 525 k elements, iii) 525 k elements vs. 626 k elements and iv) 626 k elements vs. 775 k elements. (B) Transverse wall shear stress along the inner wall of the RCA, computed with different time steps (∆t).

A similar analysis showed that a time step size of 0.005 s and two cardiac cycles were sufficient; further halving of the time step or increasing the number of cycles did not change the results. For example, Figure 3B shows transWSS along the inner wall of the RCA for four different time step sizes. The change is minimal between 0.005 and 0.0025 s. WSS metrics were calculated over the two cycles.

Influence of RCA motion and flow pulsatility on WSS metrics

Figure 4 shows maps of the four shear metrics for the three cases. Little change is visible in the patterns of TAWSS and OSI between cases, whereas substantial differences are evident in the patterns of CFI and transWSS. Table 1 gives the Spearman rank correlations between the three cases for each metric. Correlation coefficients show that wall motion did not alter TAWSS by much—the correlation coefficient between static pulsatile and dynamic pulsatile cases was 0.93—and that it had by far the largest effect on transWSS (coefficient of 0.19). Ignoring wall motion could lead to intermediate but still substantial differences for the other two metrics, as shown by their correlation coefficients of 0.58 (CFI) and 0.75 (OSI). The fact that transWSS and CFI gave lower coefficients than the OSI shows that wall motion affects truly multidirectional flow more than reversing flow. On the other hand, the consistently high coefficients (all ≥0.94) for the dynamic cases with and without flow pulsatility show that the unsteady effects of flow pulsatility are essentially negligible. Similar observations were made for maps of instantaneous WSS throughout the cardiac cycle: wall motion had a greater effect than pulsatility. Supplementary Table S2 summarises the Spearman rank correlation coefficients at different time points for WSS maps under the three conditions. Note that using rank correlation coefficients means that the patterns are compared but the magnitudes are not.

FIGURE 4

FIGURE 4. Maps of TAWSS, OSI, CFI and transWSS for the three cases. Mean RCA flow is from top to bottom.

TABLE 1

TABLE 1. Spearman rank correlation coefficient between different maps.

Influence of non-Newtonian rheology on RCA hemodynamics

Supplementary Figure S4 shows the evolution of the global importance factor over the cardiac cycle for the four published combinations of Carreau-Yasuda parameters in simulations with pulsatile flow and dynamic geometry. The threshold line identifies the level above which non-Newtonian effects are considered important. Non-Newtonian effects were important for two of the four models.

Figure 5 shows maps of the four shear metrics obtained from simulations of the dynamic pulsatile case using the Newtonian and non-Newtonian parameters of Cho and Kensey (Cho and Kensey, 1991) and Gijsen et al. (Gijsen et al., 1998). These are the most extreme cases: the global importance factor is consistently below the threshold line for the former whereas the latter has a global importance factor well over the threshold line throughout the cardiac cycle. As the parameters were determined empirically, the Newtonian viscosity is different for the two cases.

FIGURE 5

FIGURE 5. Maps of TAWSS, OSI, CFI and transWSS in the Dynamic Pulsatile case for Newtonian and non-Newtonian rheology using the models of Cho and Kensey (Cho and Kensey, 1991), and Gijsen et al. (Gijsen et al., 1998)

While the parameters of Cho and Kensey give Newtonian and non-Newtonian maps that are essentially identical, distinct differences are visible between Newtonian and non-Newtonian maps obtained with the Gijsen et al. parameters. For all four metrics, these differences are seen primarily in the proximal and distal segments of the vessel, labelled in the third row of Figure 5. Additionally, the maps of OSI and CFI show clear differences in the region labelled Segment 1. Spearman rank correlation coefficients between the Newtonian and non-Newtonian versions of the Cho and Kensey model were 0.99 for all metrics. Although the equivalent correlation coefficients were lower for the Gijsen et al. model they were still high, the lowest value being 0.95 for OSI. Despite these high rank correlations, introducing non-Newtonian rheology produced changes in pattern (as opposed to simple scaling effects) that are easily visible in all the maps in Figure 7.

Discussion

The main aim of this study was to examine the effect of wall motion on a range of computed WSS metrics, including transWSS, in the RCA using fully subject-specific geometry, motion and inflow waveform. A secondary aim was to investigate the effect of non-Newtonian blood rheology in the model. It was predicted that metrics capturing near-wall flow components perpendicular to the mean vector – i.e, the CFI and transWSS – would be affected more than the traditional TAWSS and OSI by movement and non-Newtonian rheology as a result of their dependence on time-varying secondary flows. The geometric and flow boundary conditions were obtained in a young, healthy subject; the influence of stenotic disease on WSS metrics, and on the sensitivity of those metrics to motion and rheology, was therefore not considered. Preliminary work optimised the mesh resolution, time step and number of cardiac cycles that were modelled, and showed that the method gave results similar to those previously obtained in an idealised geometry comprising a regular tube with dynamic curvature.

The maps of WSS metrics and the Spearman rank correlation coefficients computed across the different cases showed that wall motion had minimal effect on TAWSS, which is consistent with earlier studies (Santamarina et al., 1998; Zeng et al., 2003; Torii et al., 2010) and theoretical expectations (Lynch et al., 1996). The OSI was found to be more sensitive than TAWSS to changing geometry during the cardiac cycle, as reported by Torii et al. (Torii et al., 2010). A novel finding was that still larger effects of wall motion were seen for the CFI and, even more so, for transWSS. Hence coronary motion should not be ignored when mapping these two metrics. The pulsatile nature of the inflow waveform had much smaller effects.

Polar plots were generated to further understand changes to the multidirectional nature of the flow. Such plots track the direction and magnitude of the instantaneous WSS vector over the cardiac cycle at a single location on the wall. Each point in the plot represents the position of the tip of a WSS vector with its origin at 0,0. The X and Y axes in the plot represent WSS components oriented along the circumference and axis of the vessel, respectively. For truly multidirectional flow of constant magnitude, the tip of the WSS vector would trace out a perfect circle. Uniaxial flow would appear as a straight line through the origin.

Figure 6A shows the polar plot and the mean vector under static pulsatile and dynamic pulsatile conditions for a distal location in the RCA (Region 1 in Figure 4) that is far from boundaries and where transWSS differed substantially between the two cases. For the static pulsatile case, Figure 6A shows that instantaneous WSS vectors are aligned with the mean WSS vector, in the forward going axial direction, throughout most of the cardiac cycle. In the dynamic pulsatile case, on the other hand, the mean WSS vector – although of similar magnitude – has a significant circumferential component (exaggerated by the different x- and y-axis scaling in the figure) and the instantaneous WSS vectors show larger deviations from this mean direction, including in the reverse direction; the magnitudes are again not hugely different from the first case.

FIGURE 6

FIGURE 6. (A) The polar plot and mean vector for static pulsatile and dynamic pulsatile cases for a distal location in the RCA (Region 1 in Figure 4), where transWSS differed substantially between the two cases. Negative values of axial WSS indicate flow travelling down the vessel. (B) Axial vorticity for the static pulsatile case (top) and dynamic pulsatile case (bottom) across the RCA in region one at different times in the cardiac cycle. Positive and negative vorticity values beyond the range of the colour bar scale are shown as white.

Figure 6B compares axial vorticity across the distal RCA, in Region 1, for the static pulsatile and dynamic pulsatile cases. Changes in vortical strength through the cardiac cycle are greater for the dynamic case than the static case. This is consistent with a previous analytical study showing that time-varying curvature can even reverse the direction of Dean vortices (Waters and Pedley, 1999). Furthermore, the static case consistently has a positive dominant vortex throughout the cardiac cycle, presumably due to the constant out-of-plane deformation of the RCA geometry (Caro et al., 1995), whereas two distinct vortices with more equal (albeit changing) strength are seen in the dynamic case, perhaps because the degree of non-planarity is changing. The relation of vorticity to CFI and transWSS is complex because the mean axial velocity also changes over the cardiac cycle. Nevertheless, the large differences in vortical strength and structure between the two cases presumably account for the particularly large differences in CFI and transWSS.

Simulations carried out to study the effect of non-Newtonian blood rheology revealed a dependence on the choice of parameter values. These are determined empirically and may vary between subjects, implying that the effect of non-Newtonian rheology may also vary between subjects. There was a high Spearman rank correlation coefficient for all comparisons. Note, however, that magnitudes can be scaled in a non-linear way and still give a high correlation coefficient based on rank. To provide further information, Figure 7 shows maps of the percentage difference in WSS metrics between Newtonian and non-Newtonian rheology, for both the Cho and Kensey model (Cho and Kensey, 1991) and the model of Gijsen et al. (Gijsen et al., 1998). With the Cho and Kensey model, differences are minimal for all four metrics, as also seen in Figure 5, but that is not the case for the model of Gijsen et al. Not only are there large differences in magnitude between the Newtonian and non-Newtonian assumptions but there are large variations in the difference from point to point in the map, ranging from nearly 25% to approximately -250% for OSI, 15% to -80% for CFI and 60% to -5% for transWSS. Discrepancies for TAWSS are also high, but more uniform. Thus, this non-Newtonian model does appear to have an influence other than simple scaling on parameters related to the directionality of flow.

FIGURE 7

FIGURE 7. Percentage difference in values of the four WSS metrics at each location caused by a switch from a Newtonian to a non-Newtonian model. Note that the colour scale is the same for each vertical pair of maps.

Supplementary Figure S5 shows how the non-Newtonian viscosity varies with shear rate for the two models. The non-Newtonian viscosity converges to Newtonian viscosity (viscosity at infinite shear) relatively quickly for the Cho and Kensey model. This means that the non-Newtonian effect comes into play only for a limited spatio-temporal window during the cardiac cycle, where and when shear rates are low. In the Gijsen et al. model, however, the non-Newtonian viscosity does not converge to the Newtonian value even for high values of shear rate, which means that the non-Newtonian effect is significant at most locations throughout the whole cardiac cycle. According to these parameter values, the Newtonian approximation underestimates the magnitude of the dimensional metrics, TAWSS and transWSS, most likely through a direct influence of the lower viscosity. However, the lower Newtonian viscosity overestimates the non-dimensional metrics, OSI and CFI, presumably because the lower viscosity means less damping of flow disturbances. This effect will also influence transWSS.

Finally, we separate effects of fluid moving past the wall from those of the wall moving past the fluid. As the vessel deforms, distances between points on the wall are altered. In the simplest case, if a cylinder is bent into a curved tube having the same centreline length, the inner curvature must get shorter and the outer curvature must get longer. Even in the absence of any imposed flow, there will be relative motion between the wall and fluid inside the cylinder, leading to the generation of WSS. This type of effect will occur in the RCA in vivo. We have investigated its magnitude by examining WSS metrics in a model with dynamic geometry but no inflow. The resulting map of TAWSS is given in Supplementary Figure S6. Averaged over the map, TAWSS was 28% of the value seen with the same dynamic geometry but a pulsatile inflow, in which shear is generated by wall movement and fluid movement. Because of the way the problem is posed (zero flow at the inlet, free flow at the outlet, inlet anchored in space), the TAWSS generated by wall motion increases with distance down the artery.

WSS caused by wall motion is expected to be dominantly axial, and this was confirmed by creating polar plots of instantaneous WSS vectors over the cardiac cycle for regions in proximal and distal parts of the vessel (Supplementary Figure S7). The magnitude of the vectors increased between proximal and distal parts of the vessel, consistent with the trend in Supplementary Figure S6. The total forward wall displacement should equal the total backward displacement at any particular location since the net deformation of the vessel over a complete cardiac cycle is zero at all locations. WSS generated by this displacement is expected to follow the same trend and indeed the polar plots do show roughly equal excursions in the forward and backward direction. This had two consequences: the mean WSS vector was very low (on the order of 10^–2 Pa), and values of the OSI were everywhere close to its theoretical maximum of 0.5 (Supplementary Figure S6).

CFI and transWSS for the zero-inflow condition are indeterminate – the mean WSS vector, which is required for their calculation, is theoretically 0,0 when flow is purely oscillatory, and is small in magnitude and arbitrary in direction when flow is nearly so (Peiffer et al., 2013). To estimate the contribution of wall motion to these shear metrics, and to TAWSS and OSI, we subtracted the WSS vectors for the zero-inflow case from the vectors for the case of dynamic geometry with pulsatile flow. Maps of the recomputed metrics are shown in Supplementary Figure S8 and means for the entire maps are given in Supplementary Table S3. Decreases in all metrics are discernible when the zero-inflow vectors are subtracted from the physiological case (cf dynamic geometry with pulsatile inflow sub-parts of Figure 4), but only on close inspection; the overall trends are unchanged.

Limitations

The RCA does not have major bifurcations, unlike the left coronary artery system, but it does have smaller branches. These were not captured in the model since the diameter of the branch arteries is smaller than the MRI slice thickness. This means that the loss of blood to these branches is not included in the model and hence the mean flow rate will be increasingly too high on progressing down the vessel. TAWSS and transWSS magnitude will show the same trend, and hence are not suitable at this stage for comparison with, say, patterns of lesions; the data show only how these metrics are affected locally by coronary motion and blood viscosity. An estimate of the size of this effect was obtained by removing taper from the RCA, a method previously employed by Zeng et al. (Zeng et al., 2003) which assumes that mean axial velocity is constant at all points along the vessel. Substantial changes are seen (Supplementary Figure S9): in particular, the largest values of the dimensional TAWSS and transWSS metrics occur in the proximal rather than distal parts of the vessel. Smaller changes occur in the non-dimensional OSI and CFI metrics; these must reflect changes in flow structure as well as rate.

The subject-specific waveform, obtained by phase-contrast MRI at the proximal end of the vessel, indicates that cross-sectionally averaged blood flow was in the proximal-to-distal direction at all times during the cardiac cycle. Some published waveforms show periods of reverse flow at other locations. If such a waveform were imposed, it would at least have the effect of increasing OSI values, even in the absence of secondary flow.

Future work should incorporate branches, add geometries from other subjects, include a parametric study of inflow waveforms, and compare patterns of WSS metrics with the pattern of disease.

Conclusion

This study confirms previous observations that the TAWSS is not substantially affected by coronary artery motion, and that the OSI is affected more. A novel observation is that metrics capturing transverse flows – the CFI and transWSS – are affected even more than the OSI; this effect can be attributed to changing vortical structures resulting from the dynamic curvature. An additional finding was that models of non-Newtonian blood rheology can also have a significant influence if they increase viscosity at most locations and times; that leads to increased magnitude of dimensional shear metrics but also to changes in the relative magnitudes of shear metrics from location to location. Hence studies of transWSS in coronary arteries should take vessel motion into account and, depending on the characteristics of the subject, might also need to account for non-Newtonian blood rheology. This is necessary for valid comparisons with the distribution of coronary disease.

Data availability statement

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

Author contributions

PK, RT, SS, and PW contributed to the conceptualisation and writing. PK conducted the simulations. RT and PDG obtained the geometries.

Funding

Funded by the EPSRC Centre for Doctoral Training in Fluid Dynamics Across Scales and by British Heart Foundation grant RG/11/5/28743.

Conflict of interest

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

Publisher’s note

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

Supplementary material

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

References

Apostolidis, A. J., Moyer, A. P., and Beris, A. N. (2016). Non-newtonian effects in simulations of coronary arterial blood flow. J. Newt. Fluid Mech. 233, 155–165. doi:10.1016/j.jnnfm.2016.03.008

CrossRef Full Text | Google Scholar

Arshad, M., Ghim, M., Mohamied, Y., Sherwin, S. J., and Weinberg, P. D. (2021). Endothelial cells do not align with the mean wall shear stress vector. J. R. Soc. Interface 18, 20200772. doi:10.1098/rsif.2020.0772

PubMed Abstract | CrossRef Full Text | Google Scholar

Ballyk, P. D., Steinman, D. A., and Ethier, C. R. (1994). Simulation of non-Newtonian blood flow in an end-to-side anastomosis. Biorheology 31, 565–586. doi:10.3233/bir-1994-31505

PubMed Abstract | CrossRef Full Text | Google Scholar

Baran, I., and Popović, J. (2007). Automatic rigging and animation of 3d characters. ACM Trans. Graph. 26, 72. doi:10.1145/1276377.1276467

CrossRef Full Text | Google Scholar

Barnes, S. E., and Weinberg, P. D. (1998). Contrasting patterns of spontaneous aortic disease in young and old rabbits. Arterioscler. Thromb. Vasc. Biol. 18 (2), 300–308. doi:10.1161/01.atv.18.2.300

PubMed Abstract | CrossRef Full Text | Google Scholar

Caro, C. G., Fitz-Gerald, J. M., and Schroter, R. C. (1971). Atheroma and arterial wall shear. Observation, correlation and proposal of a shear dependent mass transfer mechanism for atherogenesis. Proc. R. Soc. Lond. B Biol. Sci. 177 (1046), 109–159. doi:10.1098/rspb.1971.0019

PubMed Abstract | CrossRef Full Text | Google Scholar

Caro, C. G., Parker, K. H., and Doorly, D. J. (1995). Essentials of blood flow. Perfusion 10, 131–134. doi:10.1177/026765919501000302

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Lu, X. Y., and Wang, W. (2006). Non-newtonian effects of blood flow on hemodynamics in distal vascular graft anastomoses. J. Biomech. 39, 1983–1995. doi:10.1016/j.jbiomech.2005.06.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Cherry, E. M., and Eaton, J. K. (2013). Shear thinning effects on blood flow in straight and curved tubes. Phys. Fluids 25, 073104. doi:10.1063/1.4816369

CrossRef Full Text | Google Scholar

Cho, Y. I., and Kensey, K. R. (1991). Effects of the non-newtonian viscosity of blood on flows in a diseased arterial vessel. Part 1: steady flows. Biorheology 28, 241–262. doi:10.3233/bir-1991-283-415

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman, M. H., Brinkman, A. M., Qin, J. J., and Seed, W. A. (1993). Relation between coronary artery geometry and the distribution of early sudanophilic lesions. Atherosclerosis 98, 193–199. doi:10.1016/0021-9150(93)90128-h

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghim, M., Alpresa, P., Yang, S. W., Braakman, S. T., Gray, S. G., Sherwin, S. J., et al. (2017). Visualization of three pathways for macromolecule transport across cultured endothelium and their modification by flow. Am. J. Physiol. Heart Circulatory Physiol. 313, H959–H973. doi:10.1152/ajpheart.00218.2017

CrossRef Full Text | Google Scholar

Ghim, M., Pang, K. T., Arshad, M., Wang, X., and Weinberg, P. D. (2018). A novel method for segmenting growth of cells in sheared endothelial culture reveals the secretion of an anti-inflammatory mediator. J. Biol. Eng. 12, 15. doi:10.1186/s13036-018-0107-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghim, M., Pang, K. T., Burnap, S. A., Baig, F., Yin, X., Arshad, M., et al. (2021). Endothelial cells exposed to atheroprotective flow secrete follistatin-like 1 protein which reduces transcytosis and inflammation. Atherosclerosis 333, 56–66. doi:10.1016/j.atherosclerosis.2021.08.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Gibson, C. M., Diaz, L., Kandarpa, K., Sacks, F. M., Pasternak, R. C., Sandor, T., et al. (1993). Relation of vessel wall shear stress to atherosclerosis progression in human coronary arteries. Arterioscler. Thromb. 13, 310–315. doi:10.1161/01.atv.13.2.310

PubMed Abstract | CrossRef Full Text | Google Scholar

Gijsen, F. J., van de Vosse, F. N., and Janssen, J. D. (1998). Wall shear stress in backward-facing step flow of a red blood cell suspension. Biorheology 35, 263–279. doi:10.1016/s0006-355x(99)80010-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Grøttum, P., Svindland, A., and Walløe, L. (1983). Localization of atherosclerotic lesions in the bifurcation of the main left coronary artery. Atherosclerosis 47, 55–62. doi:10.1016/0021-9150(83)90071-0

PubMed Abstract | CrossRef Full Text | Google Scholar

He, X., and Ku, D. N. (1996). Pulsatile flow in the human left coronary artery bifurcation: average conditions. J. Biomech. Eng. 118, 74–82. doi:10.1115/1.2795948

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoogendoorn, A., Kok, A. M., Hartman, E. M. J., de Nisco, G., Casadonte, L., Chiastra, C., et al. (2020). Multidirectional wall shear stress promotes advanced coronary plaque development: comparing five shear stress metrics. Cardiovasc. Res. 116, 1136–1146. doi:10.1093/cvr/cvz212

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnston, B. M., Johnston, P. R., Corney, S., and Kilpatrick, D. (2004). Non-newtonian blood flow in human right coronary arteries: steady state simulations. J. Biomech. 37, 709–720. doi:10.1016/j.jbiomech.2003.09.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnston, B. M., Johnston, P. R., Corney, S., and Kilpatrick, D. (2006). Non-newtonian blood flow in human right coronary arteries: transient simulations. J. Biomech. 39, 1116–1128. doi:10.1016/j.jbiomech.2005.01.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Keegan, J., Gatehouse, P. D., Yang, G. Z., and Firmin, D. N. (2004). Spiral phase velocity mapping of left and right coronary artery blood flow: correction for through-plane motion using selective fat-only excitation. J. Magn. Reson. Imaging 20, 953–960. doi:10.1002/jmri.20208

PubMed Abstract | CrossRef Full Text | Google Scholar

Kok, A. M., Molony, D. S., Timmins, L. H., Ko, Y. A., Boersma, E., Eshtehardi, P., et al. (2019). The influence of multidirectional shear stress on plaque progression and composition changes in human coronary arteries. EuroIntervention 15, 692–699. doi:10.4244/eij-d-18-00529

PubMed Abstract | CrossRef Full Text | Google Scholar

Ku, D. N., Giddens, D. P., Zarins, C. K., and Glagov, S. (1985). Pulsatile flow and atherosclerosis in the human carotid bifurcation. Positive correlation between plaque location and low oscillating shear stress. Arteriosclerosis. 53, 293–302. doi:10.1161/01.atv.5.3.293

CrossRef Full Text | Google Scholar

Lewis, J. P., Cordner, M., and Fong, N. (2000). “Pose space deformation: a unified approach to shape interpolation and skeleton-driven deformation,” in Proc 27th ann conf computer graphics and interactive techniques, New Orleans, LA, July 23–28, 2000, 165–172.

Google Scholar

Lynch, D. G., Waters, S. L., and Pedley, T. J. (1996). Flow in a tube with non-uniform, time-dependent curvature: governing equations and simple examples. J. Fluid Mech. 323, 237–265. doi:10.1017/s0022112096000900

CrossRef Full Text | Google Scholar

Mahmood, R., Bilal, S., Khan, I., Kousar, N., Seikh, A. H., and Sherif, E. S. (2020). A comprehensive finite element examination of Carreau Yasuda fluid model in a lid driven cavity and channel with obstacle by way of kinetic energy and drag and lift coefficient measurements. J. Mater. Res. Technol. 9, 1785–1800. doi:10.1016/j.jmrt.2019.12.010

CrossRef Full Text | Google Scholar

Mitchell, J. R. A., and Schwartz, C. J. (1965). Arterial disease. Oxford: Blackwell Scientific.

Google Scholar

Mohamied, Y., Rowland, E. M., Bailey, E. L., Sherwin, S. J., Schwartz, M. A., and Weinberg, P. D. (2015). Change of direction in the biomechanics of atherosclerosis. Ann. Biomed. Eng. 43, 16–25. doi:10.1007/s10439-014-1095-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohamied, Y., Sherwin, S. J., and Weinberg, P. D. (2017). Understanding the fluid mechanics behind transverse wall shear stress. J. Biomech. 50, 102–109. doi:10.1016/j.jbiomech.2016.11.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, J. E., Weydahl, E. S., and Santamarina, A. (2001). Frequency dependence of dynamic curvature effects on flow through coronary arteries. J. Biomech. Eng. 123, 129–133. doi:10.1115/1.1351806

PubMed Abstract | CrossRef Full Text | Google Scholar

Myers, J. G., Moore, J. A., Ojha, M., Johnston, K. W., and Ethier, C. R. (2001). Factors influencing blood flow patterns in the human right coronary artery. Ann. Biomed. Eng. 29, 109–120. doi:10.1114/1.1349703

PubMed Abstract | CrossRef Full Text | Google Scholar

Peiffer, V., Sherwin, S. J., and Weinberg, P. D. (2013). Computation in the rabbit aorta of a new metric–the transverse wall shear stress–to quantify the multidirectional character of disturbed blood flow. J. biomechanics 46, 2651–2658. doi:10.1016/j.jbiomech.2013.08.003

CrossRef Full Text | Google Scholar

Perktold, K., Hofer, M., Rappitsch, G., Loew, M., Kuban, B. D., and Friedman, M. H. (1998). Validated computation of physiologic flow in a realistic coronary artery branch. J. Biomech. 31, 217–228. doi:10.1016/s0021-9290(97)00118-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramaswamy, S. D., Vigmostad, S. C., Wahle, A., Lai, Y. G., Olszewski, M. E., Braddy, K. C., et al. (2004). Fluid dynamic analysis in a human left anterior descending coronary artery with arterial motion. Ann. Biomed. Eng. 32, 1628–1641. doi:10.1007/s10439-004-7816-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Santamarina, A., Weydahl, E., Siegel, J. M., and Moore, J. E. (1998). Computational analysis of flow in a curved tube model of the coronary arteries: effects of time-varying curvature. Ann. Biomed. Eng. 26, 944–954. doi:10.1114/1.113

PubMed Abstract | CrossRef Full Text | Google Scholar

Sebkhi, A., and Weinberg, P. D. (1994). Age-related variations in transport properties of the rabbit arterial wall near branches. Atherosclerosis 106 (1), 1–8. doi:10.1016/0021-9150(94)90077-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Sinzinger, H., Silberbauer, K., and Auerswald, W. (1980). Quantitative investigation of sudanophilic lesions around the aortic ostia of human fetuses, newborn and children. J. Vasc. Res. 17, 44–52. doi:10.1159/000158233

CrossRef Full Text | Google Scholar

Soulis, J. V., Giannoglou, G. D., Chatzizisis, Y. S., Farmakis, T. M., Giannakoulas, G. A., Parcharidis, G. E., et al. (2006). Spatial and phasic oscillation of non-Newtonian wall shear stress in human left coronary artery bifurcation: an insight to atherogenesis. Coron. Artery Dis. 17, 351–358. doi:10.1097/00019501-200606000-00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Soulis, J. V., Giannoglou, G. D., Chatzizisis, Y. S., Seralidou, K. V., Parcharidis, G. E., and Louridas, G. E. (2008). Non-newtonian models for molecular viscosity and wall shear stress in a 3D reconstructed human left coronary artery. Med. Eng. Phys. 30, 9–19. doi:10.1016/j.medengphy.2007.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Soulis, J. V., Seralidou, K. V., Chatzizisis, Y. S., and Giannoglou, G. D. (2011). “Severity parameter and global importance factor of Non-newtonian models in 3D reconstructed human Left Coronary Artery,” in 3rd micro and nano flows conference (MNF2011), Thessaloniki, August 22–24, 2011. Available at: http://bura.brunel.ac.uk/handle/2438/6835.

Google Scholar

Spiegel, M., Redel, T., Zhang, Y. J., Struffert, T., Hornegger, J., Grossman, R. G., et al. (2011). Tetrahedral vs. polyhedral mesh size evaluation on flow velocity and wall shear stress for cerebral hemodynamic simulation. Comput. Methods Biomech. Biomed. Engin. 14, 9–22. doi:10.1080/10255842.2010.518565

PubMed Abstract | CrossRef Full Text | Google Scholar

Torii, R., Keegan, J., Wood, N. B., Dowsey, A. W., Hughes, A. D., Yang, G. Z., et al. (2010). MR image-based geometric and hemodynamic investigation of the right coronary artery with dynamic vessel motion. Ann. Biomed. Eng. 38, 2606–2620. doi:10.1007/s10439-010-0008-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Van de Vosse, F. N., De Hart, J., Van Oijen, C. H., Bessems, D., Gunther, T. W., Segal, A., et al. (2003). Finite-element-based computational methods for cardiovascular fluid-structure interaction. J. Eng. Math. 47, 335–368. doi:10.1023/b:engi.0000007985.17625.43

CrossRef Full Text | Google Scholar

van Wyk, S., Prahl Wittberg, L., Bulusu, K. V., Fuchs, L., and Plesniak, M. W. (2015). Non-Newtonian perspectives on pulsatile blood-analog flows in a 180 curved artery model. Phys. Fluids 27, 071901. doi:10.1063/1.4923311

CrossRef Full Text | Google Scholar

Wang, C., Baker, B. M., Chen, C. S., and Schwartz, M. A. (2013). Endothelial cell sensing of flow direction. Arterioscler. Thromb. Vasc. Biol. 33, 2130–2136. doi:10.1161/atvbaha.113.301826

PubMed Abstract | CrossRef Full Text | Google Scholar

Waters, S. L., and Pedley, T. J. (1999). Oscillatory flow in a tube of time-dependent curvature: part I. Perturbation to flow in a stationary curved tube. J. Fluid Mech. 383, 327–352. doi:10.1017/s0022112099004085

CrossRef Full Text | Google Scholar

Weinberg, P. D. (2004). Rate-limiting steps in the development of atherosclerosis: the response-to-influx theory. J. Vasc. Res. 41 (1), 1–17. doi:10.1159/000076124

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, D., Ding, Z., Friedman, M. H., and Ethier, C. R. (2003). Effects of cardiac motion on right coronary artery hemodynamics. Ann. Biomed. Eng. 31, 420–429. doi:10.1114/1.1560631

PubMed Abstract | CrossRef Full Text | Google Scholar