Pro-arrhythmic Effects of Hydrogen Sulfide in Healthy and Ischemic Cardiac Tissues: Insight From a Simulation Study

Hydrogen sulfide (H2S), an ambient air pollutant, has been reported to increase cardiac events in patients with cardiovascular diseases, but the underlying mechanisms remain not elucidated. This study investigated the pro-arrhythmic effects of H2S in healthy and ischemic conditions. Experimental data of H2S effects on ionic channels (including the L-type Ca2+ channel and ATP-sensitive K+ channel) were incorporated into a virtual heart model to evaluate their integral action on cardiac arrhythmogenesis. It was shown that H2S depressed cellular excitability, abbreviated action potential duration, and augmented tissue’s transmural dispersion of repolarization, resulting in an increase in tissue susceptibility to initiation and maintenance of reentry. The observed effects of H2S on cardiac excitation are more remarkable in the ischemic condition than in the healthy condition. This study provides mechanistic insights into the pro-arrhythmic effects of air pollution (H2S), especially in the case with extant ischemic conditions.


INTRODUCTION
Pro-arrhythmia risks of ambient air pollution have been long established by epidemiological evidences. Inflammation, autonomic nervous system (ANS) changes, and pulmonary or systemic oxidative stress are currently recognized as potential factors responsible for arrhythmogenesis of air pollution (Link and Dockery, 2010;Martinelli et al., 2013); however, the specific effect of air pollution on cardiac functions remains incompletely elucidated owing to the complicated composition of air pollutants (Schlesinger, 2007;Mills et al., 2009;Fiordelisi et al., 2017). Hydrogen sulfide (H 2 S) is a common air pollutant characterized by an odor of rotten eggs and is largely produced in industrial activities, such as food processing, coke ovens, and petroleum refineries. Although epidemiological research has shown association between ambient H 2 S and cardiovascular disease hospitalization and mortality (Godleski et al., 2000;Bates et al., 2002;Finnbjornsdottir et al., 2015Finnbjornsdottir et al., , 2016, how it affects the functions of the cardiovascular system is unclear. Recent studies at the cellular level suggested that NaHS (H 2 S donor) can significantly shorten the action potential duration (APD) of cardiomyocytes, possibly by affecting the ATP-sensitive K + channel (K-ATP) (Bian et al., 2006;Johansen et al., 2006;Zhang et al., 2007;Abramochkin et al., 2009;Chan et al., 2018) and the L-type Ca 2+ channel (Ca-L) (Sun et al., 2008;Zhang et al., 2012). Experiments on rats have characterized the dosedependent effects of NaHS on ionic currents of the two channels, I CaL and I K,ATP (Sun et al., 2008;Zhong et al., 2010), and suggested that low doses of NaHS (i.e., 100 and 150 µmol/L) abbreviated APD by either an enhanced I K,ATP (Zhong et al., 2010) or a reduced I CaL (Sun et al., 2008). Although it is possible that other channel currents may be involved in their experiments owing to the use of a high concentration of intracellular ATP or channel blockers, it is unclear if the observed increase in I K,ATP or decrease in I CaL is sufficient to account for the abbreviation of the APD by NaHS.
It is unclear either how NaHS affects the conduction of action potential (AP) through the transmural ventricular wall. It is known that regional difference of cellular property exists in transmural ventricle cells, and such regional difference in cellular property may be altered by NaHS by modulating the intrinsic transmural dispersion of APD and the effective refractory period (ERP) of cardiac tissues, leading to increased tissue susceptibility to arrhythmogenesis. In addition, in ischemic condition, cardiac electrophysiology is remodeled (Shaw and Rudy, 1997;Trénor et al., 2007;De Diego et al., 2008). It is unclear how NaHS affects the conduction of electrical excitation waves in ischemic tissues and modulates the transmural ERP dispersion, leading to an increased pro-arrhythmic effect.
This study aimed to use a mathematical model of the heart (virtual heart) to evaluate the functional impact of H 2 S/NaHSinduced changes in I CaL and I K,ATP on the electrical properties of rat ventricular myocardium. Specifically, we modified the Pandit et al. (2001) model of healthy rat ventricular myocytes by incorporating I K,ATP and actions of NaHS on I CaL and I K,ATP . Using the model, we simulated the effect of NaHS on cardiac excitation in both control and ischemic conditions with considerations of ischemia-induced alterations to cellular and tissue properties based on available experimental data. Using cellular models, effects of NaHS on cellular AP, excitability, APD, ERP, and transmural dispersion of ERP ( ERP) were computed. Using one-dimensional (1D) model of a transmural ventricular tissue strand, effects of NaHS on ventricular conduction and susceptibility to arrhythmogenesis were also investigated. Finally, effects of NaHS on the genesis and maintenance of reentrant arrhythmias in a two-dimensional (2D) transmural ventricular sheet model were quantified.

Effects of NaHS on Action Potentials via Changes in Individual Ion Channels
The effects of NaHS by its action on I CaL alone are shown in Figure 1. In simulations, the inhibition action of 100 µmol/L NaHS on I CaL was simulated by Eq. 3, whereas the I K,ATP was simulated based on Eq. 2 with the a ratio of intracellular ATP:ADP being kept normoxic as in Sun et al. (2008) (see section "Materials and Methods"). It was shown that a decreased I CaL by NaHS abbreviated the APD for both endocardial ( Figure 1A) and epicardial ( Figure 1B) cells, with no marked effect on their amplitude or the resting potential. Such simulated changes qualitatively matched with experimental data of Sun et al. (2008), although there are some quantitative discrepancies between simulation and experimental data (Figure 1C), possibly owing to a greater I CaL current density in the Pandit model (around −11 pA/pF) as compared with the experimental data (only −3.21 ± 0.13 pA/pF), which produced a greater APD abbreviation than the experimental data of Sun et al. (2008).
The effects of NaHS on I K,ATP alone are illustrated in Figure 2. In this case, the action of NaHS on I K,ATP was calculated from Eq. 2 (see section "Materials and Methods") with a ratio of intracellular ATP:ADP as 200:4.5 for anoxic myocytes to mimic the low ATP experimental setting. It was shown that the simulated hypoxic condition abbreviated APD in both the endocardial (Figure 2A) and epicardial ( Figure 2B) cells with the ATP:ADP setting as stated above. The simulated APD abbreviation of both cell types was consistent with experimental data of Zhong et al. (2010) (Figure 2Cii), although no specific cell type was identified in their experimental study.

Effects of NaHS on APs
The effects of NaHS on the electrical APs of cardiac cells by integral actions of a reduced I CaL and an enhanced I K,ATP are presented in Figure 3. NaHS shortened the APD for both endocardial and epicardial cells (Figures 3A,B). Note that the results presented in Figures 3A,B are similar to those shown in Figure 1, where the I K,ATP was blocked by a high concentration of the intracellular ATP, suggesting that the observed APD shortening by NaHS was mainly attributable to the reduced I CaL .

Effects of NaHS in Ischemic Myocytes
NaHS may modulate differently the electrical activity of cardiac cells between ischemic and healthy conditions. To test this hypothesis, further simulations were conducted to evaluate the effect of NaHS on APs of ischemic cells, and results are shown in Figures 3C,D. Similar to the healthy condition, NaHS abbreviated APD for both endocardial and epicardial cells; however, it showed a greater APD abbreviation in endocardial than in epicardial cells. Such an observation reflects different electrophysiological responses of endocardial and epicardial cells to the same dose of NaHS. The use of 100 µmol/L NaHS barely affected the resting potential but decreased the AP amplitude (APA) of ischemic cells, which was more marked in the endocardial cell model.

Effects of NaHS on Cellular Refractoriness and Excitability
The effects of NaHS on cellular APD and ERP, and the differences between the endocardial and epicardial cells are shown respectively in Figures 3E-G. NaHS abbreviated APD for both endocardial and epicardial cells in healthy and ischemic conditions; however, it had a different effect on ERP between ischemic and healthy conditions -NaHS decreased the ERP in healthy myocytes but increased it in ischemic myocytes. Such different effects of NaHS on modulating APD and ERP were attributed to the involvement of NaHS-I K,ATP , which suppressed the upstroke and amplitude of the AP, leading to a prolonged ERP in the ischemic condition. A non-uniformly modulated ERP between endocardial and epicardial cells produced an augmented ERP dispersion across the transmural ventricular wall, as illustrated in Figure 3G, with the measured transmural  Zhong et al. (2010) without specifying cell types. In the simulations, I to , I CaL , and I K 1 were blocked in accordance with experimental settings. AP, action potential; APD, action potential duration. gradient of ERP increased by 3.1% and more significantly in the ischemic ventricular wall with ERP increased by 32.5%.
Further analyses were conducted to study the effect of NaHS on cellular excitability, measured reciprocally by the threshold of an external stimulus to evoke an AP. Results are shown in Figures 3I,J for the control and ischemic conditions with a low dose of NaHS at 100 µmol/L. In the healthy condition, NaHS did not affect the measured excitation threshold for either the endocardial or epicardial cell model but increased it in the ischemic condition, suggesting that NaHS reduced cellular excitability in the ischemic condition.
Note that the measured excitation threshold of the endocardial and epicardial cell models was affected by the ischemic condition, which is dependent on the S1-S2 intervals. In the ischemic condition, compared with the healthy condition, the measured threshold was lower at long intervals whereas greater at short intervals. Such a difference in the modulated cellular excitability by the ischemic condition was attributable to the combined effect of a prolonged ERP and a more depolarized resting potential as seen in the ischemic condition. Whereas the former increased the excitation threshold at short intervals, causing a rightward shift of the measured threshold curve (Figures 3I,J), the latter reduced the excitation threshold owing to a closer distance between the depolarized resting potential by ischemia and the threshold potential, which is easier to evoke an AP. However, if the resting potential was further elevated by more severe hyperkalemia and finally exceeded the threshold potential, it would always be harder to evoke an AP (i.e., a bigger excitation  Arrows at the top represent the epicardial or endocardial direction toward which the impulse could be conducted. VWs that are less than 0.01 ms were ignored, as they might be incorrectly calculated because of the time precision in the model. The average width of VW in all cases is compared in panel (E). In addition, the measured conduction velocity in 1D strand is shown in panel (F). threshold) in ischemia than in healthy myocyte regardless of the S1-S2 interval.

Temporal Vulnerability to Unidirectional Conduction Block
The effect of NaHS on tissue's vulnerability to genesis of unidirectional conduction block in response to a premature stimulus was also evaluated. Figures 4A,B presented examples of the simulated unidirectional conduction block in healthy and ischemic tissues. Figure 4A clearly shows that owing to the shorter APD in the epicardial region as compared with the endocardial portion, the S2 stimulus applied to the epicardial region (1.0-1.6 mm from the endocardial portion) produced an anterograde conduction. However, the same premature stimulus applied to the same location produced unidirectional conduction in the retrograde direction in ischemic tissues, owing to a reduced conduction velocity (CV) and the altered effective refractory properties. Unidirectional conduction block can only be induced at a certain time period (termed vulnerable window [VW]), or bidirectional block/propagation will happen. The computed overall distribution of VW across the tissue strand is shown in Figures 4C,D, where the light gray bars and red bars are the VW distribution with and without NaHS. The black arrows in these two panels indicated the direction of the unidirectional conduction evoked by the premature S2 stimulus, showing a change in the conduction direction. In healthy tissues, the direction of conduction changed at around 0.7 mm from the endocardial portion, and it happened in a more forward position around 1.9 mm in the ischemia case ( Figure 4D). The average width of VW in four conditions is plotted in Figure 4E, in which it can be seen that ischemia greatly increased the width of the VW by nearly 16 times from 3.2 to 50.5 ms. Marked width of VW across the whole strand suggested a highly increased tissue susceptibility to unidirectional block in the ischemic condition.
When applying 100 µmol/L of NaHS, the width of VW was increased in both healthy and ischemic conditions, and the phenomenon was more significant in the ischemic tissues. Specifically, NaHS further increased the average width of VW from 50.5 to 68.6 ms in the ischemia condition, whereas it slightly increased the average VW width from 3.2 to 4.6 ms in healthy tissues (see Figure 4). The above observation suggested that NaHS could contribute to even higher risks in the ischemic tissues that are already fragile and susceptible to unidirectional block. The adverse effect of NaHS was observed in all locations throughout the strand notwithstanding the direction of conduction. Figure 4F presents the comparison of CV, where NaHS barely affected the CV in healthy tissues with only 0.1% decrement from 19.84 to 19.82 cm/s. In ischemic tissues, NaHS decreased the CV by 2.1% from 9.16 to 8.96 cm/s.

Spatial Vulnerability of Tissue to Sustain Reentry
The effects of NaHS on the dynamic behavior of reentrant excitation waves were also investigated. Premature stimulus during the VW in 2D tissues led to a unidirectional propagation of the impulse toward the anterograde direction, and this impulse further evolved into a spiral wave, as shown in Figure 5. After initiation of the spiral wave, it required a minimum tissue's substrate to sustain, the critical length of which can be used as an index to quantify tissue susceptibility for maintaining reentry. In the healthy condition, the measured critical length was 7.7 mm (Figures 5Ai-Aviii), which decreased by almost 30% to 5.4 mm in the early phase of ischemic tissues (Figures 5Bi-Bviii). This was attributed to the reduced CV and APD in ischemic tissues, producing a short wavelength and thus facilitating the survival of the paired spiral waves.
NaHS measuring 100 µmol/L decreased the critical length in both healthy tissues and ischemic tissues (Figure 5C), which was also attributable to the shortening of APD and wavelength of excitation.

Dynamics of Spiral Waves
Spiral waves in the 2D model of transmural ventricular tissues were unstable in both healthy and ischemic tissues, leading them to self-termination with a finite life span. As shown in Figure 5D, in the healthy tissues, the life span of spiral wave was about 220 ms, which was prolonged to 232 ms by NaHS.
In the tissue model with early phase of ischemia, the life span of spiral wave was about 151 ms, which changed to about 139 ms in the NaHS condition.

DISCUSSION
Although the common industrial air pollutant H 2 S is hazardous to human cardiovascular system, its underlying mechanisms have not been elucidated yet. In this study, we computationally investigated the functional influence of H 2 S on cardiac electrical excitations in healthy and ischemic conditions on the basis of experimental data of NaHS (a donor of H 2 S). Our major findings are as follows: (i) at the cellular level, NaHS shortens APD, which was more significant in the ischemia condition, whereas it decreased ERP in healthy but increased it in the ischemic condition; (ii) NaHS augmented the ERP difference between epicardial and endocardial myocytes (i.e., ERP) in both healthy and ischemia conditions and augmented the ERP difference between ischemic and non-ischemic cells; (iii) at the tissue level, it increased tissue vulnerability for initiation and maintenance of reentry, which is more significant in the ischemic condition. It also prolonged the life span of reentry. All taken together, these findings provide explanations to the pro-arrhythmic effects of H 2 S, especially in the ischemic conditions.

Action Potential Duration and Effective Refractory Period Abbreviation
NaHS shortens APD in both healthy and ischemic myocytes through an integrated action of depressed I CaL and enhanced I K,ATP . In simulation, APD abbreviation was more obvious in ischemia owing to more opened K-ATP channels. However, changes of ERP were different in these two situations, whereas NaHS decreased the ERP in healthy myocytes but increased it in ischemic myocytes, which is also attributable to the involvement of I K,ATP in ischemia (but not in health) that enhanced the repolarization current and reduced the amplitude of APs. NaHS thus prolonged the ERP of ischemic myocytes through its action on increasing I K,ATP .

NaHS Augmented the Effective Refractory Period Dispersion
The transmural dispersion of ERP, measured by the ERP difference (i.e., ERP) between epicardial and endocardial myocytes, was augmented by NaHS in both healthy and ischemia conditions. In the healthy condition, NaHS augmented ERP as it decreased ERP more in the epicardial cells (by 9 ms) than in endocardial cells (by 8 ms) (see Figure 3F), resulting in a 1ms increment in transmural ERP (Figure 3G). In the ischemia condition, NaHS also increased ERP but not in the same way as in healthy tissues. Instead, NaHS increased the ERP in both cell types as the I K,ATP was involved in the ischemia condition. However, such ERP increasing effect was more significant in epicardial myocytes (up to 29 ms) than in endocardial myocytes (only 3 ms), as illustrated in the two data groups in the right side of Figure 3F. The unevenly altered ERP in epicardial and endocardial cells contributed to a great enhancement of transmural ERP and further increased it by 26 ms. In addition, NaHS also increased the ERP difference between ischemic and non-ischemic epicardial myocytes ( Figure 3H). The observation suggested another potential pro-arrhythmic influence of NaHS when local ischemia happened, as the ERP dispersion between the local ischemic areas and the other healthy areas in the epicardium would possibly lead to arrhythmogenesis.

Increased Tissue Vulnerability for Arrhythmogenesis
The augmented repolarization dispersion by NaHS led to the increased temporal vulnerability of tissue to unidirectional conduction block in response to a premature stimulus (Figure 4). In addition, NaHS decreased the CV of excitation waves, which also contributed to the increased width of the measured VW that was exacerbated by the ischemic condition.

Decreased Minimal Substrate Size to Sustain Reentry
NaHS decreased the critical length in healthy and ischemic tissues for sustaining reentry, suggesting an important role in facilitating reentrant arrhythmia. This can be attributable to the abbreviated APD and ERP, resulting in a shorter excitation wavelength that reciprocally characterizes the spatial vulnerability of tissues (Rensma et al., 1988;Zhang et al., 2008).

Increased Life Span of Reentry
NaHS prolonged the life span of reentry in healthy tissues but not in ischemia. Owing to a limited tissue size (only 3 mm in transmural width) in this study, the initiated spiral wave always self-terminated with finite life span, which was determined by the interplay between the refractory distance of reentry and the transmural ERP gradient of the tissues. In one way, a short refractory distance (can be simply defined as the product of CV and ERP) of an excitation wave needs a small tissue size for accommodating the leading circle of reentry (Comtois et al., 2005), which helps to sustain the reentry. In the other way, the transmural gradient of ERP induces the tip of reentry to drifta larger ERP gradient produces a more meandering of the tip to run out of the boundary, leading to self-termination. In our simulation, the spiral wave self-terminated in the control case owing to the intrinsic ventricular heterogeneity. NaHS decreased the refractory distance but did not increase the ERP gradient too much, and the two factors together resulted in a prolonged life span of reentry.
In the ischemic condition, the depressed excitability led to a prolonged ERP, and the refractory distance was enlarged despite a decreased CV. In addition, the augmented ERP dispersion also contributed to the termination of the spiral wave by accelerating the drifts. In this case, these factors interplayed, resulting in a shorter life span of spiral wave in ischemic tissues. NaHS further prolonged ERP by enhancing I K,ATP , and it also increased the ERP dispersion; these effects together terminated the spiral wave earlier.

Potential Limitations of This Study
The limitations of the Pandit et al. (2001) model were discussed elsewhere (Terkildsen et al., 2008;Gattoni et al., 2016). In the simulations, the effects of H 2 S on cardiac ion channels were based on experimental data of a H 2 S donor, NaHS, as NaHS is superior to the direct administration of H 2 S for dose control, data reproducibility, and negligible effects on the cellular Na + concentration and pH (Zhong et al., 2003(Zhong et al., , 2010. In single-cell model development, data on the effects of NaHS on I CaL and I K,ATP were based on extant data from rats as much as possible to avoid potential species dependence (Furukawa et al., 1991). However, it should be of concern that the interspecies differences regarding the transmural I K,ATP may exist. For example, experiments conducted on rabbit (Qi et al., 2000) or cat (Furukawa et al., 1991) showed a greater I K,ATP activation in epicardial cells than in endocardial cells, but this is not applicable to rats according to the rat experiments (Shimokawa et al., 2007), which reported smaller half-maximal channel inhibition (IC 50 ) of I K,ATP -ATP curve using either inside-out patch recordings or open cell-attached patch recordings, suggesting less I K,ATP activation in epicardial cells than endocardial cells when facing the same decrement of intracellular ATP. Owing to the uncertainty of interspecies difference as well as the insufficient experimental data regarding transmural I K,ATP difference in rats, the transmural heterogeneity of I K,ATP was not incorporated in this study, and the conclusion should be carefully interpreted in clinical study considering the significant electrophysiological differences between human and rats. In 2D simulations, we only considered an earlier phase of ischemia as it was difficult to initiate the spiral wave in a tissue with conditions of 10-min ischemia. Specifically, the largely increased ERP led to a rather large leading circle that required large tissue size for spiral wave tip's meandering, and the limited size of a rat's transmural tissues will no longer be able to accommodate the spiral wave in the late phase of ischemia; thus, we did not conduct the simulation on tissues with conditions of 10-min ischemia.
Note that ischemia will also create the ERP dispersion between ischemic and non-ischemic areas, and NaHS could further enhance the dispersion in the epicardium. The role of such augmented ERP dispersion around the "border zone" of ischemic region may be pro-arrhythmic (Picard et al., 1999;Bernus et al., 2005;Trénor et al., 2007). In the presence of NaHS, how it increases the vulnerability of tissue to cardiac arrhythmias warrants future study.

CONCLUSION
The pro-arrhythmic effects of H 2 S arise as a consequence of different alterations of electrical properties in epicardial and endocardial myocytes through modest changes in ion channel conductance, resulting in abbreviated ERP, augmented transmural ERP dispersion, and depressed cell excitability, leading to increased tissue susceptibility for initiation and maintenance of cardiac arrhythmia. The pro-arrhythmic effect of H 2 S was more significant in ischemic tissues, providing a novel insight into higher risks among cardiac ischemia population to H 2 S-induced ventricular arrhythmias.

Modeling K-ATP Channel in Epicardial and Endocardial Myocytes
We used a detailed mathematical model of K-ATP channel developed by Ferrero et al. (1996). More details of the model can be found in the Supplementary Material. Briefly, I K,ATP is expressed as follows: where g 0 represents the maximum channel conductance (nS); [K + ] o is the extracellular concentration of K + (mmol/L); factors f M , f N , f T respectively denote the effects of intracellular Mg 2+ , Na + , and temperature. V and E K represent the membrane potential and reverse potential of K + , respectively, (mV); and f ATP is the fraction of opened channels. The model was further modified to incorporate species-and cell-type dependence and validated on available experimental data (see Supplementary Figure 1).

Modeling Effects of Exogenous H 2 S on Ion Channels
In experiments, NaHS was used in most cases as H 2 S donor rather than the direct administration of H 2 S. Therefore, our simulations were based on experimental data of NaHS. Experimental data have shown that NaHS increased K-ATP, whereas it reduced the Ca-L channel current in rat cardiac cells. Therefore, the dosedependent effect of NaHS on I K,ATP and I CaL was modeled by the following: I NaHS K,ATP = 1 + 0.71 1 + (K m,KATP /NaHS) 0.95 I K,ATP (2) These equations were fitted to experimental data (Sun et al., 2008;Zhong et al., 2010;Zhang et al., 2012), resulting in K m,KATP as 28.57 µmol/L and K m,CaL as 376.66 µmol/L. Fitting results and experimental data are available in Supplementary Figure 2.

Single-Cell Model
Healthy Rat Ventricular Myocyte Model The above-described I K,ATP was incorporated into the Pandit et al. (2001) model [with corrected equations by Terkildsen et al. (2008)] for the AP of rat ventricular myocytes.
[ATP] i and [ADP] i were set to 6,800 and 15 µmol/L for healthy tissues (i.e., normoxic condition) based on measured nucleotide data as reported in a previous study (Weiss et al., 1992). APs were evoked in a same way as Pandit et al. (2001), that is, by a series of supra-threshold stimuli (S1; 1 Hz) with an amplitude of 0.6 nA and duration of 5 ms. Detailed methods on measuring APD, ERP, and excitability using the model are described in the Supplementary Material.

Ischemic Rat Ventricular Myocyte Model
For simulating the effect of ischemia on cellular APs, three major conditions were considered (i.e., hyperkalemia, anoxia, and acidosis), as considered in Shaw and Rudy's study (Shaw and Rudy, 1997). Simulations of the effect of 10-min ischemia on nucleotide concentrations, ion concentrations, and channel current conductance are detailed as follows: The elevated extracellular K + were based on rat experimental data as reported in Fiolet et al. (1984), Knopf et al. (1990) and Wilde et al. (1990).

Anoxia
[ATP] i = 4, 600 µmol/L (7) [ADP] i = 99 µmol/L (8) Data of nucleotide levels in ischemia were obtained from Weiss et al. (1992). P CaL,ATP represents the influence of anoxia on the open fraction of Ca-L channel; and K m = 1,400 µmol/L and H = 2.6. The effects of the above modifications on AP configurations were validated for epicardial and endocardial myocytes using available experimental data. Detailed validation results are presented later in this section.
In addition to the 10-min ischemia model, an earlier phase of ischemia was also considered for a transient stage in 2D simulation. The parameters for the ischemic condition during the transient stage are as follows: g Na,acid = 0.85 · g Na (11) g CaL,acid = 0.85 · g CaL (12) [ATP] i = 5500 µmol/L (13) [ADP] i = 64.6 µmol/L The relationship between [ATP] i and [ADP] i was calculated by fitting the calculation function proposed by Bao et al. (2013).

Validation of Single-Cell Model for Control and Ischemic Conditions
Validations of the developed single-cell model of the rat ventricular APs are shown in Figure 6 for the control and ischemic conditions. The simulated ischemic condition as stated above based on the work of Shaw and Rudy (1997) abbreviated APs in both epicardial and endocardial myocytes; decreased their APA, maximum upstroke velocity (dV/dt max ), and CV (measured in 1D tissue strand model); but elevated the resting potential from −80 to around −71 mV. These simulated effects of ischemic condition matched well to experimental observations from rat and other species (Bélichard et al., 1991;De Diego et al., 2008), validating the developed rat ventricle cell models. Validation results are illustrated in Figure 6. All the data sources are available in Supplementary Table 1.

Transmural Ventricular Tissue Model
A multicellular tissue model for simulating the AP propagation was constructed using the well-established mono-domain equation (Clayton et al., 2011): where D is the diffusion coefficient matrix describing the intercellular electrical coupling via gap junctions. When considering the 1D tissue strand, Eq. 15 can be further organized as follows: where D is a scalar value that is set to 1.8 × 10 −2 mm 2 /ms in the 1D tissue model, giving a transmural CV of 19.84 cm/s, within the measured range of 14.6 to 102 cm/s in rat ventricular tissues (Sedmera et al., 2016). As the ventricular wall in rat heart is rather thin, a length of 3 mm of transmural 1D strand was constructed with 2 mm designated as endocardial and 1 mm epicardial tissues, which was inconsistent with our previous study (Zhang et al., 2009). The spatial resolution dx was set to 0.1 mm, close to the length of single ventricular myocyte (80-150 µm); thus, the constructed strand contained 20 nodes for endocardial and 10 nodes for epicardial tissues.
The 2D transmural tissue model was constructed by expanding the 1D strand in the y direction, thus forming a 3 × 10 mm 2 tissue sheet. The spatial resolution in the y direction dy is the same as dx. Anisotropic cell-to-cell electrical coupling CV was assumed in this study. Specifically, the fiber orientation was set to vertical, and the diffusion coefficients were set to 1.8 × 10 −2 mm 2 /ms in the fiber orientation and 7.2 × 10 −2 mm 2 /ms, which were perpendicular to the fiber orientation [i.e., coefficient ratio D ⊥ :D || = 1:4 according to Benson et al. (2008)]. The above settings gave a transmural CV of 19.84 cm/s and a longitudinal CV of 41.97 cm/s, approximately in a ratio of 1:2 in healthy tissues. In the ischemic condition, the coefficient was decreased by 20% to mimic the reduced intercellular coupling.

Measurement of Vulnerability to Unidirectional Conduction Block -Temporal Vulnerability
Temporal vulnerability to reentrant arrhythmia was quantitatively measured as the width of VW in 1D tissue model. Standard S1-S2 protocol (Zhang et al., 2009) was used. Specifically, supra-threshold stimuli (S1) were applied to the 0.3-mm segment at the endocardial end for five times at 1 Hz. Then a premature stimulus S2 was applied to a 0.6-mm segment with a time delay t after S1. At a certain period that just after the refractory tail of S1 induced AP wavefront, S2 could produce a solitary wave that propagated in the anterograde but not retrograde direction, or vice versa. This was due to the different refractory durations between two directions. If S2 was applied too early, propagation failed in both directions and unidirectional conduction block could not be formed. On the other hand, if it was applied too late, bidirectional conduction could occur. The time period during which a S2 stimulus evoked unidirectional conduction block is referred to as VW.
Vulnerable window at each point in the 1D strand (except the 0.3-mm segments at the two ends of strand owing to the 0.6mm stimulating length) was measured. In our simulation, the S1-S2 protocol was iteratively performed for every point until VWs across all strand were determined. The S1 stimuli were always applied to the first 0.3-mm endocardial segment in every measurement. As for the premature stimulus S2, the stimulating length was 0.6 mm, but its location was not fixed because every point on the strand had its own VW distribution that is to be determined in the simulation. For any given point, noted x, S2 was given to the segment from (x−0.3) mm to (x+0.3) mm after the S1 stimulus was applied. The interval between S2 and the last S1 (termed S1-S2 interval) was also not fixed, and different S1-S2 intervals would produce bidirectional conduction block (S2 too early), bidirectional conduction block (S2 too late), or unidirectional conduction block (S2 in VW). The time duration that S2 produced unidirectional conduction block was recorded as the width of VW for the point x.

Measurement of Critical Size of Tissue to Support Reentry -Spatial Vulnerability
Using a similar S1-S2 protocol in 2D tissues may initiate the reentrant excitation wave. S1 stimuli were applied to the first 0.3 mm of cells in endocardial side (left side) for five times at 1 Hz. Then the premature stimulus S2 was applied during the VW to a local vertical area at the bottom of the tissues, with fixed width of 0.6 mm and variable length (or height) from 3.0 to 8.0 mm. The center of stimulated area was 1.3 mm away from the endocardial side (left side) and was consistent in all cases (i.e., control, NaHS, ischemia, and ischemia + NaHS) to avoid location-dependent variation.
During VW, excitation wave evoked by S2 stimulus could be unidirectional conduction, which evolves into spiral wave. However, the spiral wave could sustain only when a sufficient length of S2 was given, as the spiral wave required space to survive. Or if the length of S2 was too small, the tip of the spiral wave would hit the bottom boundary of the tissues and would self-terminate. Therefore, a certain length (termed critical length) could represent the spatial vulnerability of tissues to reentry.
In 2D simulations, an earlier phase of ischemia was also considered, with mild hyperkalemia from 5.4 to 6.4 mmol/L, 15% decrease of calcium and sodium channel conductance, and a reduction of intracellular ATP from 6,800 to 5,500 µmol/L. Equations in this condition are listed in Ischemic Rat Ventricular Myocyte Model, that is, Eqs 10-14.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
HZ and ZW conceptualized the idea of the manuscript. HZ provided the basic mythological structure. SGZ, SZZ, and XF developed the relevant codes in the single-cell model. SGZ, XF, and WW developed the relevant codes in 1D and 2D tissue models. HZ, SZZ, ZL, and DJ analyzed the simulation results, including the validation of cell model output and the spiral wave results. SGZ drafted the original manuscript. SGZ, HZ, and ZW wrote the manuscript.