A patient-specific lumped-parameter model of coronary circulation

A new lumped-parameter model for coronary hemodynamics is developed. This model is developed for the whole coronary network based on CT scans of a patient-specific geometry including the right coronary tree, which is absent in many previous mathematical models. The model adopts the structured tree model boundary conditions similar to the work of Olufsen et al., thus avoiding the necessity of invasive perfusion measurements. In addition, we also incorporated the effects of the head loss at the two inlets of the large coronary arteries for the first time. The head loss could explain the phenomenon of a sudden increase of the resistance at the inlet of coronary vessel. The estimated blood pressure and flow rate results from the model agree well with the clinical measurements. The computed impedances also match the experimental perfusion measurement. The effects of coronary arterial stenosis are considered and the fractional flow reserve and relative flow in the coronary vessels for a stenotic vessel computed in this model show good agreement with published experimental data. It is believed that the approach could be readily translated to clinical practice to facilitate real time clinical diagnosis.

wave intensity analysis and a 1D computational model. It was found that the geometric properties of the pulmonary musculature play an important role in the system. Hence, patient-specific geometry should be included in the modelling of the coronary network.
Detailed arterial geometry and the fluid dynamics are best captured using the finite element (FE) approach and computational fluid dynamics (CFD). Liu et al. 19 carried out an 3D FE simulation of the main circle of Willis combined with a lumped parameter model boundary condition, and the result showed good agreement with the Doppler ultrasound measurements. Using a FE approach, Al-Hassan et al. 20 and Kim et al. 21 studied the three-dimensional whole coronary system using more detailed geometries reconstructed from computed tomography (CT). Each coronary outlet boundary condition was coupled to a lumped-parameter coronary bed model and the inlet was assigned to the systemic circulation. However, such a 3D approach requires a complex model reconstruction and the geometry smoothing from CT scans can be extremely time consuming.
Stenosis (a reduced cross-sectional area) in the circulation system can result in many arterial diseases. Zhang et al. 22 studied the human common stenotic carotid arteries bifurcation. Sui et al. 23 used clinical measurement and CFD study to estimate the trend of hemodynamic parameters near atherosclerotic plaques stenosis in the carotid arteries. A potential threshold value of the CFD-based pressure gradient could serve as an additional factor for improving the accuracy of grading the severity of the stenosis 24 . The stenosis of a 0D coronary model was presented by changing the parameter to mimic the reduction of cross-section area 5 . The Fractional Flow Reserve (FFR), defined as the ratio of pressures distal and proximal to the lesion, is an important indicator in clinical diagnosis of stenosis 25 .
In this paper, we aim to develop an efficient lumped-parameter model. Our model includes the whole coronary network based on CT scans of a patient-specific geometry including the right coronary tree, which is absent in many previous mathematical models 5,7 . The model adopts the structured tree model boundary conditions based on the work of Olufsen et al. 11 , thus avoiding invasive perfusion measurements. We also incorporate the effects of the head loss at the two inlets of the large coronary arteries. In addition, effects of stenosis of different degrees in coronary vessels are modelled and analyzed using FFR.

Results
The pressure and the flow rate in coronaries. We first run a simulation based on our geometry derived from the CT scans (we refer to this as Simulation 1). The flow rates along different sections of the left coronary arteries are plotted against time in Fig. 1. The results are compared with the previous result by Pietrabissa et al. 7 . To separate the model differences and the patient-specific geometry difference, we also use their geometric parameters in our model and the results are also shown in Fig. 1 (Simulation 2). Note that Pietrabissa et al. 7 added an arbitrary resistance in their model, while our resistance is based on a rational estimation of the flow-dependent head loss. It is interesting to see that the difference in patient geometries only affects the flow rate in larger arteries (i.e., LAD, LCX, and LMCA). However, the model difference (using different resistance) is significant in most of the arteries. The diastolic dominated flow in the left coronary tree shows a typical decline during early systole. The wave forms changes slightly with locations. Due to the myocardial contraction, the blood flow rate drops at the beginning of systole and increases during diastole. We observe that the variation of coronary flow is location dependent due to the vascular pressure change.
The computed flow waveform in LCX is also compared directly to our own in vivo measurement using a Doppler catheter (with the patient's permission), as well as with the published in vivo measurements by Mynard et al. 26 , shown in Fig. 2. There is an overall agreement between our simulations and the in vivo data. Compared to the in vivo data, the simulated coronary hemodynamics is less wavy, presumably due to the fact that we ignored wave reflections in the system. Figure 2 also shows that patient-specific modelling is important in LCX.
Effect of head losses. We now study the effect of the head loss. The pressure and flow in both the RCA and PLA (with and without the local head losses) are shown in Fig. 3. We can see that the difference between the right and left coronary flow demonstrates the influence of the intra-myocardial pressure as a feedback to the cardiac system. We remark that the flow waveforms in PLA and RCA are much smaller than those without the head loss. Figure 3 shows the pressure and the corresponding flow waveforms in these vessels. We note that the pressure with head loss is in the physiological range. The phenomenon of the sudden increase of the resistance at the inlets of LAD and RCA was observed previously but not explained 5,7 . Our model suggests that this sharp increase of the resistance is due to the head losses.
Effect of coronary stenosis and FFR. The effect of a stenosis is investigated in the coronary system using our model. We define the stenosis degree α using a cross-sectional area reduction: α = A S /A 0 , where A 0 is the normal cross-sectional area and A S is the stenotic one. Then the resistance, capacitance, and the inductance of a stenotic vessel become 5 , where R 0 , C 0 , L 0 denote the resistance, capacitance, and the inductance of the non-stenotic vessel, respectively.
Here we assume that α = 100 − p, where p is the percentage of the area reduction. Figure 4 shows the details of the flow and pressure in the right coronary PDA, with the percentage of stenosis changing from 20% to 90%. The waveform of the flow in PDA decreases with the percentage of the stenosis in early systole. However, in diastole there is no significant flow drop for a stenosis below 50%. In cases of 70% or larger stenosis, the flow is greatly reduced in the late systole cardiac cycle. The results suggest that the pressure is perhaps a better indicator of the stenosis compared with the flow rate. This is because the pressure in PDA reduces monotonically with the narrowing rate, but the flow rate changes more erratically.
FFR can guide the percutaneous coronary intervention (PCI) with drug-eluting stents in patients with coronary artery stenosis. As can be seen in Fig. 5(a), FFR can be expressed as the ratio of the pressure distal to the lesion (p d ) and the pressure (p a ) proximal to the lesion 25 :  When FFR is smaller than 0.8 and the stenosis is less than 50%, stenting is recommended 20 . FFR at different levels of stenosis in PDA is shown in Fig. 5(b). The relation of FFR and stenosis was studied by Pijls & Bruyne 27 with clinical statistic data obtained from 221 patients before angioplasty or during diagnostic coronary angiography. Figure 5(b) shows that the present area reduction model gives a general agreement with the clinical results. Note that the clinical definition of stenosis is based on the radius ratio (i.e., α = r s /r, where r and r s are the normal and stenotic radii) instead of the area ratio. Therefore, the simulation result with the radius reduction is also shown for comparison. Our model can also predict the relative flow against degrees of stenosis, shown in Fig. 5(c). Again, the result agrees well with the experimental measurements by Mates et al. 28 .

Discussion
In this study, we propose a new model of simulating the coronary hemodynamics. Three new features are included in this model. Firstly, the whole coronary model geometry is reconstructed from CT scans of a real patient. Secondly, our work demonstrated that head loss at the inlets of the coronary vessels, where sudden diameter change occurs, can significantly improve the accuracy of the results. Finally, the model parameters of the root  impedance are obtained non-invasively based on a structured tree model. Using this model, the coronary blood flow rate and pressure of each vessel can be assessed quickly without time-consuming computations, and the results agree well with clinical observations and experiment data. In addition, effects of stenosis on the coronary flow are considered.
Our model yields reasonable results by considering the head loss. The FFR and relative flow in the coronary vessels for a stenotic vessel computed in this model also show good agreement with the experimental data published previously, although it is noted that the experiment data has a high level of noise.
Finally, limitations of the model are: Firstly, the model does not include specific dynamic input of the heart function. Additionally, we assumed that the vessel segments are rigid tubes, thus the model cannot predict wave propagations in the system. We also mention the limitation that measurements of coronary flow were not performed in some specific situation including an artery with large side branch. More sophisticated models can be developed to remove these limitations, although the challenges are to obtain patient-specific material properties and the spatial deformation of the whole heart.

Methods
Patient-specific coronary tree. The geometries of the coronary tree were obtained using the 64 slices CT scans of an adult male aged 65 with hypertension in the Second Affiliated Hospital of Nanjing Medical University, with the permission of the patient. We divide the coronary tree into several segments according to the American Heart Association standard 29 and the clinical data 30 .
As shown in Fig. 6, all the vessel terminals are connected to a structured tree model. The left and right coronary trees are linked to the aorta and the head loss resistance (R h ) at the junctions is included. The long right coronary artery (RCA) is further divided into three segments as in clinical practice. In this model, we consider both the left coronary artery (LCA) and the right coronary artery (RCA). This is because in addition to supplying the right ventricle with blood, RCA also supplies 25% to 35% of the blood flow to the left ventricle (LV) 31 . In particular, for right dominant patients, the RCA plays a key role in heart diseases.
The analogous circuit network for the whole coronary tree is illustrated in Fig. 7. In the network, all the terminals of the coronary arteries are assigned to a root impedance computed from the structured tree model. Similar  Mathematics of the lumped-parameter model. We assume the blood to be Newtonian fluids and the vessels to be cylindrical tubes 32 . The governing equations of the flow in each vessel are 33 where p is the pressure, Δp is the pressure drop, and q is the flow rate. The other parameters are: hydraulic resistance R, the elastic capacitance of each vessel C and the inertia of the blood flow L. Thus, analogous to resistance, capacitance, inductance in a circuit, we have 5,7 μ π where μ is the blood viscosity, ρ is the density, G is the Young's modulus of the vessels, l is the segment length, D is the diameter, and h is the wall thickness of the vessel. The parameters of all the vessels estimated from the CT data are listed in Table 1, where the heart rate of 75 beats per minute (cardiac cycle of 0.8 s) was measured. We also chose μ = 4 × 10 −3 kg · m −1 s −1 , ρ = 1.04 × 10 3 kg/m 3 and G = 2 × 10 5 Pa. The wall thickness h is 0.08D, with D being the vessel diameter.
As the blood supplies are from heart ventricles, the left ventricle (LV) and right ventricular (RV) pumping functions are chosen, following 34,35 to be where P v is the ventricular pressure, U 0 is a peak isovolumetric pressure at volume V 0 , E is the time-varying elastance, the activation function σ(t) = [1 − cos(2πt/t s )], t s is the systolic period, V is the reference volume,  V is the rate of volume change, R v is the resistance of ventricular myocardium, and E d is the passive elastance. All the pumping function model parameters are listed in Table 2. The terminal impedance. The terminal resistance of each vessel is linked to the capillaries downstream. In many previous models, the terminal resistance is estimated using a three-element Windkessel model 36 , where the total resistance and the compliance of the vascular bed need to be measured, which is often difficult to do. Hence, in this paper, we seek an alternative approach.
To estimate the terminal resistance, compliance and the inertia of the whole coronary tree, we made use of the structured tree model by Olufsen et al. 11 , as well as the perfusion measurements 7 . For simplicity, we only consider steady flow and the bifurcations of a symmetric structured tree when the diameter of vessels is less than 0.01 mm. These small arteries are treated as a structured network of straight vessels. The impedance at a vessel bifurcation is: (7) can be solved recursively to obtain the root impedance of the terminal vessel of the vascular beds 11 . Hence, the terminal impedance for a symmetric tree of N generations is where λ = 50 is a constant defining the length to radius ratio in small arteries 10,37 , r 0 is the original root radius, and α = r d /r p = 0.778 is the ratio of the daughter vessel radius (r d ) and the parent vessel radius (r p ) for a symmetric tree network when the arterial taper exponent rate is chosen to be 2.76 11,38,39 .
The values of Z at different ends of the coronary vessels (Z1-Z10) are listed in Table 3. These are compared with the corresponding perfusion resistance measurements 5,7 . A similar trend is seen between our estimations and the measurements.
Evaluation of the head loss. Energy losses associated with vessel structures and the loss coefficient can be incorporated into Bernoulli's equations 9 . To model the head loss (a sudden pressure-drop due to cross-sectional area change) at the inlets of the coronary arteries, as shown in Fig. 8 2 where ξ is the local resistance, estimated to be 2.88 from the empirical formula ξ = 288/Re by Fester et al. 40 , for a vessel constriction diameter ratio of 0.5 and Re = 100, and v is the mean velocity. The resistance R h due to the head loss is then h 2 (10) shows that the local head loss in the coronary circulation system is equivalent to a variable resistance that is proportional to the flow rate q.
The geometry measurements and flow experiments. The right and left coronary trees are reconstructed from the in vivo data from a male patient (age 61, 80 kgs). The GE 16-slice CT angiography was performed and the medical image software (GE,AW) was used for visualization and measuring vessel length and radius. The blood flow experimental data of left circumflex branch coronary vessels was measured by an intra-coronary Doppler guide wire (FloWire, 1400-Floppy). All the experiments were performed by The Second Affiliated Hospital of Nanjing Medical University, China. All the procedures were conducted according to guidelines established by The Second Affiliated Hospital of Nanjing medical University, and every effort was made to minimize patient suffering. This study was approved by the Human Experiment Committee of The Second Affiliated Hospital of Nanjing Medical University. The patient provided written consent and it was informed consent that was obtained in this paper.

Conclusion
In this work, a new lumped-parameter model based on the CT scans of coronary circulation is proposed. This model incorporates the effect of head loss at the junctions of inlets for the first time. Furthermore, the terminal resistances of the coronary arteries are computed by using structured tree models without invasive perfusion measurements. The results generated by our model are in good agreement with the previous studies and experimental data. In addition, the result shows that the estimated pressure and flow of different stenosis provide a  Table 3. Comparison of the terminal impedance estimated from the structured tree model (8) and published perfusion data. No published data are found for Z1 and Z2. useful non-invasive method for cardiac disease diagnosis. Our model can be used to provide detailed flow information of a patient-specific whole coronary circulation system due to various diseases, such as stenosis and local stiffness change, and thus can be readily translated into clinical practice and facilitate real-time clinical diagnoses.