Computational Analysis of Coronary Blood Flow: The Role of Asynchronous Pacing and Arrhythmias

: In this work we present a one-dimensional (1D) mathematical model of the coronary circulation and use it to study the effects of arrhythmias on coronary blood ﬂow (CBF). Hydrodynamical models are rarely used to study arrhythmias’ effects on CBF. Our model accounts for action potential duration, which updates the length of systole depending on the heart rate. It also includes dependency of stroke volume on heart rate, which is based on clinical data. We apply the new methodology to the computational evaluation of CBF during interventricular asynchrony due to cardiac pacing and some types of arrhythmias including tachycardia, bradycardia, long QT syndrome and premature ventricular contraction (bigeminy, trigeminy, quadrigeminy). We ﬁnd that CBF can be signiﬁcantly affected by arrhythmias. CBF at rest (60 bpm) is 26% lower in LCA and 22% lower in RCA for long QT syndrome. During bigeminy, trigeminy and quadrigeminy, respectively, CBF decreases by 28%, 19% and 14% with respect to a healthy case.


Introduction
Coronary blood flow (CBF) supplies the myocardium tissue with oxygen and other essential nutrients. The distal coronary vessels are immersed in the myocardium. Myocardium contractions produce external pressure to the immersed blood vessels and substantially elevate the terminal hydraulic resistance. Thus, the dependence of the blood flow in coronary arteries (CA) on the characteristics of the heart cycle is an essential feature of CBF.
The characteristics of the heart cycle are controlled by the electrical activity of the sinoatrial node (SAN). The SAN activity is modified by a variety of factors: signals from the sympathetic and the parasympathetic nervous system and humoral factors. Some pathological processes and artificial electric stimulation of the myocardium are also among the factors, which can modify the heart activity.
Pathological change of the myocardium contractions due to asynchronous cardiac pacing and arrhythmias produces violations of the CBF, which, in turn, results in decreased myocardium supply with nutrients and possible ischemic events. In this work, we use computational modelling to study changes in CBF, which are produced by several types of pathological heart rhythms including interventricular asynchrony due to inappropriate pacing, several types of arrhythmia including bradycardia, tachycardia, long QT syndrome and premature ventricular contraction (bigeminy, trigeminy, quadrigeminy).
Interventricular asynchrony refers to the decoordination between RV and LV contraction due to the failure in the electrical conducting system of the heart or asynchronous cardiac pacing. Cardiac pacing is performed by the pacemakers. The pacemakers stimulate the heart by electrical impulses and prompt heart beating at a regular rate. Asynchronous electrical activation of the ventricles during ventricular pacing produces irregular patterns of the mechanical stress [1][2][3] which, in turn, results in abnormal ventricle contractions and violation of CBF. Finally, it causes deficiencies in perfusion and glucose uptake even in the absence of CA disease [4,5]. In this work, we perform numerical simulations to study the effects of the asynchronous performance caused by pacemakers on the CBF. These results are also valid for interventricular asynchrony caused by other factors. We simulate asynchrony by modifying terminal hydraulic resistance in coronary vessels according to their spatial position relative to the pacemaker. We add shifts in time-dependent hydraulic resistance functions according to the asynchronous pulse generation.
Abnormal tachycardia is a stable, permanent increase of atria and/or ventricular contractions at rest above 85-100 beats per minute (bpm). It decreases the atria and ventricular filling and, thus, decreases the heart output. Various abnormalities in action potential initiation and propagation are the common reasons of tachycardia [6,7]. Increased physical load and emotional stress are the possible reasons of tachycardia in elderly people. Myocardial infarction often produces ventricular tachycardia. Tachycardia is a significant reason for morbidity and mortality in patients with ischemic heart disease. Tachycardia itself may be a reason for ischemic events and disease due to the changes in CBF. Bradycardia is a stable, permanent decrease of atria and/or ventricular contractions at rest below 55 bpm. It decreases the heart activity. Analysis of CBF during both tachycardia and bradycardia is rarely addressed in the literature.
Long-QT syndrome is a life-threatening cardiac arrhythmia syndrome that may cause sudden cardiac death [8,9]. The problems in myocardium repolarisation after a heart contraction produce increased QT intervals on the electrocardiogram (ECG) by more than 480 ms. Long-QT syndrome is associated with high duration of systole and decreased length of diastole. It may decrease CBF with an associated decrease of nutrient delivery.
Premature ventricular contractions (PVC) produce extra, abnormal heartbeats that disrupt regular heart rhythm, sometimes causing a skipped beat or palpitations [10]. They can originate from dysfunctional Purkinje fibres, ventricular or atrial tissue. In this work, we focus on PVCs with a full compensatory pause: the following SAN impulse occurs on time based on the sinus rate. The heart rhythms with one, two, or three regular heartbeats between each PVC are called bigeminy, trigeminy, or quadrigeminy [11]. PVC may increase the risk of developing arrhythmias and lead to chaotic heart rhythms and cardiomyopathy. Similar to other cases of heart rhythm failure, PVC affects myocardium contractions and, thus, the average CBF.
Mathematical models of haemodynamics utilise a variety of approaches to simulate blood flow in a network of vessels [12,13]. Mathematical modelling of the impact of tachycardia and bradycardia on blood flow [1] and electrophysiology [14] has been a subject of many studies. Coronary circulation and related haemodynamic indices in the presence of atherosclerosis are also popular topics of mathematical modelling [2,[15][16][17][18][19]. To the best of our knowledge, the blood flow in CA during abnormal heart rhythms is rarely studied in clinics and by mathematical modelling. In this work, we present a modification of the previously developed one-dimensional (1D) mathematical model of the coronary circulation, which was applied to the simulations of blood flow in CA during cardiac pacing and tachycardia [20]. We apply the new model to the computational evaluation of changes in CBF during asynchronous cardiac pacing and some types of arrhythmia. Due to the lack of appropriate patient-specific data on coronary vessels' structure and properties, we use anatomically correct arterial networks [21] and typical values of cardiovascular characteristics (vessels' elasticity, cardiac output, etc.). A similar dataset can be collected or estimated by using non-invasive data such as CT scans, arterial pressure, lifestyle conditions (sex, age, body mass index, smoking, sport) of the patient and others [22][23][24]. In our previous work, we considered the stroke volume (SV) and the ratio between systolic and diastolic phases of the heart cycle as predefined constants based on known values from the literature. Our present model updates the length of systole depending on the heart rate (HR). This model is based on the simulated duration of an action potential. We also include a regression relating SV to HR, which is based on clinical data. We show that these novel features provide substantially different numerical results on the variations of CBF during tachycardia and bradycardia. We also apply the new model to simulations of CBF during interventricular asynchrony, some types of arrhythmias, including long QT syndrome and premature ventricular contraction (bigeminy, trigeminy, quadrigeminy).

1D Mathematical Model of Blood Flow in the Coronary Vascular Network
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. Reviews and details of the 1D models of haemodynamics can be found in [13,22,25]. The 1D models were adapted and applied to the coronary circulation in [20,26,27]. In this section, we briefly present this approach and describe some novel features of the model to account for interactions with myocardium contractions. The flow in every vessel is described by mass and momentum conservation in the form

∂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-sectional 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 blood. Blood density ρ = 1 g/cm 3 and blood viscosity µ = 4 cP. Properties and characteristics of the blood are considered to be constant throughout the computational domain. This assumption is considered to be accurate in arteries with diameters larger than 1 mm under physiological conditions [22,25]. The elasticity of the vessel's wall material is characterised by where c is the velocity of small disturbances propagation in the material of the vessel wall, and f (A) is the monotone S-like function (see [28] for a review) whereÃ is the cross-sectional area of the unstressed vessel.
The boundary conditions at the vessel's junctions include the mass conservation condition and the 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), 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 by the hydraulic resistance R k . It is described by Poiseuille's pressure drop condition The hyperbolic system (1) is numerically solved within every vessel by the second-order grid-characteristic method [29]. The systems of nonlinear algebraic equations, which represent boundary conditions at the vessel's junctions (5) and (6), aortic root (7) and at the end points of terminal arteries (8) are numerically solved by Newton's method. All formulations of boundary conditions include a compatibility condition along the characteristic curve of the hyperbolic system (1), which extends outside the integration domain for every incoming and/or outgoing vessel.
The computational domain is the network of vessels including the aortic root, aorta, left and right coronary arteries and their branches ( Figure A2). Parameters of the aorta, aortic root and coronary arteries correspond to the physiological values for an adult human (Table A1). The structure of the coronary vascular network ( Figure A1) was derived from a general anatomical model [21].

Effects of the Heart Rhythm on the Coronary Circulation
The distal coronary arteries are immersed in the myocardium. Thus, myocardium contractions cause a substantial effect on the blood flow in these vessels. The variations of the heart rhythm change how the myocardium works. There are three essential features of the heart function, which cause a substantial effect on the CBF: dependency of the length of systole on the HR; dependency of the SV on the HR and dramatic increase in peripheral resistance due to compression of the terminal CAs by the myocardium during systole.
The dependency of systole duration on HR is complex. In this work, we use the action potential duration at 80% repolarisation (APD80) in a single cardiac cell as an estimation of systole duration. The action potential waveform and APD of human cardiac cells were simulated with the O'Hara-Rudy model [30]. We use the implementation of this model with a custom C++ code using the Rush-Larsen integration technique [31] with an adaptive step [32]. The minimal time step was set to 5 × 10 −3 ms. The model was paced at every HR to steady-state (100 beats). The simulated restitution curve (i.e., APD80 dependence on HR) is shown in Figure 1. It gives the approximation where τ is the length of systole and T is the period of the cardiac cycle.
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. Thus, we have The heart outflow is a function of two parameters: the SV and the length of the ventricular systole τ. When studying cardiovascular events with different or variable heart rhythms, it is especially important that both SV and τ depend on HR. In the wide range of HR values, cardiac output (Q CO ) remains constant under pacing conditions [33,34]; i.e., Thus, SV should be inversely proportional to HR. However, some experimental and clinical studies on the hearts of conscious dogs, which were paced by an implanted right atrial electrode [35], and the hearts of normal human foetus, which were auditory stimulated by a sound emitter placed on the mother's abdomen [36], report a linear relationship between SV and HR. We use the data from clinical studies on eighteen vascular surgery patients having general anaesthesia, mechanical ventilation with tidal volume 6 mL/kg and a transesophageal atrial pacemaker [37]. In this study, increasing HR from 80 to 110 bpm caused a reduction in SV from 72 to 57 mL. We use these values to obtain the linear relationship In this work, we consider the range of HR from 40 to 160 bpm where the linear regression (13) is close to the inverse relationship. Remarkable nonlinear behaviour is observed for the HR values above 200 bpm [38].
Compression of the terminal CAs during systole by the myocardium is an essential feature of coronary haemodynamics. To account for this compression, we set R k = R k (t) for the boundary condition (8) in the terminal CAs. Similar to our previous works [20,27,39,40] we assume that the dimensionless time profile of R k (t) is the same as the dimensionless time profile of cardiac output (10).
The peak value of the peripheral resistance during systole is set to R max k = 3R k , where R k is the terminal resistance during diastole [41]. It is sufficient for the complete blockage of the flow in terminal CAs during systole. The values of R k are set by the following algorithm. We assume that the total arterio-venous resistance R total of the systemic circulation produces the pressure drop ∆P = 100 mmHg [42], thus Then we perform three steps. First, we split R total between the terminal resistance of the aorta R a = R total /0.9 and the total terminal resistance of the CAs R cor = R total /0.1. These values produce the ratio of CBF to CO of about 3-6%. Second, we split R cor between the total effective resistances of major CAs and their branches: right coronary artery (RCA) and branches of left coronary artery (LCA)(circumflex artery and left anterior descending artery). R cor is divided depending on the diameters of the major CAs according to Murray's law with a power of 2.27 [41]. Third, we divide the total effective resistances of major CAs between their terminal branches according to Murray's law with a power of 2.27. We note that the error in measuring diameters of the terminal CAs may be substantial due to the quality of CT data. The major CAs generally have good visualisation on CT. The error in measuring their diameters is relatively small. Thus, the second step allows us to decrease the total error of R k assessment.

Interventricular Dyssynchrony Due to Asynchronous Cardiac Pacing
Failures in the electrical conducting system of the heart or asynchronous cardiac pacing produce dyssynchrony between RV and LV contraction. We simulate interventricular dyssynchrony by introducing a time difference of 30 and 60 ms in resistance functions (14) for the terminal branches of RCA. The disagreements among regional pre-ejection periods greater than 50 ms are considered as significant [43]. We associate the early right ventricle contraction with a pacemaker location in the right ventricle and the late right ventricle contraction with a pacemaker location in the left ventricle.
We calculate the ratios of the average blood flow in the LCA and RCA in the asynchronous and normal (synchronous) pacing conditions at different HRs. The relative change in the average blood flow in LCA is less than 0.5% in all cases because resistance functions (14) for the left CAs are always synchronous to the heart outflow (7) and, therefore, blood flow in LCA is synchronised to the contractions of the left ventricle.
From Figure 2 we observe that interventricular asynchrony causes significant changes in the average blood flow in RCA at the normal values of HR. The 60 ms early contraction of the right ventricular produces 12% increase, while 60 ms late contraction produces approximately the same decrease. The elevation of HR leads to a decrease in the relative changes in blood flow in RCA, which tends to the value corresponding to the normal pacing conditions.
Since both RCA and LCA are the branches of the aorta, the left ventricle contractions supply both RCA and LCA with the blood. It explains the asymmetric behaviour of the relative change of the blood flow in RCA and LCA.

Tachycardia and Bradycardia
In this section, we present results of CBF simulations at constant HR in the range from 40 to 160 bpm. The values of HR lower than 55 bpm are generally associated with bradycardia, and the values of HR higher than 85 bpm are related to tachycardia. We take into account the changes in systole duration (9) and stroke volume (13) with the variations in HR. Figure 3 shows the dependency between the ratio of average CBF and the average cardiac output (12) in the cases of variable (9) and constant (35%) ratio of the systole duration to the heart period from HR. From Figure 3 we observe a substantial difference in the fraction of CBF depending on the model assumptions. A constant fraction of systole duration leads to the elevation of the fraction of CBF with the increase in HR. It contradicts the clinical data [44]. The model with a variable fraction of systole results in the decrease in CBF fraction with the increase in HR. According to (9) the length of systole increases with the decreasing period of the cardiac cycle (and increasing HR). It leads to the increase in the relative length of systole and decrease in the relative length of diastole within a cardiac cycle. Thus, the relative time period with increased peripheral resistance (14) increases.
We also note that the absolute value of the average CO achieves its maximum at 100-120 bpm (see Figure 4). Figure 5 presents the corresponding values of the average CBF during the constant and variable fraction of systole. From Figure 5 we observe that average CBF achieves its maximum at 100 bpm for the variable and at 100-120 bpm for the constant fraction of systole. The early decrease in the average CBF, which is simulated by the model with a variable fraction of systole, accounts for the integrated effect of the decrease in both SV and diastolic phase of the cardiac cycle. The sensitivity of the average CBF to the fraction of systole is low for low values of HR. For HR values above 80 bpm, the assumption of a constant fraction of systole produces overestimated values of the average CBF. The difference of the absolute values of average CBF at 50, 100 and 160 bpm from the standard value at 60 bpm (see Figure 4) is less than 15% despite the permanent decrease in the fraction of average CBF to average CO (see Figure 3A), which produces 20% difference at 160 bpm.

Long QT Syndrome
Long QT interval on ECG, more than 480 ms for 60 bpm, reflects the delay in myocardium repolarisation and myocardium relaxation after a heart contraction. It is unclear how the length of systole changes with HR in the case of a long QT syndrome. We assume that this relation is similar to (9): where T LongQT is a period of the heart cycle during long QT syndrome, and a, b are constants. We estimate constants a and b from clinical data: for 60 bpm the length of systole τ LongQT = 480 ms [8]; for 113 bpm the length of systole τ LongQT = 353 ms [9]. As a result, we derive a = 624 ms and b = −143774 ms 2 . Figure 6 shows the results of blood flow simulations in RCA and LCA during long QT syndrome for values of HR from 60 to 120 bpm. From Figure 6 we observe that long QT syndrome produces a significant decrease in CBF in both RCA and LCA (Figure 6), which is more than 25%. It accounts for the decrease in the diastolic period, which is accompanied by low values of peripheral resistance according to (14). Elevation of HR results in an increase in blood flow in both RCA and LCA, which tends to the baseline values. The baseline value is the blood flow simulated by the model with a variable fraction of systole according to (9).

Premature Ventricular Contraction
PVC is an abnormal heartbeat, which occurs earlier than a scheduled normal heartbeat, which should be initiated by an action potential. The normal heartbeat is missed because the ventricles have been emptied by PVC. PVC also causes an effect on the subsequent normal pulses due to the increased ventricular filling after PVC. A regular PVC pattern includes PVC every second (bigeminy), third (trigeminy) or fourth (quadrigeminy) cardiac cycle.
We simulate PVC by modifying heart outflow function (7) in the following way. We suppose that systole starts 0.25T earlier than a periodic schedule; SV of the PVC is 71% less than SV of a normal beat due to insufficient ventricle filling [11]; the beat after PVC occurs on a schedule; the total length of PVC heartbeat and the following beat is 2T; SV of the heartbeat after PVC is 18% more than SV of a normal beat due to increased time of ventricle filling [11]. Figure 7 shows an example of modified heart outflow in the case of quadrigeminy and a HR value of 120 bpm.
We simulate average CBF in the cases of bigeminy, trigeminy and quadrigeminy for HR values from 40 to 120 bpm. Figure 8 presents the results, which include the fraction of the average CBF overall cardiac cycles of PVC pattern to the value of average CBF without PVC. From Figure 8, we observe a substantial decrease in the relative average CBF at low and normal HR. This value decreases with the increasing HR. The most pronounced effect (more, than 25% decrease at 40 bpm and 30% decrease at 120 bpm) is observed in the case of bigeminy as it produces the most frequent occurrence of PVC.

Discussion
In this work, we have presented a modified model of 1D haemodynamics in the network of coronary vessels, which includes dependencies of systole duration and stroke volume on the length of the cardiac cycle (or the heart rate). We have simulated interventricular dyssynchrony, which is caused by asynchronous cardiac pacing of the left and right ventricles. We have analysed coronary blood flow in the cases of several abnormal heart rhythms, including tachycardia, bradycardia, long QT syndrome and premature ventricular contraction.
The variable length of systole and SV produces results that are in agreement with clinical data [44]. In contrast to our previous works [20,27,39] we suggest that only terminal resistance rises due to myocardium contractions during systole, which is different from the models of general muscle contractions [26,45] and better corresponds to the anatomical features of major CAs laying outside the myocardium.
Based on experimental and clinical observations [35,36] we used linear regression for the dependency of SV on HR (13), which works well for the arrhythmia and pacing conditions. We note that this relationship becomes nonlinear for the values of HR above 160 bpm [38], and our model is not valid in this case. Some other factors (e.g., regulatory processes) may produce nonlinear behaviour even in the range from 40 to 160 bpm (e.g., physical exercise conditions).
We implicitly assumed that the length of the heart systole is approximately equal to the ventricular systole. We neglected atria systole as all applications in this work relate to the activity of the ventricles.
This assumption may produce some systematic error to the results. More detailed models are needed for the accurate simulations of the heart rhythm. We also use the single-cell model [30] for simulating action potential duration restitution curves in a cardiac cell. We assumed that the dependency of action potential duration at 80% repolarisation (APD80) on the duration of a cardiac cycle correlates with the similar dependence of the length of the ventricular systole (9). A more accurate approach should simulate the propagation of the action potential over myocardium and mechanical contractions. We imply that this process is instantaneous. It also may be a source of systematic error. We note that approximation (9) correlates with the results of QT interval measurements for different values of HR derived from ECG data. The length of the QT interval is associated with the length of the electrical and mechanical systole [46].
The cardiac muscle mechanical contraction is a result of Ca 2+ concentration increase in the cytoplasm, which, in turn, is induced by sarcolemma depolarisation and consequent calcium-induced-calcium release. Thus, despite quantitative differences in duration, all stages of the mechanical contraction are closely interrelated. Moreover, previous experiments on non-failing human hearts demonstrated that the difference between calcium transients measured at 80% recovery (CaTD80) and APD80 are approximately the same at every pacing cycle length (approximately 80 ± 20 ms for subendocardial tissue) [47]. A similar correlation was observed between twitch force duration and calcium transients in rat papillary muscle [48]. On the other hand, cardiac contraction is a complex, spatially distributed process that could be affected by numerous factors. For example, the blood vessel pressure changes and consequent baroreflex affect the sympathetic system, which has unequal effects on action potential, calcium transients and mechanical contraction [49]. However, for the purpose of this study, which is to qualitatively investigate the principal effects of HR on the CBF, APD80 provides a reasonable estimate of the mechanical contraction duration.
The developed model of CBF under various heart rhythms allows simulating the effects of medical treatment. We mention anticoagulant therapy as the primary strategy to prevent atrial fibrillation (AF) complications, including thromboembolism and stroke. Oral anticoagulants (OACs) such as vitamin K antagonists, direct thrombin inhibitors, or factor Xa inhibitors are administered routinely in patients with AF [50,51]. Anticoagulation therapy changes rheological properties of the blood and modifies oxygen and nutrients delivery to the myocardium by CBF. Adherence to OACs is the extent to which a patient takes his medication as prescribed. It is estimated as appropriate in 76.6% of patients with hypertension and other cardiovascular diseases [52]. Persistence is the act of taking drugs for the prescribed treatment duration. Persistence with OACs therapy is about 50% in a two-year perspective [53]. Thus, non-adherent and non-persistent patients make a tangible contribution to negative treatment outcomes. Multiple patient-level factors, as well as social, economic, health system causes contribute to poor adherence and persistence to OACs making this phenomenon a challenge for population studies [54]. Mathematical modelling of the CBF during arrhythmia episodes controlling for OACs therapy adherence and persistence may clarify their effects and give predictions with respect to possible negative outcomes. Some approaches to mathematical modelling of medication therapy effects on the CBF have been proposed in [55]. Anticoagulant therapy is not considered yet. We anticipate further elaboration of our model with respect to anticoagulation therapy as a possible ongoing endeavour.

Conflicts of Interest:
The authors declare no conflicts of interest.

Abbreviations
The following abbreviations are used in this manuscript: Appendix A Figure A1. 3D anatomical model of the coronary arteries [21].  Table A1. Table A1. Parameters of the arterial network in Figure A2. k is the index of the vessels, l k is length, d k is diameter, c k is pulse wave velocity index (see (3)).