Computational analysis of haemodynamic indices in synthetic atherosclerotic coronary netwroks

: Haemodynamic indices are widely used in clinical practice for deciding on a particular type of treatment. Low quality of the CT data and tachycardia complicate interpretation of the measured or simulated values. In this work, we present a novel approach for evaluating resistances in terminal coronary arteries. Using 14 measurements from 10 patients, we show that this algorithm retains the accuracy of 1D haemodynamic simulations in less detailed (truncated) geometric models of coronary networks. We also apply the variable systole fraction model to study the effect of elevated heart rate on the values of FFR, CFR and iFR. We conclude that tachycardia may produce both overestimation or underestimation of coronary stenosis signiﬁcance.


Introduction
Coronary artery disease (CAD) implicates the decrease of blood flow in the heart arteries due to the growth of the atherosclerotic plaque. CAD remains one of the leading reasons for disability or death in the world. The severity of the disease dictates a choice between noninvasive treatment (e.g. drugs administration) and invasive surgical procedures invoking stent installation. In modern medicine, the decision strongly depends on the analysis of the coefficients characterizing the haemodynamic conditions. Medical doctors widely use fractional flow reserve (FFR), coronary flow reserve (CFR), and instantaneous wave-free ratio (iFR) for evaluating the functional severity of epicardial coronary stenosis [1,2]. The measurement of the haemodynamic indices in clinics requires invasive, expensive procedures with possible side effects. Computational simulations basing on noninvasive data become an excellent alternative [3,4]. The works [5,6] present mathematical models for FFR evaluation. The iFR simulations are addressed in [7,8]. The computations of CFR are performed in [9,10]. The work [11] presents a benchmark comparison of four numerical approaches with clinical data on FFR.
Essential factors for correct clinical interpretation and decision making are the accuracy and sensitivity of the mathematical model to the quality of the input data and diagnostic implications of various coefficients. These aspects were partly studied in [9,[12][13][14][15][16]. E.g. ventricular pacing and tachycardia may potentially disturb coronary haemodynamics and change the values of FFR, CFR and iFR. Establishment and strict adherence to CCTA imaging protocols can significantly improve results of FFR estimation with mathematical models [17]. The introduction of numerical approaches into clinical practice means that FFR estimation technology will have to be implemented in various clinical centres with different CCTA imaging protocols and CT scan machines. The low quality of CT data makes small coronary arteries unrecognizable for segmentation algorithms. It leads to decreased geometric details of the coronary network and may affect calculated values of FFR, CFR and iFR.
In this work, we consider two important issues of the numerical evaluation of FFR, CFR and iFR, by 1D haemodynamic model of coronary flow. The first is the sensitivity of simulations to the quality of CT data. The second is the effect of the heart rate on the FFR, CFR and iFR variability. In Section 2.1 we briefly present 1D haemodynamic model of coronary flow. We extend it with the model of variable stroke volume (SV) and length of the systole depending on a heart rate (HR). This model was previously applied to the analysis of average coronary blood flow (CBF) during asynchronous pacing, and arrhythmias [18]. The input and modified patient data are described in Section 2.2. We assume that low-quality CT data contain less information on smaller vessels in the coronary network and remove most of them. Section 2.3 presents a novel recursive algorithm of the peripheral resistance distribution from the aortic root to distal vessels basing on Murray's law at every level. This algorithm is applied both to the input and reduced networks. Section 2.4 contains definitions of haemodynamic indices (FFR, CFR, iFR) used in this work. In Section 3.1 we compare haemodynamic indices in the input and reduced networks. In Section 3.2 we study the effect of HR elevation on the FFR, CFR and iFR. The results are discussed in Section 4.
We conclude that the new algorithm of peripheral resistance distribution throughout a network makes the mathematical model almost insensitive to the absence of the smallest arteries. It produces less or similar error than the distribution algorithm among terminal vessels regardless of the parent's vessels. It also does not require information on the patient's coronary dominance type, unlike other approaches [11]. We also conclude that tachycardia can lead to overestimating or underestimating stenosis significance due to the complex relationship between haemodynamic indices and heart rate.

1D mathematical model of the blood flow in the network of coronary vessels
The blood flow in the coronary vascular network and the aorta is simulated by a 1D reduced-order model of unsteady flow of viscous incompressible fluid through the network of elastic tubes. The details of this approach can be found in [19][20][21] The 1D models were adapted and applied to the coronary circulation in [14,18,22]. In this section brief description of the model is presented. The flow in every vessel is described by mass and momentum balances ∂V ∂t where t is the time, x is the distance along the vessel counted from the vessel's junction point, ρ is the blood density (constant), A(t, x) is the vessel cross-section area, p is the blood pressure, u(t, x) is the linear velocity averaged over the cross-section, ψ is the friction force µ is the dynamic viscosity of the blood. Relationship between pressure and cross-section is defined by wall-state equation: where ρ w is vessel wall density (constant), c is the velocity of small disturbances propagation in the material of the vessel wall, F(A) is monotone S-like function (see [23] for a review) whereÃ is the cross-sectional area of the unstressed vessel. At the vessel's junction points we impose mass conservation condition and continuity of the total pressure: where k is the index of the vessel, M is the number of the connected vessels, {k 1 , . . . , k M } is the range of the indices of the connected vessels, ε = 1,x k = L k for incoming vessels, ε = −1,x k = 0 for outgoing vessels.
The boundary conditions at the aortic root include the blood flow from the heart, which is set as a predefined time function Q H (t) In this work, we use a simple approximation of the heart outflow Q H (t) in the time domain. We define it as a half-sine function during ventricular systole and set it to zero otherwise where SV is the stroke volume of the left ventricle, T is the period of the cardiac cycle, τ is the duration of the systole. Thus, the stroke volume is and cardiac output (CO) by definition [24] equals The dependency of τ from HR is nonlinear. Previously we derived the function τ(T) as a regression basing on the results of simulations of action potential dynamics in human cardiac cells with the O'Hara-Rudy model [18] Experimental and clinical studies [25,26] report a linear relationship between SV and HR. Thus we use linear regression obtained in [18] We adjust SV 0 according to the patient's data on HR and SV. According to clinical studies [25,26] α = 2[bpm/ml] and this value vary only sligtly among patients. Finally, we use (11) and (12) for parametrisation of Q H (t) in (8) for various heart rates. The outflow boundary conditions assume that a terminal artery with index k is connected to the venous pressure reservoir with the pressure p veins = 8 mmHg through a hydraulic resistance R k . It is described by Poiseuille pressure drop condition where p k , A k , u k are blood pressure, cross-sectional area and blood velocity at the terminal point of k-th vessel. We use the same outflow condition for the aorta and coronary arteries.
The hyperbolic system (1) inside every vessel is numerically solved by the grid-characteristic method of the second order [27]. The systems of nonlinear algebraic equations in vessel's junctions (5),(6), aortic root (7) and at the end points of terminal arteries (13) are numerically solved by the Newton's method. To close the system of nonlinear equations, we add the mass conservation condition and second-order compatibility conditions of hyperbolic set (1) (see [28]).
The computational domain is the network of vessels, including the aortic root, aorta, left and right coronary arteries and their branches. Aorta and other systemic arteries are simulated as a single vessel with the length set to 80 cm and diameter set to 2.17 cm. We refer to this vessel as the 'aorta' since most of its properties correspond to the patient's aorta (diameter, elasticity). Aorta's parameters adjusted to represent compliance of a systemic circle and to get an adequate arteriovenous pressure drop (see Section 2.3). On the inflow of the aortic root, we set boundary condition (7). The aortic root splits into three branches: aorta left coronary artery (LCA) and right coronary artery (RCA). The sum of average blood flow through LCA and RCA is a coronary blood flow (Q cor in Section 2.3). The geometrical properties of the aortic root and coronary arteries are extracted from patients' CT scans.

Patient data
We investigate a cohort of 10 patients with 14 sites of FFR measurements. Patient's data are freely available and presented in detail in [11]. The data includes 3D surface meshes, 1D meshes, patient clinical data with FFR and cFFR values, and sketches from physicians indicating the approximate location of FFR measurement. Patient clinical data include blood pressure measurements (systolic, diastolic and average), body mass index (BMI), age, heart rate during CT acquisition. Table 1 briefly summarizes clinical data used in this work.
We use the patient's age to estimate pulse wave velocity [29] that equals to the velocity of small disturbances propagation c in (3). Stroke volume is estimated based on BMI values according to [30]. We use blood pressure measurements to estimate resistance for aortic outflow boundary condition (13). Resistances for terminal coronary arteries are set according to a recursive algorithm (see Section 2.3).
Patient-specific geometries are obtained from the CT images by the method of automatic CT scans processing [16]. We refer to these geometries as 'full' coronary trees. The method includes four stages: aorta segmentation; computation of Frangi vesselness; ostia points detection and coronary vessel segmentation; skeletonization of segmented vessels, and graph construction.
In this work, we investigate the sensitivity of haemodynamic indices to the degree of coronary network details. We synthetically generate 'cut' low-detailed coronary trees by cutting full trees. We remove all branches distal to the location of FFR measurement (provided there are no other sites of FFR measurements downstream). All junctions proximal to the location of the FFR-measurement are preserved to ensure adequate flow distribution. All junctions without distal FFR-measurement sites are removed. An example is shown in Figure 1. Cut trees represent synthetically generated coronary trees with low segmentation quality. We cut full trees as much as possible without affecting the FFR-measurement site. This cutting algorithm can also save computational time, provided it does not significantly affect the values of haemodynamic indices.

Terminal resistance
Myocardium compression is an essential feature of coronary haemodynamics. To account for the compression, we set R k = R k (t) for the boundary condition (13) in the terminal CAs. Similar to our previous works [6,14,16] we set time profile of R k (t) similar to the time profile of a cardiac output (8).
The peak value of the terminal resistance during systole is set to R max k = 3R k , where R k is the terminal resistance during diastole [11]. It is sufficient for the complete blockage of the flow in terminal CAs during systole.
The values of R k are set according to the following algorithm. First, we estimate the outflow resistance of an aortic part (R a ) and the total effective resistance of a coronary part (R cor ). Aortic 6 of 16 resistance is prescribed at the end of the virtual vessel representing aorta and systemic circulation. We assume that R a produces the pressure drop ∆P = P ave − p veins mmHg, where p veins = 8 mmHg is venous pressure [24], P ave is measured average blood pressure in systemic circulation and the ratio of coronary blood flow (CBF) to minute cardiac output (CO) is β = 0.05 (Q cor = βQ CO ). Thus we have where R cor is total resistance of coronary microcirculation and Q cor is average coronary blood flow (CBF). As a result, we have R cor = 19R a . In our simulations, these values produce the calculated ratio of β between 0.03 and 0.06, which belongs to the well known physiological range [24,31,32]. The value of β varies depending on the particular resistance of an arterial network. Second, we assign terminal resistances R k for each terminal coronary artery. We propose a recursive algorithm that splits resistances in each branching point according to Murray's law with a power of 2.27 [33]. Let us consider parent branch and N child branches with diameters d i , i = 1, 2, . . . , N (Fig. 2). Parent branch supplies a region of the myocardium with resistance R. Child branch with index i supplies a subregion with resistance R i . Our goal is to calculate resistances R i , i = 1, 2, . . . , N. We assume, that where p = 2.27 and Q i is an average blood flow through i-th branch. Using (16) we derive R i as a function of R 1 Finally, from (16) and (17) we get We use expressions (17) and (18) to calculate R i , i = 1, 2, . . . , N. The recursive algorithm starts by applying procedure (17), (18) to the roots of LCA and RCA. We consider R cor as a parent branch resistance (R in (18)) and LCA and RCA as two child branches (N = 2). We calculate resistances corresponding to LCA and RCA roots. After that, we treat them as parent branches and traverse all coronary networks with the depth-first search until we calculate resistance R k for each terminal artery. We will refer to this algorithm as a 'recursive algorithm'. Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 2 August 2021 doi:10.20944/preprints202108.0026.v1 A general method to distribute terminal resistances involves dividing R cor directly between terminal points of an arterial network, ignoring diameters of all branches in-between. We will refer to this approach as the 'terminal points algorithm'. This approach is typically modified to consider coronary dominance type [11] or arterial network resistance [34]. Our calculations use a basic version of the 'terminal points algorithm' since most of the modifications can be improved by implementing a recursive version.

Haemodynamic indices
We simulate stenosis as a separate part of the vessel with reduced diameter and increased velocity of small disturbances propagation c (3) (by a factor of 2). We perform two calculations for stenosed tree -under normal conditions and with vasodilatation (hyperemia). Vasodilatation is simulated with the same model with reduced terminal resistances R k of coronary arteries for the boundary condition (13). We reduce resistances by 70% [11,35].
Fractional flow reserve (FFR) is calculated as the ratio of average pressure in coronary artery distal to stenosis (P h dist ) to average aortic pressure (P h aortic ) during hyperemia.
Coronary flow reserve (CFR) is calculated as the ratio of average blood flow through stenosed vessel during hyperemia (Q h ) to average blood flow through stenosed vessel under nonhyperemic normal condition (Q n ).
Instantaneous wave-free ratio (iFR) is defined as the ratio between average pressure in coronary artery distal to stenosis (P w dist ) and average aortic blood pressure (P w aortic ) during the diastolic wave-free period (WFP) under the nonhyperemic normal condition. WFP begins in 25% of the way into diastole and ends 5 ms before the end of diastole in accordance with the general definition of iFR [36].

Haemodynamic indices sensitivity to coronary tree detalization
In this section, we present the results of calculating haemodynamic indices for two resistance distribution algorithms: recursive algorithm and terminal points algorithm (see Section 2.3). We compare the computed values for the cases of full trees and cut trees.
We start with evaluating blood flow distribution between coronary arteries. Figure 3 shows CBF distribution between LCA and RCA calculated with the recursive algorithm on full coronary trees. Patients 1, 4 and 5 have the highest ratio of blood flow in LCA to blood flow in RCA. These patients were marked as left-dominant type by a cardiologist. The rest of the patients were marked as right-dominant type. The mean value of LCA blood flow proportion is 79.6% among left-dominant patients (patients 1, 4, 5) and 58.5% among right-dominant patients (patients 2, 3, 6-10). These values are in good agreement with clinical data [31]: 76.1% and 57.8%, respectively. We conclude that the recursive algorithm automatically considers coronary dominance type and distributes blood flows accordingly. It reduces the need for cardiologist expertise during cardiovascular simulations.  Table A1.
To investigate the effect of coronary tree details on CBF distribution, we compare blood flows for full trees and cut trees (see Section 2.2). The average difference in blood flow through LCA between full trees and cut trees for the recursive algorithm of resistances distribution is 4%. For RCA, the average difference is 3%. For the terminal point algorithm of resistances distribution, the average difference in blood flow through LCA between full trees and cut trees is 10%. For RCA, the average difference is 18%. It demonstrates that the terminal points algorithm is quite sensitive to the degree of details of an arterial tree. This algorithm may produce significant errors when calculating haemodynamic indices.
Next, we investigate the effect of coronary tree details on haemodynamic indices. We compare sensitivity to the degree of details of a tree in two resistance distribution algorithms: recursive algorithm and terminal points algorithm. FFR, CFR and iFR values for full trees are considered to be the base values. We calculate haemodynamic indices for cut trees and present root-mean-square error (RMSE) relative to the base values. We also present normalized RMSE as a ratio between RMSE and the mean base index value (sum of all indices divided by the number of stenoses). Table 2. Sensitivity of FFR, CFR, iFR to coronary tree detalization for two cases: 1) resistances distributed with recursive algorithm, 2) resistances distributed with terminal points algorithm. RMSE -root-mean-square error, NRMSE -normalized root-mean-square error. RMSE Table 2 shows mean calculated values of haemodynamic indices, RMSE and NRMSE values between indices calculated on full trees and cut trees. For all patients haemodynamic indices were less sensitive to coronary tree detalization in case of recursive algorithm (Figures 4, 5, 6).
The most significant deviations are achieved for patient 1 (for terminal points algorithm). The distribution of flows between LCA and RCA for this patient changed significantly after cutting the full tree. Change in LCA blood flow for patient 1 equals 22% (most significant among all patients), which caused deviations in calculated indices. Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 2 August 2021 doi:10.20944/preprints202108.0026.v1

Haemodynamic indices sensitivity to heart rate
In this section, we calculate values of FFR, CFR and iFR for various heart rates from 60 bpm to 120 bpm. All presented results correspond to full trees and recursive algorithm of resistances distribution. We start with calculating FFR values for baseline values of HR and SV (table 1). Results are presented in Figure 7. RMSE between calculated and measured FFR is 0.04 or 5%, with sensitivity and specificity equal to 1.0. Figure 7. Difference between calculated (full trees, recursive algorithm) and measured FFR Next, we investigate sensitivity of haemodynamic indices to heart rate. We adjust stroke volume (12) and systolic duration (11) for each value of heart rate. . CFR for various heart rates (full tree, recursive algorithm) Figure 10. iFR for various heart rates (full tree, recursive algorithm) Figures 8,9, 10 demonstrate sensitivity of haemodynamic indices to heart rate. Exact values are presented in tables A2, A3, A4. In our previous works [6] we concluded that FFR decreases with an increase of HR or SV. It means that the haemodynamic significance of stenosis increases with cardiac output. However, in this work, we consider the relation between HR, SV and systolic duration. In most cases, haemodynamic indices do not change significantly within 60-120 bpm range. In some cases index decreases initially and then increases (FFR, patient 7, LAD, table A2). In other cases index steadily decreases (FFR, patient 4, table A2) or increases (CFR, patient 7, LCX, table A3). We assume that FFR, CFR and iFR behaviour depends on various factors: stenosis degree, artery diameter, patient's cardiac output, etc. Additional studies should be performed to investigate these factors further. Clinical studies [37] also show that the impact of the HR on haemodynamic indices can vary and depends on the artery and lesion properties.

Discussion
In this work, we studied the sensitivity of calculated haemodynamic indices to the degree of coronary tree and heart rate details. We presented a new recursive algorithm of terminal resistances distribution. This algorithm significantly reduces the model's sensitivity to the degree of details of a coronary tree and automatically considers coronary dominance type. In other works, this problem is partially solved by introducing additional coefficients [34], changing algorithm based on the coronary dominance type [11] or performing preliminary calculations [38]. All these approaches are valid and effective, but most of them can be further improved by implementing the recursive algorithm.
Decreased sensitivity to the details of an arterial tree means that calculated haemodynamic indices are less susceptible to low-quality CT images. It also means that we can synthetically reduce details of an arterial network to reduce calculation time without significant error.
Although the significant influence of HR on coronary haemodynamic is a known phenomenon [39], the literature examining the role of HR on FFR, CFR and iFR is sparse. Results of our study imply that tachycardia might be responsible for an overestimation or underestimation of haemodynamic indices. We did not find a clear pattern to describe a relation between HR and haemodynamic indices. Blood flow in the LAD occurs mainly during diastole and drops to zero (or is even reversed) during systole. However, systolic reduction of blood flow in the RCA is much less evident due to lower right ventricular pressure. It means that changes in systolic duration due to increased HR have a different effect on various arteries. iFR is measured during the diastolic wave-free period. Relative diastole duration shortens significantly at high heart rates, leaving less time for a blood flow to occur. As a result, iFR behaviour with increased HR can be different from FFR or CFR. We hypothesize that the increase in haemodynamic indices at higher heart rates can be attributed to decreased cardiac output during significant tachycardia. Lower values of blood flow tend to decrease the haemodynamic significance of a stenosis.
Lack of CFR and iFR measurements (as well as measurements of FFR for various heart rates) constrains our study. We assumed similar relation between SV and HR for all patients. This approach may produce increased error in FFR, CFR, iFR estimations at higher heart rates. Some further studies should be performed that include measurements of haemodynamic indices or at least measurements of stroke volume and systolic duration for various heart rates. The latter could identify the relation between SV, systolic duration and HR for each patient and improve our estimate of tachycardia effect on coronary circulation. These potential results can substantially affect clinical practice, especially regarding clinical decision-making when FFR, CFR, iFR values are borderline.

Abbreviations
The following abbreviations are used in this manuscript: