Personalized coronary and myocardial blood flow models incorporating CT perfusion imaging and synthetic vascular trees

i. Patient population

We performed simulations on a subset of 4 patients (figure 1) from a cohort who were scheduled to undergo coronary artery bypass graft (CABG) surgery, and also underwent coronary CT angiography (CCTA) and dynamic CT myocardial perfusion (MPI_CT) imaging at Stanford University School of Medicine, Stanford CA, USA as part of an ongoing study (NCT03894423). This study was approved by the Institutional Review Board at Stanford University School of Medicine. Written informed consent was received prior to patient participation. The exclusion criteria for the cardiac CT exam were: (1) Age<40 years; (2) left ventricular ejection fraction <30%; (3) repeat CABG; (4) contra-indications to iodine contrast medium, including eGFR<45 ml/min, known contrast allergy; (5) contra-indications to vasodilators, including bronchial asthma, advanced atrioventricular block, sinus node disease, clinically significant carotid artery narrowing, severe aortic stenosis or left ventricle outflow tract narrowing, blood pressure <90mmHg, use of dipyridamole or aminophylline; (6) conditions technically challenging the examination, including severe valvular disease, atrial fibrillation and body weight >100kg; (7) other conditions of that would render the examination unsafe, including pregnancy, unstable ischemia, severe arrhythmia, heart failure. From the patient cohort that underwent a cardiac CT exam, the inclusion criteria for this computational modeling study were the availability of CCTA and MPI_CT images of sufficient quality to re-construct the 3D coronary anatomy and myocardial blood flow (MBF), as well as clinical measurements of heart rate, blood pressure, cardiac output and ejection fraction measured at the time of the CCTA.

• Download figure
• Open in new tab

Figure 1:

Coronary and LV anatomical models for the four patients in this study. The colors on the coronary arteries qualitatively represent simulated average pressure over a cardiac cycle and the colors on the LV represent MBF from MPI_CT. Figures (a), (b), (c) and (d) show patients 1, 2, 3 and 4, respectively.

ii. Clinical imaging protocols and processing

MPI_CT followed by CCTA were performed on a third-generation dual-source CT scanner (SOMATOM Force, Siemens Healthineers). Hyperemia was induced by slow bolus injection of 0.4 mg regadenoson. A 40-ml bolus of iopamidol (Isovue, 370 mg/ml, Bracco, Italy) was injected at 5.5 ml/s followed by 40 ml saline. The dynamic perfusion protocol requires that the table position is repeatedly alternated (shuttle mode) to completely cover the myocardium and acquire up to 15 data samples during a 30s breath hold. Acquisition parameters were 80-kV tube voltage, automated exposure control (300 mAs/rotation at 80 kV as reference), 96×0.6-mm detector collimation, 105-mm z-axis coverage, 250-ms gantry rotation time, and 66-ms temporal resolution. Reconstructed images with a 3.0-mm slice thickness, and 2.0-mm reconstruction interval were processed on a dedicated workstation (Syngo.CT Myocardial Perfusion, Siemens Healthineers). Serial samples were aligned using a motion correction algorithm. The arterial input function was sampled in the descending aorta. Time-attenuation curves were created for each voxel within the segmented left ventricle (LV) myocardium. Using a dedicated parametric deconvolution algorithm based on a 2-compartment model of intra- and extravascular space, a fit curve was created. MBF was calculated as the ratio between the maximum slope of the fit curve and the peak arterial input function, and the results were exported as three-dimensional MBF maps (figure 2a). Comprehensive descriptions of the myocardial perfusion protocol and the calculation of MBF can be found in [35].

• Download figure
• Open in new tab

Figure 2:

Overview of the computational modeling workflow. (a) Raw volume mesh of the LV, with color contours showing MBF from MPI_CT. (b) Segmented image-based coronary tree and the LV models after co-registration and processing. (c) LV perfusion volumes corresponding to each coronary artery. (d) The flow simulation setup, consisting of 3D models of the coronary arteries and LV, synthetic vascular trees perfusing the LV, and closed-loop boundary conditions. RA, RV, LA, LV and C denote the right atrium, right ventricle, left atrium, left ventricle, and coronary outlet LPN models respectively. (e) Simulated aortic pressure, aortic flow and left coronary flow over a cardiac cycle.

A standard CCTA was performed 10 minutes after MPI_CT by prospective ECG-triggered axial scan mode. The tube voltage ranged between 100 and 120 kV. Most patients received 50mg aminophylline after the perfusion scan, as well as intravenous metoprolol if the heart rate exceeded 70 bpm. All patients received sublingual nitroglycerin just before the contrast-enhanced CT exam.

Since the CCTA and MPI_CT were performed at mid-diastole and systole respectively, they were co-registered using the affine registration procedure implemented in 3D Slicer (www.slicer.org). The LV myocardial volume was then segmented from the co-registered MPI_CT image using thresholding in ParaView (www.paraview.org). It was smoothed and meshed with tetrahedral elements using Meshmixer (www.meshmixer.com) and TetGen [70], respectively. Finally, the MBF values at each voxel of the MPI_CT scan were interpolated on to the nodes of the LV mesh. Note that all the software used for the processing of these images are freely available, and all except Meshmixer are open-source. Figure 2(b) shows the 3D volumes of the co-registered LV and coronary vasculature for patient 1.

Following the construction of co-registered 3D models of the LV myocardium and coronary arteries, the LV volume was divided into non-overlapping sub-volumes corresponding to the perfusion territory of each coronary artery (figure 2c). This was done using Voronoi tessellation, where each node in the tetrahedral mesh for the LV volume was assigned to its closest coronary artery outlet. Therefore, each coronary artery outlet was associated with a distinct region of the LV myocardium. Subsequently, the flow-rate associated with each coronary outlet was estimated by integrating the MBF obtained from MPI_CT within each corresponding LV perfusion territory. We note that the LV perfusion territories corresponding to each coronary artery can also be computed based on the distance of each point in the LV to the closest coronary artery (instead of coronary artery outlet) [62, 86] as well as weighted measures of distance [13, 49]. This approach was chosen for simplicity, however, the framework presented in this work is equally applicable to any method for computing perfusion regions in the LV.

iii. Image-based modeling and coronary flow simulations

We used the open-source SimVascular software [77] for segmentation of coronary arteries from CT angiography images, construction of 3D anatomical models, and coronary flow simulations. Tetrahedral finite element meshes were created using the open-source TetGen package [70], which is included with SimVascular. Three dimensional flow simulations were performed using the svSolver flow solver within SimVascular, which uses a stabilized finite element method with linear tetrahedral elements for spatial discretization and generalized-α time-stepping with second-order time accuracy [23, 83].

The governing equations for coronary blood flow were the three-dimesional incompressible Navier-Stokes equations, where ⃗u and p are the blood flow velocity and pressure, respectively. Blood was assumed to be a Newtonian fluid with viscosity µ = 0.04 dynes/cm² and density ρ = 1.06 g/cm³. The Newtonian fluid assumption is valid because non-Newtonian effects begin to be apparent only for blood vessels with diameters below 300 µm [65], which was approximately the limit of the smallest vessels modeled in this study.

Coronary artery walls were treated as deformable, and fluid-structure interaction arising from the coupling between blood flow and the artery walls was handled using the coupled momentum method [16]. Material properties for the artery walls were selected based on literature data following previous work [20, 39]. The elastic modulus for the aorta and coronary arteries were set to 0.25 MPa and 1.15 MPa, respectively [19, 8, 59]. Wall thicknesses were based on published radius-thickness ratios for the aorta [8] and morphometric data for the coronary arteries [53].

iv. Synthetic vascular trees

We augmented our image-based coronary artery anatomical models with synthetic vascular trees to model coronary vasculature with diameters below the imaging resolution of the CCTA. Constrained Constructive Optimization (CCO) [63] was used to generate synthetic vascular trees that are constrained by vascular scaling laws and hemodynamic factors while minimizing vascular volume. CCO has recently been adapted to generate multiple non-intersecting trees within non-convex ventricle-like tissue volumes [24], although this process is computationally expensive. In this work we used an improved open-source implementation of Sexton et al. [68] to more efficiently generate synthetic vascular trees at the outlet of each coronary artery that perfused the LV. While the CCO implementation used here allows the generation of several competing vascular trees within the same tissue volume, we chose to generate one independent vascular tree within the LV volume associated with each coronary artery. In our experience, this approach resulted in the LV being more uniformly vascularized and with significantly lower computational cost by avoiding the generation of unphysiological competing tree morphologies from approximately 15 coronary artery outlets.

For each case simulated in this work, we generated synthetic vasculature with a total of 500 outlet branches perfusing the LV (compared to approx. 16 branches perfusing the LV in the purely image-based anatomical models). The total number of synthetic outlet branches was divided amongst the vascular trees appended to each coronary artery outlet based on the size of their corresponding LV perfusion volumes. The smallest vessels produced were approximately 300 µm in diameter. Flow in the synthetic vasculature was modeled using a lumped-parameter/0D representation of the vessels, which included resistance, capacitance and inductance to model viscus losses, wall dilation, and fluid inertia respectively. The lumped-parameter resistance, capacitance and inductance of each vessel were given by linearization and Poiseuille flow simplifications of the Navier Stokes equations [50]. The lumped-parameter flow in these synthetic vascular trees was simulated using the open-source svZeroDPlus solver, which is part of the SimVascular software suite, and was coupled to the hemodynamics in the 3D image-based anatomical model using a modular implicit coupling method [15].

v. Coronary flow boundary conditions

We used a closed-loop lumped parameter network (LPN) model to enforce boundary conditions for the coronary flow simulations at the aortic inlet and the aortic and coronary outlets of the 3D model [61, 74, 39] (figure 2d). The LPN included specific features to model distal systemic resistances using Windkessel models [81] and the effect of the four heart chambers. It also captured the distal resistance and intra-myocardial pressure experienced by coronary arteries, which produces out-of-phase coronary flow with respect to the cardiac cycle [28]. Figure 2(e) shows representative pressure and flow waveforms for one case. As seen in figure 2(d), the LPN boundary conditions at the coronary artery outlets consisted of three resistance elements – R_a, R_µ and R_v, to model the resistance of small arteries, microvasculature, and the venous system respectively. These resistance and capacitance elements distal to each coronary artery in the model, as well as other parameters of the closed-loop LPN, such as the heart chamber elastances and capacitances, and systemic resistances, were estimated to match clinically measured data (see section vii). As with the lumped parameter representation of the synthetic vascular trees discussed in section iv, the closed-loop LPN boundary conditions were simulated using svZeroDPlus and coupled with the 3D flow solver using a modular implicit coupling method [15].

We also simulated coronary flow at hyperemia to recapitulate common clinical measurements that assess coronary hemodynamics under stress (including FFR and MPI_CT). We simulated hyperemia by scaling the resistance boundary conditions distal to coronary arteries by a factor of 0.24 compared to baseline resting values [84], as commonly done in FFR_CT [72].

vi. Myocardial blood flow simulations

Myocardial perfusion in the arterioles and capillaries within the LV myocardium was modeled as flow through a porous material due to the prohibitive computational cost of resolving discrete microvascular blood vessels [5, 9, 22, 13, 49, 29]. Flow through porous materials is described by Darcy’s Law, and in this work we used a single-compartment Darcy model governed by the equations, where w∼ and p are the Darcy flow velocity and pressure, respectively. The permeability of the LV myocardial tissue is given by K. We assumed constant and isotropic permeability for simplicity, as in prior Darcy models for myocardial perfusion [5, 13, 14, 49], i.e. K(⃗x) = kI, where I is the identity tensor and k is the scalar, spatially-constant permeability. We fixed k = 2 10⁻⁵ cm² Pa⁻¹ s⁻¹, which is informed by previous work [5, 49]. Myocardial perfusion was driven by the coronary flow in the upstream vasculature via pressure source terms, which are denoted by p_src in equation 2. These pressure sources, corresponding to the outlet of each coronary artery that is perfusing the LV, were specified as the microvascular pressure in the coronary outlet boundary condition model (see section v and P_µ in figure 2d). Therefore, each perfusion territory associated with a coronary artery outlet, discussed in section ii, was assigned a corresponding pressure source. The pressure of the venous return system is given by p_snk in equation 2, and represents a sink for myocardial perfusion. We assumed p_snk = 0 mmHg in this work [49]. Lastly, the terms β_src and β_snk represent the conductance of flow entering and leaving the tissue, respectively, and were parameterized as follows,

In equation 3, Q_LV and V_LV are the total flow perfusing the LV and the volume of the LV respectively. The spatially-averaged pressure source (p_src) is denoted as p^-_src and p_cap is the target capillary pressure assumed to be p_cap = 15 mmHg [5, 49]. Note that for a given tissue volume V_i V_LV, such as one vessel-specific perfusion territory, the MBF within that volume is given by,

We simulated perfusion in patient-specific LV volumes segmented from MPI_CT, as described in section ii. Unstructured tetrahedral meshes were generated using the open-source TetGen package [70], and pressure sources (p_src) were imposed from flow simulations in the coronary arteries for each patient-specific case. A finite element solver for equation 2 employing linear tetrahedral elements was implemented in the open-source svFSI software, which is part of the SimVascular software suite. We used zero-flux boundary conditions on the surface of the LV mesh.

vii. Parameter estimation for personalized boundary conditions

We employed a multi-stage surrogate-based optimization procedure to estimate the parameters of the closed-loop LPN boundary condition model (discussed in section v) to recapitulate measured clinical data for each patient. This was performed for purely image-based models as well as those augmented with synthetic vascular trees. The clinical targets we used were the vessel-specific flows from MPI_CT; systolic and diastolic aortic pressure from blood pressure cuff measurements; left ventricular ejection fraction, stroke volume, cardiac output and systolic/diastolic ventricular volumes from echocardiography; and physiological targets derived from literature including pulmonary pressure and coronary flow waveforms [31, 48, 74]. We note here that the vessel-specific flow and MBF targets from MPI_CT were calculated relative to the total LV flow because MPI_CT is known to underestimate the absolute flow [46, 86]. As a basis for comparison, we also tuned each model to match all the above clinical targets except vessel-specific flows. In these latter models, the distribution of flow amongst the coronary arteries was prescribed based on their lumen diameters using Murray’s law [43, 88], as is more conventionally done.

The surrogate model we used for the optimization was a lumped parameter representation of each patient-specific coronary artery anatomical model. This was generated using an automated method [50] that extracts the lumen centerlines of all vessels in the anatomy, computes vessel radii by traversing along the centerlines, and distills this information into a network of resistance-capacitance-inductance lumped parameter models for each vessel in the anatomy. As in the lumped-parameter representation of the synthetic vascular trees discussed in section iv, the resistance, capacitance and inductance of each vessel were given by linearization and Poiseuille flow simplifications of the Navier Stokes equations [50]. We henceforth refer to this as a 0D surrogate model.

vii.1 Optimizing for cardiac function

In the first stage we optimized the parameters of the closed-loop system to match all the above clinical targets except the vessel-specific flows. This consisted of 36 parameters governing the hemodynamics and function of the four heart chambers, the pulmonary circulation, the intramyocardial pressure experienced by coronary arteries, and the total distal resistance and capacitance at all the aortic and coronary outlets. For details about the parameters optimized in this stage, see reference [39]. This optimization was performed using the Nelder-Mead method, a deterministic gradient-free optimization technique [45]. The convergence tolerance for the optimization was 0.1, maximum number of iterations was 200, the parameter update step size was 0.1 of the range for each parameter, and the number of restarts was 50.

The initial guess for the total vascular resistance for each patient was determined from the clinically measured average aortic pressure divided by cardiac output. The total coronary flow was given by the morphometric scaling with the LV volume, Q_LV = 3.41 V^0.75 mL/min for branches perfusing the LV and 0.2 Q_LV for other branches [6, 10, 49, 86]. For the models tuned to match the flow distribution measured by MPI_CT, the total coronary resistance was distributed amongst the branches of the coronary tree based on the ratio of total MBF to vessel-specific MBF for each vessel. On the other hand, for the models tuned to match flow distributions resulting from Murray’s law, the total coronary resistance was distributed amongst all coronary branches based on their diameters scaled to an exponent of 2.6 [88]. The total capacitance was set as 0.28 10⁻³ cm⁵/dynes for aortic branches and 10⁻⁶cm⁵/dynes for coronary branches based on previous work [61]. The total capacitance was distributed amongst branches to be proportional to their outlet areas [6]. When tuning all the models, the permissible range for each parameter was assigned based on physiological values and the initial guess was based on the maximum posterior of the distribution of parameters computed using Markov Chain Monte Carlo techniques in our previous work using similar closed-loop LPN models for coronary hemodynamics [74].

vii.2 Optimizing for blood flow distribution: Without synthetic vasculature

The second stage of optimization focused on vessel-specific outlet boundary conditions. In this stage, we kept the 36 parameters that were optimized in the first stage fixed, and tuned the outlet resistances (represented by R_a, R_µ and R_v in figure 2(d)) of all the branches in the coronary tree to match clinically measured flow distributions from MPI_CT.

For the models that were not augmented with synthetic vascular trees, this consisted of two parameters for each of the N_outlets^CT coronary outlets perfusing the LV – the total outlet resistance for each branch (Rⁱ_a + Rⁱ_µ + Rⁱ_v, where i = 1,…, N^CT_outlets is the index of the outlet), and the ratio of proximal vessel to microvascular resistance at each outlet (Rⁱ_a/Rⁱ_µ). The former determines the flow distribution in the coronary tree while the latter dictates the distribution of microvascular pressure that drives myocardial perfusion. We first used Nelder-Mead optimization to tune the relative distribution of total resistance amongst all the coronary artery outlets so that the resulting distribution of flow in the coronary tree matched that measured from MPI_CT. This was done using a N^CT_outlets-dimensional optimization procedure that estimated a scaling factor constrained to the range [0.5, 2.0] for the resistance at each outlet, such that the total outlet resistance of the model was kept equal to the value obtained from the first stage of optimization. The convergence tolerance for the error was set to 0.1, the maximum number of iterations was 200, the parameter update step size was 0.1 of the range, and the number of restarts was 50.

We then fixed the total resistance of the complete tree as well as the total resistance at each outlet, and estimated the ratio of proximal to microvascular resistance that matched the required distribution of microvascular pressure. Note that the required microvascular pressure source (p_src in equation 2) in each vessel-specific perfusion territory can be computed from equation 4. This is given by pⁱ_src = MBF_Vi/(β_srcV_i) + p_cap, where a guess for all the pressure sources is required to compute p^-_srcin β_src. In addition, the simulated microvascular pressure at each outlet i in the flow simulation is given by Pⁱ_μ = Pⁱ − Qⁱ Rⁱ_a, where Pⁱ and Qⁱ are the pressure and flow at outlet i (see figure 2d). We therefore ran a 3D simulation with deformable walls to estimate the outlet pressure Pⁱ and microvascular pressure Pⁱ at each coronary outlet. For this initial simulation, we assumed the ratio of proximal to microvascular resistance was Rⁱ_a/Rⁱ_μ = 0.38 at each outlet, following previous work [61, 39]. This computation of the pressure loss along each branch, combined with the required pⁱ for each branch, informed a direct update of Rⁱ_a/Rⁱ_μ at each outlet to obtain the required microvascular pressure. This updated ratio of proximal to microvascular resistance, along with the optimized distribution of outlet resistances, was then used to run a 3D simulation with deformable walls at hyperemia. The microvascular pressure obtained at each outlet from this simulation, along with the required microvascular pressure source for each outlet from equation 4, was then used to inform a second update of the ratio between the proximal and microvascular resistance at each coronary outlet.

vii.3 Optimizing for blood flow distribution: With synthetic vasculature

For the models augmented with synthetic vascular trees, a different approach was employed for tuning the boundary conditions at the outlets, since these models consisted of N^syn_outlets = 500 coronary outlets (compared to 16 outlets without synthetic vasculature), rendering the approach above computationally intractable. Moreover, since the synthetic trees generated by CCO are not patient-specific, we noticed that the vascular resistance along some paths from the aorta to outlet vessels led to unphysiological outlet boundary conditions (such as negative resistance) in order to match clinically measured flow and pressure distributions. Therefore, we developed a novel iterative approach to personalizing synthetic vasculature that estimates not only the parameters at the coronary outlets, but also dilates or constricts the synthetic vessels to match patient-specific flow and pressure targets.

Each iteration of the procedure consisted of: (1) sub-iterations to match the required flow distributions in the 0D surrogate to within 5% error, which involved updating the distal resistance boundary condition at each outlet vessel and dilating the synthetic vascular trees if unphysiological negative resistances were required to match the flow targets; (2) a simulation of the updated surrogate model and comparison of the outlet pressure distribution with that required to match clinically imaged MBF; (3) the dilation or constriction of all the synthetic vascular trees to obtain the required pressure distribution if the average pressure error was greater than 5%. The above steps were iterated upon until both the flow and pressure errors (based on the 0D surrogate model) were below 5%. Each of these steps is described in detail below.

With the parameters derived from the first stage of optimization, we first ran a 0D simulation at hyperemia to compute the pressure at all outlets. Using the aortic pressure and outlet pressures, we estimated the linear vascular resistance along each path from the aortic inlet to a coronary outlet, given by Rⁱ_path = (P_aorta − Pⁱ)/Qⁱ, where Pⁱ and Qⁱ are the pressure and flow at an outlet synthetic vessel indexed by i. We also computed the resistance of the paths from each outlet of the image-based epicardial vessels to all outlets of the synthetic tree that was appended to that epicardial vessel. This is given by Rⁱ_path,syn = (Pⁱ_epi − Pⁱ)/Qⁱ where P^iepi is the pressure at the outlet of the image-based epicardial vessel which is upstream of the synthetic vessel outlet i.

Given the resistance Rⁱ_path along each aorta-to-outlet path and the total resistance of the model estimated from the first stage of optimization (R_tot), we calculated the required resistance distal to each outlet as Rⁱ_out,reqd = R_totQ_LV/Qⁱ_reqd − Rⁱ_path, where Qⁱ_reqd is the target flow through outlet i, according to MPI_CT. This informed an updated resistance distal to each outlet. For some outlets in the model, this estimated outlet resistance was negative. This was due to excessive pressure drop along the path from the aorta to outlet i, i.e. Rⁱ_path was too large. Since each Rⁱ_path consisted of the sum of resistance from the upstream image-based epicardial vessels and downstream synthetic vessels, we addressed this by dilating the synthetic vascular trees. This achieved the required resistance along each path while maintaining the anatomy of the upstream branches that were based on patient-specific imaging. For each Rⁱ_out,reqd that was negative, we estimated the required resistance of the path from the outlet of the corresponding image-based vessel to the synthetic outlet i, which we denote as R^{∼ i}_path,syn, such that Rⁱ_out,reqd > 0. Note that R^{∼ i}, since we are aiming to decrease Rⁱ_path. Each of these synthetic vessel outlets was then assigned a radius scaling factor to achieve this reduced resistance, given by α = ⁱ_path,syn/(Rⁱ_path,syn)^1/4 > 1. The synthetic trees were then dilated using the method described below. The above procedure was iteratively performed until physiological outlet resistances were attained for each synthetic vessel outlet and the average error in the outlet flow fractions, ɛ_flow = ∑_i(|Qⁱ_reqd − Qⁱ_reqd)N^syn_outlets, was less than 5%.

We scaled the radii of synthetic vessels after each outlet vessel was assigned a radius scaling factor, αⁱ, by traversing up each synthetic tree from the outlet vessels to the root of the tree. As we traversed up the tree, each parent vessel was assigned a scaling factor that was the maximum of the scaling factors of its two daughter vessels. Moreover, the maximum scaling factor for each vessel was constrained by the radius of its parent vessel. This ensured the preservation of the morphometric structure of the tree. In this way, each vessel in every vascular tree was assigned a scaling factor αⁱby propagating the outlet vessel scaling factors up the tree. We then updated the resistance, inductance and capacitance values of all the synthetic vessels based on the new radius.

Once the outlet boundary conditions were updated through the above procedure, we tuned the vascular resistance of the synthetic vasculature to obtain the required microvascular pressure distribution at the coronary outlets based on the MBF distribution from MPI_CT. As discussed in section vii.2, the target microvascular pressure source at each coronary outlet (p_src in equation 2) can be estimated from equation 4. This is given by = MBF_Vi /(β_srcV_i) + p_cap. Therefore, the target pressure at every synthetic vessel outlet with index i can be calculated as Pⁱ_reqd = pⁱ_src + Qⁱ_reqdRⁱ_a. Furthermore, the required vascular resistance along each path from the aortic inlet to each synthetic vessel outlet can be calculated from the required outlet pressure and required outlet. Therefore, we ran a second 0D simulation at hyperemia to estimate the pressure at all the synthetic vessel outlets. We computed the average relative pressure error over all the outlets as ∈_pressure = Σ_i(|Pⁱ − Pⁱ_reqd~/Pⁱ_reqd)N^syn_outlets. If this error was greater than 5%, we proceeded as follows to update the resistance of synthetic vascular trees. The 0D simulation again allowed us to calculate the path-specific vascular resistance from the aortic inlet to each synthetic tree outlet (Rⁱ_path). Similarly, we calculated the resistance from each image-based epicardial coronary outlet to the synthetic vessel outlets that were appended to that image-based vessel. Given the required resistance along each path computed above (Rⁱ_path,reqd, and Rⁱ_path,syn) the separate contributions to this resistance from the image-based vessels and synthetic vessels (Rⁱ_path), we computed the required vascular resistance of the synthetic vessels along the path from the outlet of the corresponding upstream image-based vessel to the synthetic outlet i. This resistance, denoted by R^{∼ i}_path,syn, then informs a radius scaling factor) α_i = (Rⁱ_path,syn/Rⁱ_path,syn)^1/4 for the each synthetic vessel outlet i. We constrict/dilate the synthetic vascular trees using the procedure described above, and then return to the start of the iteration to estimate new outlet resistances and the flow distribution resulting from this updated anatomy.

i. Impact of MPI_CT-informed coronary flow distributions

A major take-away from this study is that MPI_CT-based models resulted in significantly different flow distributions than Murray’s law-based models. The largest differences in coronary artery flow between occurred distal to lesions. Moreover, our findings suggest that there is no consistent relationship across patients between the diameter of a coronary artery and its corresponding flow as measured by MPI_CT, albeit in a small cohort. This implies that the errors in coronary flow distribution arising from Murray’s law based boundary conditions are not predictable a priori. However, this should be confirmed with larger patient cohorts and may differ in healthy patients, according to the work of van der Giessen et al. [78]. Crucially, the differences in flow within epicardial vessels between the models informed by MPI_CT and Murray’s Law affected the computation of FFR, particularly in vessels with borderline lesions of 50% stenosis. Risk stratification based on hemodynamics is particularly valuable in these lesions [73, 44], hence, the implications of these findings for the clinical use of computational models warrants more attention.

The difference in flow distributions resulting from Murray’s law and MPI_CT likely stems from inter-patient variability arising from coronary regulation, metabolic demand, etc., which are especially important in diseased coronary arteries. Our findings agree with previous studies that have shown significant hemodynamic differences, especially in wall shear stress, between models that used boundary conditions based on Murray’s law versus in vivo flow measurements [78, 62, 32]. We note, however, that these studies only included models of individual coronary arteries or single bifurcations, rather than the entire coronary artery tree simulated here.

Another reason for the inaccurate flow distributions resulting from Murray’s law-based boundary conditions is its strong dependence on the outlet diameter of arteries. In the context of image-based computational models, these diameters are derived from segmentations that are error-prone. In contrast, one can view the use of MPI_CT-based boundary conditions as a lumped-parameter correction for inconsistencies that arise in the reconstruction of vessel lengths and lumen diameters in image-based models. For example, the resolution of clinical imaging dictates the length of each vessel visible in the image. However, this is not reflective of the actual length of the vessel, and its resultant hydrodynamic resistance or pressure drop. This therefore can be corrected for by optimizing lumped-parameter boundary condition resistances to compensate for inaccuracies in image-based segmentation. Indeed, recent work by Xue et al. [86], who demonstrated a framework to compute FFR from computational models informed by MPI_CT, showed the robustness of such computational models to inaccuracies in segmenting coronary epicardial vessels.

In the context of clinical translation, a particular advantage of MPI_CT-informed models is the ability to provide high-resolution, full coverage quantification of MBF. The ability to resolve vessel-specific MBF distributions (in contrast to lower-resolution perfusion imaging methods such as PET/SPECT) is especially valuable in the context of the framework demonstrated here. Moreover, there are clear advantages to using a single imaging modality to produce both CCTA and MPI_CT so that the entire framework can be integrated into a single non-invasive imaging protocol. The prognostic value of MPI_CT for detection of functionally significant CAD was previously demonstrated in several studies [1, 18, 33, 47], and its incorporation into computational models has substantial potential for added clinical value.

ii. Synthetic vascular trees and multi-scale modeling

Another contribution of this work is the incorporation of the automated model personalization with multi-scale models of coronary heomdyamics. Our framework included flow in vessels smaller than the CCTA imaging resolution, microvascular perfusion in the myocardium, and LPN models of cardiac function and systemic circulation.

Our use of synthetic vascular trees appended to image-based coronary vessels is motivated by studies that have highlighted the effect of inaccurate coronary artery branching on modeled hemodynamics [82, 79]. Synthetic vascular trees provide a viable method to address this challenge even when using imaging of insufficient resolution, especially when combined with published morphometric data [26] to create realistic vascular anatomies. In addition, we demonstrated that models of myocardial perfusion that include more uniform and realistic vasculature can more accurately recapitulate clinically imaged MBF distributions. This also translates to the assimilation of data from MPI_CT data into computational models, i.e. higher resolution coverage of the perfusion volume by synthetic vasculature allows higher-fidelity assimilation of MPI_CT data into computational models.

Furthermore, our incorporation of closed-loop LPN models representing cardiac function allows us to further personalize models based on clinical measurements, rather than using literature-based boundary conditions. Such LPN models have been previously used for model-based treatment planning and to model patient-specific changes in cardiac function [64]. Finally, the coupling of the upstream hemodynamics in the coronary artery tree with the microvascular blood flow in the myocardium could enable models of coronary revascularizaton on myocardial tissue. This is infeasible in prior models of coronary hemodynamics that focus exclusively on epicardial hemodynamics [72]. In addition, this multi-scale framework opens up the possibility of studying microvascular disease using computational models.

iii. Limitations and future work

We recognize several limitations in this study. The framework presented was deterministic, utilizing single measurements of the clinical targets for each patient to optimize the computational models. However, each clinical measurement is associated with various sources of uncertainty from the measurement protocols, patient physiology, as well as post-processing. Future work should account for uncertainties in the estimation of patient-specific parameters as well as in the model predictions; this can be performed following recent work in uncertainty quantification by our group and others [66, 67, 17]. Similarly, the sensitivity of the model to different parameters also needs to be assessed. Additionally, the dataset used here did not include clinically measured intra-coronary flow to compare with the estimated flows from MPI_CT. While the goal is to build a non-invasive framework based entirely on MPI_CT, comparisons with intra-coronary measurements would be valuable to increase confidence.

Another area for improvement is the complexity of the modeling utilized. While we used 3D simulations of coronary hemodynamics to compute FFR, the parameter estimation used 0D surrogate models that deviated from the true hemodynamics. In particular, 0D models do not accurately capture the flow within stenoses and at junctions between arteries. This contributed to larger errors in the flow distributions predicted from 3D models than the thresholds specified in the surrogate-based parameter estimation. Improved surrogate models will also be valuable for applying this method within time-frames that are feasible for clinical use. Another simplification was the single-compartment Darcy model we used to simulate myocardial perfusion. Previous work [22] has shown heterogeneity and anisotropy in the myocardial microvasculature, which are ignored by the assumption of constant, isotropic permeability. Moreover, the permeability was not tuned to patient-specific measurements.

This framework would also benefit from a more physiological partitioning of the myocardium into regions perfused by each coronary artery. Although there is evidence that each part of the coronary artery tree perfuses a distinct volume of the myocardium [71, 9], there is no consensus method to perform this partitioning in previous studies [62, 13, 49, 86]. However, these perfusion regions affect the simulated flow through each coronary artery. They also influence the structure of the synthetic vascular trees that extend into the myocardium from epicardial coronary arteries. While we generated vascular trees within pre-defined LV perfusion regions, future work will include methods to grow larger synthetic vascular trees that compete to perfuse physiological regions of the myocardium.

As mentioned earlier, the main goal of this work was to develop models for personalized CAD treatment planning. To that end, future work will explore applying these tools to predict the hemodynamics of patients after CABG surgery. This will involve studies in larger patient cohorts and methods to account for post-surgical vascular and cardiac adaptation.