Myocardial Perfusion Simulation for Coronary Artery Disease: A Coupled Patient-Specific Multiscale Model

Patient-specific models of blood flow are being used clinically to diagnose and plan treatment for coronary artery disease. A remaining challenge is bridging scales from flow in arteries to the micro-circulation supplying the myocardium. Previously proposed models are descriptive rather than predictive and have not been applied to human data. The goal here is to develop a multiscale patient-specific model enabling blood flow simulation from large coronary arteries to myocardial tissue. Patient vasculatures are segmented from coronary computed tomography angiography data and extended from the image-based model down to the arteriole level using a space-filling forest of synthetic trees. Blood flow is modeled by coupling a 1D model of the coronary arteries to a single-compartment Darcy myocardium model. Simulated results on five patients with non-obstructive coronary artery disease compare overall well to [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{15}$$\end{document}15O]\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {H}_{{2}}$$\end{document}H2O PET exam data for both resting and hyperemic conditions. Results on a patient with severe obstructive disease link coronary artery narrowing with impaired myocardial blood flow, demonstrating the model’s ability to predict myocardial regions with perfusion deficit. This is the first report of a computational model for simulating blood flow from the epicardial coronary arteries to the left ventricle myocardium applied to and validated on human data. Electronic supplementary material The online version of this article (10.1007/s10439-020-02681-z) contains supplementary material, which is available to authorized users.

rotation time of 270 ms. To visualize the coronary artery lumen a bolus of 100 mL iobitidol (Xenetix 350) was injected intravenously 5.7 mL s −1 , with immediately after a 50 mL saline chaser. The scan was triggered with an automatic bolus tracking technique. The region of interest was placed in the descending thoracic aorta with a threshold of 150 HU. Metoprolol 50 to 150 mg was administered orally if patients had a prescan HR ≥ 65 beats per minute (bpm) one hour before the start of the CT protocol. If necessary, 5 to 25 mg metoprolol was given intravenously during the scan to achieve a heart rate < 65 bpm. All patients received 800 µg of sublingual nitroglycerine immediately before cCTA.

[ 15 O]H 2 O Positron Emission Tomography (PET)
Patients had to refrain from taking products containing caffeine or xanthine 24 hours before imaging. Patients fasted for at least 4 hours before the scan protocol. All patients were  Using the initialized flows in eq. 3 (main paper), we calculate ideal baseline resistance at each terminal segment: For each terminal segment, an expected minimum resistance is estimated as follows: Here, we assume a uniform factor but this can vary for each patient. Anatomically, the factor represents a maximal radius dilation capacity of 40% based on uniform dilation of a tree to achieve a 4-fold reduction in resistance 5 . Note that minimum resistances are calculated once and fixed for the whole parameterization process.
Using the terminal segment flows as boundary conditions for the first iteration, flow and pressure are computed for the entire tree by solving the 1D equations. The computed terminal segment pressure p T,i and the assigned terminal segment flow q T,i rest determine the simulated terminal segment resistance R T,i sim : If the simulated terminal segment resistance is lower than the minimum resistance value R T,i min , we update both the geometry and the flow. We dilate the radius of the terminal segment and its upstream synthetic segments by multiplying them with a uniform factor where R T,i p = If the simulated terminal segment resistance is greater or equal to R T,i min , we update the geometry by multiplying by a factor of R T,i base R T,i sim 0.25 but maintain the same flow q T,i rest . Note that each terminal segment can have a different radius dilation factor as the resistance of the terminal segment will depend on the pressure loss accumulated along the path. For the upstream segments branching to child segments with different dilation factors, bigger dilation factors are chosen to dilate them. Several iterations are necessary to obtain convergence as both the geometry and boundary conditions are updated in each iteration.
The convergence criteria used: with n the iteration counter. Note that for synthetic trees arising from the segmented vessel outlets, the dilation propagation stops before reaching the root segment since this was directly observed in the CT data. For synthetic trees arising along a segmented vessel, the dilation propagation is applied up to and including the root segment. This ensures that only arteriole and small artery vessels are dilated since they contribute most to coronary resistance modulation and are dilated most under hyperemic conditions. Further, the CT data is obtained with the use of sublingual nitrates so the segmented tree is assumed to be fully dilated at rest.

Hyperemic conditions
The minimum resistance at each terminal segment is computed using the ideal terminal segment flow: Since the whole system is already dilated to the maximum radius dilation capacity, only terminal flows can be changed. The system is first solved using q T,i stress as a boundary condition at each terminal segment outlet. With the resulting simulated pressures p T,i stress and set flows, we calculate the simulated resistance: In each iteration, the terminal segment flow is incremented by δq by solving the following equation: Consequently when the simulated resistance R T,i sim,stress is less than the minimum resistance R T,i min,stress , the flow is decreased. Otherwise flow is increased. Simulation is repeated with the updated terminal flow values. Upon iteration, the computed terminal resistance values converge to the minimum resistance values. The convergence criteria is: 3 Coupling method In the coupling iterations, the resulting pressure of the coronary model at each terminal segment end p T,i k is used as input to the porous model, defining the source pressure in each tessellation territory Ω i . The porous model is then solved for flow values for each perfusion territory, that are used as new inputs for the coronary model, q T,i k+1 . For iterations with k > 0, the coronary model solves for pressure with given flows while keeping the same geometry and boundary conditions. The coupling loop is illustrated in Fig. 1d (main paper). The coupling convergence is established considering the terminal segment flow values between iteration k and iteration k + 1: In practice for a few cases, in particular the patient with significant obstructive disease, a relaxation scheme was necessary.
Note that in principle, each perfusion territory is assigned to one and only one segment.
However in some cases, the scale of the segment outlets is much smaller than the resolution of the myocardial mesh. As a result, during the Voronoi tessellation computation, two segments can be assigned to the same perfusion territory, which can cause problems to the coupling algorithm. In order to avoid this, we locally refine the myocardial mesh and rerun the tessellation. If the issue persists as the segments are too close to each other and an excessively refined mesh resolution is required, we modify the vasculature by trimming the smallest terminal segment and recompute the tessellation.

Sensitivity to the synthetic vascular network
In this section, we assess the impact of the vascular network on the results by generating different synthetic vasculatures for the same patient (Patient 1) as in section 3.1 (main paper).

Varying number of terminal segments
We generated three distinct synthetic vasculatures with 3000, 6000 and 12000 terminal segments, respectively.
Vasculature analysis We evaluate the vascular depth reached in each vasculature using the Strahler order system 4 . Overall, increasing the number of terminal segments leads to a greater percentage of synthetic segments with lower Strahler orders (5-6) and a lower percentage of segments with higher Strahler orders (7-9) (Fig. 4.1a). This implies that vasculatures with higher number of terminals achieved greater depth. In particular, close to 20% of synthetic segments in the 12000-terminals vasculature were extended down to Strahler order 5, while < 5% and < 1% of segments reached that order in the 3000 and 6000-terminal vasculatures, respectively.
Further investigating the diameter distribution of terminal segments reveals that both smaller diameters and greater terminal diameter uniformity is achieved with an increase in number of terminal vessels (Fig. 4.1b). Mean and standard deviation of terminal diameters are 195 ± 48 µm, 155 ± 38 µm and 123 ± 30 µm for the 3000, 6000 and 12000 terminal vascular trees, respectively. Hemodynamic results Blood flow in the three generated vasculatures is summarized in Table 3 (main paper). In terms of blood flow in the coronary arteries, increasing the number of terminals leads to a more homogeneous flow distribution at the outlet level for both coupled and coronary models, as suggested by the lower standard deviation of terminal flows (Table 3, main paper). This behaviour reflects the more homogeneous diameter distribution.
Naturally, the mean terminal segment flow is also lower, as the same amount of total flow is distributed in more segments.
In terms of blood flow in the myocardium, mean MBF remains almost identical regardless of the number of terminals and the model used. While the coronary model MBF values remain extremely heterogeneous, a slight reduction is observed when increasing the terminal segments (see reduced SD, Table 3  To conclude, the 12000-terminals vasculature demonstrates preferable morphological features (greater vascular depth, more uniform outlet diameter size) leading to a minor, albeit noticeable, improvement in the coronary model results. In order to explicitly model as many vessels as possible while maintaining a manageable additional computational cost, we consider 12000 terminal segments as the default choice for this paper.

Varying randomness of the vasculature
Here we generate 5 vasculatures comprising 12000 terminal segments. All synthetic trees start from the same roots on the segmented vessels, but grow differently based on the inherent randomness in the tree generation as described in 3 . Thus these 5 "seeds" lead to different but similar-looking synthetic networks.
MBF results for both models across different seeds are presented in Fig. 4 The associated perfusion maps display some regional differences, but their spatial distribution of flow has similar main characteristics (Fig. 4.3c

Patient selection and characteristics
For this study, a total of 6 patients were selected to have high cCTA image quality and contain a wide distribution of heart dominance (3 right, 2 codominant, 1 left). They are grouped based on non-obstructive CAD (5 patients) or obstructive CAD (1 patient), with obstructive CAD defined as > 30% reduction in diameter in at least one location on the larger segmented vessels.   Vasculatures comprise 12000 terminal segments.