Cardiac Alternans Occurs through the Synergy of Voltage- and Calcium-Dependent Mechanisms

Cardiac alternans is characterized by alternating weak and strong beats of the heart. This signaling at the cellular level may appear as alternating long and short action potentials (APs) that occur in synchrony with alternating large and small calcium transients, respectively. Previous studies have suggested that alternans manifests itself through either a voltage dependent mechanism based upon action potential restitution or as a calcium dependent mechanism based on refractoriness of calcium release. We use a novel model of cardiac excitation-contraction (EC) coupling in the rat ventricular myocyte that includes 20,000 calcium release units (CRU) each with 49 ryanodine receptors (RyR2s) and 7 L-type calcium channels that are all stochastically gated. The model suggests that at the cellular level in the case of alternans produced by rapid pacing, the mechanism requires a synergy of voltage- and calcium-dependent mechanisms. The rapid pacing reduces AP duration and magnitude reducing the number of L-type calcium channels activating individual CRUs during each AP and thus increases the population of CRUs that can be recruited stochastically. Elevated myoplasmic and sarcoplasmic reticulum (SR) calcium, [Ca2+]myo and [Ca2+]SR respectively, increases ryanodine receptor open probability (Po) according to our model used in this simulation and this increased the probability of activating additional CRUs. A CRU that opens in one beat is less likely to open the subsequent beat due to refractoriness caused by incomplete refilling of the junctional sarcoplasmic reticulum (jSR). Furthermore, the model includes estimates of changes in Na+ fluxes and [Na+]i and thus provides insight into how changes in electrical activity, [Na+]i and sodium-calcium exchanger activity can modulate alternans. The model thus tracks critical elements that can account for rate-dependent changes in [Na+]i and [Ca2+]myo and how they contribute to the generation of Ca2+ signaling alternans in the heart.


Introduction
The alternating strong and weak beats in the left ventricle are known as pulsus alternans or mechanical alternans which was first described in the 19th century by Traube [1]. Another type is electrical alternans (or T-wave alternans) which describe the beat-to-beat variation in direction, amplitude, and duration of any components in an ECG waveform [2]. The two are distinguished, yet they may both coexist [3,4].
Pulsus alternans is associated with different pathophysiological conditions, e.g., aortic stenosis, tachycardia, ischemia, acidosis and hypertrophic cardiomyopathy [5]. Other tested conditions that may lead to cellular or subcellular alternans include, but are not limited to: decrease RyR2 open probability (P o ) to increase the variability of Ca 2+ transient between the different areas of the cell, metabolic deficiencies, e.g., acidosis, and/or abnormal calcium handling [4,6,7]. Pulsus alternans may lead to pulseless activity, i.e., when there is an electrical activity, but the heart either does not contract, or the contraction is not strong enough to produce a sufficient cardiac output to generate a pulse [8]. Figure 1. Schematic diagram of a ventricular cell that shows only 2 T-tubule branches with 4 or the 20,000 Ca 2+ release units. Abbreviations: jSR-junctional sarcoplasmic reticulum; nSR-network sarcoplasmic reticulum; ryr-ryanodine receptor; T-tubule-transverse tubule; NCX-sodium-calcium exchanger; Na/K-sodium-potassium ATP ase; PMCA-plasmalemmal calcium ATPase; SERCA-sarcoplasmic and endoplasmic reticulum calcium ATPase; INa-sodium current; IKss-steady-state (sustained) potassium current; IKtof-fast transient outward potassium current; Itos-slow transient outward potassium current; IK1-inward rectifier potassium current; IbCa = background calcium current; IbK-background potassium current; IbNa-background sodium current.

Calcium Release Site (CRU)
A calcium release site is formed by the dyadic subspace and contains a cluster of 49 RyR2 channels and 7 LCCs channels. At each calcium release site, dynamic calcium buffering is implemented for three different endogenous buffers: calmodulin (CaM), sarcolemmal (SL) buffer and sarcoplasmic reticulum (SR) buffer rather than using the rapid buffering approximation or fixed buffering [27,30].

Ryanodine Receptor Type-2 Model
The 2-state ryanodine receptor model incorporates cytosolic calcium-dependent and luminal calcium dependent gating as described previously with only a small modification to the luminal dependence function to match spark characteristics [27]. The states are with the transition rate being increased by subspace Ca 2+ ([Ca ] ) and junctional SR Ca 2+ ([Ca ] ) and being a constant. The RyR2s are arranges in a cluster and display coupled gating as using our previous formulation [27].

Calcium Release Site (CRU)
A calcium release site is formed by the dyadic subspace and contains a cluster of 49 RyR2 channels and 7 LCCs channels. At each calcium release site, dynamic calcium buffering is implemented for three different endogenous buffers: calmodulin (CaM), sarcolemmal (SL) buffer and sarcoplasmic reticulum (SR) buffer rather than using the rapid buffering approximation or fixed buffering [27,30].

Ryanodine Receptor Type-2 Model
The 2-state ryanodine receptor model incorporates cytosolic calcium-dependent and luminal calcium dependent gating as described previously with only a small modification to the luminal dependence function to match spark characteristics [27]. The states are with the transition rate k + ryr being increased by subspace Ca 2+ ( Ca 2+ i ds ) and junctional SR Ca 2+ ( Ca 2+ i jsr ) and k + ryr being a constant. The RyR2s are arranges in a cluster and display coupled gating as using our previous formulation [27]. (1)

L-Type Ca 2+ Channel Model
The 6-state L-type Ca 2+ channel (LCC) model, Figure 2, is derived from 5-state LCC model for rat ventricle myocyte from Sun and co-workers with parameters adjusted for the new spark model [31,32]. The Sun model was developed to work in the range of −30 mV + 30 mV of the transmembrane potential reproducing single channel dwell times and currents from populations of channels matched to experiments. As Sun and coworkers suggested in their original paper, the 6-th state C6 was added to work with stronger depolarization (<−40 mV), so that all the channels stay in this state when the cell is at rest. The 6-state L-type Ca 2+ channel (LCC) model, Figure 2, is derived from 5-state LCC model for rat ventricle myocyte from Sun and co-workers with parameters adjusted for the new spark model [31,32]. The Sun model was developed to work in the range of -30 mV + 30 mV of the transmembrane potential reproducing single channel dwell times and currents from populations of channels matched to experiments. As Sun and co-workers suggested in their original paper, the 6-th state C6 was added to work with stronger depolarization (< -40 mV), so that all the channels stay in this state when the cell is at rest.

Figure 2.
Schematic diagram of the model (top) for the L-type calcium channel (DHPR channel) from [29]. C1 and C6 are closed states with C6 the resting state. O2 and O3 are the open states that conduct Ca 2+ ions with equal conductance. C4 is Ca 2+ -dependent inactivated state. C5 is the voltagedependent inactivated state. The fraction of LCC channels in different states during the voltageclamp of Vstim = +10 mV for 200 ms (bottom).
The inactivation of the LCC was modelled via 2 separate pathways: VDI (voltagedependent inactivation (O2C5)) and CDI (Ca 2+ -dependent inactivation (O2C4)). At each release site, the Ca 2+ level that controls the inactivation is the subspace calcium. The source of contributing calcium to this microdomain is the influx of calcium via LCC and the release of calcium from SR via RyR. The much higher level of calcium in this subspace, i.e., ~100-fold compared to cytosolic bulk calcium, enhances the rate of inactivation of LCC and thereby preventing calcium overload.
The endogenous buffer CaM bound to Ca 2+ (CaCaM) is the effector for Ca 2+ -dependent inactivation (CDI) of L-type Ca 2+ channel. The calcium-free CaM is called apo calmodulin (apo-CaM) with two homologous domains, known as lobes [31]. For each apo-CaM molecule, there are four different calcium binding sites: two at the E-F hand motifs of the N-terminal (N-lobe) and two at the E-F hand motifs of the C-terminal (C-lobe) of CaM (for review see [32]). In L-type Ca 2+ channels (Cav1.2), CDI is triggered by the binding of two Ca 2+ ions to the C-lobe of CaM [33][34][35]. It means that to model CDI, two Ca 2+ binding at C- Schematic diagram of the model (top) for the L-type calcium channel (DHPR channel) from [29]. C1 and C6 are closed states with C6 the resting state. O2 and O3 are the open states that conduct Ca 2+ ions with equal conductance. C4 is Ca 2+ -dependent inactivated state. C5 is the voltage-dependent inactivated state. The fraction of LCC channels in different states during the voltage-clamp of V stim = +10 mV for 200 ms (bottom).
The inactivation of the LCC was modelled via 2 separate pathways: VDI (voltagedependent inactivation (O2→C5)) and CDI (Ca 2+ -dependent inactivation (O2→C4)). At each release site, the Ca 2+ level that controls the inactivation is the subspace calcium. The source of contributing calcium to this microdomain is the influx of calcium via LCC and the release of calcium from SR via RyR. The much higher level of calcium in this subspace, i.e.,~100-fold compared to cytosolic bulk calcium, enhances the rate of inactivation of LCC and thereby preventing calcium overload.
The endogenous buffer CaM bound to Ca 2+ (CaCaM) is the effector for Ca 2+ -dependent inactivation (CDI) of L-type Ca 2+ channel. The calcium-free CaM is called apo calmodulin (apo-CaM) with two homologous domains, known as lobes [31]. For each apo-CaM molecule, there are four different calcium binding sites: two at the E-F hand motifs of the N-terminal (N-lobe) and two at the E-F hand motifs of the C-terminal (C-lobe) of CaM (for review see [32]). In L-type Ca 2+ channels (Cav1.2), CDI is triggered by the binding of two Ca 2+ ions to the C-lobe of CaM [33][34][35]. It means that to model CDI, two Ca 2+ binding at C-lobe is necessary. This differs from the original Sun model, which used the Hill coefficient as 3 based upon earlier data [36]. In addition, the original model didn't incorporate the loss of calcium in the subspace due to binding to calmodulin (CaM). This is important as the level of free [Ca 2+ ] ds control the gating of RyR2 channels. This was also corrected in our modified model. Experiments suggest the existence of non-junctional DHPR (10-20%), located on the external sarcolemma and not forming the release sites with RyR2 [37,38]. Hence the contribution of calcium from a small fraction of DHPR I dhpr,nj (15%) was also added.

Na + Channel Model
The rapid inward Na + current in the cardiac cells exhibits a bi-phasic time courses for recovery during inactivation and one activation gating variable [39]. We have modified the Pandit model so that the model more closely simulates data for the adult rat ventricular myocyte than does the more commonly used Luo-Rudy model as shown in Figure 3 [40][41][42][43][44]. The new model better simulates the correct upstroke velocity for rat as 213.8 ± 6.6 V/s with the overshoot is 48.76 ± 1.09 (mV). Experiments estimate the velocity range as 150-190 V/s at 23 • C V/s in Wistar rat heart which is calculated to be 184-233 V/s using a Q10 of 1.23 [42,45,46]. The new model also better estimates the value of h∞ to be 0.812 at −80 mV at body temperature. This value of h ∞ has been experimentally estimated as 0.7-0.8 at −80 mV in rat and mice [39,47]. In contrast, the value of h∞ in Pandit is 0.6; and the value in Luo-Rudy is 0.97. Our estimate takes into consideration the temperature dependence for recovery kinetics with the Q10 ranging from 1.7 to 2.3 for membrane potential in the range −76 mV to −62 mV, and 1.5 to 1.8 for the membrane potential in the range −84 mV to −65 mV [43]. The maximum current conductance is 8 mS/cm 2 which is in the range 3-25 mS/cm 2 [39,48].

K + Channel Models
There are four different K + channel currents (I tof , I, I K1 , and I Kss ). The formula for fast and slow transient outward currents (I Ktof , I Ktos , respectively) are based on the model based on experimental data observed in mice [49].

Sarcoplasmic Reticulum Ion Pumps
The sarco (endo)plasmic reticulum Ca 2+ -ATPase (SERCA) pump re-sequesters Ca 2+ back to the SR/ER during each excitation-contraction cycle to facilitate muscle relaxation by pumping two calcium ions per ATP molecule hydrolyzed [50]. We used the 2-state formulation by Tran and co-workers developed because it is constrained both by the thermodynamic and kinetic data for the SERCA pump [51].

Calcium Buffers
The three endogenous buffers of calmodulin (CaM), troponin (Trpn), and the phospholipids of the SR membrane (SRbuf) are used for the bulk myoplasm. The troponin complex consists of three different subunits. The troponin complex as modeled includes the binding of calcium (troponin C), the inhibition of actomyosin interaction (troponin I), and the binding to tropomyosin (troponin T).

Membrane Potential
During AP, the dynamics of the membrane voltage V m are governed by the ionic currents described above.

K + Channel Models
There are four different K + channel currents (Itof, I, IK1, and IKss). The formula for fast and slow transient outward currents (IKtof, IKtos, respectively) are based on the model based on experimental data observed in mice [49].

Sarcoplasmic Reticulum Ion Pumps
The sarco (endo)plasmic reticulum Ca 2+ -ATPase (SERCA) pump re-sequesters Ca 2+ back to the SR/ER during each excitation-contraction cycle to facilitate muscle relaxation by pumping two calcium ions per ATP molecule hydrolyzed [50]. We used the 2-state formulation by Tran and co-workers developed because it is constrained both by the thermodynamic and kinetic data for the SERCA pump [51].

Calcium Buffers
The three endogenous buffers of calmodulin (CaM), troponin (Trpn), and the phospholipids of the SR membrane (SRbuf) are used for the bulk myoplasm. The troponin complex consists of three different subunits. The troponin complex as modeled includes the binding of calcium (troponin C), the inhibition of actomyosin interaction (troponin I), and the binding to tropomyosin (troponin T).

Numerical Methods
We used our patented Ultra-fast Monte Carlo Simulation Method to solve for the states of the RyR2 and L-type Ca 2+ channels in each release unit [52]. Using this method the RyR2 channels at each release site are modeled as a single stochastic cluster, and the L-type channels are modeled similarly. Rather than keeping track of individual channel's state, we used a mean-field approach in which we assume all channels in the cluster see the same local calcium concentration in the dyadic subspace [53,54]. Thus, the individual channel's states are ignored, and only the number of channels in each state is important. Each release site is fed with a different sequence of pseudo-random numbers. These Monte-Carlo simulations are computed on Fermi-GPU cards, with pseudo-random numbers were derived from the Saru PRNG algorithm implemented on GPU provided by Steve Worley (Private communication at GTC'12) [55]. Instead of using a fixed timestep, an adaptive time-step strategy is used. 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 time [56,57]. This limits Type II error with the hypothesis that there is only channel state transition in the cluster per time step. In fact, when a full Monte Carlo Simulation is performed there are two channels undergoing state transitions in each timestep <0.6% of the time.
The system of ordinary differential equations comprising the model is solved using the explicit Euler method. The small and adaptive timestep (10-100 ns) which is required to simulate the fast and stochastic gating of DHPR and RyR2 channels is sufficient to ensure numerical stability.

Results
The model integrates the complex mechanisms involved in excitation-contraction coupling by describing the 20,000 stochastic calcium release units. In the model components were validated in the model described above and the model dynamics below in the results section. For example, the model demonstrates the same mechanism of release as our previous work and fully accounts for the SR Ca 2+ visible and invisible leak by flux through the RyR2 channels in the forms of Ca 2+ sparks and non-spark openings, respectively ( Figure A1) [27,58,59]. Details of the ionic currents are shown in Figure A2. Figure 4 shows for 1 Hz pacing the time courses for a train of action potentials, myoplasmic calcium transients, network, and SR calcium transients. In our model, the ratio of SR calcium release over the influx of calcium during a twitch is 10.0 ± 0.3. It means that, on average, the SR-release contributes about 90.07% and calcium influx contributes 9.03%. This approximates the value 92% of SR contribution estimated for rat ventricular myocytes [9,60]. During one cycle, our model produced a total extrusion of calcium of about 11.5 µM/cell/s. Among this,~1 µM is extruded by the PMCA and~10.5 µM is extruded by the NCX. This is in agreement with the value given above by [61]. The ionic currents during a twitch-relaxation cycle are given in Figure A2. The model also simulates pacing, with APD20~10.69 ms, APD50~22.57 ms and APD90~78.3 ms as is evident in Figure 4.

Mechanisms of Alternans
The model produces stable action potential and calcium transient trains below 8 Hz, however, above 8 Hz pacing frequency, periods of alternans are observed. One such period of AC concordant alternans is shown in Figure 5

Mechanisms of Alternans
The model produces stable action potential and calcium transient trains below 8 Hz, however, above 8 Hz pacing frequency, periods of alternans are observed. One such period of AC concordant alternans is shown in Figure 5 with both [Ca 2+ ] myo and V m alternating in phase. There are also alternating small and large amplitude and duration of other observable quantities such as [Ca 2+ ] SR , RyR P O , and I Ca . It is also important to note that the amplitude of I Na and I CaL are greatly reduced. Note that both the [Ca 2+ ]myo and [Ca 2+ ]SR fail to recover back to the normal diastolic level between systoles.
Many studies have suggested that alternans is a defect of the calcium subsystem. The conceptual framework of refractoriness, randomness, and recruitment has been used to understand alternans [62]. To study the refractoriness of calcium release beat-to-beat changes of the CRU's states were analyzed. A CRU is activated if a Ca 2+ spark occurs. The CRUs are grouped into either activated or inactivated at each beat. Therefore, when considering a pair of beats, there are 4 groups: act-act, act-inact, inact-act, and inact-inact.
In the transition from beat-1 to beat-2, there is a large fraction of CRUs changing from Inactivated to Activated (green), and also a large fraction of CRUs which are Activated in beat-1 continue to Activate in beat-2 (red)-the beat that has a strong contraction. For example, the number of activated release sites increases from Beat1-2 as indicated by the increased green and red. Corresponding to this, Figures 6B shows that beat 2 has a higher RyR open probability than beat-1. This suggests that alternans might originate at the level of the CRU. Figure 6C tracks the [Ca 2+ ]jsr local depletion for individual CRUs as shows that the CRUs (red, blue, green, and black) activate at alternating beats. Note that both the [Ca 2+ ] myo and [Ca 2+ ] SR fail to recover back to the normal diastolic level between systoles.
Many studies have suggested that alternans is a defect of the calcium subsystem. The conceptual framework of refractoriness, randomness, and recruitment has been used to understand alternans [62]. To study the refractoriness of calcium release beat-to-beat changes of the CRU's states were analyzed. A CRU is activated if a Ca 2+ spark occurs. The CRUs are grouped into either activated or inactivated at each beat. Therefore, when considering a pair of beats, there are 4 groups: act-act, act-inact, inact-act, and inact-inact.
In the transition from beat-1 to beat-2, there is a large fraction of CRUs changing from Inactivated to Activated (green), and also a large fraction of CRUs which are Activated in beat-1 continue to Activate in beat-2 (red)-the beat that has a strong contraction. For example, the number of activated release sites increases from Beat1-2 as indicated by the increased green and red. Corresponding to this, Figure 6B shows that beat 2 has a higher RyR open probability than beat-1. This suggests that alternans might originate at the level of the CRU. Figure 6C tracks the [Ca 2+ ] jsr local depletion for individual CRUs as shows that the CRUs (red, blue, green, and black) activate at alternating beats.
Membranes 2021, 11, x FOR PEER REVIEW 10 of 37 Figure 6. Alternans occurs at 8 Hz pacing rate. (A) An example that show a CRU that fires at two contiguous beats, e.g., (X) marks an individual release site that is Act-Act. (B) Alternans in the calcium release seen by the changing amplitude and duration of the [Ca 2+ ]myo. Beat-to-beat variation in CRU's states where we examine act-act-the fraction of CRU that activate in beat (i) (red) and continue to activate in beat (i + 1), similarly with act-inact (blue), inact-inact (yellow), and inact-act (dark green). However, a Chi-squared test with the 2 × 2 contingency table shown in Figure 7A did not support this hypothesis (p < 0.61). Even though there are more L-type Ca2+ channels assuming the Ca2+ dependent inactivated state during 8 Hz pacing, their number does not exceed 2 or 3 per CRU out of the 7 available channels. This leaves sufficient number of channels in a state to be activated during the subsequent beat. The next hypothesis is that after a CRU has been triggered it is less likely to activate again in the following beat. To quantify this hypothesis, whether or not a spark occurred in the first beat versus whether or not a spark occurred in the second beat, we used a Chi-squared test with a 2 × 2 contingency table ( Figure 7B). With a Chi-squared statistic of 55.6 and a p-value of p < 0.001, this showed that if a spark occurred in the first beat for a given CRU, there was a slightly lower probability that it would occur in the next. Figure 6B depicts the explanation for this. The jSR depletes after a CRU triggers and does not fully recover until the next beat. This lowers the likelihood of CRU activation both through the effects of jSR calcium Figure 6. Alternans occurs at 8 Hz pacing rate. (A) An example that show a CRU that fires at two contiguous beats, e.g., (X) marks an individual release site that is Act-Act. (B) Alternans in the calcium release seen by the changing amplitude and duration of the [Ca 2+ ] myo . Beat-to-beat variation in CRU's states where we examine act-act-the fraction of CRU that activate in beat (i) (red) and continue to activate in beat (i + 1), similarly with act-inact (blue), inact-inact (yellow), and inact-act (dark green). However, a Chi-squared test with the 2 × 2 contingency table shown in Figure 7A did not support this hypothesis (p < 0.61). Even though there are more L-type Ca 2+ channels assuming the Ca 2+ dependent inactivated state during 8 Hz pacing, their number does not exceed 2 or 3 per CRU out of the 7 available channels. This leaves sufficient number of channels in a state to be activated during the subsequent beat. The next hypothesis is that after a CRU has been triggered it is less likely to activate again in the following beat. To quantify this hypothesis, whether or not a spark occurred in the first beat versus whether or not a spark occurred in the second beat, we used a Chi-squared test with a 2 × 2 contingency table ( Figure 7B). With a Chi-squared statistic of 55.6 and a p-value of p < 0.001, this showed that if a spark occurred in the first beat for a given CRU, there was a slightly lower probability that it would occur in the next. Figure 6B depicts the explanation for this. The jSR depletes after a CRU triggers and does not fully recover until the next beat. This lowers the likelihood of CRU activation both through the effects of jSR calcium on RyR open probability and in the event if a RyR does open, the reduced flux will be insufficient to trigger adjacent RyRs. As a result, the model suggests that refractoriness in jSR is an important component of alternans.  These results suggest that the recovery of junctional SR Ca 2+ after a release event can contribute to the development of alternans. To assess the SR Ca 2+ content, we split the CRUs at each beat in two groups: those that are activated (Act) and those are not activated (Inact) and investigated two things: (1) the values of [Ca 2+ ]jSR at these release sites, right at before the next beat to start, (2) the nadir of [Ca 2+ ]SR during calcium release ( Figure 8). Larger SR [Ca 2+ ] comes with larger Ca 2+ release and the beats with larger release have more activated CRUs. We can see that the ensemble initial [Ca 2+ ]SR is not an indicator of alternans formation. If we look at the nadir [Ca 2+ ]SR, the result shows a similar small beat-tobeat alternation, in agreement with the study of [16]. These results suggest that the recovery of junctional SR Ca 2+ after a release event can contribute to the development of alternans. To assess the SR Ca 2+ content, we split the CRUs at each beat in two groups: those that are activated (Act) and those are not activated (Inact) and investigated two things: (1) the values of [Ca 2+ ] jSR at these release sites, right at before the next beat to start, (2) the nadir of [Ca 2+ ] SR during calcium release ( Figure 8). Larger SR [Ca 2+ ] comes with larger Ca 2+ release and the beats with larger release have more activated CRUs. We can see that the ensemble initial [Ca 2+ ] SR is not an indicator of alternans formation. If we look at the nadir [Ca 2+ ] SR , the result shows a similar small beat-to-beat alternation, in agreement with the study of [16].
The mechanism by which randomness of calcium release arises in the system is also explained by the model. The reduction in I Na current which reduced the action potential magnitude contributed a huge effect in reducing the number of L-type channels opening that opens the gateway for alternans ( Figure 5). When there is less I Na ( Figure 5E), there is a smaller AP ( Figure 5C) and reduced L-type current ( Figure 5D). These beat also have a smaller RyR open probability ( Figure 5B) and Ca 2+ transient ( Figure 5A) This is in agreement with the result in which a small depolarizing potential under voltage clamp can induce alternans [16], and confirms the hypothesis that SR Ca 2+ release is graded with the beat-to-beat alternation in I CaL [63,64]. The reduction in L-type opening reduces the synchronization of activation of the CRUs. Furthermore, the elevated myoplasmic [Ca 2+ ] increases the stochastic activation of Ca 2+ sparks. At the normal quiescent condition, as shown in our previous study, the activation of a CRU is the result of 8-10 RyR2 opening [27]. During field-stimulus simulation with high pacing, this number reduces to 5, where Ca 2+ alternans is observed (Figure 9). During the large beat, there can be up to 80% of CRUs activated, though the time to peak of activation is slower, e.g., 60%. The mechanism by which randomness of calcium release arises in the system explained by the model. The reduction in INa current which reduced the action po magnitude contributed a huge effect in reducing the number of L-type channels o that opens the gateway for alternans ( Figure 5). When there is less INa ( Figure 5E), a smaller AP ( Figure 5C) and reduced L-type current ( Figure 5D). These beat also smaller RyR open probability ( Figure 5B) and Ca 2+ transient ( Figure 5A) This is in ment with the result in which a small depolarizing potential under voltage clam induce alternans [16], and confirms the hypothesis that SR Ca 2+ release is graded w beat-to-beat alternation in ICaL [63,64]. The reduction in L-type opening reduces th chronization of activation of the CRUs. Furthermore, the elevated myoplasmic [C creases the stochastic activation of Ca2+ sparks. At the normal quiescent condit shown in our previous study, the activation of a CRU is the result of 8-10 RyR2 o [27]. During field-stimulus simulation with high pacing, this number reduces to 5, Ca 2+ alternans is observed (Figure 9). During the large beat, there can be up to CRUs activated, though the time to peak of activation is slower, e.g., 60%.  The mechanism by which randomness of calcium release arises in the system is also explained by the model. The reduction in INa current which reduced the action potential magnitude contributed a huge effect in reducing the number of L-type channels opening that opens the gateway for alternans ( Figure 5). When there is less INa (Figure 5E), there is a smaller AP ( Figure 5C) and reduced L-type current ( Figure 5D). These beat also have a smaller RyR open probability ( Figure 5B) and Ca 2+ transient ( Figure 5A) This is in agreement with the result in which a small depolarizing potential under voltage clamp can induce alternans [16], and confirms the hypothesis that SR Ca 2+ release is graded with the beat-to-beat alternation in ICaL [63,64]. The reduction in L-type opening reduces the synchronization of activation of the CRUs. Furthermore, the elevated myoplasmic [Ca 2+ ] increases the stochastic activation of Ca2+ sparks. At the normal quiescent condition, as shown in our previous study, the activation of a CRU is the result of 8-10 RyR2 opening [27]. During field-stimulus simulation with high pacing, this number reduces to 5, where Ca 2+ alternans is observed (Figure 9). During the large beat, there can be up to 80% of CRUs activated, though the time to peak of activation is slower, e.g., 60%. Figure 9. The fraction of CRUs that has 1 or more RyR openings during each beat (black = beat-1, red = beat-2, blue = beat-3, yellow = beat-4, green = beat-5). Figure 9. The fraction of CRUs that has 1 or more RyR openings during each beat (black = beat-1, red = beat-2, blue = beat-3, yellow = beat-4, green = beat-5).
We then explored the hypotheses that a high level of [Na + ] i increases the likelihood of alternans through the modulation of Na + /Ca 2+ exchanger activity (NCX). The normal function of NCX is to extrude one calcium ion in exchange for three sodium ions. It has been hypothesized that at high level of [Na + ] i , it may fail to extrude calcium, thus maintaining a high diastolic cytosolic [Ca 2+ ]. In other words, increasing Na + lead to more reverse mode NCX, leading to increases [Ca 2+ ] myo and [Ca 2+ ] SR . With 8 Hz pacing simulation with our model, alternans were observed at [Na + ] i at~10.2 mM. We tested the effect of high sodium concentration by fixing the sodium level with [Na + ] i = 15 mM. Under this condition, alternans were also observed yet the levels of cytosolic calcium, at both basal and peak, are higher ( Figure 10). with our model, alternans were observed at [Na ]i at ~ 10.2 mM high sodium concentration by fixing the sodium level with [Na condition, alternans were also observed yet the levels of cytosol and peak, are higher ( Figure 10). Simulations were performed with [Na + ]i held at 1 mM, 12 mM increase the diastolic (peak) [Ca 2+ ]nsr and the [Ca 2+ ]myo both increas pected, with increasing [Na + ]i the reverse NCX current increases the cell ( Figure 11C). The forward NCX (Ca 2+ extrusion) also inc [Ca 2+ ]myo. When the NCX rate is increased intracellular [Ca 2+ ] fall trusion ( Figure 11D). Increasing [Na + ]i also attenuates alternans and A5). Simulations were performed with [Na + ] i held at 1 mM, 12 mM, and 20 mM. As [Na + ] i increase the diastolic (peak) [Ca 2+ ] nsr and the [Ca 2+ ] myo both increase ( Figure 11A,B). As expected, with increasing [Na + ] i the reverse NCX current increases bringing more Ca 2+ into the cell ( Figure 11C). The forward NCX (Ca 2+ extrusion) also increases due to the higher [Ca 2+ ] myo . When the NCX rate is increased intracellular [Ca 2+ ] falls due to the increase extrusion ( Figure 11D). Increasing [Na + ] i also attenuates alternans ( Figures A4 and A5).
If reduction of NCX activity by high Na + is crucial for the formation of alternans, NCX upregulation should reduce alternans. Consistent with this hypothesis, high expression levels of NCX reduces SR calcium, along with cytosolic [Ca 2+ ] and thus reduce alternans amplitude ( Figure 12C). Large reduction of NCX (to 50% control) still generates alternans, however, the calcium alternations are attenuated. It again confirms that there is a range of changes that can be critical for the arrhythmogenesis to occur. When the temperature is reduced (23 • C), the alternans become more prominent (Figure 12). In this test case, the simulation was done using a stochastic model for the Na + current developed by [40]. This is in agreement with many other experiments that low-temperature increases the likelihood of alternans development [65][66][67], and we anticipated that Na + current plays an important role here. If reduction of NCX activity by high Na + is crucial for the formation of alternans, N upregulation should reduce alternans. Consistent with this hypothesis, high express levels of NCX reduces SR calcium, along with cytosolic [Ca 2+ ] and thus reduce alterna amplitude ( Figure 12C). Large reduction of NCX (to 50% control) still generates alterna however, the calcium alternations are attenuated. It again confirms that there is a range changes that can be critical for the arrhythmogenesis to occur. When the temperature reduced (23 o C), the alternans become more prominent (Figure 12). In this test case, simulation was done using a stochastic model for the Na + current developed by [40]. T is in agreement with many other experiments that low-temperature increases the lik hood of alternans development [65][66][67], and we anticipated that Na + current plays an i portant role here.

Discussion
Cellular alternans in cardiac myocytes have been shown in experiments and modeling to have a mechanism that depends both on the membrane currents and on the Ca 2+ subsystem. This modeling study demonstrates that for alternans produced at high pacing rates, both mechanism act synergistically to produce alternans. Under the hypothesis that relies on membrane currents, action potential restitution is the underlying cause of cardiac alternans [9][10][11][12][13]. Studies that have suggested that modified intracellular calcium cycling plays a role in occurrence of mechanical and electrical alternans have found that it is possible to have alternans in calcium release without requiring action potential alternans [14][15][16][17][18][19][20][21]. Current computational models are unable to recreate this phenomenon; unless certain

Discussion
Cellular alternans in cardiac myocytes have been shown in experiments and modeling to have a mechanism that depends both on the membrane currents and on the Ca 2+ subsystem. This modeling study demonstrates that for alternans produced at high pacing rates, both mechanism act synergistically to produce alternans. Under the hypothesis that relies on membrane currents, action potential restitution is the underlying cause of cardiac alternans [9][10][11][12][13]. Studies that have suggested that modified intracellular calcium cycling plays a role in occurrence of mechanical and electrical alternans have found that it is possible to have alternans in calcium release without requiring action potential alternans [14][15][16][17][18][19][20][21]. Current computational models are unable to recreate this phenomenon; unless certain modifications to the ionic currents were made [9,23]. Newer studies have suggested that both mechanisms are possible under different conditions [24,25].
Computational studies have suggested that the randomness of Ca 2+ sparks; recruitment of a Ca 2+ spark by neighboring Ca 2+ sparks; and refractoriness of Ca 2+ release units are the important factors required for cardiac alternans [68]. The work here is consistent with this idea. While we take into account the stochastic nature of ionic channel gating at individual release site, the subspace area is treated as a single compartment and thus all ion channels (RyRs and LCCs) in the same calcium release site sense the same level of calcium. The potential spatial distribution of the receptors that may govern the allosteric interaction between them, underlying the long-or short-range correlation mechanism between receptors is approximated using a mean-field approach which governs the non-linearity via a non-linear coupling function as shown in Appendix B.2. One question that can be probed using the model is how Ca 2+ -dependent inactivation of the L-type Ca 2+ channel varies beat to beat during alternans. Figure 6D shows more L-type Ca 2+ during large beats than small beats during alternans. If this is caused by increased Ca 2+ dependent inactivation, it would be expected that there be a faster rate of decay of the L-type current during the large currents. The rate of L-type current decay is shown in the Figure A6 in Appendix A by overlapping the currents at a big and small beat during alternans. Fitting an monoexponential curve (f(x) = a × exp(−x/b)) to the decay yielded a = −9.37 pA and b = 0.21 s for the large beat and a = −55.8 pA and b = 0.46 s for the small beat with a R2 = 0.977 and R2 = 0.980, respectively. The large beat where Ca 2+ is more elevated displays a faster decay than the small beat.
Refractoriness in Ca 2+ release is demonstrated to be crucial to alternans by our statistical analysis of sparks at each release site (Ca 2+ release unit). Our studies also show that an uncoupling of depolarization to Ca 2+ release occurs as the L-type current is reduced both through Ca 2+ -dependent inactivation and smaller action potential amplitude due to reduced Na + current. Our model is not a spatial model and hence does not model neighboring sparks. However, we do see recruitment as elevated bulk myoplasmic Ca 2+ during rapid pacing, in part produced by more sparks, promotes Ca 2+ sparks at additional Ca 2+ release units.
Previous studies have shown that block of NCX by ORM-10962 has been shown to attenuate alternans in experiments [69]. In our model with decreasing NCX activity the NSR Ca 2+ load increases and NCX attenuates alternans. The model suggests that large reduction of NCX (to 50% control) still generates alternans, however, the calcium alternations are attenuated. This suggests that there is a range of conditions that can be critical for the occurrence of alternans. Other studies have shown that blocking late Na + current by ranolazine attenuates alternans, presumably by reduction of reverse NCX although this has not been proven [70,71]. While reduction of NCX, can reduce Ca 2+ entry through reverse-NCX, it will suppress the greater role of Ca 2+ extrusion leading to more calcium overload. The rat ventricular myocyte does not have a late component of the Na + current so we cannot test this directly. In our simulations, increasing intracellular Na + to 12, 15, and 20 mM, the NSR Ca 2+ load increases, yet alternans are attenuated. There is a decrease in the action potential duration with increasing [Na + ] i as shown in the Figures A3-A5 allowing for recovery of the Na + current and production of a regularization of the action potential. It suggests that blocking late Na + current allows the Na + channels to recover to open in the next beat resulting in attenuated alternans.
Calcium sensitive K + channels (SK channels) have been found in rat ventricular myocytes. These channels activate and allow hyperpolarizing outward current when [Ca 2+ ] myo is elevated with a K 0.5 = 0.5 µM [72]. The SK channel conductance in rat peaks at 1-2 pA/pF [73]. This Ca 2+ sensitivity is conferred though the binding of Ca 2+ with calmodulin with experimental studies suggesting calmodulin variants can alter SK channel function and potentially lead to arrhythmia [74]. Furthermore, SK channel expression has been observed to increase during heart failure or after myocardial infarction and have been suggested to contribute to arrhythmia [75,76]. Computational rabbit ventricular myocyte models have suggested that when these channels are blocked, a mild to moderate action potential prolongation occurs depending on the conductance level assigned to the channel [73,76]. However, at the physiological conductance levels, these models predict a very small effect. Additionally, experimental studies have und AP clamp suggested that blocking the K + using Cs under AP clamp, resulted in no change in the difference current between beats in alternans suggesting little role for the SK current alternans. The same study however found that the Ca 2+ -activated Cl − current did play a role [77]. For this reason, we did not include them in the current model and left them for future studies of arrhythmia and disease.
Calmodulin has high and low affinity Ca 2+ binding sites. Experiments have measured high affinity sites in the C-lobe to have a K d = 1 µM and the low affinity sites in the N-lobe to have a K d = 12 µM [78,79]. We have chosen to only include the low affinity N-lobe site because in the microdomain near the L-type channel upon L-type opening that high affinity site will saturate rapidly leaving the low affinity site to play the regulatory role. This is supported by the role the N-lobe has been found to play in Ca 2+ dependent inactivation of the L-type Ca 2+ channel [80]. Furthermore, at high pacing rates where Ca 2+ is elevated, it even more likely that the C-lobe is saturated and plays less of a regulatory role. We did not choose to use the model by Limpitikul et al., which uses a C-lobe K d = 1.15 µM and a N-lobe K d = 0.9 µM as it would not be appropriate for a model which explicitly models the dyadic subspace [81]. According to Uniprot, the D96V, D130G and F142L variants are implicated in LQTS. While located in the C-terminus, molecular dynamics simulations show that they affect the positional relation between the lobes such as the linker distance and dihedral angles between the lobes, so the effect of the variants is not simply an effect on Ca 2+ binding affinity [82]. Including the low affinity sites would not affect the model results presented here. However, to model LQTS all four sites would need to be considered. This is left for future work.
The model is a set non-linear differential equations with stochastic elements. The conversion to alternans is period doubling behavior and has been observed both in experiments and in other models of the cardiac action potential [83][84][85]. For example, with rapid pacing the dog heart can develop alternans. With increasing pacing rate there is period doubling, a repeating sequence of four beat amplitudes. With further increases, fibrillation or chaos occurs. We have also observed this behavior in our previous Guinea pig model and this model (at 12 Hz) [86].
Ca 2+ oscillations have been observed in cardiac ventricular myocytes under Ca 2+ overload conditions. Similar to deterministic systems such as our previous model, this model is capable of Ca 2+ oscillations due to the dynamics of the Ca 2+ subsystem because it is a non-linear excitable system [86]. However, the cardiac ventricular myocyte is a driven system with a periodic applied current to trigger action potentials similar to experiment that mimic the periodic excitation of a ventricular myocytes by adjacent cells during the heartbeat. If the rapid pacing is abruptly ceases, there will be some spontaneous Ca 2+ release events (calcium oscillations) similar to experiment [87]. Furthermore, we have developed a spatial model of the rat ventricular myocyte that displays oscillatory Ca 2+ waves under calcium overload conditions with no depolarization stimulus similar to the phenomena observed in isolated rat cardiac ventricular myocytes [88,89].
Anomalous diffusion has been proposed to explain the spatial diffusion of Ca 2+ . Anomalous diffusion by definition is a diffusion where the mean squared displacement does not depend linearly on time [90]. This is different from the Fickian diffusion that causes Brownian motion which has been observed in cells. Other modeling studies have invoked anomalous diffusion to create Ca 2+ spark models because existing models using Fickian diffusion have only achieved FWHM of 1.0 µm compared to the experimentally observed FWHM of 1.8 µm [91]. With their anomalous diffusion model (subdiffusion) it is possible to get a FWHM of 2.0 µm. However, we have shown previously that with our model that includes Fickian diffusion we can get a FWHM of 1.8 µm [89]. In our modeling efforts there has been no need to invoke novel diffusion formulations to reproduce experimentally observed phenomena.

Conclusions
In conclusion, we presented here a study using a stochastic computational model to study the cellular mechanism underlying cardiac alternans in rat ventricular myocyte. The model helps to explain a modest role of [Ca 2+ ] jSR in forming alternans, while it's suggested that disturbing I Na , I CaL and membrane potential plays a dominant role in the forming of pulsus alternans. In addition to this, the model was able to reproduce results at conditions that have been known for alternans like lowering the temperature, high [Na + ] i or reducing alternans amplitude by up/down regulation of NCX. The limitation of the model is the inability to investigate the spatial effect on the generation of cardiac alternans, i.e., subcellular Ca 2+ alternans or calcium waves in alternans. This will be the next step in our study with a full-scale spatio-temporal model of the cardiac ventricular myocyte to investigate subcellular Ca 2+ alternans [92][93][94].

Conflicts of Interest:
The authors declare that there is no conflict of interest.

Appendix A Appendix A.1. Spark Termination
The termination of Ca 2+ sparks involves the combination of luminal calcium depletion, which reduce the driving force for the Ca 2+ flux; the stochastic closing of the channels and their failure to re-open due to less subspace and jSR Ca 2+ , and the effect of the increasing number of closed channels on the open ones through allosteric coupling, as in our previous model ( Figure 8) [27,58]. Clearly, this predates the concept of pernicious attrition and induction decay in which the smaller Ca 2+ flux at each RyR2 due to local luminal calcium depletion devastates inter-RyR2 CICR [95]. In fact, the idea the reduced Ca 2+ flux through the RyR2s in the dyad are insufficient to sustain CICR as suggested in pernicious attrition has been suggested previously [53,86].

Appendix A.2. Dynamics of Calcium Sparks and Calcium Leak
The model predicts three possible mechanisms of the SR Ca 2+ leak that are sufficient to account for balancing Ca 2+ uptake via SERCA2a during the diastolic phase, similar to our previous work [27]. They are visible Ca 2+ sparks (~120 sparks/cell/s under quiescent), invisible Ca 2+ leaks known as Ca 2+ quarks (~4300 quarks/cell/s under quiescent) and calcium leak via a small population of RyRs at a distance from the release sites, known as non-junctional RyR or 'rogue' RyR ( Figure 8) [96]. The so-called non-spark events or invisible Ca 2+ quarks are the result of one or a few RyR channels opening whose signals cannot be detected using standard confocal imaging microscopies [97].
The visible leak i.e., 120 Ca 2+ sparks per cell per second is in the physiological range (100-200 sparks/cell/s) and accounts for 75% of the leak. Part of the invisible leak results from the stochastic gating of one or a few RyR2 channels, accounting for 25% of the leak. This yields a spark fidelity of~2.8% which is the ability a single calcium release to trigger the Ca 2+ spark. The total quiescent leak rate was~1.2 µmol/(L cyt.)/s which is in the range 0.5-5 µM [98]. The Ca 2+ spark peak was~100.1 µM (Figure A1), which is in agreement with the estimate that due to the small subspace, the elevation of free calcium can be two orders of magnitude larger than the resting value [99]. The plateau in the spark profile of free calcium explains the energy coupling effect between the closed and open channels, before the Ca 2+ spark is completely terminated. Similar to the study of [27], the effect of IP 3 R has not been incorporated in this model study. Thus, the involvement of other leak pathways like via IP 3 -receptors (IP3R) has not been excluded, though the functional role of IP 3 R at rest has not been confirmed. As shown in Figure A1D, the higher frequency of [Ca 2+ /CaM] elevation in the quiescent cell model provided the evidence of non-spark events where calcium release via the opening of one or a few RyRs bind to calmodulin buffer, as well as SR and SL buffers.

Appendix B-Model Equations and Parameters
Calcium Release Site (CRU) The differential equations describing Ca 2+ at the release site are as follows: [ where (i) is an index indicating the i th specific CRU, myo ds ds V V / = λ is the volume fraction that scale the fluxes, defined based on myoplasmic volume, to the subspace volume compartment. The membrane buffers used here are similar to those used previously [58,100].

Appendix B. Model Equations and Parameters
Appendix B.1. Calcium Release Site (CRU) The differential equations describing Ca 2+ at the release site are as follows: where (i) is an index indicating the ith specific CRU, λ ds = V ds /V myo is the volume fraction that scale the fluxes, defined based on myoplasmic volume, to the subspace volume compartment. The membrane buffers used here are similar to those used previously [58,100].

Appendix B.2. Ryanodine Receptor Type-2 Model
The states of the RyR2 are Coupling of the RyR2 channels in the cluster uses allosteric coupling energies: ε OO , ε CC as before [27]. At each release site, the transition between two cluster states Si and Sj is represented as where N O→C (N C→O ) represents the number of channels in cluster state Si (Sj) that can be switched from Open to Closed (Closed to Open). The coupling gating in the cluster is given by the formula with a j is the average allosteric connectivity for each RyR based on mean-field approach [54]. The variables N closed and N open represent the current number of non-conducting, and conducting channels in the cluster state Si, respectively. The non-junctional or 'rogue' RyR is modelled deterministically using the formula

Appendix B.3. L-Type Ca 2+ Channel Model
The 6-state L-type Ca 2+ channel (LCC) model, Figure 2, is derived from 5-state LCC model for rat ventricle myocyte from Sun and co-workers with parameters adjusted for the new spark model [29]. Permeation of the L-type Ca 2+ channel used the Goldman-Hodgkin-Katz (GHK) formalism was used [101] with (i) represents the index of a release site N ds is the calcium concentration in a dyadic subspace), [Ca 2+ ] o is the extracellular calcium concentration; P dhpr is the single channel permeability; z Ca = 2 is the valence of Ca 2+ ion; R is the universal gas constant, T is the temperature and F is the Faraday constant. The partition coefficients are β 0 = 0.341, β 1 = 1 in the case of Ca 2+ channel [102,103]. NOTE: If V m = 0, this formula is used instead Even though the role of extracellular calcium was not investigated in this study, the model also considered the effect of extracellular calcium [Ca 2+ ] o on the model by modulating single channel current and shifting the half-activation voltage V1/2.
It has been suggested that CaM tethers to a region of C-terminal between an EF and IQ domain of LCC [104,105], with 2 different segments (A and C sequences) within the pre-IQ motif of LCC have been indicated to be critical for CDI [106], and this tethering occurs at very low level of calcium concentration (20-100 nM) [35]. This is the range of resting calcium level; thus, it is assumed that all CaM tethers to the channels all the time. So, the total [CaM] is assumed fixed in the subspace, and calcium-dependent inactivation of LCC (k2→4) is modelled via calcium-Calmodulin (CaCaM) complex. In the model, similar to the experimental data, it is assumed that a single bound Ca 2+ /CaM complex is both necessary and sufficient for CDI [107]. Hence, the transition rate for CDI is formulated as k24 = c24 [Ca 2+ /CaM].
Not only Ca 2+ , but also other ions can also permeate via LCC [108,109]. However, due to the large permeability of Ca 2+ compared to other ions (e.g., PCa/PNa > 1000), in this study, only Ca 2+ current is modelled. Due to the nonlinearity in the I-V curve and based on the assumption of independent permeation between ion species and constant-field theory, the Goldman-Hodgkin-Katz (GHK) formalism was used [101].
The equations for the LCC model are given in Table A1. The fraction of LCC channels in each state during a +10 (mV) voltage clamp is shown in Figure 2B. State 2 and 3 are open states; while State 4 is calcium-dependent inactivation and State 5 is voltage-dependent inactivation.

Appendix B.4. Na + Channel Model
The whole-cell Na + current was derived from the formula with the ODE for gating variable follows the first-order differential equation [108] dx Na where x Na represents the dimensionless gating variables, either m, h, or j, in the range between 0 and 1.0 and x ∞ is the steady-state value for the gating variable. The unit for the time constants is in seconds.
The initial values are given in Table A2.
with A m is the SL surface area and V myo is the myoplasmic volume The NCX is given by the formula with I ncx is the maximal NCX current. The parameters are given in Table A3.
In the case of sodium and potassium pumps, we have The three endogenous buffers of calmodulin (CaM), troponin (Trpn), and the phospholipids of the SR membrane (SRbuf) are used for the bulk myoplasm. The buffering parameters are listed in Table A4 I background = I bCa + I bNa + I bK (A56) NOTE: The units: membrane potential and reversal potential are in mV, while other parameters are given in the form of unit-density, e.g., the specific membrane capacitance C sc = 1 µF/cm 2 , the ionic currents are given in (µA/cm 2 ), except the L-type currents are measured at whole-cell level (µA) which needs to be converted to current density by dividing the membrane surface area A m (cm 2 ).

Appendix B.11. Parameters
The rat ventricular cell volume was chosen 25 pL to match to the number 20,000 CRUs per cell with Cm/Vcell = 7.12 pF/pL [99,100]. The schematic diagram of the compartmental model of the rat ventricular myocyte with all the ionic currents is given in Figure A1. Table A1. Parameters for L-type Ca 2+ channel model.