Critical Requirements for the Initiation of a Cardiac Arrhythmia in Rat Ventricle: How Many Myocytes?

Cardiovascular disease is the leading cause of death worldwide due in a large part to arrhythmia. In order to understand how calcium dynamics play a role in arrhythmogenesis, normal and dysfunctional Ca2+ signaling in a subcellular, cellular, and tissued level is examined using cardiac ventricular myocytes at a high temporal and spatial resolution using multiscale computational modeling. Ca2+ sparks underlie normal excitation–contraction coupling. However, under pathological conditions, Ca2+ sparks can combine to form Ca2+ waves. These propagating elevations of (Ca2+)i can activate an inward Na+–Ca2+ exchanger current (INCX) that contributes to early after-depolarization (EADs) and delayed after-depolarizations (DADs). However, how cellular currents lead to full depolarization of the myocardium and how they initiate extra systoles is still not fully understood. This study explores how many myocytes must be entrained to initiate arrhythmogenic depolarizations in biophysically detailed computational models. The model presented here suggests that only a small number of myocytes must activate in order to trigger an arrhythmogenic propagating action potential. These conditions were examined in 1-D, 2-D, and 3-D considering heart geometry. The depolarization of only a few hundred ventricular myocytes is required to trigger an ectopic depolarization. The number decreases under disease conditions such as heart failure. Furthermore, in geometrically restricted parts of the heart such as the thin muscle strands found in the trabeculae and papillary muscle, the number of cells needed to trigger a propagating depolarization falls even further to less than ten myocytes.


Introduction
Ill-timed propagated ectopic beats can initiate a potentially fatal arrhythmia, with the most serious being ventricular fibrillation [1,2]. Both atrial and ventricular fibrillation are a rapidly growing global health problem mainly due to aging of the human population, diseases such as heart failure, the adoption of unhealthy lifestyles, and channelopathies and other genetic factors in younger patient populations [3,4]. Despite efforts by clinicians and scientists, cardiac arrhythmia remains a significant health problem accounting for 15-20% of all deaths [5]. Previous studies have demonstrated that aberrations in Ca 2+ dynamics can play a role in cardiac arrhythmia [6][7][8][9][10]. These focal excitations in an individual myocyte can be caused by spontaneous Ca 2+ release resulting in Ca 2+ waves. A Ca 2+ wave in a ventricular myocyte is initiated when spontaneous Ca 2+ is released from the sarcoplasmic reticulum (SR) via ryanodine receptors ((RyR2)/Ca 2+ release channels) and diffuses to adjacent release sites (dyads) where they trigger additional Ca 2+ release. The spread of this process to different release sites results in a Ca 2+ wave. The rapid change of Ca 2+ results in a depolarization inward current through the Na + -Ca 2+ exchanger (I NCX ), which depolarizes the cell membrane sufficiently to trigger a beat, as demonstrated by previous experiments and computational modeling [6,[11][12][13].
When an action potential (AP) travels through cardiac tissue, the wavefront acts as a source of depolarizing current for the repolarized myocytes nearby (the sink). The source current density must be high enough to activate the sink, and if the source-sink mismatch is too significant, propagation will fail [14]. The propagating wavefronts of an AP have been studied extensively previously. Some earlier studies that have investigated this question suggest that as many as approximately~700,000 myocytes must undergo such behavior to initiate a propagating action potential or an arrhythmia in a block of heart tissue [14,15]. The development of successful, empirical methods for the treatment of cardiac arrhythmias is possibly only through knowledge of the underlying biology [16]. This paper addresses the basic questions of how many myocytes are needed to trigger cardiac arrhythmia in 1-D, 2-D, and 3-D tissue. Simulations in 1-D, 2-D, and 3-D represent physical structures in the heart (see [17] for more details). Simulations in 1D can represent fiber strands such as Purkinje or cardiac trabeculae [18]. Simulations in 2-D represent myocardia sheets reflecting that the heart is arranged in fiber sheets with the conductivity in the sheet being greater than between the sheets. The 3-D simulations can represent blocks of heart tissue such as the ventricular wall.
In order to understand how Ca 2+ dynamics plays a role in arrhythmogenesis, model simulations explore normal and dysfunctional Ca 2+ signaling in cardiac ventricular myocytes at a high temporal and spatial resolution. These findings suggest processes that lead to the initiation of an arrhythmia by current injection and spontaneous Ca 2+ release during heart failure. The goal of this study is to achieve a realistic and adequate quantitative understanding of the cardiac arrhythmias, i.e., the number of myocytes needed to trigger a sustained depolarization in a tissue. The model describes the most elementary event of cardiac Ca 2+ release, the calcium spark, with a stochastic model that explains the mechanisms of Ca 2+ release termination, graded Ca 2+ release, Ca 2+ homeostasis, and the sarcoplasmic reticulum calcium leak, and the generation of arrhythmias from defects in Ca 2+ signaling.
The computational modeling of cardiac arrhythmias presented here provided a new understanding of Ca 2+ -entrained arrhythmias and suggest that the aberrant activation of only a small number of myocytes is required to generate an arrhythmia. The model also suggests the role of recruitment of Ca 2+ release sites play in the initiation of arrhythmias.

Model
This study integrated our previous local-control model of excitation-contraction coupling rat ventricular myocytes into a network of myocytes to describe ventricular tissue [11,19]. Local-control is the feature that excitation-contraction coupling is modulated locally at the release site level. This model features the fully stochastic gating of ion channels from 20,000 calcium release sites using the algorithm developed in the previous paper [20]. Each release site contains a cluster of 49 RyR2 channels and a smaller cluster of 7 L-type calcium channels in closed proximity to each other, forming the dyadic subspace known as calcium release units (CRUs).
The model used in this paper has several important features, some of which were introduced in [11,[19][20][21], such as a spark model well-constrained by experimental data and a thermodynamically and biochemically constrained model for SERCA2a (sarco-and endoplasmic reticulum Ca 2+ ATPase isoform 2a) pump [22]. Other changes include using (1) explicit buffer dynamics in the subspace, (2) a different L-type Ca 2+ channel that incorporate Ca 2+ -bound calmodulin-dependent inactivation and updated parameters based on newer experimental data from rat ventricular myocyte [23]. In addition, by using the novel Cells 2022, 11, 1878 3 of 21 computational method and GPGPU (general purpose graphical processing unit) technology, it now allows us to do larger scale simulations that provide insights into calcium dynamics.
The current Ca 2+ spark model was developed using the mean-field energetic coupling between RyR2s distributed in an array that has been observed in biological systems [24,25]. The RyR2 model, which features both calcium-dependent activation from dyadic cleft Ca 2+ and RyR2 open probability dependence on SR lumen Ca 2+ , can explain the robust mechanism of spark initiation and spark termination [26]. The model has been used to explain the diastolic Ca 2+ leak with how stochastic gating of RyR2s and a thermodynamic SERCA2a model can balance SERCA2a-pump activity at the quiescent condition [20]. Recently, experiments have shown that dyad models with RyR2 gating without Ca 2+ -dependent inactivation produce more realistic Ca 2+ spark results than those that assumed [27]. By increasing SR buffer capacity, it increases Ca 2+ spark, which supports the idea that depletion of SR release plays a more dominant role in spark termination than Ca 2+ -dependent inactivation [28,29].
The whole-cell compartmental model employed the Ultra-fast Monte Carlo algorithm introduced in the (Williams, Chikando et al. 2011) for channel gating, with 20,000 CRUs per cell, and was simulated using Fermi CUDA-capable GPGPU [20,30]. The single cell model is incorporated into a tissue model in which the individual myocytes are coupled to the four neighboring myocytes using electrical coupling. The membrane voltage of each myocyte is governed by where V m is the transmembrane membrane potential, C sc is specific membrane capacitance. I Na is the whole cell Na + current and the two extraction pathways for the calcium via SL (sarcolemmal) are the plasma-membrane Ca 2+ /ATPase (PMCA) and the Na + /Ca 2+ exchanger (NCX). I background represents three background currents, Ca 2+ , Na + , and K + that follow the linear Ohm's Law. Aside from the background current and the Na + /K + exchanger, there are four other K + currents. The formulas for fast and slow transient outward currents (K tof and K tos , respectively) were created for mice by [31] for the detailed equations of each individual ionic currents, readers are advised to check our recently published articles [11,19]. The term n is the number of coupled neighbors and G k gap is the gap junction conductance between the myocytes and its kth neighbor. The tissue model has provided some insights into the number of myocytes, under different conditions, required to generate an ectopic heartbeat or sustained spontaneous AP propagation across the tissue.

Numerical Methods
The system of ordinary differential equations comprising the model is solved using the explicit Euler method. The simulation uses no-flux boundary conditions based upon the assumption that there was not a significant gradient across cells at the border of the simulation. The adaptive timestep (10-100 nanoseconds) is required for numerical stability and is also necessary to capture the fast and stochastic gating of DHPR (dihydropyridine receptor) and RyR2 channels. Each release site uses a different sequence of pseudo-random numbers to control channel gating. These Monte Carlo simulations are computed on Fermi-GPGPU cards, with pseudo-random numbers derived from the Saru pseudo-random number generator on GPGPU provided by Asfar et al. [32]. When the channel fires, a smaller time-step is selected; first to ensure numerical stability, second to limit maximum 10% of the CRUs having state changes to occur at a specified time [33,34]. Channel gating depends upon the local dyadic subspace (Ca 2+ ). The myoplasmic calcium indirectly links all the release sites through modulation of the subspace (Ca 2+ ). Thus, the fraction 10% was selected so that the amount of calcium release does not change much in the bulk Cells 2022, 11, 1878 4 of 21 myoplasm of the system so that these CRUs are assumed to experience the same bulk myoplasmic calcium.

Results
The modeling study considers two possible mechanisms for the initiation of any arrhythmia: (1) current injection to replicate the membrane channel activation of the action potential, and (2) spontaneous Ca 2+ release to model a Ca 2+ -entrained arrhythmia.

Current Injection
The first set of simulations explores how the depolarization of cardiac ventricular myocytes by current injection can initiate an arrhythmia. The compartmental rat ventricular myocytes model developed previously is integrated into a network of myocytes to simulate one-dimension (1-D), two-dimension (2-D), and three-dimension (3-D) tissues [20]. Simulations in 1-D, 2-D, and 3-D represent physical structures in the heart (see [17] for more details). Simulations in 1-D can represent fiber strands such as Purkinje or cardiac trabeculae [18]. Simulations in 2-D represent myocardia sheets reflecting that the heart is arranged in fiber sheets with the conductivity in the sheet being greater than between sheets. The 3-D simulations can represent blocks of heart tissue such as the ventricular wall. The introduction of a "coarse grained model", that describes a cardiac myocyte as a 1-D lattice of 200 independent Ca 2+ release sites (each of these sites represents the behavior of 100 identical Ca 2+ release sites to achieve a model with 20,000 total Ca 2+ release sites), increases computational efficiency. Here the depolarization during an AP propagates from a cardiac myocyte to an adjacent no-refractory myocyte due to a fire-diffuse-fire mechanism. The 1-D chain of myocytes is stimulated by current injection that is applied at one end of the cable (Figure 1, red region), all myocytes in the cable are connected by gap junctions with conductance G k gap , i.e., gap junctions enable propagation of the cardiac action potential from cell to cell. Figure 1A shows the schematic of 1-D tissue. In order to generate a sustained propagating action potential in 1-D tissue, current needs to be injected into only seven myocytes (magenta), as shown in Figure 1B. In these simulations, the gap junction conductance is estimated at 2500 nS for the AP conductance velocity V c to be 0.36 m/s, which is in the range of experimentally reported values [35,36]. The number of cells needed to initiate an arrhythmia will increase as the value of the gap junction conductance increases (Appendix B Figure A1). There is a close relationship between intracellular Ca 2+ and the membrane potential in cardiac myocytes. When the cell depolarizes, Ca 2+ influx through open L-type Ca 2+ channels in the cell membrane triggers the release of Ca 2+ by RyR2s on the sarcoplasmic reticulum into the dyad resulting in the myoplasmic Ca 2+ transient. A simulated propagating action potential and the accompanying triggered Ca 2+ wave is shown in Figure 1C. The excitation of seven myocytes is needed to initiate a propagating action potential. Appendix B Figure A2 shows that if five cells are activated, an action potential that can propagate to adjacent cells is not produced.
The simulations were extended into 2-D and 3-D tissue. In the 2-D homogeneous sheet, a 100 × 100 myocyte grid with anisotropic gap conductance is used. In this case the gap junction conductance is about three-fold (800 nS) greater along the longitudinal direction compared to the transverse direction. The anisotropy of conductance is visible in a 2-D tissue image in the propagation of AP shown in Figure 2A and the Ca 2+ wave propagation in Figure 2B. Due to the anisotropy of gap junction conductance, propagation is faster (elliptical shape) along the x-axis compared to the y-axis. The propagation of AP and Ca 2+ waves are shown in the Supplementary Videos S1 and S2. For the propagation of AP in 2-D tissue the number of myocytes increased to 49 myocytes and for the 3-D simulation (100 × 100 × 12 cells) with the anisotropy gap conductance the number of myocytes needed for the propagating AP was only 294. The 2-D simulations of AP and Ca 2+ release propagation illustrated in Figure 2A used a network of 100 × 100 rat ventricular myocytes with stochastic Ca 2+ dynamics. When the myocytes receiving the simulation The simulations were extended into 2-D and 3-D tissue. In the 2-D homogeneous sheet, a 100 × 100 myocyte grid with anisotropic gap conductance is used. In this case the gap junction conductance is about three-fold (800 nS) greater along the longitudinal direction compared to the transverse direction. The anisotropy of conductance is visible in a 2-D tissue image in the propagation of AP shown in Figure 2A and the Ca 2+ wave propagation in Figure 2B. Due to the anisotropy of gap junction conductance, propagation is faster (elliptical shape) along the x-axis compared to the y-axis. The propagation of AP and Ca 2+ waves are shown in the supplemental videos S1 and S2. For the propagation of AP in 2-D tissue the number of myocytes increased to 49 myocytes and for the 3-D simulation (100 × 100 × 12 cells) with the anisotropy gap conductance the number of myocytes needed for the propagating AP was only 294. The 2-D simulations of AP and Ca 2+ release propagation illustrated in Figure 2A used a network of 100 × 100 rat ventricular myocytes with stochastic Ca 2+ dynamics. When the myocytes receiving the simulation are not adjacent (spatially random stimulation pattern) the number of myocytes needed increases significantly. Simulation studies also investigated the response of the myocytes in tissue when an instantaneous and gradual current injection was applied by keeping the mass of the stimulus constant. In case of a gradual stimulus, a current of −6.52 µA/cm 2 was injected for 5 ms for a total current mass of the stimulus to be −32.60 µA/cm 2 ·mS. For the instantaneous case, a stimulus of −16.30 µA/cm 2 was applied for only 2 ms to yield the same total stimulus current mass. In the instantaneous case, the number of myocytes to produce a sustained propagating AP was reduced. Since myocytes in the tissue are coupled by gap junction conductance G gap k , the gradual current injection dissipates to the neighboring myocytes making it hard for the simulated myocytes to depolarize. On the other hand, with a large duration current injection of a short duration there is sufficient current to depolarize in the tissue. The numbers of cells and percentage of cells needed for generating a successful propagation of depolarization when the current is injected instantaneously and gradually are given in Table 1. The numbers are comparable in 1-D to experimental observations in Purkinje fibers [37].

Spontaneous Ca 2+ Release during Heart Failure (HF)
Spontaneous Ca 2+ release is another mechanism for generating depolarization in cells by the activation of NCX. Previous studies have used the whole-cell model to examine the changes in local calcium release in HF [38]. Here, simulations are performed to suggest that the mechanism behind, and the number of myocytes needed to trigger a propagating AP in a tissue by implementing previously established changes in ion transport proteins: fast and slow K + currents (I to and I K1 ), were reduced by 20%; NCX protein expression was increased by 100%; SERCA2a protein was decreased by 30%; and RyR2 Ca 2+ sensitivity was increased by 50% to mimic the increased activity from chronic hyperphosphorylation [7,[39][40][41]. Our previous studies have reported the experimentally observed transverse tubule (TT) changes during HF [38]. In this study a computational model was developed to capture these changes. This included increased spacing between TTs and RyR2 clusters at CRUs of only Cells 2022, 11, 1878 6 of 21 25% of TT junctions (whereas 75% of CRUs nanodomains remained unchanged), consistent with earlier studies [42]. This orphaning of RyR2 clusters during HF due to T-tubule remodeling was implemented for 25% CRUs through 30-fold increases in the subspace volume, whereas 75% of CRUs remains unchanged. The increase in subspace volume emulates the uncoupling between the L-type Ca 2+ channel and RyR2 during an AP. The parameter values for normal and HF myocytes are described in the Appendix A. Simulation studies also investigated the response of the myocytes in tissue whe instantaneous and gradual current injection was applied by keeping the mass o stimulus constant. In case of a gradual stimulus, a current of −6.52 μA/cm 2 was inje for 5 ms for a total current mass of the stimulus to be −32.60 μA/cm 2 ·mS. For instantaneous case, a stimulus of −16.30 μA/cm 2 was applied for only 2 ms to yield same total stimulus current mass. In the instantaneous case, the number of myocyt produce a sustained propagating AP was reduced. Since myocytes in the tissue coupled by gap junction conductance Ggap k , the gradual current injection dissipates to neighboring myocytes making it hard for the simulated myocytes to depolarize. On other hand, with a large duration current injection of a short duration there is suffi current to depolarize in the tissue. The numbers of cells and percentage of cells nee for generating a successful propagation of depolarization when the current is inje instantaneously and gradually are given in Table 1. The numbers are comparable in to experimental observations in Purkinje fibers [37].  In failing hearts, the increase in spontaneous Ca 2+ leak is simulated by an increase in the opening rate constants. Figure 3 shows an AP generated by the rat ventricular myocytes under the control and the heart failure condition. The AP is slightly prolonged in simulated heart failure compared to the control ( Figure 3A). The (Ca 2+ ) myo transient is slightly smaller and the diastolic level decreased during heart failure similar to the experimental observation in rat ventricular myocytes ( Figure 3B) [43]. Increased RyR2 phosphorylation restores the RyR2 open probability in heart failure levels nearly similar to the control during the Ca 2+ transient generated by the AP (Figure 3C). The diastolic RyR2 open probability is increased from 0.00051 for the control to 0.000204 for heart failure. The SR Ca 2+ content is decreased in heart failure similar to the experiment ( Figure 3D) due to the decreased SERCA expression and the increased diastolic Ca 2+ leak. Note that this increased RyR2 open probability during heart failure helps to restore the (Ca 2+ ) myo transient to near-control levels. The delayed subcellular Ca 2+ release and AP prolongation showed in the simulation results confirmed that the RyR2 cluster reorganization may contribute to dysfunctional AP and Ca 2+ dynamics (Supplementary Videos S3 and S4). Ca 2+ release is observed in cardiac myocytes by random opening of RyR2 channels located on the SR membrane. There is a possibility that Ca 2+ waves can be triggered via the Ca 2+ -induced Ca 2+ release phenomenon (CICR). CICR is the positive feedback property of the RyR2 that occurs when a rise calcium activates the RyR2 to release additional Ca 2+ . The rise in calcium can be from an influx from an L-type calcium channel, from spontaneous opening of the RyR2, or from the release of caged Ca 2+ . The abnormal rise of Ca 2+ induced by these mechanisms will be refered to as a spontaneous Ca 2+ wave. The rise of Ca 2+ in the myoplasm caused by the Ca 2+ wave activates the inward NCX current (forward mode I NCX ), which triggers an extrasystole and hence depolarizes the myocytes to produce the sustained propagating AP in the myocardial tissue. In the simulations shown in Figure 4, 10 channels in a RyR2 cluster of about 49 RyR2s are opened to release Ca 2+ into the dyadic cleft in 50% of the Ca 2+ release units in HF phenotype myocytes. The intial Ca 2+ release from the 10 activated channels in these HF myocytes will activate Ca 2+ release in the remaining channels within the same RyR cluster. In the simulations shown in Figure 4, 10 channels in a RyR2 cluster of about 49 RyR2s are opened to release Ca 2+ into the dyadic cleft in 50% of the Ca 2+ release units in HF phenotype myocytes. The intial Ca 2+ release from the 10 activated channels in these HF myocytes will activate Ca 2+ release in the remaining channels within the same RyR cluster. Figure 4A shows the NCX current in the tissue containing the failing myocytes (magenta). With the opening of the 50% of release a large forward mode NCX current is created because 20% of the NCX are in the dyadic space and thereore experience very high Ca 2+ levels. This large NCX current facilitates ectopic depolarization caused by spontaneous Ca 2+ release for failing myocytes. The AP produced in these failing myocytes propagates to the normal myocytes by the gap junction conductance G gap as shown in Figure 4B. Only seven myocytes (magenta Figure 4B) were needed in 1-D to produce a propagating AP.  Table 2 shows the number of cell needed to trigger an AP by current injection and spontaneous Ca 2+ release. Increasing the fraction of Ca 2+ release units from 50% to 80% in simulated 1-D, 2-D, and 3-D tissues with the HF phenotype decreases the number of myocytes required to depolarize the tissue from 512 to 100 in 3-D tissue, from 70 to 25 myocytes in 2-D tissue, and from 7 to 4 in 1-D tissue, as shown in Figure 4D.

Influence of Cardiac Geometry
The heart has many fine structures such as Purkinje fibers and trabeculae (lining the chambers of the heart) that are in effect a 1-D structure as shown in Figure 5A. Studies have suggested that an arrhythmia originates in trabeculae or Purkinje fibers [37]. The simulations of 3-D blocks of tissues such as the ventricular tissues require hundreds of myocytes, which is potentially too large a number to generate an arrhythmia as all these myocytes would have to simultaneously depolarize. Cardiac trabeculae that are found on the inner wall of the ventricular lumen were simulated as 1-D structures connected to a 3-D mass simulating the heart wall. Figure 5B shows a schematic diagram of this structure with three rectangles (1), injected current at the top in 2 × 2 × 8 myocytes show that a total of 32 myocytes is sufficient to depolarize the myocytes to generate an AP propagating into the 3-D tissue (ventricular wall). The number of myocytes required for the case with instantaneous current injection is reduced to only 2 × 2 × 5 = 20 myocytes. The simulations suggest that initiation site geometry plays an important role in triggering AP. The 20-32 myocytes required for triggering arrhythmia by gradual current injection compared to 292 myocytes in the ventricular wall suggests that structures such as the trabeculae or Purkinje system are likely sites for the initiation of arrhythmia. In the case of failing heart tissue, the activation of 20 myocytes by gradual current injection in the trabecula is sufficient to initiate depolarization in the ventricular wall. If current is injected instantaneously then the number of myocytes is reduced to 20 from 100. In case of the spontaneous Ca 2+ release in a failing heart, the number of myocytes is reduced to only 32 in the 1-D trabecula compared to 512 in the 3-D ventricular wall. The number is significantly smaller because the small size of the trabecular tissues reduced the electrotonic load that must be overcome to activate adjacent tissue.
(A) (B)  The depolarization of this 3-D structure by current injection and spontaneous Ca 2+ release was simulated. Figure 6 shows the triggered AP at different times after injecting the current at the top of the trabecula as it propagates into the 3-D tissue. The gradually injected current at the top in 2 × 2 × 8 myocytes show that a total of 32 myocytes is sufficient to depolarize the myocytes to generate an AP propagating into the 3-D tissue (ventricular wall). The number of myocytes required for the case with instantaneous current injection is reduced to only 2 × 2 × 5 = 20 myocytes. The simulations suggest that initiation site geometry plays an important role in triggering AP. The 20-32 myocytes required for triggering arrhythmia by gradual current injection compared to 292 myocytes in the ventricular wall suggests that structures such as the trabeculae or Purkinje system are likely sites for the initiation of arrhythmia. In the case of failing heart tissue, the activation of 20 myocytes by gradual current injection in the trabecula is sufficient to initiate depolarization in the ventricular wall. If current is injected instantaneously then the number of myocytes is reduced to 20 from 100. In case of the spontaneous Ca 2+ release in a failing heart, the number of myocytes is reduced to only 32 in the 1-D trabecula compared to 512 in the 3-D ventricular wall. The number is significantly smaller because the small size of the trabecular tissues reduced the electrotonic load that must be overcome to activate adjacent tissue.

Distributed Activation
The final set of simulations tested the effect of current injection in a non-contiguous mass of myocytes rather than an adjacent cell mass. The first current was injected into alternate myocytes. In this simulation successful depolarization of these myocytes was not possible even if current was injected into many myocytes because of the continuous leak to the unstimulated myocytes. If a non-contiguous mass of cell as shown in Figure 7A is stimulated (red myocytes stimulated), 80% more myocytes (nine myocytes compared to five) are needed to produce depolarization in the 2-D tissue. The numbers of myocytes not only increased, but also needed to have placed at least three myocytes at the end of the cable for the depolarization of the myocytes.
In 2-D tissue the number of myocytes increased almost three-fold to produce successful depolarization of the tissue ( Figure 7B). The red represents the stimulated myocytes, and the black represents the unstimulated myocytes in the tissue in the simulation. This simulation suggests that a large number of myocytes is needed to depolarize the tissue if the stimulated myocytes are not in a contiguous cell mass and suggests that even if a large number of uncoupled myocytes display an aberrant behavior, arrhythmia is unlikely the stability of normal heart rhythm. In order to generate arrhythmia, many myocytes in the close proximity need to display coordinated arrhythmogenic behavior.

Distributed Activation
The final set of simulations tested the effect of current injection in a non-contiguous mass of myocytes rather than an adjacent cell mass. The first current was injected into alternate myocytes. In this simulation successful depolarization of these myocytes was not possible even if current was injected into many myocytes because of the continuous leak to the unstimulated myocytes. If a non-contiguous mass of cell as shown in Figure 7A is stimulated (red myocytes stimulated), 80% more myocytes (nine myocytes compared to five) are needed to produce depolarization in the 2-D tissue. The numbers of myocytes not only increased, but also needed to have placed at least three myocytes at the end of the cable for the depolarization of the myocytes. myocytes, and the black represents the unstimulated myocytes in the tissue in the simulation. This simulation suggests that a large number of myocytes is needed to depolarize the tissue if the stimulated myocytes are not in a contiguous cell mass and suggests that even if a large number of uncoupled myocytes display an aberrant behavior, arrhythmia is unlikely the stability of normal heart rhythm. In order to generate arrhythmia, many myocytes in the close proximity need to display coordinated arrhythmogenic behavior. In 2-D tissue the number of myocytes increased almost three-fold to produce successful depolarization of the tissue ( Figure 7B). The red represents the stimulated myocytes, and the black represents the unstimulated myocytes in the tissue in the simulation. This simulation suggests that a large number of myocytes is needed to depolarize the tissue if the stimulated myocytes are not in a contiguous cell mass and suggests that even if a large number of uncoupled myocytes display an aberrant behavior, arrhythmia is unlikely the stability of normal heart rhythm. In order to generate arrhythmia, many myocytes in the close proximity need to display coordinated arrhythmogenic behavior.

Discussion
Cardiac arrhythmias require the aberrant behavior of a sufficient number of myocytes so that they can disrupt the normal depolarization pattern of the heart. This is true whether the arrhythmia is initiated by an ectopic focus or by the excitability/action potential dynamics of the cardiac myocyte. When an ectopic focus elicits a propagating action potential in front of the normal propagating depolarization wave front it can break the wave front, which can lead to spiral waves that are responsible for ventricular tachycardia, polymorphic ventricular tachycardia, and fibrillation [44]. In order for an ectopic focus to exist, it has to be protected from resetting by the normal depolarizing wave, which can occur due to a non-conductive zone (ischemia, fibrosis, and necrosis) or have a spontaneous release triggering a depolarization after the normal wave has passed [44][45][46][47]. Earlier modeling studies have suggested that a large number of myocytes (>700,000 myocytes) were required to initiate an arrhythmia [14]. Experimental studies that motivated their work indicated that an implanted mass of myocytes in the heart wall to take over heart pacing was approximately this size [15]. For such an implanted structure, the connections to the rest of the heart wall were most likely less well developed than connections between myocytes in the heart wall. Other simulation studies indicated that for the sinoatrial (SA) node, the depolarization of 1000 myocytes were needed to serve as an ectopic focus to initiate an arrhythmia [48]. Furthermore, the fine strands of atrial tissue that penetrate into and interdigitate with the SA node help to overcome the electrotonic loading effects [49].
The contribution of the Na + -Ca 2+ exchange to Ca 2+ removal varies across species. In rat ventricular myocytes, 92% of the Ca 2+ transient is removed by SERCA2a and 7% by the NCX. In humans and rabbits, 70-75% of the Ca 2+ transient is removed by SERCA2a, respectively, and 20-25% by the NCX [50,51]. Ca 2+ sparks in paced rat and human ventricular myocytes are similar. Ca 2+ sparks in human ventricular myocytes have been observed to have a F/F0 of 1.3 in normal heart and 1.19 in a failing heart [52]. Ca 2+ sparks in rat ventricular myocytes have been reported to have a F/F0 of 1.2-1.6 in paced ventricular myocytes [53,54]. However, in resting rat ventricular myocytes the F/F0 has been measured as high as 3.5-4 in the control and 3.0 post-myocardial infarction ventricular myocytes due to loading of the SR. These two observations suggest that there would be more depolarizing current through Na + -Ca 2+ exchange in human and rabbit ventricular myocytes and hence be easier to generate an arrhythmia.
The changes to the NCX expression during heart failure varies between species [55]. In humans, some studies have failed to find a correlation between NCX expression levels and the different stages of heart failure [56]. Simulations performed at different levels of expression indicated that with an increase NCX expression, there is a greater likelihood of generating an action potential that can propagate. This is consistent with the experimental findings that in heterozygous NCX knock-out mouse ventricular myocytes, the occurrence of delayed afterdepolarizations and early afterdepolarizations is reduced compared to the wild type [57]. However, it should be noted that if the aberrant Ca 2+ release is larger, a smaller amount of the NCX overexpression is needed to generate a propagating arrhythmia. Appendix B Figure A3 shows how the number of cells needed to initiate a propagating AP decreases with increasing the NCX maximum current density. Studies in humans have observed that NCX protein levels were inversely correlated with the frequency-dependent rise in the diastolic force [39]. Interestingly, the NCX protein levels have been positively correlated with the plasma norepinephrine levels in heart transplant patients [56].
Similarly, the changes to Ca 2+ dynamics will also vary across species and with the extent of progression of the heart failure. In compensated heart failure, beta adrenergic stimulation is able to rescue the myoplasmic Ca 2+ transient to be close to normal. In uncompensated heart failure the Ca 2+ transient is decreased in amplitude and of a prolonged duration. In the case of decompensated heart failure, the AP can also be very prolonged. In the model presented here, the Ca 2+ transient amplitude is only slightly decreased from the control because the beta-adrenergic stimulation has increased the RyR opening. The diastolic myoplasmic Ca 2+ level can also vary being elevated, decrease, or unchanged depending upon the species, pacing rate, and extent of heart failure in human patients even with depleted SR Ca 2+ [58]. The model suggests that this is likely related to the changes in NCX expression since the NCX controls the total intracellular Ca 2+ . In cases where there is upregulation of the NCX the diastolic myoplasmic Ca 2+ will be lower [58,59].
Previous studies have shown that changes to the T-tubular network also influence the propensity for spontaneous Ca 2+ release and waves [60,61]. Factors such as the increased distance between sites can limit Ca 2+ wave propagation and the increased amount of RyR2 phosphorylation can enhance wave propagation. Our earlier studies have reported the experimentally observed transverse tubule changes during HF [38]. In this study a computational model was developed to capture these changes. This included increased spacing between TTs and RyR2 clusters at CRUs of only 25% of TT junctions (whereas 75% of the CRUs nanodomains remained unchanged), consistent with earlier studies [42]. This orphaning of RyR2 clusters during HF due to T-tubule remodeling was implemented for 25% of CRUs through 30-fold increases in the subspace volume, whereas 75% of the CRUs remain unchanged. The increase in subspace volume emulates the uncoupling between L-type Ca 2+ calcium channel and RyR 2 during an AP.
The model simulated the major changes associated with heart failure. There are additional changes that could be studied in the future. For example, the scientific literature has reports that phospholamban changes during heart failure and other reports where it does not [62]. A more recent study claims that there is little evidence of phospholamban down-regulation in heart failure [63]. Due to the controversy over these changes, we have not included the effect of changes in phospholamban expression.
A simulation series was performed to understand the importance of synchrony of excitation in the excited cell mass. When the release in different myocytes in a cluster are temporally distributed, that is, the triggering Ca 2+ release events do not happen synchronously, it becomes more difficult to generate a propagating depolarization as the temporal spread of myocyte depolarizations gets larger. However, as long as the depolarizations overlap such as would occur if the depolarization spread from cell to cell, a small number of cells is still needed. Further studies are needed to explore how a Ca 2+ wave can spread between cells to form the mass of cells needed to trigger an arrhythmia.
The simulations presented here indicate that in 1-D-like structures, the excitation of as little as four myocytes is required to initiate a propagating action potential to the rest of the structure in agreement with experimental findings [37]. In fact, if the arrhythmia is initiated in a1-D-like structure such as a cardiac trabecula, Purkinje fiber, or the atrial penetrations into the SA node, the number of myocytes to propagate excitement to a mass of tissue such as the ventricular wall is small, as little as 20 myocytes. This can explain the experimental finding that arrhythmias originate in the trabeculae or Purkinje fibers. It can also explain how in damaged/fibrosed tissue the non-conductive scar tissue can isolate lower dimensional structures that lower the threshold for initiating an arrhythmia. The simulations suggest that in a 2-D sheet of myocytes. As few as 36 cells are needed for an ectopic focus to initiate a cardiac arrhythmia. In experiments with a cultured 2-D sheet of ventricular myocytes, they found that 8-50 myocytes were sufficient to initiate an arrhythmia at the edge of an ischemic region consistent with simulations of a similar condition [46,47].
The ventricular wall has different cell types in the epicardium, midmyocardium, and endocardium. There is also a transmural gradient moving across the wall transmurally and from apex to base. In the 3-D simulations of a tissue, there was no transmural gradient included in the simulation because in the block of ventricular tissue the number of cells in the simulations were small compared to the wall thickness. For simplicity a single ventricular cell type was used. This is justified because during the activation of the action potential they all have a fast upstroke velocity varying by less than 10% [64]. Furthermore, studies have suggested that a difference in the action potential duration rather than the initiation are important for the generation of arrhythmia [65,66].
In this study, compartmental models of the rat ventricular myocyte were used instead of a spatially resolved model. The simulation of spatially resolved ventricular myocytes can simulate spontaneous Ca 2+ waves and the mechanisms that contribute to their initiation and propagation [19]. Experimental studies have shown that a spontaneous Ca 2+ wave can be accompanied by the depolarization of the myocyte [67]. The depolarization will activate more Ca 2+ release units. Hence, activation of 50% of the release units was used to simulate spontaneous Ca 2+ release in the ventricular myocytes. Because depolarization of a single ventricular myocyte is uniform across the cell, compartmental models of the ventricular myocyte are sufficient for this study. The computationally expensive stochastic local control model used has advantages over the computationally fast common pool model. The local control model produces graded release and captures the important interactions between SR Ca 2+ release and the Ca 2+ -dependent inactivation of the L-type channel, which the common pool model fails to accomplish [11,20,38,68].
The origination of ectopic foci in both atrial and ventricular tissue has been studied previously. Atrial arrhythmia can originate at the pulmonary veins [69]. Other veinous structures have also been found to serve as the initiation sites for arrhythmias [70]. The simulations shown here demonstrate how a reduction in the dimensionality of the system can lower the number of myocytes needed for the generation of an arrhythmia. Structures such as the pulmonary veins effectively reduce the dimensionality from 3-D to 2-D. Furthermore, the ventricular wall structure contains sheets of cells that also could have a reduced dimensionality. The Purkinje fiber system and trabeculae have been identified as a source for arrhythmia [71][72][73][74][75]. These arrhythmias can arise as a result of ion channel dysfunction, for example through mutation. Alternately, aberrant Ca 2+ dynamics can also be a cause of these arrhythmia [71]. The Purkinje system and trabeculae are effectively a 1-D structure that serves as the conduction system to excite the ventricular tissue. In the Purkinje system, experiments show that as little as four myocytes need to depolarize together to create a propagation depolarization [37]. The simulations shown in the current study agree with this experimental finding that only a small number of myocytes are needed to trigger an arrhythmia given the appropriate geometry. Other studies have shown that arrhythmia can initiate in the trabeculae [18,75]. The trabeculae are columnar strands of muscle that are found on the interior of the ventricular walls. They effectively form a 1-D structure from which an arrhythmia can initiate with a small number of myocytes similar as demonstrated in the model simulations. Arrhythmia occurring during ischemia and reperfusion have been shown to involve calcium overload, particularly in the border zone between the ischemic and perfused regions [47,76]. In this border zone, gap junctions have been observed to decrease by 40% [77]. The combination of these two factors results in an increase in spontaneous depolarization and a reduction in dimensionality that can increase the possibility of the initiation of an arrhythmia by a mechanism described by the simulations in this paper.
Experimental studies have considered the remuscularization of the heart as a means of repairing damaged heart muscle [78,79]. Multiscale modeling studies have demonstrated that for implanted tissues, arrhythmias tend to initiate at the border of the implanted tissue [80]. This implanted tissue has been suggested to have a lower oxygenation that can impact graft survival [81]. This can lead to an ischemic region that can serve as an initiation site for arrhythmia, as described above.
One important question not addressed by the current modeling study is how a mass of myocytes synchronize to serve at the ectopic focus for initiating an arrhythmia. This has been studied in other modeling studies. Previous modeling studies have considered the coupling of an SA node cell using a 7 × 7 grid of myocytes [82]. These studies demonstrated how the coupling between myocytes plays a role in the synchronous firing of the SA node. A more recent modeling study demonstrated using a simplified model that at low pacing frequencies Ca 2+ waves do not contribute to synchronization [83]. However, at faster pacing frequencies the Ca 2+ wave can synchronize over millions of myocytes and lead to dramatic action potential fluctuations. The modeling in this paper does not spatially resolve individual myocytes and instead treats the myoplasm as one compartment as described previously [11]. The contribution of Ca 2+ waves can be studied using our stochastic spatial rat ventricular myocyte model and is left for future work [19].
The model predictions of the critical number of cells needed to initiate an ectopic propagating action potential suggest that a small number of cells is needed. These results might vary with other computational models for rat or for other species. For example, Xie and co-workers found very different numbers using a different model with different assumptions about the system being modelled [84]. For example, Wilhelms and co-workers discovered differences when using alternate atrial myocyte models on action potential properties in 1-D, 2-D, and 3-D [85]. Furthermore, modeling studies in ventricular myocytes have shown that the Na + channel dynamics can change model behaviors during arrhythmia [86,87]. This is a limitation of the current study. Future work can explore the effect of different models for rat and other species on the number of cells. Alternately, a perturbation analysis of the current model to explore how changes to the model affect simulation predictions similar to the approach used by Sutanto can be applied [88].

Conclusions
The experimentally constrained computational model developed here permits an understanding of how Ca 2+ sparks trigger spontaneous Ca 2+ release in the rat ventricular myocyte. The detailed 1-D, 2-D, and 3-D modeling of myocytes enables the investigation of how Ca 2+ instability arises, and how it may lead to novel electrical activity and arrhythmia. The model presented here supports the idea that spontaneous calcium release during conditions such as heart failure or in an EAD can result in an action potential through the activation of NCX. This is a Ca 2+ -activated electrogenic transport protein and is one of the critical components of the heart cell that links (Ca 2+ ) i changes with membrane current. Thus, when there is a large increase in (Ca 2+ ) i , as there is with every heartbeat, there is a parallel inward (depolarizing) current. Abnormal Ca 2+ release events activate the NCX current (I NCX ), which may trigger an extrasystole and hence produce an arrhythmia. The subcellular and nanoscopic characterization of (Ca 2+ ) i , which activates I NCX , has also been investigated. This requires high temporal and spatial resolution modeling of local Ca 2+ signals originating from the junctional sarcoplasmic reticulum (JSR). Both aspects of this behavior have been investigated.
Simulations suggest that in order to generate depolarization in the heart, the myocytes need to be in close proximity, and only then a small number of myocytes will be able generate arrhythmia in heart. If these stimulated myocytes were not close to each other, then number of the myocytes to produce successful depolarization increased a lot. We conclude that this could be the reason that makes arrhythmia a rare event. Finally, the simulations suggest that a relatively small number of myocytes are needed to trigger a propagating action potential with the number increasing as the dimensionality increases from 1-D to 2-D to 3-D. Fine structures in such as Purkinje fibers and trabeculae provide a virtual 1-D media in which an arrhythmia can be initiated. This greatly reduces the number or myocytes needed to initiate an arrhythmia in the 3-D heart.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cells11121878/s1, Video S1: Simulated action potential propagation across normal ventricular myocytes in a 100 × 100 network initiated by current injection; Video S2: Simulated action potential triggered calcium release across normal ventricular myocytes in a 100 × 100 network initiated by current injection; Video S3: Simulated action potential propagation across failing ventricular myocytes in a 100 × 100 network initiated by current injection; Video S4: Simulated action potential triggered calcium release across failing ventricular myocytes in a 100 × 100 network initiated by current injection; Video S5: Simulated action potential propagation across normal ventricular myocytes in a simulated trabeculae connected to the ventricular wall initiated by current injection; Video S6: Simulated action potential propagation across failing ventricular myocytes in a simulated trabeculae connected to the ventricular wall initiated by current injection. Model Parameters are available in Appendix A below. The rat ventricular myocyte model was recently published, and the parameters described below are for the normal and HF setting. Please check [11,19] for the detail and description of the parameters of the model and the Data Availability Statement.  Data Availability Statement: Model code is complex and might need explanation. The codes will be made available upon request to aullah3@gmu.edu.

Conflicts of Interest:
The authors declare no conflict of interest.  Figure A1. The number of cells needed to initiate a propagating AP increases as the gap junction conductance increases. The base gap junction conductance used in the simulations is 2400 ns. The formula and R 2 for trendline is included in the inset. Figure A1. The number of cells needed to initiate a propagating AP increases as the gap junction conductance increases. The base gap junction conductance used in the simulations is 2400 ns. The formula and R 2 for trendline is included in the inset. Figure A1. The number of cells needed to initiate a propagating AP increases as the gap junction conductance increases. The base gap junction conductance used in the simulations is 2400 ns. The formula and R 2 for trendline is included in the inset.   Figure A3. The number of cells needed to initiate a propagating AP increases as the maximal Na + -Ca 2+ exchanger maximal current density increases (equivalent to increase expression. The base value of vncx used in the simulation is 750 s −1 . The formula and R 2 for trendline is included in the inset.