Impact of Arrhythmia on Myocardial Perfusion: A Computational Model-Based Study

: (1) Background: Arrhythmia, which is an umbrella term for various types of abnormal rhythms of heartbeat, has a high prevalence in both the general population and patients with coronary artery disease. So far, it remains unclear how different types of arrhythmia would affect myocardial perfusion and the risk/severity of myocardial ischemia. (2) Methods: A computational model of the coronary circulation coupled to the global cardiovascular system was employed to quantify the impacts of arrhythmia and its combination with coronary artery disease on myocardial perfusion. Furthermore, a myocardial supply–demand balance index (MSDBx) was proposed to quantitatively evaluate the severity of myocardial ischemia under various arrhythmic conditions. (3) Results: Tachycardia and severe irregularity of heart rates (HRs) depressed myocardial perfusion and increased the risk of subendocardial ischemia (evaluated by MSDBx), whereas lowering HR improved myocardial perfusion. The presence of a moderate to severe coronary artery stenosis considerably augmented the sensitivity of MSDBx to arrhythmia. Further data analyses revealed that arrhythmia induced myocardial ischemia mainly via reducing the amount of coronary artery blood ﬂow in each individual cardiac cycle rather than increasing the metabolic demand of the myocardium (measured by the left ventricular pressure-volume area). (4) Conclusions: Both tachycardia and irregular heartbeat tend to increase the risk of myocardial ischemia, especially in the subendocardium, and the effects can be further enhanced by concomitant existence of coronary artery disease. In contrast, properly lowering HR using drugs like β -blockers may improve myocardial perfusion, thereby preventing or relieving myocardial ischemia in patients with coronary artery disease. R-R interval irregularity severities coronary artery


Introduction
Arrhythmia is a common cardiovascular disease featured by abnormally low/high or irregular heartbeat. The prevalence of arrhythmia is high, especially in patients with basic cardiovascular disease [1][2][3]. Arrhythmia may cause various symptoms, such as palpitations, dyspnea, heart failure, and confusion, or even lead to sudden death [4][5][6]. Clinical studies have demonstrated that a physiologically low heart rate is associated with lower cardiovascular mortality [7][8][9], and that in patients with stable ischemic heart disease and left ventricular systolic dysfunction, reducing heart rate helps to improve prognosis [8]. On the other hand, tachycardia has been found to be an independent risk factor for poor outcome after coronary revascularization [9]. A theoretical study has demonstrated that reducing heart rate is beneficial to the distribution of blood toward the subendocardium, which is more susceptible than the subepicardium to ischemia in the presence of coronary artery disease [10]; this may partly explain the aforementioned clinical findings. However, it remains unclear how different types of arrhythmia combined with coronary artery disease would affect myocardial perfusion, especially the severity of myocardial ischemia.
Quantitative measurement of myocardial perfusion parameters for assessing the severity of ischemia during the onset of arrhythmia is usually hard to implement in clinical settings, which is often replaced by diagnosing the absence or presence of myocardial ischemia through analyzing the characteristics of electrocardiograms [11]. Relatively, computational modeling provides a more convenient means of quantifying the impact of arrhythmia on myocardial perfusion and uncovering underlying mechanisms. Computational models of the coronary circulation have been widely applied to study the characteristics of coronary blood flow in the context of arrhythmia or coronary artery disease [10,12,13]. However, few studies have been dedicated to addressing the impact of arrhythmia on myocardial perfusion. In the present study, a computational model of the coronary circulation coupled to the global cardiovascular system was employed to address the issue. In addition, we proposed an index for assessing the severity of transient myocardial ischemia based on the relationship between the amount of coronary blood supply and the oxygen demand of myocardium.

Computational Model of the Coronary Circulation
The computational model was adapted from the models developed in our previous studies where the modeling methods and corresponding numerical schemes have been described in detail [10,14]. In brief, a zero-one-dimensional (0-1-D) multi-scale modeling method was employed to represent the coronary circulation coupled to the global cardiovascular system. The 1-D modeling method was used to represent large epicardial coronary arteries and main systemic arteries, while the 0-D (i.e., lumped parameter) modeling method was applied to represent other portions of the cardiovascular system, including intramyocardial vessels, systemic microvasculature, veins, the pulmonary circulation, and the heart (see Figure 1). Coupling the 0-D and 1-D models yielded a closed-loop model capable of describing both coronary and systemic hemodynamics. In particular, the model allowed for incorporation of varying heart rate and coronary artery disease, thereby providing a flexible platform for investigating the characteristics of coronary hemodynamics and myocardial perfusion under various arrhythmic conditions combined with various severities of coronary artery disease. Hereafter, we will give some details on mathematical modeling of hemodynamics in the coronary circulation. The governing equations for hemodynamics in each coronary artery were expressed in the form of 1-D mass and moment conservation equations that describe the time-dependent changes in vascular cross-sectional area (A) and volumetric blood flow rate (Q) along the axis of the artery [10].

∂A ∂t
∂Q ∂t Here, t, x, and ρ denote the time, axial coordinate and density of blood (which herein was set to be 1060 kg/m 3 ), respectively; and A, Q and P represent the cross-sectional area of the artery, volumetric flow rate and blood pressure, respectively. α and F r refer respectively to the momentum flux correction coefficient and friction force per unit length, which were herein set to 4/3 and −8πυ (υ is the kinematic viscosity of blood) based on the assumption that the cross-sectional velocity profile is parabolic [15].
ner, was introduced to the model, the automatic regulatory mechanism of coronary blood flow was incorporated by means of tuning the resistances of coronary arterioles (i.e., R1 in the left upper panel of Figure 1) via a proportional-integral-derivative (PID) feedback control loop to preserve coronary blood flow against the decrease in perfusion pressure caused by coronary artery disease [10]. Incorporating the regulatory mechanism enabled the model to represent the long-term influence of coronary artery disease on coronary hemodynamics and myocardial perfusion.

Modeling of Arrhythmias
Under resting conditions, healthy subjects usually have a heart rate (HR) of 50~90 bpm (beats per minute). Bradycardia and tachycardia are two common types of arrhythmia, which are defined as resting HRs of less than 50~60 bpm [4] and HRs of over 100 bpm [5], respectively. There is another common type of arrhythmia characterized by irregular patterns of electrical activation, i.e., R-R intervals do not follow a repetitive pattern [18]. Based on the ranges of HRs reported in the literature [19,20], we set the baseline resting HR at 80 bpm, and the lower and upper limits of HR to be 36 bpm and 184 bpm, respectively. In this study, we assumed that resting HRs below 50 bpm represent bradycardia, while HRs above 100 bpm represent tachycardia. Given the existence of three unknowns (i.e., A, Q and P) in the system of Equations (1) and (2), an equation describing the relationship between P and A was further introduced so that the equation system could be closed [10,15]: where A 0 and r 0 represent the cross-sectional area and radius of artery at the reference pressure (P 0 = 80 mmHg), respectively. E is the elastic modulus of arterial wall, which together with wall thickness (h) and Poisson's ratio (σ) of wall materials determines the stiffness of arterial wall. τ σ and τ ε are viscoelastic parameters representing the relaxation times for constant stress and strain of arterial wall, respectively. At the junctions of coronary arteries, the conditions of mass conservation and total pressure continuity were imposed so that hemodynamic variables in the entire coronary arterial tree could be solved as a whole, and such conditions were also adopted to link hemodynamic variables at the distal ends of coronary arteries to those in the 0-D model of downstream intramyocardial vessels.
To simplify the modeling work, intramyocardial vessels downstream from each distal coronary artery were allocated to three myocardial layers, namely the subepicardium, midwall and subendocardium. Applying the 0-D modeling method to the intramyocardial vessels yielded a network of parameters that represented the lumped properties (i.e., compliance (C), resistance (R) and inertance (L)) of various types of vessels in different myocardial layers (see the left upper panel of Figure 1). Accordingly, the governing equations for hemodynamics were defined at the C and L elements following the principle of mass and momentum conservation [16]: where P j and P j−1 denote the blood pressures at the C j and C j−1 elements, respectively, which are calculated as P j = V j /C j + P im and P j−1 = V j−1 /C j−1 + P im , where V j and V j−1 are the volumes of blood contained by the C j and C j−1 elements, and P im the intramyocardial tissue pressure which changes with cardiac contraction/relaxation and differs among myocardial layers [10]. Q j and Q j+1 represent, respectively, the inflow to and outflow from the C j element; R j and L j are the resistance and inertance that link the C j and C j−1 elements. It is noted that the values of C and R were not kept constant but varied with V to represent the effects of the large deformation of intramyocardial vessels induced by the changes in P im during each cardiac cycle (please see reference [10] for more details).
The system of the 1-D partial differential equations governing hemodynamics in the coronary and systemic arteries and the system of 0-D ordinary differential equations describing hemodynamics in the intramyocardial vessels and other portions of the cardiovascular system were solved with the two-step Lax-Wendroff method and the fourth-order Runge-Kutta method, respectively. The solutions of the two equation systems were communicated at the interfaces between the 1-D and 0-D models where an iterative numerical scheme was implemented to enforce conservation of mass flux and momentum across the interfaces. More details on the numerical algorithms have been described elsewhere [17].
The parameters of the model were initially assigned based on relevant data reported in the literature so that the model-predicted values of major coronary/systemic hemodynamic variables fell in the population-averaged ranges of in vivo data acquired from subjects free of arrhythmia or coronary artery disease (please refer to [10] for more details), thereby obtaining a baseline model to which various types of arrhythmia or coronary artery disease were later introduced. It is noted that when simulating arrhythmia by varying heart rate (HR), other model parameters were held at their default values so that the model-simulated results could reflect the pure effects of short-term HR variations; whereas, as coronary artery disease, which generally develops in a chronic manner, was introduced to the model, the automatic regulatory mechanism of coronary blood flow was incorporated by means of tuning the resistances of coronary arterioles (i.e., R 1 in the left upper panel of Figure 1) via a proportional-integral-derivative (PID) feedback control loop to preserve coronary blood flow against the decrease in perfusion pressure caused by coronary artery disease [10]. Incorporating the regulatory mechanism enabled the model to represent the long-term influence of coronary artery disease on coronary hemodynamics and myocardial perfusion.

Modeling of Arrhythmias
Under resting conditions, healthy subjects usually have a heart rate (HR) of 50~90 bpm (beats per minute). Bradycardia and tachycardia are two common types of arrhythmia, which are defined as resting HRs of less than 50~60 bpm [4] and HRs of over 100 bpm [5], respectively. There is another common type of arrhythmia characterized by irregular patterns of electrical activation, i.e., R-R intervals do not follow a repetitive pattern [18]. Based on the ranges of HRs reported in the literature [19,20], we set the baseline resting HR at 80 bpm, and the lower and upper limits of HR to be 36 bpm and 184 bpm, respectively. In this study, we assumed that resting HRs below 50 bpm represent bradycardia, while HRs above 100 bpm represent tachycardia.
The third type of arrhythmia was represented by imposing irregular cardiac R-R intervals (herein a R-R interval = 60/HR) during a period of successive cardiac cycles. In the present study, we constructed a probability density function to control the distribution of random R-R intervals.
Here, x denotes random cardiac R-R intervals with a mean value of µ (herein set to 750 ms, which corresponds to the baseline HR of 80 bpm), and standard deviation of σ. By changing the value of σ, one can adjust the range of x distribution so as to control the severity of R-R interval irregularity. In this study, the value of σ was varied from 100 ms to 500 ms at an interval of 100 ms to represent the increasing severity of R-R interval irregularity. In addition, in consideration of the fact that all R-R intervals should be physiologically reasonable (e.g., negative R-R intervals are physiologically nonexistent), we only applied Equation (6) to calculate the probabilities of R-R intervals falling in the range of µ − σ and µ + σ. With the constriction, the largest range of R-R intervals was from 250 ms to 1250 ms when σ = 500 ms, which was roughly in line with the range of HRs measured in clinical settings [19]. Panels (A) to (E) of Figure 2 show the probabilities of R-R intervals between µ − σ and µ + σ. In order to generate a sequence of random R-R intervals for 100 cardiac cycles, we first generated 100 R-R intervals based on the probabilities of R-R intervals, and then randomly sequenced them using a random operator embedded in MATLAB. Figure 2F shows the obtained random R-R intervals for 100 cardiac cycles given different values of σ.
Mathematics 2021, 9, x FOR PEER REVIEW 5 of 16 The third type of arrhythmia was represented by imposing irregular cardiac R-R intervals (herein a R-R interval = 60/HR) during a period of successive cardiac cycles. In the present study, we constructed a probability density function to control the distribution of random R-R intervals.
Here, x denotes random cardiac R-R intervals with a mean value of μ (herein set to 750 ms, which corresponds to the baseline HR of 80 bpm), and standard deviation of σ. By changing the value of σ, one can adjust the range of x distribution so as to control the severity of R-R interval irregularity. In this study, the value of σ was varied from 100 ms to 500 ms at an interval of 100 ms to represent the increasing severity of R-R interval irregularity. In addition, in consideration of the fact that all R-R intervals should be physiologically reasonable (e.g., negative R-R intervals are physiologically nonexistent), we only applied Equation (6) to calculate the probabilities of R-R intervals falling in the range of μ − σ and μ + σ. With the constriction, the largest range of R-R intervals was from 250 ms to 1250 ms when σ = 500 ms, which was roughly in line with the range of HRs measured in clinical settings [19]. Panels (A) to (E) of Figure 2 show the probabilities of R-R intervals between μ − σ and μ + σ. In order to generate a sequence of random R-R intervals for 100 cardiac cycles, we first generated 100 R-R intervals based on the probabilities of R-R intervals, and then randomly sequenced them using a random operator embedded in MATLAB. Figure 2F shows the obtained random R-R intervals for 100 cardiac cycles given different values of σ.

Indices for Evaulating Myocardial Perfusion
In this study, the subendocardial-to-subepicardial flow ratio (Q endo /Q epi ) was used to evaluate the transmural distribution of blood flow in the myocardium [10]. A low value of Q endo /Q epi indicates the redistribution of coronary blood flow from the subendocardium toward the subepicardium, which may increase the risk of ischemia in the subendocardium. To more clearly evaluate the severity of myocardial ischemia, we further Here, Q rs,0 , T 0 and PVA rs,0 represent, respectively, the reference values of the mean blood flow rate in a specific coronary artery during a cardiac cycle, the time period of a cardiac cycle and the left ventricular pressure-volume area (PVA) under normal resting conditions, with their counterparts under other physiological or pathological conditions being denoted by Q rs , T and PVA rs , respectively. Left ventricular PVA is the total area circumscribed by the end-systolic pressure-volume (P-V) line, end-diastolic P-V curve and systolic P-V trajectory, which is often taken as a measure of the total workload of the left ventricle, and has been found to correlate linearly with myocardial oxygen consumption [21][22][23]. Since oxygen supply to the myocardium in a cardiac cycle is proportional to the amount of coronary blood supply (i.e., Q rs ·T), the ratio of Q rs ·T to PVA rs can be taken as a measure of the balance between coronary blood supply and myocardial perfusion demand. If we assume that Q rs,0 ·T 0 /PVA re,0 represents the optimal state, when Q rs ·T/PVA rs = Q rs,0 ·T 0 /PVA rs,0 , the value of MSDBx is zero according to Equation (7), which represents the optimal (or baseline) state of myocardial perfusion; accordingly, a decrease in Q rs ·T/PVA rs relative to Q rs,0 ·T 0 /PVA re,0 will lead MSDBx to increase above zero, indicating that myocardium is in an ischemic state. By definition, the larger the value of MSDBx is, the more severe myocardial ischemia is.

Results
In our numerical experiments, apart from varying HR, we also considered coronary artery disease by introducing a stenosis to a branch of the LAD (see Figure 1 where the artery is named as LAD b1 ). The severity of the stenosis, herein quantified by diameter stenosis rate (SR), was varied from 20% to 70%. Accordingly, the values of Q endo /Q epi and MSDBx used for evaluating the condition of myocardial perfusion were calculated based on the model-simulated hemodynamic data in LAD b1 and intramyocardial vessels distal to it. The aforementioned data treatments are applicable to all the relevant results to be presented later, unless stated otherwise.

Hemodynamic Effects of Variations in HR
The panels A and C of Figure 3 show the model-simulated blood flow waveforms in LAD b1 and left ventricular P-V loops under different HR conditions. The corresponding amount of coronary blood flow (Q T , which is obtained by integrating the transient blood flow rate over a cardiac cycle) and left ventricular PVA are shown in Figure 3B,D, respectively. Variations in HR induced significant changes in both the coronary blood flow waveform and P-V loop. Quantitatively, although Q T decreased monotonously following the increase in HR, PVA showed more complex responses, for instance, PVA increased with HR when HR was lower than 60 bpm, but decreased with HR when HR was higher than 60 bpm. In addition, the range of changes in Q T was overall larger than that of PVA, indicating that Q T is more sensitive to variations in HR.

Effects of Variations in HR on Myocardial Perfusion
Myocardial perfusion was firstly evaluated with respect to the sensitivity of Qendo/Qepi to variations in HR. From the results presented in Figure 4A, Qendo/Qepi decreased evidently with the increase in HR when HR was higher than 60 bpm; relatively, Qendo/Qepi was insensitive to variations in HR when HR was low (e.g., < 60 bpm). In addition, the presence of a coronary artery stenosis of moderate to high severity (i.e., SR = 50~70%) considerably reduced Qendo/Qepi and augmented its sensitivity to variations in HR. These results indicate that tachycardia and coronary artery disease play additive roles in depressing the supply of blood to the subendocardium, thereby increasing the risk of subendocardial ischemia. Figure 4B shows the relationships between MSDBx and HR under various coronary artery stenotic conditions. MSDBx exhibited an overall trend of rising with the increase in HR, with the trend being most remarkable when HR was lower than 100 bpm. Recalling the results presented in Figure 3 and the definition of MSDBx (i.e., Equation (7)), the increase of MSDBx with HR is caused mainly by the reduction in the amount of coronary artery blood flow supplied to the myocardium in each cardiac cycle following the increase in HR. Relatively, varying the severity (i.e., SR) of coronary artery stenosis had little influence on MSDBx when the stenosis severity was at the mild to moderate level, because under such conditions the stenosis-caused reduction in coronary artery flow could be well compensated for by the coronary automatic regulatory mechanism [24]. When SR reached up to 70%, there was a marked increase in MSDBx under all HR conditions since the automatic regulatory mechanism is no longer

Effects of Variations in HR on Myocardial Perfusion
Myocardial perfusion was firstly evaluated with respect to the sensitivity of Q endo /Q epi to variations in HR. From the results presented in Figure 4A, Q endo /Q epi decreased evidently with the increase in HR when HR was higher than 60 bpm; relatively, Q endo /Q epi was insensitive to variations in HR when HR was low (e.g., <60 bpm). In addition, the presence of a coronary artery stenosis of moderate to high severity (i.e., SR = 50~70%) considerably reduced Q endo /Q epi and augmented its sensitivity to variations in HR. These results indicate that tachycardia and coronary artery disease play additive roles in depressing the supply of blood to the subendocardium, thereby increasing the risk of subendocardial ischemia. Figure 4B shows the relationships between MSDBx and HR under various coronary artery stenotic conditions. MSDBx exhibited an overall trend of rising with the increase in HR, with the trend being most remarkable when HR was lower than 100 bpm. Recalling the results presented in Figure 3 and the definition of MSDBx (i.e., Equation (7)), the increase of MSDBx with HR is caused mainly by the reduction in the amount of coronary artery blood flow supplied to the myocardium in each cardiac cycle following the increase in HR. Relatively, varying the severity (i.e., SR) of coronary artery stenosis had little influence on MSDBx when the stenosis severity was at the mild to moderate level, because under such conditions the stenosis-caused reduction in coronary artery flow could be well compensated for by the coronary automatic regulatory mechanism [24]. When SR reached up to 70%, there was a marked increase in MSDBx under all HR conditions since the automatic regulatory mechanism is no longer sufficient to preserve coronary blood flow. With respect to myocardial ischemia, increasing HR to a high level could lead to myocardial ischemia (i.e., MSDBx > 0) even in the absence of coronary artery stenosis; on the other hand, reducing HR could convert the myocardial ischemic state (i.e., MSDBx > 0) into a non-ischemic state (i.e., MSDBx < 0), even in the presence of severe coronary artery stenosis. Figure 5 further shows the relationships between MSDBx and HR in the subendocardium and subepicardium, respectively. While MSDBxes in the subendocardium and subepicardium both showed a trend of increasing with HR, the sensitivity of MSDBx to HR was much higher in the subendocardium than in the subepicardium. In addition, the effects of coronary artery disease on MSDBx were more evident in the subendocardium than in the subepicardium when the severity of stenosis was at a moderate to high level (SR ≥ 50%), especially under the conditions of high HR (>100 bpm). These results imply that the perfusion state of the subendocardium is much more sensitive than that of the subepicardium to variations in HR, and that tachycardia is more likely to induce ischemia in the subendocardium even in the absence of coronary artery disease or in the presence of moderate coronary artery disease.
Mathematics 2021, 9, x FOR PEER REVIEW 8 of 16 sufficient to preserve coronary blood flow. With respect to myocardial ischemia, increasing HR to a high level could lead to myocardial ischemia (i.e., MSDBx > 0) even in the absence of coronary artery stenosis; on the other hand, reducing HR could convert the myocardial ischemic state (i.e., MSDBx > 0) into a non-ischemic state (i.e., MSDBx < 0), even in the presence of severe coronary artery stenosis. Figure 5 further shows the relationships between MSDBx and HR in the subendocardium and subepicardium, respectively. While MSDBxes in the subendocardium and subepicardium both showed a trend of increasing with HR, the sensitivity of MSDBx to HR was much higher in the subendocardium than in the subepicardium. In addition, the effects of coronary artery disease on MSDBx were more evident in the subendocardium than in the subepicardium when the severity of stenosis was at a moderate to high level (SR ≥ 50%), especially under the conditions of high HR (> 100 bpm). These results imply that the perfusion state of the subendocardium is much more sensitive than that of the subepicardium to variations in HR, and that tachycardia is more likely to induce ischemia in the subendocardium even in the absence of coronary artery disease or in the presence of moderate coronary artery disease.

Effects of Irregular HRs on Coronary Artery Flow and Myocardial Perfusion
Each set of 100 irregular R-R intervals generated by assigning different values to σ in Equation (6) (shown in panel F of Figure 2) was imposed to the model to simulate the

Effects of Irregular HRs on Coronary Artery Flow and Myocardial Perfusion
Each set of 100 irregular R-R intervals generated by assigning different values to σ in Equation (6) (shown in panel F of Figure 2) was imposed to the model to simulate the changes in coronary hemodynamics and myocardial perfusion. For the coronary artery (i.e., LAD b1 ), in addition to the normal state (i.e., no stenosis), moderate (SR = 50%) and severe (SR = 70%) stenotic diseases were considered. The cardiac cycle-specific mean coronary blood flow rates and MSDBxes are presented in the form of scatter plots in panels (A)-(F) of Figure 6 and in the form of statistical box plots in panels (G) and (H) of Figure 6. From the scatter plots, it is evident that the stochastic distributions of both mean coronary blood flow rate and MSDBx were enhanced following the increase in σ, which corresponds to the enlargement of the error bars in the box plots. Increasing the severity of coronary artery stenosis led to an overall decrease in blood flow rate, while at the same time an increase in MSDBx. In the meantime, the ranges of stochastic distributions of mean coronary blood flow rate and MSDBx were significantly narrowed when the severity of coronary artery stenosis was high (i.e., SR = 70%). The cardiac cycle-specific mean coronary blood flow rates and MSDBxes were further analyzed over each set of 100 consecutive cardiac cycles to calculate the total amount of blood flow (Q T100 ) and the mean value of MSDBxes larger than zero, with the results being presented in Figure 7 in the form of histograms. It was observed that Q T100 , as expected, decreased significantly with the increase in the severity of coronary artery stenosis, but only slightly decreased following the increase in σ. In contrast, the mean value of MSDBxes > 0 increased considerably with the increase in σ, with the degree of increase being largest in the absence of coronary artery stenosis, while smallest in the presence of severe coronary artery stenosis. We further performed a stratified data analysis by dividing positive MSDBxes into three groups, i.e., 0 < MSDBx ≤ 0.25, 0.25 < MSDBx ≤ 0.5, and MSDBx > 0.5, which were herein defined to represent mild myocardial ischemia, moderate myocardial ischemia and severe myocardial ischemia, respectively. The grouping of MSDBx allowed us to evaluate the overall risk of myocardial ischemia of certain severity by calculating the percentage of cardiac cycles, with the values of MSDBxes falling in each group among the 100 cardiac cycles. From the results shown in Figure 8, when the coronary artery was free of stenosis, although the total proportion of cardiac cycles with MSDBx > 0 was almost not affected by σ, the proportions of cardiac cycles with 0.25 < MSDBx ≤ 0.5 and MSDBx > 0.5 increased considerably with the increase in σ, indicating that severe irregularity of HRs may result in the occurrence of moderate to severe transient myocardial ischemia even in the absence of coronary artery disease. When there was a moderate or severe coronary artery stenosis (SR = 50% or 70%), the proportion of cardiac cycles with MSDBx > 0.5 increased significantly, although the total proportion of cardiac cycles with MSDBx > 0 and the proportion of cardiac cycles with 0.25 < MSDBx ≤ 0.5 both decreased.
In order to uncover the mechanisms underlying the changes in MSDBx, we further statistically analyzed the cardiac cycle-specific values of Q T (i.e., total amount of coronary artery blood flow in each individual cardiac cycle) and PVA in the groups of cardiac cycles with mild myocardial ischemia (0 < MSDBx ≤ 0.25), moderate myocardial ischemia (0.25 < MSDBx ≤ 0.5) and severe myocardial ischemia (MSDBx > 0.5), respectively. From the results presented in Figure 9, under all conditions characterized by different values of σ or severities of coronary artery stenosis, the mean value of Q T exhibited more evident inter-group differences than that of PVA, with the group of cardiac cycles with severe myocardial ischemia having the smallest value, while the group with mild myocardial ischemia had the largest value. In addition, the inter-group differences in Q T augmented considerably following the increase in σ. From the definition of MSDBx and the results presented in Figure 9, it is easy to conclude that the changes in Q T , rather than the changes in PVA, upon irregular variations in R-R interval (or HR) dominate the changes in MSDBx; in other words, the marked decreases in Q T are a major cause of the increased values of MSDBx in the groups with severe myocardial ischemia. to severe transient myocardial ischemia even in the absence of coronary artery disease. When there was a moderate or severe coronary artery stenosis (SR = 50% or 70%), the proportion of cardiac cycles with MSDBx > 0.5 increased significantly, although the total proportion of cardiac cycles with MSDBx > 0 and the proportion of cardiac cycles with 0.25 < MSDBx ≤ 0.5 both decreased.  R-R interval irregularity when the coronary artery is normal (no stenosis) or has a 50% or 70% stenosis. The statistical box plots of mean flow rates and MSDBxes computed for each set of 100 cardiac cycles are presented in panels (G,H), respectively.   R-R interval irregularity when the coronary artery is normal (no stenosis) or has a 50% or 70% stenosis. The statistical box plots of mean flow rates and MSDBxes computed for each set of 100 cardiac cycles are presented in panels (G,H), respectively.   vere myocardial ischemia having the smallest value, while the group with mild myocardial ischemia had the largest value. In addition, the inter-group differences in QT augmented considerably following the increase in σ. From the definition of MSDBx and the results presented in Figure 9, it is easy to conclude that the changes in QT, rather than the changes in PVA, upon irregular variations in R-R interval (or HR) dominate the changes in MSDBx; in other words, the marked decreases in QT are a major cause of the increased values of MSDBx in the groups with severe myocardial ischemia.

Discussion
In the present study, computational model-based numerical experiments were carried out to investigate the influences on myocardial perfusion of various types of arrhythmia combined with coronary artery disease. In particular, we introduced two indices (i.e., Q endo /Q epi and MSDBx) in addition to generally considered coronary hemodynamic parameters for the purpose of quantitatively evaluating myocardial perfusion and exploring the mechanisms underlying the occurrence of myocardial ischemia. The main findings are as follows: (1) there was an overall decrease in Q endo /Q epi , while an increase in MSDBx, following the increase in HR; (2) MSDBx in the subendocardium was more sensitive than that in the subepicardium to variations in HR; (3) the presence of coronary artery disease of moderate severity could considerably augment the sensitivities of Q endo /Q epi and MSDBx to variations in HR, especially under tachycardic conditions (HR > 100 bpm); (4) severe irregularity of HRs (or R-R intervals) alone could induce the occurrence of transient myocardial ischemia, and could increase the risk of severe myocardial ischemia when combined with coronary artery disease; and (5) the amount of coronary blood flow in each cardiac cycle was a major determinant factor for the condition of myocardial perfusion under various conditions of arrhythmia.
The first and second findings indicate that ischemia is more likely to occur in the subendocardium under tachycardic conditions due to the redistribution of coronary blood flow from the subendocardium toward the subepicardium. Such a trend will be further enhanced in the presence of moderate coronary artery disease according to the third finding. Although moderate coronary artery disease is generally considered to have a relatively low risk of inducing severe myocardial ischemia [25], our findings suggest that the combination of moderate coronary artery disease with tachycardia can significantly increase the risk and severity of ischemia in the subendocardium. Therefore, for patients with moderate coronary artery disease, the level of resting HR should be taken as an important physiological parameter when assessing the condition of myocardial perfusion, especially for the subendocardium. The fourth finding suggests that severe irregularity of HRs can either be an independent risk factor or be an additive factor to coronary artery disease for myocardial ischemia. With regard to the mechanisms underlying the influences of arrhythmia on myocardial perfusion, our numerical results reveal that the sharp decrease in the amount of coronary artery blood flow in each cardiac cycle with the increase in HR is a main cause of myocardial ischemia, although both coronary artery blood flow and PVA of the left ventricle decrease following the increase in HR.
In the treatment of arrhythmia, β-blockers, which have a role of reducing HR and blood pressure through blocking the access of catecholamines to their receptors, are commonly used [26]. However, whether lowering HR with β-blockers can effectively improve the prognosis of myocardial ischemia remains unclear. In the guidelines of the ACC/AHA, the use of β-blockers is recommended in the management of patients with stable ischemic heart disease or patients with normal left ventricular function after myocardial infarction or acute coronary syndrome [27]; however, this is not recommended in the ESC guidelines due to the lack of randomized trials on the use of β-blockers in patients with stable coronary artery disease [28]. At this point, the use of β-blockers may be supported by the findings of our study regarding the beneficial role of lowering HR in improving myocardial perfusion (especially in the subendocardium) in the presence of moderate to severe coronary artery disease. In addition, clinical studies have shown that arrhythmia is significantly associated with the severity and progression of coronary artery disease [29,30]. Our study demonstrates that tachycardia or severe HR irregularity can considerably increase the risk of severe myocardial ischemia in the presence of moderate coronary artery stenosis, which implies that arrhythmia might be involved in the pathological changes of myocardium through its influence on myocardial perfusion. Given the fact that coexistence of coronary artery disease and arrhythmia is not uncommon in the population [31], arrhythmia should be well considered in order to identify patients with clinically insignificant coronary artery stenoses while severe myocardial ischemia.
Although the findings of our study provide theoretical insights for understanding the characteristics of myocardial perfusion under various arrhythmic conditions combined with coronary artery disease, they must be considered in the context of certain limitations. A major limitation arises from the disregard of the regulatory responses of systemic hemodynamics to HR variations. For example, lowering HR can not only alter hemodynamic conditions in the coronary circulation, but can also induce decreases in arterial blood pressure [10] and cerebral blood flow [4]. Under in vivo conditions, a decrease in arterial blood pressure will alter the reflex responses of the baroreceptors located at the carotid sinus and aortic arch, which will further trigger the regulatory responses of the autonomic nervous system aimed to recover blood pressure via modulating cardiac contractility and the tones of small resistance arteries, and such systemic regulatory responses will act along with the local regulatory responses (e.g., vasodilation or remodeling) of cerebral vessels to preserve cerebral blood flow [32]. These regulatory mechanisms play an important role in maintaining the homeostasis of hemodynamics, but were not incorporated in our model. In particular, arterial blood pressure has been demonstrated in previous studies [10,33] to have moderate influence on the distribution of blood flow in the myocardium; therefore, the regulation of arterial blood pressure may be a factor worthy of consideration when HR varies slowly or is stabilized at a certain level for long time. In this sense, our study only revealed the transient influences of HR variations on myocardial perfusion. In addition, the model parameters have been assigned based on population-averaged clinical data, which renders the model generic and unable to represent the characteristics of coronary hemodynamics and myocardial perfusion in a specific patient. Accordingly, the presented findings would contribute as a theoretical reference for interpreting population-averaged rather than patient-specific clinical observations. In future studies, the model may be calibrated to the clinical data of individual patients [34,35] to predict the risk or severity of myocardial ischemia upon the onset of arrhythmia, or be utilized as a mathematical tool to generate in-silico data representing a large number of virtual patients so that statistical methods or machine learning techniques could be applied to identify specific pathophysiological conditions under which the arrhythmia-induced risk of myocardial ischemia is high.

Conclusions
Model-based numerical experiments have been performed to systematically study the influences of different types of arrhythmia combined with coronary artery disease on coronary hemodynamics and myocardial perfusion. The obtained results showed that tachycardia and severe irregularity of HRs could induce myocardial ischemia even in the absence of coronary artery disease. In contrast, reducing HR was found to have a beneficial role in improving myocardial perfusion, especially in the subendocardium. In addition, the influences of arrhythmia on myocardial perfusion were augmented in the presence of coronary artery disease. These findings suggest that arrhythmia is a factor worthy of note in the diagnosis and management of patients with coronary artery disease.

Data Availability Statement:
The data presented in this study are available from the corresponding author on reasonable request. For readers who are developing similar models, technical support from the corresponding author are also available.

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

Abbreviations
MSDBx, myocardial supply-demand balance index; HR, heart rate; bpm, beats per minute; LM, main trunk of the left coronary artery; LAD, left anterior descending artery; LCx, left circumflex coronary artery; RCA, right coronary artery; Q endo /Q epi , subendocardial-to-subepicardial flow ratio, PVA, pressure-volume area; SR, stenosis rate; ACC, American College of Cardiology; AHA, American Heart Association; ESC, European Society of Cardiology.