Received 1 April 2016; accepted 27 May 2016; published 30 May 2016
Cardiovascular disease is a leading cause of mortality and morbidity in developed countries  . Narrowing and obstruction of blood vessels resulting from thrombus or plaque deposit, called stenosis, lead to a decrease in blood flow, complex flow patterns, and abnormal wall motions  . A coronary artery surrounding the heart muscle plays a significant role in providing the heart with oxygen and nutrients. The dysfunction of the coronary artery due to the stenosis is a major cause of serious cardiovascular diseases, such as ischemic stroke or cardiac infarction  .
So far, there have been numerous studies to investigate the hemodynamics of the stenotic arteries for the treatment of the related diseases. Significant advances in the understanding of the relationship between hemodynamics, mechanical factors, and atherosclerotic changes have been made over the past few decades. Many of the studies have revealed that the wall pressure, wall shear stress (WSS), or pressure drop of the stenotic blood vessel (between the inlet and outlet) become much larger than the normal blood vessel   -  . In contrast, some studies found that low WSS and wall pressure is a significant indicator for predicting plaque-prone site or atherosclerosis   . One of the most effective ways to treat blood clots or plaques in the stenotic area is systemic administration of thrombolytic drugs or clot-lysing agents into the blood vessels. Korin et al. developed a therapeutic strategy for the stenosis treatment by employing micro-sized particles composed of aggregates of nanoparticles (PLGA) through a microfluidic device  . They adopted the basic principle that high shear stress in the stenotic area can disassemble the microparticles into nanoparticles, which subsequently are attached to the target stenotic area. Vijayaratnam et al. investigated hemodynamic flow and drug transport of Newtonian and non-Newtonian models using Paclitaxel, an antiproliferative drug, in stented arteries  . Stangeby et al. utilized low density lipoprotein (LDL), which is one of the major causes for stenosis, to examine the mass transfer in a stenotic artery  .
Still, there are several challenges in the blood flow and drug delivery of the blood vessels. First, the dose of drugs as a therapeutic strategy should be performed with caution. This is because long-term infusion of drugs can cause the potential risk of bleeding  . In addition, the blood flow is pulsatile flow with periodic variations, so it is not easy to accurately measure the transport of the drug and estimate the resulting efficacy in vivo with the varying time  . To address these difficulties, this study tested a hypothesis with fluid structure interaction (FSI) numerical simulation that if we are able to find the critical injecting time point for the best performance of the drugs where they deposit most in the stenotic area, the cost, time, and potential risk for the clinical therapy can be reduced.
Motivated by these facts, the aim of this study was to: 1) develop a computational model of blood flow for the stenotic right coronary artery (RCA), 2) investigate the effect of wall compliance on the drug transport in the vessel using fluid-structure interaction (FSI), and 3) determine the injecting time point for the best effectiveness of the drug delivery to the stenotic sites. First, the effect of Newtonian and non-Newtonian viscosity models on velocity profiles in the five major locations was examined. Second, the area-averaged pressure and velocity profiles between rigid and flexible RCAs at the throat of the stenosis were compared. Third, the particle transport was simulated using rigid and flexible RCAs for a range of times, and the effective injecting time point was discussed. Finally, limitations of this study were addressed.
2. Computational Methods
2.1. Reconstruction of Right Coronary Artery
A three-dimensional (3D) stenotic RCA was reconstructed based on a 62-year-old patient as reported as shown in Figure 1  . The normal diameter is about 3.21 mm, the diameter of stenotic area is 1.43 mm, and the wall thickness is 1 mm   .
The stenosis severity is estimated to be approximately 80% by the following form  :
where As and An are the area of the stenotic cross section (occluded at the throat of the stenosis) and normal cross section of the RCA, respectively.
The hybrids composed of hexahedral and tetrahedral mesh were generated. In addition, since the numerical results should be independent of the number of mesh, the mesh independence test was performed for five cases as shown in Table 1. The results showed that the maximal difference in averaged velocity and pressure at the throat among different mesh densities is found to be less than 3%.
Figure 1. The RCA geometric model reconstructed from clinical angiogram.
Table 1. Area-averaged velocity and pressure at the throat of the stenosis with the different mesh densities.
2.2. Mathematical Equations for Two-Way FSI
2.2.1. Fluid Model
The Newtonian fluid model was considered and the flow was assumed to be transient, isothermal, incompressible and laminar due to low Reynolds number (Re = 2678) at the peak flow. The Womersley number (α = R(w/ν)0.5) is approximately between 2 and 5 depending on vessel diameter, viscosity model, and blood density. For fluid mechanics computation, the Navier-Stokes equations with arbitrary Lagrangian-Eulerian formulation were employed as follows:
where ρf is the blood density (1050 kg/m3)  , u is the fluid velocity, p is the pressure, ug is the moving coordinate velocity and µ is the dynamic viscosity. Thus, u − ug is the relative velocity of the fluid with respect to the moving coordinatevelocity.
2.2.2. Solid Model
It has been reported that the arterial wall, composed of layered collagen, is dynamic and compliant during the cardiac cycle  . Furthermore, it is a heterogeneous and anisotropic material. Based on this fact, a nine-para- meter Mooney-Rivlin hyperelastic model was adopted to describe the behavior of the arterial wall  . The strain energy function is shown with the following form:
where c10 = 0.07 MPa, c20 = 3.2 MPa, c21 = 0.0716 MPa, and cij = 0 MPa. Ī1 and Ī1 are the first and second deviatoric strain invariants, dm (2/K, K is bulk modulus) is the material incompressibility parameter (10−5 Pa−1), and J is the determinant of the elastic deformation gradient tensor.
The momentum equation with the ignoring gravity for the wall segment of the RCA based on Lagrangian coordinate system is given by Equation (5).
σs is the solid stress tensor, ρs is the arterial wall density (1366 kg/m3)  , and Us is the solid displacement.
2.2.3. Physiological Boundary Conditions
As for the boundary conditions, the time-varying velocity and pressure wave were applied at the inlet and outlet, as previously reported  . Figure 2 shows the inlet velocity (peak flow at 0.15 s, red solid line) and outlet pressure profiles (peak pressure at 0.4 s, blue dash line).
At the FSI interface, the following boundary conditions are applied: 1) displacements of the fluid and solid domain are equal, 2) traction forces are compatible, and 3) the non-slip condition is chosen for the fluid:
whered, σ, and n are displacements, stress tensors, and normal vectors with the subscripts f and s for fluid and solid, respectively. The inlet and outlet ends were set to be fixed.
2.2.4. Particle Transport Model
To describe the drug transport in the RCA, the Lagrangian particle transport equation was employed:
where mp, dp, and CD are the particle mass, particle diameter, and drag coefficient. The used particle, composed of aggregates of multiple nanoparticles, was taken to be micro-sized PLGA with 3.8 µm of diameter and 1.34 g/cm3 of density, respectively  . The term on the left-hand side is a summation of all of the forces acting on the particle expressed in terms of the particle acceleration. The first term on the right-hand side is the drag force, the second term is the force exerting the particle due to the pressure gradient, and the third term is the force to accelerate the virtual mass of the fluid in the volume occupied by the particle. The number of injected particles is
Figure 2. Physiological boundary condition with inlet velocity and outlet pressure. Inlet peak velocity (red solid line) appears at 0.15 s in systolic phase and outlet peak pressure (the blue dash line) is shown at 0.4 s in diastolic phase.
100 at each four-time point (0 s right before acceleration, 0.15 s for inlet peak velocity in systole, 0.4 s for outlet peak pressure in mid-cycle, and 0.7 s in diastole) during the cardiac cycle.
2.3. Numerical Methods
The fully coupled solid and fluidmodels were solved interactively with 0.01 s of the time step and 1 s of the total time (i.e., one cardiac cycle) based on a heart rate of 60 beats per minute using a commercial program ANSYS for structural mechanics analysis and ANSYS-CFX for fluid mechanics analysis. The fluid pressure at the fluid- solid interface can be applied as a load for the structural analysis, and the resulting displacement, velocity, or acceleration obtained in the structural analysis can be passed back for fluid analysis. As a solver scheme, a high resolution for advection and second order backward Euler for the transientscheme were utilized, and the convergence criterion was set as 10−4. Iteration between the structural analysis solution and fluid mechanics solution continues until overall equilibrium is reached.
3. Results and Discussion
3.1. Comparison of Velocity Profiles of Newtonian and Non-Newtonian Models
In reality, the blood has been known to be a non-Newtonian fluid where the blood viscosity changes with the hematocrit and macroglobulin concentration  , but many studies have assumed the blood as a Newtonian fluid under the high shear rate (i.e., shear rate > 100 s−1)  . As mentioned previously, I assumed the blood as a Newtonian fluid for the FSI analysis. In order to justify this assumption, it should be guaranteed that physical quantities such as velocity profiles, shear rates, or wall shear stress between the Newtonian and non-Newtonian models should be similar. For this, I employed various blood viscosity models for examining their effect on the hemodynamics. Table 2 shows several blood models for the effective viscosity. The non-Newtonian models include power law, Carreau, Casson, and Walburn-Schneck. As a Newtonian model, plasma and blood with the constant viscosity were utilized. For the comparison between viscosity models, the rigid RCA wall was used.
Figure 3(a) shows five major areas of the stenotic RCA to extract velocity profiles: in the upstream (L1), in the pre-stenotic region (L2), at the throat of the stenosis (L3), in the post-stenotic regions (L4), and in the downstream (L5). The diameter of each location is 3.47 mm (L1), 3.16 mm (L2), 1.43 mm (L3), 2.61 mm (L4), and 2.50 mm (L5), respectively. Figure 3(b) shows velocity profiles at three time points (0.15 s, 0.4 s, and 0.7 s) for each location. Overall, the velocity magnitude at the throat of the stenosis (L3) is highest for all the time points, regardless of blood models. In addition, at the peak flow in systolic phase (0.15 s), the maximum velocity values (around 4 m/s) are approximately two times as high as those (around 2 m/s) of other time points (0.4 s and 0.7 s). Noticeably, at the post-stenotic region (L4) and the downstream (L5), I found that the pattern of the velocity profiles were different from other locations (L1, L2, and L3). The velocity profiles are highly skewed due to the curved shape of the RCA and the resulting higher inertial force of the fluid.
Table 2. The Newtonian and non-Newtonian blood models for the effective viscosity.
Figure 3. Comparison of Newtonian and non-Newtonian models: (a) Five different locations to extract velocity profiles: in the upstream (L1), in the pre-stenotic region (L2), at the throat of the stenosis (L3), in the post-stenotic regions (L4), and in the downstream (L5). (b) Velocity profiles with three time points for five locations.
For the comparison of the viscosity models, the Casson model predictsslightly higher velocity magnitudes than others, especially for L1, L2, and L4, whereas the plasma model predicts relatively lower values at the same locations. Among non-Newtonian models, the Walburn-Schneck model tends to underpredict the velocity magnitudes, particularly for L1, L2, and L4 compared to others. Overall, only subtle difference in velocity mag- nitudes among the viscosity models―except for the plasma model and Walburn-Schneck―was observed for all locations. It should be noted that the Newtonian blood model in this RCA model presents relatively similar velocity patterns and values for all times and locations to other non-Newtonian models. For example, for the comparison of the blood Newtonian model to non-Newtonian models except for the Walburn-Schneck model, the maximum difference is found to be within 4%. In light of the result, it is believed that applying the Newtonian blood model to FSI analysis (i.e., flexible wall) does not have a significant impact on the hemodynamic results, compared to non-Newtonian models.
3.2. Total Displacement and von Mises Stress of Flexible Wall
Figure 4 shows a total displacement of the flexible RCA wall. The extracted total displacement is a quantity of the movement at the wall interface at each time relative to the initial state. Overall, the total displacement was
Figure 4. Total displacement of the flexible RCA wall. The center of the stenotic area at the initial state is expressed as a black cross.
found to be highest at the stenotic throat. The maximum displacements for 0.15 s, 0.4 s, and 0.7 s were shown to be about 1.67 mm, 2.14 mm, and 1.46 mm, respectively. These values relatively agree with data reported in the other study by showing the same order of the magnitude  .
The maximal von Mises stress was found to be 92.06 kPa at the downstream away from the throat. This is less than 10% of the breaking strength of the RCA wall, which has been known to be approximately 1 MPa  . However, the von Mises stress at the throat yields a very low value.
3.3. Comparison of Pressure and Velocity Profiles between the Rigid and Flexible RCA Wall
Figure 5 shows the comparison of temporal averaged velocity (a) and pressure profiles (b) at the cross-section of the stenotic throat (L3). It is presented that velocity and pressure patterns between the two models are similar but there exists quantitative difference. First, the averaged velocity magnitudes for the flexible wall are smaller than the rigid wall for most of the time points. In particular, the degree of the difference in velocity is found to be largest (15%) at the inlet peak flow (0.15 s). This is because the high fluid pressure force (due to the peak flow in the flexible RCA) is transferred to the wall and subsequently causes wall movement and deformation. Correspondingly, there are more spaces where fluid can move and thus, the overall velocity magnitude might be lower. On the other hand, there is very small difference (less than 5%) in velocity at the 0.43 s. This is thought to be because the low inlet velocity conditions lead to the small fluid pressure exerting on the wall, which does not much affect the wall movement or deformation. Consequently, this might lead to no significant difference between the two models at this time point. In contrast to the case of the velocity, the pressure values for the flexible wall are slightly larger, but the difference in pressure between the two models is on average smaller.
Figure 5. Comparison of averaged velocity and pressure profiles between the rigid and flexible RCA wall at the cross- section of the stenotic throat (L3): (a) velocity and (b) pressure profiles, respectively.
3.4. Comparison of the Particle Behavior between Rigid and Flexible RCA Wall
Figure 6 shows the comparison of the particle behavior within the rigid and flexible RCAs using four different injecting time points (tin = 0 s, 0.15 s, 0.4 s and 0.7 s). The red color indicates trajectories of deposited particles and the green indicates trajectories of particles exiting from the RCA. This qualitative result clearly shows that the number of particles exiting (green) is higher than deposited particles (red) within the RCA.
Such a trend was shown in quantitative results (Figure 7) as well. Figure 7(a) shows one of the sample images of the deposited particles in stenotic (yellow oval) and non-stenotic area. Figure 7(b) shows quantitative results of particle deposition in the RCAs. As in the case of the qualitative results in Figure 6, it was clearly observed that the overall deposited particles were less than 50% of all injected particles, regardless of rigid and flexible wall. In other words, more particles tend to escape the RCAs than deposit. The ratio of total particle deposition (blue color) in the rigid RCA wall is 43%, 37%, 25%, and 23% for 0 s, 0.15 s, 0.4 s, and 0.7 s, respectively. In the case of the flexible RCA wall, the total particle deposition (red color) of 41%, 29%, 24%, and 25% for the same time points above was shown. The result clearly shows that the number of particles deposited in the rigid wall is generally higher than the flexible wall, except for 0.7 s. In particular, at 0.15 s (inlet peak flow), there was the most notable difference in particle deposition between the rigid and flexible wall showing 8% discrepancy. On the other hand, the difference was not more than 2% for other cases. When considering the number of particles deposited only in stenotic areas, the difference between rigid and flexible models became substantially lower. Such a small difference might be associated with the fact that dynamic vessel motion on flow behaviors is less significant than other factors such as arterial geometry or pulsatile flow condition  . Similarly, the flexible wall might not affect the particle behaviors that much either. This result implies that for the prediction of drug transport within the RCAs, rigid wall model as well as flexible wall model can be one of the options.
Interestingly, for the time period up to 0.15 s as shown in Figure 7(b), the number of particles deposited in the stenotic area was larger than the non-stenotic area for both rigid and flexible wall. On the contrary, from the time point of 0.4 s, the trend became inverse. The most notable feature is that the injecting time of 0 s was most effective for the particle deposition for both two models. During this injecting time, the particle deposition in the stenotic area is found to be highest presenting approximately 24% for both the rigid wall and flexible wall. This is presumed due to the fact that, from the initial state (0 s) to peak flow (0.15 s), fluid flow tends to be accelerated, whereas fluid decelerates after then as shown in the inlet velocity condition of Figure 2. During this time period, the maximum particle traveling time is shown to be slightly less than 0.2 s. Thus, the injected particles at 0 s might be strongly affected by higher inertia force due to high velocity magnitudes and thus deposit more, compared to other time points.
Figure 6. Comparison of particle deposition rate of the rigid and flexible RCA wall. The red color indicates trajectories of deposited particles and the green is trajectories of particles exiting from the RCA.
Figure 7. A sample image (a) and quantitative results (b) of particle deposition of the rigid and flexible wall: stenosis, non- stenosis, and total indicate particles deposited in the stenotic area, non-stenotic area, and both in the stenotic and non-stenotic area.
In this study, the numerical simulations of the curved stenotic RCA were quantitatively performed in terms of velocity profiles, wall displacement, maximal von Mises stress, and particle deposition in the rigid and flexible RCAs. I employed the real pulse wave for the inlet and outlet condition to make the simulation closer to the actual physiological phenomena. In addition, I justified the utilization of the Newtonian blood viscosity model by comparing velocity profiles in five major areas of the RCAs with non-Newtonian viscosity models.
Since the characteristics of blood vessel are highly complicated, modeling a 3D realistic blood vessel is essential for accurate surgical prediction. Complex flow such as turbulent flow cannot be precisely described in two- dimension (2D). In contrast, our 3D model and simulation of deformable blood vessel and drug transport enable to describe more realistic phenomena, and thereby can give exact information and feedback to the patient. In addition, it can reduce the amount of time for model validation.
Nevertheless, several limitations still exist that need to be addressed. First, only a specific particle size (3.8 µm, microparticle composed of aggregates of nanoparticles) was adopted as reported  . However, in reality the size of the drug varies, ranging from submicron to micron, and thus different-sized particles should be employed to investigate the size effect. In particular, submicron particles tend to moverandomly in a fluid, Brownian motion, so the related force (i.e., Brownian force) should be taken into account for the particle transport in this case. In addition, I do not consider surface property or signaling factor for the attachment of the drugs. Finally, the simulation was carried out based on the patient-specific size and shape of the stenotic RCA. Thus,varioussized and shaped RCA models of different patients should be tested to guarantee the reliability of the efficacy of the drug. These addressed factors will be the primary goals in our later research.
I have developed a computational model with FSI for blood flow and drug delivery of the curved RCA with stenosis. The key points in this research are as follows:
1) The Newtonian blood viscosity model used in this study shows similar flow patterns and velocity magnitudes to non-Newtonian models.
2) The maximal total displacement and von Mises stress of the RCA wall are 2.14 mm and 92.06 kPa, which are both on the same order of the magnitude reported in the other study.
3) There is no notable difference in the particle sedimentation to the stenotic areas between the rigid and flexible RCA models.
4) The effective time point for the drug injection is found to be between 0 s and 0.15 s (i.e., fluid acceleration region).
I believe that this information will significantly contribute to the design of a drug delivery system for the treatment of the stenotic arteries by targeting drugs selectively to the stenotic sites.
This work was done in the mechanical engineering computer laboratory (MECL) at Purdue University.