Patient-Specific Multiscale Modeling of Blood Flow for Coronary Artery Bypass Graft Surgery

Abstract

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.

Appendix

Here, we provide details on computing the WSS gradients used in this study. For each element e, let the normal be denoted by n ^e and the mean shear stress be denoted by $\bar{\tau}^e = \,\parallel \bar{\tau} \parallel{s^{e}_1}.$ The direction s ^e₂ is computed as $s^{e}_2=\frac{s^{e}_{1} \times n^{e}}{\parallel s^{e}_{1} \times n^{e} \parallel }.$ 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:

$$ \bigtriangledown_{s_i} \tau^e_{s_i} = s_i \cdot \bigtriangledown \tau^{e} \cdot s_i. $$

(7)

The gradient of the shear stresses are evaluated at each tetrahedral element as:

$$ \bigtriangledown \tau^{e}_{i,j} = \sum_{A=1}^{3} N_{A,j}\tau^{i}_{A}. $$

(8)

The corresponding MWSSG is represented as $\overline{\tau^{G,e}}.$ We then derive the nodal values as

$$ \int\limits_{\Upomega} \overline{\tau^{G}(x)} w(x) d\Upomega = \int\limits_{\Upomega} \overline{\tau^{G,e}} w(x) d\Upomega. $$

(9)

Discretizing this equation and choosing the test functions as shape functions, we have

$$ {\mathbf A}_{e=1}^{nel} \left( N_{A}, N_{B} \right) \overline{\tau^{G}_B} = {\mathbf A}_{e=1}^{nel} (N_{A},1) \overline{\tau^{G,e}}, $$

(10)

where A denotes the assembly matrix. Integrating this equation, we observe that:

$$ {\mathbf A}_{e=1}^{nel} \frac{ | J_e |}{24} \left( \begin{array}{lll} 2 &1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{array} \right) \left( \begin{array}{l} \overline{\tau^{G}_{1}} \\ \overline{\tau^{G}_{2}} \\ \overline{\tau^{G}_{3}} \end{array} \right) = {\mathbf A}_{e=1}^{nel} \frac{ |J_e |}{6} \left( \begin{array}{l} 1 \\ 1 \\ 1 \end{array} \right) \overline{\tau^{G,e}}. $$

(11)

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.