Abstract
A diseased coronary artery has been modeled to study the implications of plaque morphology on the fluid dynamics. In our previous study, we have successfully classified the coronary plaques of 42 patients who underwent intravascular ultrasound (IVUS) into four-types (Type I, Type II, Type III, and Type IV) based on the plaque morphology. In this study, we demonstrate that, for the same degree of stenosis (height of the plaques), hemodynamics parameters are strongly dependent on the plaque shape. This study is the first one to clearly demonstrate that in addition to wall shear stress, presence of turbulence and location of transition from laminar to turbulence state are additional hemodynamics parameters to identify plaques vulnerable to rupture.
1. Introduction
Coronary artery disease (CAD) is a progressive disease characterized by the accumulation of plaques on the artery walls. CAD is initiated by the deposition of fatty materials in the coronary artery resulting in the thickening and formation of streaks of plaque on the artery walls. During these early stages, the plaques are not of significant consequence to the flow dynamics, and the flow does not deviate from the laminar state present in a normal coronary artery [1–3]. As time progresses, these plaques start growing inwards into the lumen (channel in which the blood flows), resulting in the localized narrowing of the artery or stenosis of the artery lumen, and thus, playing a critical role in altering the flow characteristics. It is clear that once a stenosis is developed, the blood flow is further disturbed and hemodynamic parameters continue to play a crucial role as the stenosis progresses [4]. Hence, as the stenosis increases this alters the flow characteristics causing a laminar to turbulent transition in the blood flow. The evidence of turbulence in regions distal to a stenosis was first demonstrated based on detection of high frequency pressure signals for lesions with 23%–76% stenoses [5]. Based on clinical findings as murmurs [6], laboratory experiments [1, 7, 8], and medical imaging methods [9–11], it is now well established that stenosed coronary artery creates high levels of turbulence, thus significantly modifying the flow characteristics.
As of today, the degree of stenosis has been a classical metric to define the extent of the disease for medical purposes. The stenoses are commonly assessed as a percentage of obstruction in the diameter of the lumen. It is a common practice to characterize the stenosis using the percentage obstruction or the height of the obstruction criterion. However, recent studies have clearly shown that the characteristic flow dynamics due to stenosis are strongly influenced by factors such as morphological structure of the stenosis developed, the pulsatile nature of the flow, and numerous others [5, 12, 13]. Under laminar flow assumptions, Stroud et al. [14] have demonstrated that for idealized test cases with identical degree of stenosis, axisymmetric and nonaxisymmetric lesions exhibited significantly different flow fields and wall shear stress distributions, thus suggesting the importance of plaque morphology rather than a percentage based criterion to characterize the flow. Thus, this method to define the extent of the disease proves to be rudimentary and incomplete since it provides the simplest and most inaccurate representation of a patient’s actual condition. Hence, it is becoming clear that the degree of stenosis is not a true representation of the severity or characteristics of the turbulent flow generated by the stenosis.
With advances in medical imaging, it is now possible to obtain a comprehensive view of the coronary anatomy to assess plaque morphology using imaging techniques such as intravascular ultrasound (IVUS), computed tomography (CT), X-ray angiography, and magnetic resonance imaging (MRI). Hence, there is a clear need for well-defined metrics in addition to plaque height in order to identify patient-specific point of transition from laminar to turbulent flow and also to identify patient-specific poststenotic behavior such as turbulence levels and exact locations where the flow relaminarizes. This information is of clinical relevance because it provides important clinical interpretation. Turbulence is of pathological importance because localized bursts of shear stress may erode already weakened endothelial patches from the arterial wall [15] and increase the probability of rupture risk. Prediction of stent placement can be achieved if patient-specific location of transition is identified, hence preventing the flow to become turbulent. Furthermore, if patient-specific poststenotic flow behavior (i.e., extent of turbulence region) is identified, an assessment of the plaque vulnerability to rupture can be achieved since turbulence regions are associated with strong vorticity and high levels of kinetic energy which are dissipated as high viscous stresses.
The focus of the present study is to conduct a fluid dynamic simulation in diseased coronary artery to predict the transition to turbulence state, if any. Beaumont et al. [16] have conducted fuzzy logic analysis on plaques based on IVUS measurements, and their study has revealed four types of plaques (see Figure 1). In the present study the differences in the location of transition to the turbulent state are studied. The paper is organized as follows: Section 2 discusses the numerical methodology and validation tests, Section 3 discusses the results, and the conclusions are presented in Section 4.
(a)
(b)
(c)
(d)
2. Methodology
In this section we discuss the numerical methodology including assumptions, the construction of the 3D solid model, boundary conditions, validations, and verifications. The governing equations have been solved numerically using the commercially available software FLUENT 6.0. A 3D model of the artery is constructed using the plaque morphology from the IVUS. Refer to Beaumont et al. [16] and Bhaganagar et al. [17] for details on extraction of plaque morphology from IVUS. Unsteady flow simulations are conducted for a long time starting from laminar state.
2.1. Assumptions
The length of the artery is assumed to be 72 mm in length, and the outer diameter and the inner diameter dimensions are obtained from the IVUS measurements. The plaque inlet was positioned at an axial distance of 40 mm from the inlet of the artery in order to allow the flow to develop naturally before entering the stenotic region; while the plaque exit, after using the 13 different, 1 mm spaced cross sections, was position at 52 mm from the inlet of the artery. For the patient-specific flow simulations, a physiological waveform flow was specified at the inlet of artery. The waveform specified is based on Womersley’s solution [18]. The waveform profile is shown in Figure 2. The Womersley number used for all the patient-specific simulations was chosen to be 5, which is characteristic of arteries. The time duration for one cardiac cycle of the specified waveform is 0.33 seconds. The Reynolds number achieved during peak systole is approximately 2100, and the Reynolds number achieved during peak diastole is 1700. The Womersley and Reynolds number were based on the inlet diameter. Furthermore, a turbulence intensity (Tu) value of 1.5% was specified at the inlet. Justification of the Tu value of 1.5% as an inlet boundary condition is explained in detail at Section 2.4. Furthermore, an open boundary condition with a relative static pressure of 0 Pa was specified at the outlet of the diseased coronary artery for all cases. The inlet Tu value was introduced in order to represent initial disturbances in the laminar flow, which allows the transition to occur realistically. It has been shown from in vivo studies that low disturbances in a laminar flow are present just before transition [16], thus justifying the use of low inlet turbulence intensities as a realistic representation of an initial disturbance level at the inlet. Additionally, the arterial wall was treated as a rigid solid with a no-slip boundary condition where the assumption of a rigid boundary is reasonable since development of atherosclerosis in arteries causes a considerable reduction in the elastic property of its wall. Blood was assumed to be a Newtonian and incompressible fluid with a density of 998 kg/m3 and dynamic viscosity of 1.0 mPa·s. The assumption of Newtonian fluid is justified, as it is known that blood behaves as a Newtonian fluid in large arteries especially at high shear rates. The convergence criterion for all the residuals was manually defined as 10−7 in order to control the accuracy of the solution. The timestep size was chosen as 0.0002855 s with approximately 1160 timesteps per cardiac cycle. A timestep independency test is described in Section 2.4. The maximum number of iterations per timestep was defined as 200 iterations. Approximately 50 iterations were needed to reach convergence at each timestep. Five cardiac cycles were simulated for each transient calculation in order to ensure statistical convergence (cycle-to-cycle differences less than 1%).
2.2. 3D Solid Model
Initially, the data obtained from the IVUS at four circumferential locations was imported into SolidWorks. Based on the data, the plaque height was measured from the outer wall. The data between these four points has been interpolated using spline interpolation. The procedure was repeated for the 13 cross sections along the length of the artery. It was ensured that the accurate shape of the plaque is represented after the spline interpolation. After modeling the stenosis, the geometry was imported into the ANSYS software where a coarse, medium, and fine mesh were created for a mesh independency test. The independency test results are described is Section 2.4. The tetrahedral meshes generated for the lumen are shown in Figure 3.
2.3. Unsteady Flow Simulations Using Transition Turbulence Model
The flow is started from initial laminar flow conditions. At the inlet, the flow Reynolds number is 200 based on the diameter of the artery. At the output constant pressure, boundary condition has been prescribed. A transition model is applied to capture the transition, if any, to turbulent state. For this purpose, the shear stress transport (SST) transition model has been applied. The SST transition model is based on the coupling of Menter’s transport equations with two other transport equations, one for the intermittency and one for the transition onset criteria. Menter’s SST model is utilized to effectively blend the robust and accurate formulation of the model in the near-wall regions with the free-stream independence of the model in the far field. With correlation strategy based on local flow parameters but keeping the effects of flow history, the transition model [19, 20] introduces two further transport equations, intermittency and transition momentum thickness Reynolds number. The model is based on established correlations as a starting point, such as pressure gradient, external turbulence level, transition length, and the model for low turbulence intensity [21]. It also includes an innovative treatment of separation-induced transition of importance for the present study. The full model introduces new proprietary empirical correlations to cover the large range of transition mechanisms and has been validated in conjunction with the SST turbulence model [20]. Futher details of these models with validation and justication of using this model in physiological flow conditions were reported by Tan et al. [22] and are shortly explained in Section 2.4.
The experimental results of Ahmed and Giddens [23] were used for initial validation of the model.
2.4. Validations and Verification
A mesh independence test was carried out using the SST transition model with realistic boundary conditions for Type I (protrusion). The boundary conditions will be further described. Results in terms of maximum axial wall shear stresses (WSS) and peak turbulence kinetic energy (TKE) at the peak systolic phase were analyzed. Initial flow simulations results were mesh independent after using a grid of 250,000 tetrahedral elements with elements concentrated near the wall where the velocity gradients were expected to be high. A very fine near wall resolution of the wall boundary cells along the length of the stenosis was ensured in order to satisfy the requirements of the turbulence and transition models. Therefore, a value of less than 2 along the walls was ensured. The differences between a 130,000 elements mesh and a 250,000 elements mesh in terms of maximum WSS and peak TKE during the peak systolic phase of the cardiac cycle were less than 3%. Additionally, the differences between the 250,000 elements mesh and a 600,000 elements mesh were less than 2%. Details of the mesh independency test are shown in Table 1. Furthermore, a timestep independency test was carried out using the SST transition model. Transient calculations of a timestep size of 0.001142 s, 0.000571 s, and 0.0002855 s were performed. Such timestep sizes resulted in an approximation of 290, 580, and 1160 number of timesteps per cycle, respectively. Five cardiac cycles were simulated for each transient calculation in order to ensure statistical convergence or periodicity. After five complete cycles, the cycle-to-cycle differences were less than 1%. Similar to the mesh independency test, maximum WSS and peak TKE during the peak systolic phase were analyzed. The differences between the timesteps were approximately 2%. The results are tabulated in Table 2.
Additionally, a sensitivity test of the inlet turbulence intensity (Tu) was carried out by using the SST transition model. The experimental results of Ahmed and Giddens [23] were used for the initial validation of the sensitivity test. The test was carried out by using the idealized geometry and steady flow conditions () depicted in Deshpande and Giddens [24]. The form of the stenosis chosen was where , with at the throat of the stenosis, is the unobstructed tube diameter, taken as 1 cm, and is the local radius which varies with the axial direction . The constants given correspond to a stenosed tube with a 75% area reduction as those reported by Ahmed and Giddens [23]. The tube was extended upstream and downstream of the throat by and , respectively. The inlet and outlet lengths are those found in Ryval et al. [25] who found that the specified inlet and outlet lengths were sufficient to predict recirculation length.
Figure 4 shows different levels of Tu between 1.5% and 6.75% at three axial positions ( cm, cm, and cm) at the poststenotic region. In the turbulence intensity profiles, the radial position, , was normalized by the local radius where corresponds to the upper wall of the model, while corresponds to the center of the tube. It is shown that values between 1.5% and 2% are able to capture the experimental data quite well. In comparison, it is shown that values between 3% and 6.75% cause the flow to recover rapidly at the post-stenotic region. Hence, for the simulations performed in this study, a Tu value of 1.5% was specified at the inlet. The results obtained for this sensitivity test are in complete agreement with those performed by Tan et al. [22]. Furthermore, a validation test was conducted in order to confirm that the SST transition model is capable of predicting turbulence under pulsatile flow conditions. The pulsatile flow conditions were based on the experiments conducted by Ahmed and Giddens [26]. Using the geometry specified for the Tu sensitivity test, the imposed inlet pulsatile flow condition varied approximately sinusoidally with a frequency of 4 Hz, resulting in a Womersley number of 7.5. The fluid properties used were those based on the experiments by Ahmed and Giddens [23] where they used a solution of water and glycerin. The solution was 63% glycerol by weight ( stokes). The amplitude of the transient flow ranged from of 200 to 1000.
(a)
(b)
(c)
Figure 5 shows the axial velocity profiles at peak systole for different axial positions starting at the throat ( cm) up to the post-stenotic region ( cm). All velocities were normalized by , one-half of the sum of the maximum and minimum upstream centerline velocities. It can be seen from the figure that the model represents the experimental data fairly well. From the experimental data, it can be noted that there is a recirculation region from cm to cm where the recirculation length was captured fairly well by the model. The results obtained for the validation test employed are also in good agreement with those obtained by Tan et al. [22], hence validating the capability of the model to predict turbulence well under pulsatile flow condition.
(a)
(b)
(c)
(d)
(e)
(f)
3. Results
The simulations have been performed for protrusions of Type I through Type IV (see Figure 1) using the physiological waveform depicted in Figure 2. In this section, a comparison of the similarities and differences in the flow characteristics amongst these types is performed. The flow characteristics of interest are the three components of the mean velocity profiles, the wall shear stresses (WSS) due to the viscous and turbulent terms, the turbulent kinetic energy (TKE), and the turbulence intensity (Tu), where the last two characteristics are an important representation of the turbulent state of the flow.
3.1. Mean Velocity Profiles
First, the differences and similarities in the axial component of the mean velocity are analyzed during the peak systolic phase of the cardiac cycle. Figure 6 shows axial component of the mean velocity. Seven different axial locations along the length of the artery have been selected for the analysis. The locations correspond to mm, 35 mm (inlet region), mm (proximal region of stenosis), mm (peak location), mm (distal region of stenosis), and mm, 60 mm (exit region). In the mean velocity profiles, the wall normal distance () has been scaled by the lumen height () where and correspond to the lower wall and upper wall, respectively. At the inlet region ( mm, 35 mm), all types have similar mean axial velocity profiles with a maximum mean axial velocity of approximately 0.6 m/s. Such similarity is expected since the blood flow is natural and has not encountered the stenotic region. Differences in the mean axial velocity profiles start emerging in the proximal region of the stenosis ( mm). For Type I, the profile is no longer symmetric where the maximum centerline mean velocity is shifted towards the lower wall. In comparison, for Type II and Type III, the profile is no longer symmetric where the maximum centerline mean velocity is shifted towards the upper wall. For Type IV, the profile remains almost unchanged. For Type I through Type III, the maximum mean axial velocity is approximately 1.5 m/s compared to a maximum mean axial velocity of approximately 1.2 m/s exhibited by Type IV. At the peak location ( mm) of the stenosis, the mean axial velocity profile for Type I remains unchanged. Type IV exhibits minimal flow reversal at the upper wall. In contrast, the velocity profiles for Type II and Type III exhibit flow reversal near the lower wall. This trend continues for Type II and Type III for the distal portion of the stenosis ( mm) until the flow begins to stabilize at the exit region ( mm, 60 mm). For Type I, flow reversal is present at the distal portion of the stenosis ( mm), whereas at the exit region ( mm, 60 mm), the flow also begins to stabilize. In comparison, for Type IV, the flow starts to stabilize at distal portion of the stenosis ( mm) rather than at the exit region.
(a)
(b)
(c)
(d)
Similarly, the differences and similarities in the radial component of the mean velocity are analyzed during the peak systolic phase. Figure 7 shows the radial component of the mean velocity for the four types. At the inlet region, the flow is unidirectional with no radial component of velocity for the four types. As the flow approaches the proximal region of the plaque, the flow becomes three-dimensional with a nontrivial radial component of velocity. At the proximal region, for Type I, Type II, and Type IV, the flow is directed towards the wall in the upper-half region of the lumen, while the flow is directed away from the wall in the lower-half region of the lumen. In contrast, for Type III, the flow is directed away from the wall near the upper and lower wall regions. At the peak region, for Type I and Type IV, the flow is now directed away from the wall near the upper and lower wall regions, whereas for Type II, the flow is directed away from the wall near the upper wall region while remaining unidirectional on the lower wall. In comparison, the flow for Type III is directed towards the wall in the upper-half region of the lumen and directed away from the wall in the lower-half region of the lumen. The highest radial mean velocity values are observed for Type I and Type II corresponding to −0.18 m/s and −0.3 m/s, respectively. At the exit region, the flow is almost unidirectional for Type II and Type IV, whereas for Type I and Type III, the flow is still disturbed.
(a)
(b)
(c)
(d)
Furthermore, the three mean velocity components have been averaged in the plane and plotted against the wall normal distance in order to study the overall effect of the stenosis in the natural blood flow. Figure 8 shows the mean velocity profiles averaged in plane plotted against wall-normal distance during peak systole for the axial, radial, and circumferential components of the mean flow for Type I through Type IV. For the axial mean velocity component, Type I and Type IV exhibit minor disturbances in the natural blood flow. Type I and Type IV are almost symmetric with a peak velocity of around 0.6 m/s. This profile is similar to those observed in the inlet region in Figure 6. In contrast, Type II and Type III exhibit a major overall disturbance of the natural blood flow due to the shape of the stenosis. For both types, the profile is no longer symmetric with peak velocities concentrated above the centerline. The peak axial mean velocities exhibited for Type II and Type III are 0.75 m/s and 0.85 m/s, respectively. For the radial mean velocity component, Type IV exhibits minor disturbances in the overall radial flow, behaving unidirectionally. In comparison, Type I and Type II exhibit major disturbances in the radial mean flow where the flow is overall directed away from the wall in the lower-half region of the lumen and directed towards the wall in the upper-half region, whereas for Type III, the flow is overall directed away from the wall in the lower-half region of the lumen while remaining unidirectional in the upper-half region. For the circumferential mean velocity component, the circumferential flow overall remains undisturbed and unidirectional with low magnitudes for all types, whereas for Type II, the circumferential flow attains a maximum of 0.02 m/s near the upper wall.
(a)
(b)
(c)
3.2. Flow Parameters
Next, an analysis of flow parameters such as turbulence intensity (Tu), turbulence kinetic energy (TKE), and wall shear stresses (WSS) is presented. Such analysis is of importance because these parameters are true representations of the turbulent nature of the flow. First, the differences and similarities in the turbulence kinetic energy profiles for Type I through Type IV are analyzed during the peak systolic phase of the cardiac cycle.
Figure 9 shows ths TKE profiles for Type I through Type IV. Similar to the mean velocity profiles, seven different axial locations along the length of the artery have been selected, whereas the wall normal distance () has also been scaled by the lumen height (). At the inlet region, all the types show no generation of TKE. Similarly, at the proximal region of the stenosis, all the types exhibit no TKE generation. Differences start emerging at the peak location of the stenosis where for Type III, generation of TKE is notorious, whereas for Type III minimal TKE generation is present, thus suggesting a much earlier transition to turbulence for Type II and Type III in comparison to the other types. In contrast for Type I and Type IV, generation of TKE remains unseen. At the distal region of the stenosis, Type I through Type III reach a peak generation of TKE, whereas TKE generation for Type IV remains considerably low. For all types, the peak TKE is present in the lower-half region of the lumen. At the exit region, dissipation of TKE is evident for Type I through Type III. In contrast, TKE values for Type IV remain considerably low suggesting that no transition to turbulent state is present. Furthermore, at the exit region, TKE values for Type I are much smaller in comparison to Type II and Type III, suggesting that the turbulence region extents further for these two types. Additionally, the differences and similarities in the turbulence intensity (Tu) profiles for Type I through Type IV are also analyzed during the peak systolic phase of the cardiac cycle. Figure 10 shows the Tu profiles for Type I through Type IV. Similar to the other profiles, seven axial locations were also selected. No significant differences in the Tu profiles are observed in the inlet region of the stenosis. Significant differences are observed in the peak and distal region of the stenosis for Type I through Type III. The conclusions are similar to those observed in the TKE profiles.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Additionally, the viscous and turbulent wall shear stresses (WSS) at the wall along the length of the stenosis have been computed during the peak systolic phase. The viscous and turbulent components in the mean flow direction normal to the wall have been included. Figure 11 shows the viscous shear stresses at the wall along the length of the artery for the four types. For Type I through Type III, the viscous WSS values reach a peak value at their corresponding peak locations in the stenosis, whereas the viscous WSS values for Type IV are evenly distributed along the length of the stenosis. The position of the peak viscous WSS values for Type I through Type III are well correlated with the position where TKE generation starts to grow, thus suggesting that the sudden increment of viscous WSS at the wall is a true indicator of point of transition. Furthermore, the trend followed by the viscous WSS stress values is in complete agreement with the numerical simulations performed by Bhaganagar et al. [17, 27] for the four types under the laminar regime.
(a)
(b)
(c)
(d)
4. Conclusions
Patient-specific simulations with realistic physiological flow conditions are conducted to understand the effect of plaque morphology in altering the flow characteristics in diseased coronary artery. In our previous study Bhaganagar et al. [17], we have developed novel fuzzy classification to classify 42 plaque morphologies into four well-defined types. This classification based on morphology is of great significance for fluid dynamics analysis. Bhaganagar et al. [17] employed steady, laminar flow conditions to demonstrate that the differences in the intertypes flow characteristics are more significant compared to intratype flow. The previous study was restrictive as it assumed (1) laminar flow conditions and (2) steady forcing flow.
In the present study, we extend to realistic physiological flow conditions by accounting for the unsteady flow conditions (systole/diastole) as well as the transition from laminar to turbulent state. We use transition model to predict the location of flow transition to turbulent state. Under these realistic flow conditions, we demonstrate for same degree of stenosis (35%), and both the mean flow conditions as well as the disturbed flow conditions are significantly type-dependent. For the purpose of analysis, we select several axial locations along the length of the artery. The locations correspond to 30–35 mm (inlet region), m (proximal region of plaque), mm (peak location of plaque), mm (distal location of plaque), and 56–60 mm (exit region).
The mean velocity (time and phase averaged) profiles exhibits significant differences between the types. All the types except Type IV exhibit flow reversal at distal portion of the stenosis. Peak velocity for Types I, II, and III is around 1.5 m/s, and it is only 1.2 m/s for Type IV. The velocity profile loses its symmetry along the centerline for Types I, II, and III. However, the location of the peak mean velocity shifts to the upper half of the artery for Types II and III, whereas this shift is towards the lower half of the artery for Type I. Similar trend has been observed during peak systole where Type II and Type III have higher mean velocity compared to the other types. Next, we analyze the flow disturbances. In particular, we are interested in understanding the differences in the location of the transition to turbulence, total stresses, and intensities of turbulent fluctuations between the types. The flow undergoes a transition from laminar to turbulent state for Types I, II, and III, whereas the flow continues to be in laminar state for Type IV. The location of the transition varies between 47 and 51 mm depending on the type. The transition for Type III occurs at much proximal location compared to other types suggesting that flow alters at earlier location for this type.
Peak turbulence intensities are present for Type III, followed by Type II and Type I, respectively. The flow continues to be in a disturbed state close to the exit region for Type II and Type III, whereas the flow relaminarizes to laminar state beyond the distal location of the plaque for Type I. Type III also exhibits high WSS, followed by Type II and Type I, respectively.
Further, the peak viscous WSS values for Type I through Type III are well correlated with the locations where TKE generation starts to grow, thus suggesting that the sudden increment of viscous WSS at the wall is a true indicator of point of transition. This study clearly shows that, for the same degree of stenosis, (a) the presence of turbulence, (b) location of transition to turbulence, (c) turbulence intensity, and (d) region of turbulence are type-dependent. This study is of great significance to determine the risk of rupture of the plaques and thus identify the vulnerable plaques prone to rupture. It should be noted that transition to turbulence is just one of the several biomechanical factors that contribute to risk of rupture plaque (see Fukumoto et al.) [28].
Conflict of Interests
The coauthors do not have a direct financial relation with the trademarks mentioned in the paper that might lead to a conflict of interests for the coauthors.
Acknowledgment
This material is based upon work supported by the National Science Foundation under Grant number HRD 0932339.