Mechanism underlying impaired cardiac pacemaking rhythm during ischemia: A simulation study.

Ischemia in the heart impairs function of the cardiac pacemaker, the sinoatrial node (SAN). However, the ionic mechanisms underlying the ischemia-induced dysfunction of the SAN remain elusive. In order to investigate the ionic mechanisms by which ischemia causes SAN dysfunction, action potential models of rabbit SAN and atrial cells were modified to incorporate extant experimental data of ischemia-induced changes to membrane ion channels and intracellular ion homeostasis. The cell models were incorporated into an anatomically detailed 2D model of the intact SAN-atrium. Using the multi-scale models, the functional impact of ischemia-induced electrical alterations on cardiac pacemaking action potentials (APs) and their conduction was investigated. The effects of vagal tone activity on the regulation of cardiac pacemaker activity in control and ischemic conditions were also investigated. The simulation results showed that at the cellular level ischemia slowed the SAN pacemaking rate, which was mainly attributable to the altered Na+-Ca2+ exchange current and the ATP-sensitive potassium current. In the 2D SAN-atrium tissue model, ischemia slowed down both the pacemaking rate and the conduction velocity of APs into the surrounding atrial tissue. Simulated vagal nerve activity, including the actions of acetylcholine in the model, amplified the effects of ischemia, leading to possible SAN arrest and/or conduction exit block, which are major features of the sick sinus syndrome. In conclusion, this study provides novel insights into understanding the mechanisms by which ischemia alters SAN function, identifying specific conductances as contributors to bradycardia and conduction block.

Ischemia in the heart impairs function of the cardiac pacemaker, the sinoatrial node (SAN). However, the ionic mechanisms underlying the ischemia-induced dysfunction of the SAN remain elusive. In order to investigate the ionic mechanisms by which ischemia causes SAN dysfunction, action potential models of rabbit SAN and atrial cells were modified to incorporate extant experimental data of ischemia-induced changes to membrane ion channels and intracellular ion homeostasis. The cell models were incorporated into an anatomically detailed 2D model of the intact SAN-atrium. Using the multi-scale models, the functional impact of ischemia-induced electrical alterations on cardiac pacemaking action potentials (APs) and their conduction was investigated. The effects of vagal tone activity on the regulation of cardiac pacemaker activity in control and ischemic conditions were also investigated. The simulation results showed that at the cellular level ischemia slowed the SAN pacemaking rate, which was mainly attributable to the altered Na þ -Ca 2þ exchange current and the ATP-sensitive potassium current. In the 2D SAN-atrium tissue model, ischemia slowed down both the pacemaking rate and the conduction velocity of APs into the surrounding atrial tissue. Simulated vagal nerve activity, including the actions of acetylcholine in the model, amplified the effects of ischemia, leading to possible SAN arrest and/or conduction exit block, which are major features of the sick sinus syndrome. In conclusion, this study provides novel insights into understanding the mechanisms by which ischemia alters SAN function, identifying specific conductances as contributors to bradycardia and conduction block. Sinus bradycardia, a phenotype of sick sinus syndrome (SSS), may be caused by sinoatrial node (SAN) ischemia. Recent experimental studies have focused on characterizing ischemia-induced changes to membrane ion channels of nodal cells. However, it is still unclear how such changes at the cellular ion channel level may cause sinus bradycardia, compromising the ability of the SAN to pace and drive the surrounding atrium, given that the cardiac pacemaker is an integrative and complex nonlinear system. In this study, we address these issues with multi-scale models of the SAN, from single cell to 2D SAN-atrium tissue levels. The functional impact of ischemia-induced electrophysiological changes on cardiac pacemaking action potentials and their conduction was investigated. Effects of vagal tone activity on the regulation of cardiac pacemaker activity in control and ischemic conditions were also studied. The results provide new insights into understanding the mechanism by which ischemia causes cardiac pacemaker dysfunction, and why SSS patients have significantly lower heart rates (even sudden death) at night when the vagal tone is high.

I. INTRODUCTION
The sinoatrial node (SAN), the pacemaker of the mammalian heart, possesses intrinsic automaticity. Dysfunction of the SAN may lead to abnormal cardiac rhythms, manifested by intermittent sinus bradycardia, sinus arrest, sinus pause, slow SAN-atrium conduction, sinus exit block or alternating bradycardia and atrial tachycardia. 1,2 SAN dysfunction may arise from a variety of conditions, such as ageing, 3 gene mutations, 4,5 and cardiac ischemia. 6 In the latter case, the sinus node artery, which carries the main blood supply to the SAN, may suffer from hypoperfusion, causing ischemia that impairs SAN function. Experimental studies have shown that occlusion of the sinus node artery produces various dysrhythmias, including sinus slowing and a) Author to whom correspondence should be addressed: henggui.zhang@ manchester.ac.uk sinoatrial block. [7][8][9] With decreasing sinus rhythm, cardiac output and blood pressure decreases, and in some extreme cases, sinus arrest may occur causing sudden death. 10 Although ischemia-induced SAN dysfunction can be fatal, the exact underlying ionic mechanisms remain elusive.
Prior studies have been carried out to characterize the effects of ischemia on the electrophysiological properties of mammalian atrial and ventricular myocytes (e.g., . It has been found that ischemia is associated with acidosis, hyperkalemia and hypoxia, each of which produces alterations to the electrophysiological properties of cardiac cell and tissues, including changes to ion channel properties and intercellular electrical coupling. 12,14 These changes may be responsible for the genesis of cardiac arrhythmias during ischemia. 14,15 However, there is little experimental information on the effects of ischemia on the electrophysiological properties of SAN due to its complicated and heterogeneous anatomy. Experimental evidence suggests that ischemia can also modulate ion channel function in SAN cells, producing an increased L-type calcium channel current (I CaL ) 16 and hyperpolarization activated inward current (I f ), but decreased T-type calcium channel current (I CaT ), rapid and slow components of the rectifier potassium channel currents (I Kr and I Ks ), and Na þ -Ca 2þ exchanger activity (I NaCa ). 6,16,17 However, it is still unclear if the changes identified at the cellular ion channel level are sufficient to account for the ischemia-induced sinus bradycardia, a nonlinear behaviour arising from the interactions of millions of SAN cells. The relative roles of ischemia-induced changes to altered particular ion channel currents in compromising SAN pacemaking activity and the ability of the SAN to drive the surrounding atrial muscle are yet to be elucidated. Acetylcholine (ACh), the neurotransmitter released by the vagal nerves is known to slow down the pacemaking rate (PR). [18][19][20] Measurements from anaesthetized dogs have demonstrated that vagal action on the heart rate can be potentiated during SAN ischemia. 21 However, whether and how the vagal activity compromises the SAN pacemaking activity during ischemia has not yet been elucidated.
Computational models of the heart have provided powerful tools to study the functional effects of gene mutations, 4,5 heart failure, 22,23 and ageing 24,25 on cardiac pacemaker activities. Being constructed from, and validated against experimental data, they provide a means for quantitatively predicting the functional roles of altered molecular dynamics and ionic channels in a systematic way that is difficult to achieve in an experimental setting. 4,26,27 In previous computational work, Butters et al. showed that SCN5A gene mutations are causally linked to the sick sinus syndrome. 4 Modelling has also been used to investigate the mechanisms by which ankyrin-B defects increase susceptibility to SAN dysfunction and atrial fibrillation and the functional impact of heart failure on SAN. 28,29 In this study, we developed multi-scale (cell and tissue) models of the rabbit SAN-atrium in control and ischemic conditions. Using these models we investigated the functional impact of ischemia on the initiation and conduction of SAN action potentials and their conduction into the atrium.

A. Mathematical models of single cell and 2D SAN-atrium tissue
The consequences of ischemia-induced electrophysiological alterations for SAN cells were investigated by using: (1) electrophysiologically detailed central and peripheral SAN cell models; 30 (2) a 2D anatomical model of the intact SAN-atrium tissue, incorporating accurate single cell models of the SAN 30 and the right atrium (RA), 31 together with the reconstructed tissue geometry. 32 These single cell models have been validated in prior studies. 30,31 In numerical simulations, the time step was chosen to be 0.01 ms at the cellular level. In the 2D intact SAN-atrium model, we used the same time step with a space step of 0.04 mm which gave accurate numerical solutions.

B. 2D tissue slice model
The 2D anatomical model of the intact SAN-atrium tissue used in this study was reconstructed from histological and immunohistochemistry imaging data. 32 In brief, the model considered distinct regions of the SAN-atrium, with consideration of the intrinsic electrical heterogeneity in the atria and central and peripheral SAN cells. 32,33 The anatomical geometry of the model was based on the reconstructed endocardial surface of the 3D SAN-atrium tissue 32 (see Fig. S1 in the supplementary material), which was discretized by a high spatial resolution of 0.04 mm to form a regular Cartesian grid of 385 Â 250 nodes. Each node had a flag variable to identify the cell type (the SAN central, peripheral and RA cell), based on immunohistochemistry mapping data. 32 The central and peripheral SAN cells were modelled by the Zhang et al. equations, 30 and atrial cells of the RA were modelled by the Aslanidi et al. equations. 31 Similar to our prior studies, 4,34 in the 2D intact SANatrial tissue model we incorporated local variations of cellular electrophysiology as observed experimentally; these included regional differences of cell membrane capacitance (C m ), the diffusion coefficient (D; as a determinant of conduction velocity) and cellular electrophysiological properties 32 (see Fig. S2 in the supplementary material). Spatial distributions of current densities were correlated with the cell membrane capacitance, C m , which increases gradually from the centre to the periphery of the SAN. 30,34 The gradient distribution of D was introduced according to Zhang et al. 34 In order to model the effects of ACh on SAN cells, the ACh-activated K þ current, I K,ACh was incorporated, I CaL was partially inhibited, and the activation curve of I f was shifted toward more negative potentials. Detailed equations and parameters used in the model have been given in the work of Butters et al. 4 The 2D anatomical model of the intact SAN-atrium was validated by its ability to reproduce the correct sequence of action potential (AP) initiation and conduction through the rabbit intact SANatrium, as found experimentally 3,35 (see Fig. S1in the supplementary material).
C. Ischemic SAN model

Changes to ion channel currents
In order to model effects of ischemia on the SAN, the conductances of some of the main ionic channel currents responsible for cellular depolarization and repolarization were adjusted, based on experimental data on the effects of ischemia on the rabbit SAN (see Table I). 6,16 These changes include an increased L-type calcium channel conductance (g CaL ; increased by 110%), 16 a decreased T-type calcium channel conductance (g CaT ; decreased by 36.8%), 6 decreased rapid-and slowdelayed rectifier potassium channel conductance (g Kr and g Ks ; decreased by 16% and 26% respectively) 6 and decreased Na þ / Ca 2þ exchanger (g NaCa ; decreased by 42%) 6 current. It was shown in a previous experimental study that in "ischemic" Tyrode's solution the funny current (I f ) of SAN cells increased by 15.3% at À60 mV. 16 As sodium channel current (I Na ) is present in the peripheral SAN cells and does contribute to the pacemaking and driving ability of the SAN, 30,36 an ischemiainduced reduction in the channel conductance (g Na ; decreased by 25%) 14,37 as observed in ventricular cells was also considered in simulations (see Table I, the right column).

Hyperkalemia
For simulating hyperkalemia during ischemia, the extracellular potassium concentration ([K þ ] o ) was elevated by 50% as observed experimentally. 38,39 3. ATP-sensitive channel current (I KATP ) due to hypoxia Previous experimental studies suggest that the metabolic depletion-activated K-ATP channel current (I KATP ) may play an important role in the pathogenesis of SAN dysfunction during cardiac ischemia. 40,41 It has been shown that I KATP channel openers regulate the maximum diastolic potential (MDP) and pacemaking rates of rabbit SAN cells. 42,43 Similar functional effects of I KATP on regulation of cardiac pacemaking APs have also been observed in mice during hypoxia. 41,44 Therefore, in simulations we incorporated the I KATP formulation developed by Ferrero et al. 45 Details of the I KATP formulation can be found in the supplementary material.

Changes of ionic homeostasis and energy metabolism
The homeostasis of H þ in rabbit SAN cells is maintained by the Na þ /H þ exchanger, 46 which may lead to an increased intracellular Na þ concentration ([Na þ ] i ) during ischemia. In addition, some other studies from ventricular cells have found that the intracellular free Na þ concentration is increased due to the hypoxia-induced depletion of the cellular ATP content. 47,48 However, these experimental data are absent for rabbit SAN cells. In order to consider all possible actions of ischemia in simulations, some data from ventricular ischemia experiments were also considered, including changes in intracellular ATP (ATP i ), the extracellular ADP (ADP i ), and the intracellular Na þ concentration. For details, please see Table I.

A. SAN effects of ischemia at the single cell level
The effects of ischemia-induced changes on cellular ion channels and ion homeostasis (see Table I) on spontaneous action potentials of SAN cells were first investigated. The results are shown in Fig. 1, in which simulated action TABLE I. Ischemia-induced electrophysiological changes in rabbit SAN cells. All the experimental data on ischemia-induced ion channel remodeling were from rabbit SAN cells with 5 or 7 min ischemia. Other experimental data were based on rabbit ventricular cells with mild ischemia condition.

Changed factors
Percentage of change (normalized to control condition) Changed factors Percentage of change (normalized to control condition) potentials in the ischemic condition are superimposed on those in the control condition [ Fig. 1(a)]. Ischemia increased the overshoot of the APs, which corresponded to the experimental observation of Du and Nathan; 6 prolonged the action potential duration, and hyperpolarized the MDP in the central cell model, but depolarized the MDP in the peripheral cell model. In both cell models, the phase-4 depolarisation process was slowed down, leading an increased pacemaking cycle length (PCL; the time interval between two successive pacemaking APs); consequently there was a deceleration of the pacemaking rate. Ischemia increased PCL by about 41.2% in the central cell model, which was close to the experimental data of Du and Nathan 6 (PCL increased by 43%). A 23.5% increase of PCL in the peripheral cell model was observed, which was somewhat smaller than experimental data 6 [as shown in Fig. 1(b)]. Therefore, our simulation results show that the ischemia-induced changes in ion channel currents, ionic concentrations, and the activation of I KATP are sufficient to account for the bradycardia observed in ischemic rabbit SAN cells. The more depolarised MDP in the peripheral cell model with ischemia matched the experimental observation of Du and Nathan. 6 The modest hyperpolarization of MDP in the central cell model with ischemia differed from the experimental data of Du and Nathan. 6 In the study of Du and Nathan, 6 SAN cells were isolated from a rectangular piece of nodal tissue with a size of (1.0-1.5) Â (3.0-3.5) mm located at about 0.5-1.0 mm from the crista terminalis (CT), with cell capacitance ranging from 20 pF to 62 pF; this wide range of cell sizes suggests that the isolated cell population may have contained both central (small capacitance) and peripheral (large capacitance) cells. In simulations, our cell models considered the intrinsic regional heterogeneity of the SAN, 31 which assumed a cell capacitance of 20 pF for the central cell and 65 pF for the peripheral cell models. This may account for the discrepancy between the simulation and the experimental data of Du and Nathan for the modest hyperpolarized membrane potential. The roles of each individual ischemia-induced cellular change in alterations to SAN cell APs were investigated by using either an "inclusive" method or an "exclusive" method. With the inclusive method, each individual ischemic action was considered alone in simulations. With the exclusive method, the individual ischemic action of interest was excluded from the integral action of ischemia. Figure 2 shows the results obtained from the inclusive method. It illustrates that in both the central and peripheral cell models the ischemia-induced increase of I CaL played a primary role in prolonging the action potential duration (APD) (by 20.  21,25 Our results suggested that the altered I CaT , I Kr , I Ks , and I NaCa had almost no effect on PCL in the peripheral cell model [ Fig. 2(bii)], while the ischemia-induced decrease of I Na , activation of I KATP , increases of [Na þ ] i and I CaL caused an increase of PCL by 3.9%, 1.4%, 7.2%, and 17.2%, respectively. Note that in the peripheral cell model the increase of I CaL contributed the most to the PCL increase, which resulted from an increase of the APD as shown in Fig. 2(aii). Simulations from the exclusive method showed similar results with respect to the influence of each individual ischemic action on SAN APs to the inclusive method, as shown in Figs. S3, S4, and S5 in the supplementary material.
Further simulations were performed to elucidate the mechanisms by which the [Na þ ] i increase ([Na þ ] i ¼ 13.28 mM) during ischemia slowed down pacemaking. The results are shown in Fig. 3 for central (left panels) and peripheral (right panels) cells in control and [Na þ ] i . In both cell models, it was found that an increased [Na þ ] i resulted in an increased PCL and more hyperpolarized MDP  Fig.  4(d)]. In the ischemic condition, application of ACh increased the PCL of the central and peripheral cell models by 72.9% and 30.2% respectively compared to control; these changes were markedly greater than those under ischemic condition (i.e., PCLs were increased by 41.2% and 23.5% respectively as compared to control condition). As such, the negative chronotropic effects of ACh were amplified by ischemia in both the centre and the peripheral cell models.
The dose-dependent effects of ACh on APs in control and ischemic conditions were also considered. Figure 5 presents the results obtained from the central and peripheral cell models for various "physiological" [ACh]. [18][19][20] It shows that an increase of ACh concentration resulted in an increase in the PCL of both cell types. The negative chronotropic effect of all ACh concentrations was amplified by ischemia. At an ACh concentration of 8.0 Â 10 À8 mol/l, central SAN cells still exhibited pacemaking (though with a prolonged PCL) in the control condition, but became quiescent in the ischemic condition [ Fig. 5(bi)]. Similar results were seen in the ischemic peripheral cell at an ACh concentration of 4.0 Â 10 À7 mol/l [ Fig. 5(cii)]. Figure 5(d) summarizes the simulated dose-dependent effect of ACh on PCLs in control and ischemic conditions. It shows that ischemia shifted the dose-dependence of the PCL leftwards, indicating a more suppressive effect of ACh on SAN cells during ischemia.

C. Two-dimensional tissue effects
Due to the electrotonic interaction between the SAN and atrium, the functional impact of ischemia seen at the isolated single-cell level may be different to that seen in the intact SAN-atrium tissue level. Thus, in addition to single cell simulations, multicellular tissue models must be used to evaluate the functional consequences of ischemia and ACh in the setting of ischemia. Therefore, simulations of activity from intact SAN-atrium tissue under control and ischemia conditions were conducted, the results of which are shown in Fig. 6. Figure 6(ai) shows the snapshots of the initiation and conduction pattern of pacemaking APs in the 2D SANatrium model under control condition at different timings. The AP was first initiated in the centre of the SAN. Once initiated, the AP propagated preferentially from the centre towards the periphery of the SAN, and then to the crista terminalis (CT) before entering the atrium in the direction towards the right atrium. However, in the direction towards the atrial septum, the AP conduction was blocked in the conduction block zone, which was encircled by excitation waves from the superior and inferior tissues surrounding the zone. The simulated activation and AP conduction pattern in the 2D model matches prior experimental data from the rabbit heart. 31 Figure 6(aii) plots the spatial (running vertically) and temporal (running horizontally) profiles of APs [recorded from representative cells along the black line across the SAN-atrium model as shown in Fig. 6(ai)] for three consecutive cycles, demonstrating a stable excitation and conduction pattern in the intact SAN-atrium model. The measured PCL was about 367 ms, which is close to experimental data obtained from multicellular tissue recordings of the intact SAN-atrium of the rabbit heart. 31 Note that the measured PCL at the tissue level was greater than that of isolated central SAN cells, which controlled the pacemaking rhythm of the SAN. This is due to the electrotonic interaction between the SAN and the atrium, 31 which depressed the pacemaking APs of SAN cells as observed experimentally in previous studies. 31,51 SAN ischemia [ Fig. 6(bii)] or application of ACh (1.0 Â 10 À6 mol/l) [ Fig. 6(cii)] in the 2D model did not alter AP conduction patterns (the relative snapshots of the initiation and conduction pattern of the pacemaking action potentials in the 2D tissue shown in Figs. 6(bi) and 6(ci) but shifted the leading pacemaking sites towards the periphery by 0.12 mm in both conditions. An ACh-induced shift in the leading pacemaker site has been observed experimentally by optical mapping. 52 In simulations, SAN ischemia and ACh also slowed down pacemaking rate, resulting in a PCL of 396 and 560 ms, respectively. It is of interest that AP conduction was also slowed, as demonstrated by the slower spread of the excitation wavefront in the tissue, as compared to control tissue shown in Fig. 6(ai). In the normal tissue, the excitation wavefront spread over the 2D tissue model by about 80 ms, which increased to 90 and 120 ms in SAN ischemia and ACh conditions, respectively.
Combined actions of ACh (1.0 Â 10 À6 mol/l) and ischemia amplified the negative chronotropic effects of each intervention. The results are shown in Fig. 6(di) for the activation pattern and 6(dii) for the space-time plots of APs from the representative cells in the tissue. With an integrated action of both ACh and SAN ischemia, SAN pacemaking activity was further slowed, with a measured PCL of 593 ms. It is also notable that the original leading pacemaking site was shifted toward the periphery SAN, and a region in the peripheral SAN tissue close to the superior vena cava became a secondary leading pacemaking site. The APs generated by either site were unable to drive the surrounding atrial muscle, leading to conduction exit block-a phenomenon typically seen in the sick sinus syndrome. With an increased ACh concentration during ischemia, SAN exit block became more likely to occur than in the control condition with the same level of ACh (data are not shown).
Detailed analyses were also performed to investigate the mechanisms by which ischemia affects APs and conduction in the intact SAN-atrial tissue model and PCL. The results are shown in Fig. 7. Figure 7(a) plots the measured activation time for the representative cells across the tissue in control, ischemia and in ACh with and without ischemia; the corresponding conduction velocity of APs across the tissue was calculated and shown in Fig. 7(b). Both ischemia and ACh increased the activation time, but the latter had greater effects, indicating an increased time required for the AP to

093934-6
Bai et al. Chaos 27, 093934 (2017) conduct from SAN to the atrial muscle. In Fig. 7(a), the sinoatrial-atrium activation time was reduced in the direction towards the atrium septum in the ischemia and ACh condition [blue line in Fig. 7(a)]. This was due to the presence of two leading pacemaking sites in that condition, one of which was at the peripheral region of the SAN that was close to the right atrium (i.e., the leading pacemaker site shifted, see Fig.  6(di), snapshot at 40 ms), therefore, the time taken to reach the right atrium was decreased, though the sinoatrial-atrium conduction was actually slowed down. The increased activation time in the direction towards the crista terminalis corresponded to a decreased AP conduction velocity as shown in Fig. 7(b). The increased activation time and the corresponding decreased conduction velocity were attributable to a decreased upstroke velocity of cellular APs as shown in Fig.  7(c). Ischemia reduced dV/dt max in both right atrium and SAN, with a greater effect on the right atrium. Application of ACh further reduced dV/dt max in ischemic tissue, most severely in the SAN peripheral region. A decreased dV/dt max in both SAN peripheral regions indicated a reduced driving power of the SAN and an increased electrical load of the atrium to the SAN, resulting in an increased PCL and even a SAN conduction exit block, in which case the SAN failed to drive the atrial muscle.

IV. DISCUSSION
In this study, we modeled the electrophysiological activities of central and peripherial SAN cells during ischemia based on experimental data, and incorporated these models into an electrophysiologically and anatomically detailed 2D intact rabbit SAN-atrium model. We investigated the functional role of ischemia-induced alterations at the cellular level to elucidate which and how these may cause SAN bradycardia during ischemia, and whether or not these alterations play the same functional role in both central and peripheral SAN cells. With the tissue model, we further investigated the functional effects of SAN ischemia in impairing the ability of the SAN to pace and drive the surrounding atrial muscle. Our simulation results provide evidence substantiating the causative link between SAN ischemia and sinus bradycardia. The major findings of this study are as follows: (1) at the cellular level, the ischemia-induced bradycardia in both central and peripheral cells can be attributed to the intracellular Na þ increase, I KATP activation, increased I CaL , and the reverse mode of I NaCa ; (2) at the intact SAN-atrium tissue level, ischemia decreased the pacemaking rate and slowed down the AP conduction across the SAN-atrium; (3) the bradycardiac effects of SAN ischemia are likely to be amplified by vagal nerve activity: simulated addition of ACh to the tissue with ischemia not only slowed pacemaking and AP conduction, but also compromised the ability of the SAN to pace and drive the atrium, leading to a higher probability of conduction block at the SAN-atrium junction (a phenomenon known as SAN exit block); this effect was greater than that seen with ACh alone. This finding may provide further insight into the mechanism underlying high risk of cardiac arrest in ischemic patients at night, when the vagal tone is greater. 47 Previous experiments on rabbit SAN cells have demonstrated that the MDP changed and pacemaking rate slowed with the addition of ATP-sensitive K þ channel openers. 42,43 Our simulated effect of I KATP is concordant with these experimental findings. Experimental data have shown that Na þ -K þ pump activity was increasingly inhibited during ventricular ischemia leading to intracellular Na þ overload. 54 We assumed that a similar intracellular [Na þ ] i overload may also occur during SAN ischemia. The combined effect of the augmented reverse-mode of I NaCa and/or I NaK , and intracellular Na þ overload further slowed down pacemaking significantly. However, in the peripheral cell model the causative link between SAN ischemia and bradycardia is not exactly the same as in the central cell model. In the peripheral cell model, the increased I CaL and intracellular [Na þ ] i are the major causes of bradycardia, while I KATP activation has little effect. The rise in intracellular [Na þ ] i resulted in functional changes in ion channel currents (Fig. 3) similar to those in the central cell. Altogether, our simulation results suggest that ischemia-induced changes in I NaCa and I CaL , and activation of I KATP play important roles in ischemia-induced bradycardia, which could be potential targets for the treatment of SAN effects of ischemia. Reduction of I CaT had little effect on AP in our simulation results (Fig. 2) as some prior experimental/modelling studies found. [55][56][57] At the intact SAN-atrium tissue level, ischemia slowed down the pacemaking rate, as well as the AP conduction (AP) (Fig. 6). The simulation results of slowed AP conduction through the SAN-atrium during ischemia are in agreement with experimental observations. 58 We simulated the effect of ischemia by considering both ischemia-induced changes to various ionic currents by altering their maximum membrane conductances and altered [K þ ] o and [Na þ ] i . As altered [K þ ] o and [Na þ ] i might produce a secondary effect on the membrane ion channel currents of interest, further simulations were conducted to dissect such a combined process. For this, six different scenarios were considered, and the results are shown in the Fig.  S6 (in the supplementary material)

B. ACh effects on the ischemic SAN
Previous studies have confirmed the negative chronotropic effect of ACh on cardiac pacemaking cells. [18][19][20] However, until now there has been a paucity of information about effects of ACh on ischemic SAN. The present study addresses this issue. With increasing ACh concentrations, a more suppressive effect of ACh on ischemic SAN has been shown at both cellular and 2-D tissue levels. The augmented bradycardiac effect of ACh can be attributed to the reduction of I f (Fig. 4), as well as the activation of I K,ACh .
Our simulations demonstrated synergistic effects of ACh and ischemia: ACh not only further slows down the pacemaking rate in ischemic SAN cells, but also impairs the ability of the SAN to drive the surrounding atrial muscle, as revealed by the results shown in Figs. 6 and 7. It was shown that ACh caused failure of AP conduction from SAN to the atrium in the ischemic condition, which may contribute to the observed conduction failure risks, e.g., atrial fibrillation 59 or sinus arrest 7 (even sudden death) at night when the vagal tone is high.

C. Limitations
The limitations of the rabbit cell and tissue models used in our study have been described previously. 30,31 In order to produce ischemic SAN cell models, we modified the equations for I CaL , I CaT , I Kr , I Ks , I NaCa , and I f based on available experimental data from rabbit SAN cells. 6,16,17 Other changes incorporated into our models came from mammalian ventricular ischemia experiments. 14,37,38,45,49,50 The I KATP formulation here was based on that from a guinea pig ventricular cell model. 45 In cases where data from a specific species/cell type are unavailable, it is accepted practice to use data from other species or cell-types in cardiac modeling. In this regard, it is notable that our simulation results of ischemia-induced bradycardia at the single-cell level showed good agreement with experimental data from ischemic rabbit SAN cells, 6,16 suggesting that the parameters in our ischemic models are reasonable. Although our 2D anatomical tissue model considered the intrinsic electrical heterogeneity of SAN-atrial tissue, further 3-D anatomical models may be required to investigate fully AP initiation in the SAN and conduction from SAN toward atrium (e.g., to establish precise locations of the sinus exit block sites with different concentrations of ACh and/or with different extents of SAN ischemia), which may depend on the details of the tissue spatial structure, 3D heterogeneity and anisotropy. Nevertheless, although the dimensions of our tissue model have been reduced from full 3-D to 2-D, the results of our simulations provide valuable insights into mechanistic links between SAN ischemia and ischemia-induced sinus bradycardia. Moreover, gap junction uncoupling may happen during SAN ischemia as reported in ventricular ischemia, 60 which warrants further study by using bi-domain models to capture gap junction remodeling during ischemic conditions.

SUPPLEMENTARY MATERIAL
See supplementary material for the model equations, parameters and additional results.