|Home | About | Journals | Submit | Contact Us | Français|
We present a computational framework for multiscale modeling and simulation of blood flow in coronary artery bypass graft (CABG) patients. Using this framework, only CT and non-invasive clinical measurements are required without the need to assume pressure and/or flow waveforms in the coronaries and we can capture global circulatory dynamics. We demonstrate this methodology in a case study of a patient with multiple CABGs. A patient-specific model of the blood vessels is constructed from CT image data to include the aorta, aortic branch vessels (brachiocephalic artery and carotids), the coronary arteries and multiple bypass grafts. The rest of the circulatory system is modeled using a lumped parameter network (LPN) 0 dimensional (0D) system comprised of resistances, capacitors (compliance), inductors (inertance), elastance and diodes (valves) that are tuned to match patient-specific clinical data. A finite element solver is used to compute blood flow and pressure in the 3D (3 dimensional) model, and this solver is implicitly coupled to the 0D LPN code at all inlets and outlets. By systematically parameterizing the graft geometry, we evaluate the influence of graft shape on the local hemodynamics, and global circulatory dynamics. Virtual manipulation of graft geometry is automated using Bezier splines and control points along the pathlines. Using this framework, we quantify wall shear stress, wall shear stress gradients and oscillatory shear index for different surgical geometries. We also compare pressures, flow rates and ventricular pressure–volume loops pre- and post-bypass graft surgery. We observe that PV loops do not change significantly after CABG but that both coronary perfusion and local hemodynamic parameters near the anastomosis region change substantially. Implications for future patient-specific optimization of CABG are discussed.
Bypass graft (BG) surgeries require surgical construction of a conduit or graft over a blocked blood vessel. Coronary artery bypass graft (CABG) surgery is usually performed when one or more coronary arteries, that supply blood to the heart muscle, are fully or partially blocked due to the build up of atherosclerotic plaque. Approximately 500,000 CABG surgeries are performed annually in the U.S. Coronary grafts are typically harvested from native tissue such as saphenous vein from the leg, left/right internal thoracic artery or radial artery from the forearm, or harvested from cadavers, or (rarely) composed of synthetic materials such as poly tetra-fluoro ethylene (PTFE). Percutaneous coronary intervention (PCI) with angioplasty (inserting and inflating a balloon in the coronary artery) and stenting (inserting and expanding a permanent mesh structure) are alternative treatments, which may be chosen when progression of disease is less severe, the patient is not a candidate for surgery, or when indicated by other clinical factors.
The operative mortality rate for CABG is approximately 3%. Following surgery, the patency and longevity of a CABG depends on multiple factors including the type of conduit used, the local hemodynamic environment, and structural properties of the walls. It is known that arterial grafts such as internal mammary artery (IMA) and radial artery grafts provide, in general, better patency than vein grafts. Approximately 15–30% of saphenous vein grafts occlude within the first year of surgery, with the rate increasing to over 50% after 10 years.32 Therefore, typically arterial grafts are preferred. However, because multiple vessels are bypassed in most CABG surgeries, vein grafts are usually used as well.
The dominant complications post-CABG are restenosis36 and graft occlusion with associated myocardial infarction,5 and are attributed primarily to deposition of atherosclerotic plaque,11,28 which may subsequently rupture leading to thrombosis, and associated hypo-perfusion. Graft incompatibility, and hemodynamic factors such as blood recirculation, wall shear stress (WSS), and wall shear stress gradients (WSSG) play an important role in the onset and development of plaque (atherogenesis).2,22,23,26,27,38,54 Because these factors are difficult to predict a priori, surgeons use intuition and geometric guidelines to make decisions about graft design, despite the well known links between hemodynamics and restenosis. While mean flow and pressure waveforms can be obtained from invasive coronary catheterization, local hemodynamic information cannot be obtained reliably using currently available methods. Hence, patient-specific modeling and accurate numerical simulations provide a means to obtain data on hemodynamics and WSS that is currently unavailable.15,34,50
It has been shown experimentally (in vivo) that anastomosis angle effects flow fields47 and that flow patterns effect graft patency.19 Using three angles, Rittgers et al.39 did not observe a significant correlation between angle and IH on external iliac arteries for dogs. Using three in vivo patient-specific cases, Giordana et al.12 showed that 80° is better than 30° (since platelet activation is reduced with increased flow mixing) and the low angles studied are predisposed to occlude earlier. These in vivo do not present a standardized or holistic case study to recommend a specific anastomosis angle and further experiments are still needed.13 Earlier computational studies of the influence of CABG on mechanical environment have largely focussed on idealized geometrical models, simplifying physical assumptions, and/or a limited trial-and-error approach to identify BG shapes that improve outcome.1 In our previous work, we studied the influence of anastomosis angles and graft radius on flow in idealized BG models41 demonstrating the effect of geometrical parameters on local hemodynamics. Using regions of low WSS as a measure of undesirable (stagnant) flow patterns, we observed that these regions correlated with three known locations of BG failure (toe, heel, and arterial floor region) and that they could be minimized by changing the graft anastomosis angles.31,41 We also showed that uncertainty quantification42 influences the outcome of surgical optimization. Holzapfel et al.14 showed the interplay between BG shape, suture length, residual stresses and the structural stresses induced by CABG. They showed that small variations of the arterial incision have big effects on the size of the arterial opening and thereby, anastomosis angles and WSS.
Modeling flow in the coronary arteries is particularly challenging due to the influence of ventricular contraction, making traditional boundary conditions inappropriate. Sankaranarayanan et al.43 used constant pressure boundary conditions to study flow in an idealized CABG. This requires knowledge of patient-specific coronary pressure information, which is not usually available from standard clinical data. Pressure boundary conditions are also limited in their capability, since they cannot be used to predict a post-surgical outcome. Dur et al.8 investigated surgical planning for CABG through shape optimization, but their studies were limited to hemodynamic efficiency. Recently, Kim et al.20,21 developed an algorithm that models coronary circulation and microcirculation using a lumped parameter network (LPN). In this work, a set of ODEs describing the coronary LPN was solved analytically to prescribe the relation between flow and pressure at the coronary boundary faces. Inlet boundary conditions were switched from Neumann to Dirichlet when the aortic valve was closed, and the authors employed a penalty method assuming a parabolic flow profile to prevent divergence due to back flow. Here, we build on this previous work by employing an implicit numerical approach to solve the coupled PDE/ODE system. Using a finite element framework (FEM),51 we simultaneously solve the coupled boundary condition problem, modeling the effect of heart behavior and contraction in a closed loop network. This work was also motivated by the closed loop LPN models of Migliavacca and colleagues,18,24,33 which have been applied to model single ventricle physiology. This technique generalizes the computational framework for cases in which analytical solutions are not possible, and eliminates the need to switch between two different boundary conditions (i.e., Dirichlet and Neumann) for the same surface.
In this work, we propose a general cost function for CABG design, but point out that the choice of cost function can change with the acquisition of future experimental data. The use of a formal cost function is meant to set the stage for future work using optimization algorithms for shape design. A cost function can be chosen for BG design using WSS, oscillatory shear index (OSI), gradient oscillatory number (GON), aneurysmal flow index (AFI), or a combination of them. Our cost function is based on a simple linear combination of various risk factors for IH, and further experimental studies are needed to refine this choice. Though exact numerical values that correspond to diseased states are not fully understood, experimental studies can provide guidance on acceptable values.7,29 WSS values between 4 and 15 dynes/cm2 are, in general, acceptable. In coronary arteries, WSS between 10 and 25 dynes/cm2 are considered normal.40 WSSG values of 400 dynes/cm2 or higher are considered abnormal leading to excessive cell proliferation. Values of oscillation in WSS over 5 dynes/cm2 are undesirable, which roughly translates to values of OSI greater than 0.25. Acceptable ranges for other derived variables, GON and AFI, have not yet been experimentally determined. We also develop tools to automatically post-process hemodynamic quantities and test different graft designs for the CABG surgery.
The three goals of this paper are (1) to present a computational framework that could be used (with proper validation) in future work to predict coronary flow based solely on non-invasive clinical measurements, (2) to investigate the dependence of local hemodynamics on anastomosis angle, and (3) to compare the pre- and post-surgical flow conditions using multiscale modeling. The paper is organized as follows. In “Patient-Specific Simulations of Coronary Flow” section, we provide details of patient-specific geometric modeling. We then provide details of boundary condition implementation and post-processing. In “Parameterizing Surgical Geometry” section, we describe a spline-based method to parameterize surgical geometry. Finally, we present patient-specific modeling results.
In this section, we describe a general framework for hemodynamic simulations in CABG patients. The features of our modeling process, in the order performed, are listed below:
We describe each of these steps in the following sections.
We obtained post-operative anatomic data from a clinically indicated CT angiogram (CTA) of a 63 year old male patient (with IRB approval) who underwent a triple BG surgery at UCSF (Fig. 1). The patient received three BGs: two native saphenous vein grafts (one to the right coronary artery and one to an obtuse marginal branch of the left circumflex artery) and one IMA graft (to the left anterior descending artery) (Fig. 2). The saphenous vein grafts originate from two locations in the ascending aorta, distal to the aortic valves. Blood flow was reinstated distal to two stenoses in the left coronary artery (one IMA and one saphenous vein) and one stenosis in the right coronary artery (one saphenous vein).
A 3D iso-contour of the pixel intensities obtained from the image slices is shown in Fig. 2. Three steps were used to reconstruct 3D geometric models from the CTA: (i) creation of centerline paths for vessels of interest including coronary arteries, aorta, brachiocephalic artery, IMA, etc. (ii) segmentation of the vessel lumen through a combination of pixel intensity thresholding and level set method, and (iii) lofting of the segments to construct the 3D geometry. Models are created using a customized version of the open sourced Simvascular software package (simtk.org),48,51 as shown in Fig. 2. The reconstructed arteries are discretized into tetrahedral elements using the commercial software Meshsim (Simmetrix, Inc, Clifton Park, NY).
The strong form of the governing Navier–Stokes equations is given by
where Ω represents the 3D domain, ρ is density of blood, is blood velocity, p is pressure, n is the normal to the vessel wall, τ is stress on the wall, f represents body forces, t is time and x represents space. The Neumann and Dirichlet boundaries are denoted by hΩ and gΩ respectively, and values for traction and velocity are denoted by h and g, respectively.
The corresponding weak form is given by
A stabilized finite element technique using the generalized-α method16,49 is employed with linear basis elements using an in-house custom solver. Walls are assumed to be rigid, though we use capacitance elements in each of the outlets to capture the global effect of flexible walls. Newtonian constitutive behavior of the fluid is assumed, with viscosity of blood set to 0.04 g/cm and density set to 1.06 g/cm3. A fourth order Runge–Kutta time stepping method is used for the 0D model, which is coupled to the 3D domain. We use a recently developed backflow stabilization technique9 to prevent divergence of the numerical solution in the presence of backflow. At the aortic inlet and all the outlets of the model, a pressure (Neumann) BC is prescribed, from the solution to a lumped parameter heart model, as explained in the next section. Pressures in the coronary LPN model are also connected via a closed loop network to the heart model.
We adapt a recently developed coupling framework (with applications in pediatric surgery18) to numerically couple the pressure and flow resulting from finite element solutions at the inlet and outlet boundaries of the 3D model to the 0D lumped parameter model. Using Neumann boundary conditions, pressure is passed from the 0D to the 3D model and flow is passed from the 3D to the 0D model. Previous coupling methods have used an explicit approach, but the semi-implicit method that we use here10 is more robust and has a less restrictive time-step choice. In the semi-implicit method, the tangent matrix in the finite element equation is calculated just once at the beginning of the simulation, as opposed to a fully implicit method where it is updated at each iteration. Numerical stability, accuracy and further details of the coupling technique using explicit, semi-implicit and fully implicit methods are described in Esmaily Moghadam et al.10
The 0D LPN model includes blocks for the four-chambered heart model, the pulmonary arteries, the coronary model accounting for micro-coronary arteries and coronary veins, and RCR circuits that model the downstream resistances and capacitances of all other outlets. The LPN circuit model is governed by a set of ODEs that can be numerically solved, given boundary conditions imposed by the flow at the inlets and outlets of the 3D model.
The heart model circuit has four one-way valves which are modeled using an indicator function, (x) which is 0 if x ≤ 0, and 1 otherwise. The circuit also contains resistors , capacitors (P V), inductances and elastance (P f(V, t)), where P and V represent pressure and volume of blood. For the Runge–Kutta method, the derivatives of variables with respect to time are used to update the variables at each time step.
The set of ODEs governing the closed circuit heart model is presented below. The variables are denoted as Xi and their time derivatives as . Qt is the flow rate through the circuit, P is pressure, R is resistance, L is inductance and C is capacitance. Subscripts ra, rv, la, lv, and p refer to the right atrium, right ventricle, left atrium, left ventricle and pulmonary systems, respectively. Subscript pa represents pulmonary artery, pp and pd represent pulmonary proximal and distal resistances, Caorta refers to aortic capacitance, Qaorta refers to flow in the inlet to aorta, and Ei represents the elastance function which is different for each of the four chambers. The parameter values are given in Table 1.
The governing ODEs at the coronary outlets are
The governing ODEs at the RCR outlets are
A schematic of the LPN is shown in Fig. 3. The Xs refer to either flow rates, volumes, or pressures. In the figure, X2, X4, X7, X9 are the flow variables denoted by arrows, X1, X3, X6, X8 are variables corresponding to ventricular or atrial volume, and the rest (X5, X10) are pressures. Y is used for pressures at the coronary outlets and Z is used for pressures at all the other outlets.
The pressures in the atria and ventricles are computed using the relation Pi = Ei × (Vi − Vi,0) where i denotes the chamber (la, lv, ra, or rv), and Vi,0 is the stress free volume.35 The resistances of the coronary arteries act against the coronary-ventricular pressure difference, while the resistances in the RCR blocks act against the right atrial pressure (which is typically very low). Further, the atrial and ventricular pressures at each time step are calculated as a product of the elastance with the volume of the corresponding vessels.
The net flow into the right atrium Qt is computed from the set of ODEs governing the LPN circuits. The coronary bed boundary conditions are comprised of a coronary arterial circulation and microcirculation (capillaries), and coronary venus circulation and microcirculation that drains the blood into sinus of Valsalva, feeding into the right atrium (Fig. 3). The capacity of the coronary arteries to store blood depends on the ventricular pressure (this is indicated by the line connecting coronary capacitance to left ventricle). This coupling enables us to reproduce realistic coronary waveforms, in which peak flow occurs during diastole, and minimum flow occurs during systole. Traditional boundary conditions such as RCR fail to accurately capture this behavior. Pressure volume curves, pressure tracings, volume and flow rate are shown in Fig. 4.
Since post-CABG clinical data was unavailable for this patient, we assumed normal hemodynamic conditions were achieved post surgery, and fit our parameters accordingly. The elastance function for the left and right ventricles are derived from the normalized curve of Senzaki et al.20,35,44 The values for other parameters of the heart model are also determined from literature values for a normal adult patient.20 It is important to note that if data is available, our computational framework has the capacity to incorporate it, either iteratively or through a rigorous optimization procedure.
The total resistance for all the outlets is . This is split into the total resistance for the coronary outlets (cor) and the rest (br). Since the circuits are in parallel, we have . We prescribe the average flow split between the coronaries and the aorta. Let this fraction be represented as β, i.e., Qcor = βQbr. The net coronary and branching resistances are then shown to be and Rbr = (1 + β)Rtotal. In this work, β was chosen to be 0.04 based on typical literature values.4 Capacitance does not affect the mean flow and is tuned to adjust the range and variation of the pressure values. The capacitance values and flow rate are chosen to be proportional to the outlet areas,53 Ci Ai. We first fix the capacitance of the aortic outlet and then compute the proportionality constant to fix the capacitance of all other outlets. The capacitance of the aortic outlet is iterated on, using trial and error, to match the typical blood pressure of 120/80. For instance, a higher capacitance value will shift the peak towards the right.
The net resistance of the branches, Rbr, is split for each outlets as follows. Given the net resistance and area of each branch,21 we have Ri = k/Ai where k = RbrAbr,total. The resistance of each branch is then Rp,i = kp × Ri and Rd,i = kd × Ri, where subscripts p and d refer to proximal and distal respectively. The constants kp and kd are 0.09 and 0.91, respectively.20 The RCR values for each of the outlets are given in Table 2.
At the coronary outlets, we split the resistance for the left and right coronaries to match literature data.17 Given a flow split γ between left and right coronaries, we have . We choose γ = 7/3 to enforce a typical 70–30% flow split between LCA and RCA.17 Hence, . The resistance for each branch is then Rcor = Ra,cor + Ra−μ,cor + Rv−μ,cor + Rv. We choose typical values for splitting the coronary resistance,6,20 Ra,cor = 0.32Rcor, Ra−μ,cor = 0.52Rcor and Rv + Rv−μ = 0.16Rcor. The values of the coronary BC parameters are given in Table 3.
The capacitance values for each outlet are chosen to be proportional to their area. For each outlet, the capacitance value is split according to typical values Ca = 0.11Caorta and Cim = 0.89Caorta.
In this section, we explain the methodology used to compute the WSS and WSSG.
Wall shear stress is computed in a post-processing step from the velocity () and pressure () fields of the Navier–Stokes solution. Rearranging the weak form in Eq. (1), with velocity and pressure solutions and = ∑A NA, we have:
The weighting function, , is non-zero only at the surfaces where the WSS is to be computed. Choosing to be the surface shape functions, we solve for the boundary tractions. Two additional quantities of interest are derived from the WSS solution, τ(x, t). First, OSI is defined as
and is an indicator of the directionality of flow. OSI has minimum value of 0, corresponding to unidirectional traction, and maximum value of 0.5, corresponding to zero mean traction with equal time spent in both directions. Second, we compute the AFI,30 originally developed to measure likelihood of aneurysm formation, by quantifying the temporal fluctuations of the direction of WSS. The AFI is defined as the instantaneous projection of the wall shear direction on the mean wall shear direction. Previous studies support the relationship between endothelial response and temporal fluctuations in WSS.37 Since it is well known that endothelial cells orient in the direction of mean WSS,30 the AFI quantifies the time averaged tensile/compressive force that endothelial cells are subject to.
Wall shear stress gradients quantify the spatial variations in hemodynamic forces on the vessel wall.25,37 WSSG quantifies the magnitude of shear-stress gradients using the average wall shear direction and its in-plane (plane of the wall) normal. Off-diagonal coupled-derivatives, such as gradient of mean wall shear in a direction normal to the wall, are neglected since they have been shown to have little clinical relevance.25 The formulation for computation of WSSG is given below.
From the velocities, tractions are defined as ti = σijnj and shear stresses as τi (x, t) = ti − (tk nk) ni. Let s1(x) be a unit vector in the direction of mean wall shear traction and let s2(x) represent a direction orthogonal to s1 in the plane of the wall (i.e., nw is the normal to the wall, s2 = nw × s1). The WSS is decomposed into two components, (τs1, τs2) = (τ · s1, τ · s2) and the WSSG is defined as25,38:
where the individual terms are computed as Finally, the mean WSS gradient is computed as:
The GON45 is defined from the WSSG as:
where G = siτsi.
Details of how the GON is computed are given in the Appendix.
We parameterize the surgical graft shape using cubic splines. There are three BGs in the patient. Initially, we identify the center pathline of each vessel of interest and identify a set of control points along the pathline (7 for this patient). These pathlines are parameterized by a variable t which varies from 0 to 1 along the length of the path. Piecewise cubic splines were chosen to interpolate the paths, after second order polynomials were found to be inadequate. To compute the coefficients, the spline is made to pass through all the control points and the slope is made continuous at the control points. The number of equations is one less than the number of control points.
Different BG shapes are obtained by manipulating the location of the control points. Let r0 represent the location of the control point for the original BG. We specify a maximum deviation dm as well as a direction n and define the possible new shapes of the BG as r = r0 + λdmn where λ is between −1 and 1 and represents the extremes of the BG shape. The direction of motion n is chosen perpendicular to the pathlines (at the control point) in an arbitrary plane.
Let t be the parameter that represents distance along a saphenous vein graft and s be the parameter that represents distance along a coronary artery. The anastomosis angle between these arteries is given by represents the tangent to the centerline path of the saphenous vein and s′, which is defined similarly, represents tangent to the coronary artery. A different anastomosis angle can be tested by modifying the control points. The endpoints are fixed, and if the control point closest to the anastomosis is moved by a value λ normal to the spline, the rest of the control points move by where λk is a vector normal to the cubic spline (i.e.) λk = λ × nk, where nk is the normal to the spline. A plot of λ vs. anastomosis angle θ is shown in Fig. 5.The plot shows that a full range of realistic anastomosis angles can be obtained by this parameterization.
For the patient-specific model, parameterization of the BGs is shown in Fig. 5. A set of seven control points is used for each of the BGs. Different anastomosis angles are then obtained by manipulating the location of these control points. We automated the process of building a geometry given a trial anastomosis angle. Through a simple inversion process, we obtain the location of the control points for a given trial anastomosis angle. This is then converted into a new pathline, and segmentations constructed from CT are translated accordingly. These segmentations are then lofted together and blended to form the new BG geometry.
To ensure convergence, we performed mesh-independence and time-step independence tests. We gradually refined the mesh density and compared pressure, velocity, energy, and WSS. Figure 6 depicts the change in simulation results with mesh density. We conclude that the mesh density with 4 million elements is sufficient for our simulations.
A patient-specific measured heart-rate of 61 beats per minute is used for the cardiac cycle. We used a time step size of 0.5 ms for the multiscale model to ensure convergence. We used 4 Newton–Raphson linearization iterations per time step for all the models. We simulate flow in the first model geometry for 10 cardiac cycles to ensure convergence, and then perform subsequent simulations for other geometries using the first case as the initial condition. Only 3–4 cardiac cycles are subsequently needed for the results to stabilize. When we restart the simulations, we use 5 Newton–Raphson iterations for the first 100 time steps and then continue with 4 Newton–Raphson non-linear iterations, to ensure convergence.
The multiscale model results in a coronary flow peak at diastole, agreeing with clinical observation, as expected. The mean flow ratio matches the ratio of the prescribed resistance between the coronary and aortic circulation, as expected. The flows were all measured distal to the anastomosis and branching locations in five arteries in the left coronary system (one diagonal, LAD after diagonal bifurcation, two obtuse marginals, LCX after second obtuse marginal bifurcation) and one artery in the right coronary system. We iterated the total capacitance of the coronary circuit to match the position and ratio of pressure peaks to match typical clinical data. An LCA-RCA flow split of 70–30% was achieved, as shown in Fig. 7.
Using the coupled multiscale model, we compute the flow conditions pre- and post-CABG surgery and compare quantities of interests, keeping all LPN boundary condition parameters fixed. The coronary arteries have a mean and maximum Reynold’s number of ~424 and ~256 respectively, and Womersley nomber of ~4.71. In the aorta, maximum Reynold’s number is ~3200 and Womersley number is around ~12.53. As shown in Fig. 9, coronary flow is partially restored after the CABG is performed. Prior to CABG surgery, we assume that the bypassed right coronary artery and left obtuse marginal have a 75% stenosis and create the corresponding model keeping the rest of the patient-specific geometry intact. This changes the effective hemodynamic resistance of the model. We then measure how the cardiac output and other quantities of interest in the LPN and the 3D model change with inclusion of BG resistance.
Pressure-volume curves of the left and right ventricle before and after the surgery are compared in Fig. 8. As expected, the cardiac output is not significantly increased by the reduction of hemodynamic resistance induced by the BG. There are two factors that contribute to this. First, the coronary flow is a very small fraction of the full aortic flow, and second, the model resistance is dominated by the rest of the circuit, implying that addition of the BG results in insignificant changes to the net resistance.
We compare the coronary perfusion (measured distal to the BG) before and after CABG surgery in Fig. 9. After CABG surgery, blood is drawn through the BG which offers lower resistance than the stenosed coronary arteries. The figure clearly shows an improvement in the coronary flow-rate, and hence perfusion, after the BG is inserted in the two stenosed coronary arteries. We also confirm the diastolic flow peak in the coronary arteries both pre- and post-surgery. Prior to implanting the BG, the model resistance of the stenosed coronary arteries dominates the downstream boundary conditions. After the surgery, the BGs impose a lower resistance to flow and hence provides a conduit for almost all coronary flow in the stenosed arteries. The flow in the remaining coronaries is unaffected.
In addition, we changed the BG radius to half its original value, keeping the rest of the geometry intact. We observed that changing the graft radius in this range had little effect on the cardiac output or coronary flow. This is because downstream coronary artery resistance is still much greater than the BG, which has a very small contribution to the overall coronary resistance. However, dropping the radius below this value, we start to see an effect on coronary flow since the resistance varies with the radius.4
Figure 5 shows the shape of the BG for values of θ = 70°, 62°, 50°, 35°, 10°. A wide range of physically realistic anastomosis angles is achieved with this parameterization. In this section, we restrict our discussion to the parameterization of the saphenous vein BG, which is anastomosed to the left circumflex, LCx.
Volume renderings of the velocity magnitude at diastole for three BG shapes, θ = 70°, 50°, 10° are shown in Fig. 10. Coronary flow is maximum in diastole. The figure shows that there is a change in the aortic flow resulting from different bypass shapes, while the flow in the non-stenosed coronary arteries remains almost the same. Figure 11 shows the difference in flow characteristics for different graft angles at the anastomosis region. A lower anastomosis angle implies a longer suture region, and hence we observe that (a) a higher bypass angle leads to rapid flow mixing, and (b) low BG angles ensure a streamlined mixing of blood at the location of the anastomosis and larger areas of low WSS. Plots of WSS and OSI for five different BG shapes are shown in Fig. 12. The WSS plots reinforces earlier observations27,28 showing that regions of critical hemodynamic importance in BGs are the toe, heel, and arterial floor regions. High anastomosis angle results in reduced energy efficiency and larger flow separation regions. Lower anastomosis angles, however, result in higher areas of low WSS at the heel and toe regions of the BG.
Figure 13 shows the variation of different quantities of interest with the BG parameterization. All quantities are normalized to a value of 1 based on the maximum value across the five designs. Areas of low wall shear, OSI and GON have the same minima corresponding to θ = 70° (an anastomosis angle of roughly 70°). Overall, OSI and AFI are the least sensitive with respect to shape. If we consider the sum of heel, toe and anastomosis floor regions, the higher impact angle (θ = 70°) has the optimal configuration when considering area of low WSS. This is because the arterial floor region contributes most to the total area of low WSS, and low WSS in this region is minimized in designs with higher anastomosis angle. The trend in GON is the result of two competing factors: a higher anastomosis angle results in (1) a higher contribution to the GON via WSSG, and (2) a lower contribution to the GON via variance in the WSS direction. The plot shows that the latter term is dominant, resulting in the optimal GON at θ = 70°. If we consider a cost function defined as = 0.25 × (A(τ ≤ 1) + OSI + GON + AFI), we conclude from Fig. 13 that θ = 70° is the optimal configuration. The first term indicates the area of the lumen wall that is subject to low shear stress. The appropriateness of this choice of cost function and determination of whether the cost function and resulting optimal angle should be patient-specific will be the subject of future research. The optimal bypass shape ultimately ought to involve a detailed combinatorial optimization involving several different clinically motivated cost functions and constraints.
To test the sensitivity of simulation outcomes to the total coronary flow (which we assumed was 4% of the aortic flow), we performed two additional sets of simulations with total coronary flow set to 3% and 5% of the aortic flow. The resistances were adjusted accordingly, and the total flow rate was maintained. We observed that the coronary flow profile shifted in magnitude by ±25% keeping the shape of the flow waveform shape the same. However, we observed a change of less than 1% in the cost function values since the cost function values are generally normalized by their mean values.
We expanded upon an implicitly coupled multiscale framework to perform CFD simulations in a patient who underwent CABG. The described method provides a means to perform realistic simulations of flow in the coronary arteries and BGs by modeling the entire circulatory system as a closed loop. This technique is useful for virtual surgical parameterization in which the surgical geometry could affect the boundary conditions, and this is implicitly captured by the described technique. When available, clinical data was incorporated into our computational framework and nominal literature-based parameter values were assumed otherwise.
We used cubic splines with seven control points to parameterize the BG shapes, and this method could potentially be used in the parameterization of other surgical procedures, or in a formal optimization procedure. We are able to virtually manipulate surgical shape, anastomosis angles and correlate hemodynamics with different surgical shapes.
Our framework reproduces the clinically observed coronary flow peak during diastole and also incorporates data such as flow split between coronary arteries and aorta, when available. We observe 3D flow profiles at the anastomotic region for high anastomosis angles. We observed that performing the BG surgery increases coronary perfusion, as expected. The increase in flow rate in the stenosed arteries is observed to be as much as 200%. While previous studies have also examined the effect of geometry on BG hemodynamics,8,43 this was the first to do it using only CT and non-invasive clinical measurements, without the need to assume pressure and/or flow waveforms in the coronaries.
We observed that local variables such as WSS and WSSG are affected by the BG shape, and different cost functions have different optima. Local velocities in the aorta are also observed to be affected by changes in BG shape. Global quantities such as energy efficiency and PV loops are unaffected by the resistance of the inserted BG and its shape. In patients with impaired left ventricular systolic function prior to surgery, these quantities will also depend on the ventricular contractility, which often improves after BG surgery. Our computational framework has the potential to model this through a change in elastance function, whenever such data is available.
We conclude that BG shape has a strong influence on local quantities of interest that would potentially affect patient longevity, and that different hemodynamic cost functions have different optima. Normalizing each cost function to a value of 1 and assuming that each cost function is weighted equally, we observe that for the range considered, an anastomosis angle of 70° would be hemodynamically optimal for this patient. Hence, our results provide further evidence that performing patient-specific optimization on the surgical geometry should be considered and could lead to improved longevity of the graft.
Limitations of this work include the use of rigid walls, Newtonian flow, and assumptions of coronary flow distribution. In a recent study,46 the difference in WSS between FSI and rigid wall simulations in coronary arteries was found to be less than 5%, however studies in other arteries have found larger differences.3 The physical location of the surgery prevented us from exploring higher anastomosis angles without adding additional geometric parameters. It remains to be seen how the choice of cost function presented here impacts long term clinical outcomes, or if different patients need different weighting functions based on a priori clinical factors. In the future, we plan to perform a detailed study of flow in coronary arteries in multiple patients and perform statistical analyses to test the generality of our observations. We also plan to perform patient-specific optimization on different patients to test whether a one-size-fits-all solution is sufficient or if different patients should have different optimal BG shapes.52 Further validation with post surgical clinical data is needed. The effects of uncertainty in simulation inputs on their predictive capability should also be analyzed using formal uncertainty quantification tools.42 We also plan to incorporate the effects of exercise and autoregulative mechanisms of the heart.
The authors acknowledge funding from an American Heart Association postdoctoral fellowship, the Leducq Foundation, a Burroughs Wellcome Fund Career Award at the Scientific Interface, and NIH grant RHL102596A for funding this work. We thank Hyun Jin Kim for expertise on coronary modeling, and Francesco Migliavacca for expertise on multiscale modeling. Dr. Guccione acknowledges funding from the National Institute of Health, grant numbers R01HL077921 and R01HL086400.
Here, we provide details on computing the WSS gradients used in this study. For each element e, let the normal be denoted by ne and the mean shear stress be denoted by . The direction is computed as . The components of elemental shear stress in these two directions are computed by taking a dot product with the directions. We can compute the gradients as follows:
The gradient of the shear stresses are evaluated at each tetrahedral element as:
The corresponding MWSSG is represented as . We then derive the nodal values as
Discretizing this equation and choosing the test functions as shape functions, we have
where A denotes the assembly matrix. Integrating this equation, we observe that:
This equation is found from the least squares method (minimizing the difference between τG and τG,e). For the purpose of visualization, we approximate the LHS system to choose nodal values that are weighted averages of the elemental solutions.