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

Abstract

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, [Ca²⁺]_myo and [Ca²⁺]_SR respectively, increases ryanodine receptor open probability (P_o) 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 [Ca²⁺]_myo and how they contribute to the generation of Ca²⁺ signaling alternans in the heart.

1. 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²⁺ 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].

It has been suggested that a voltage-dependent mechanism underlies cardiac alternans. Under this hypothesis, action potential restitution is the underlying cause of cardiac alternans [9,10,11,12,13]. With the shorter diastolic interval with faster pacing rates, the sarcolemmal ion channels do not completely recover from one beat to the next. Therefore, after a large/long action potential, since ion channels have not fully recovered, they are not available to participate in the next action potential resulting in a small/short action potential. In these studies, a steep restitution curve is the primary requirement, and pathological conditions or experimental manipulations can increase the slope of this restitution curve making alternans more likely.

Other studies have suggested, that modified intracellular calcium cycling plays a role in occurrence of mechanical and electrical alternans [14,15,16,17,18,19]. In fact, many studies have found that it is possible to have alternans in calcium release without requiring action potential alternans [18,20,21]. Calcium alternans is a beat-to-beat variation in intracellular Ca²⁺ transient amplitude. Typically, calcium alternans occurs at high heart rates, yet the frequency threshold varies by different conditions such as ischemia or ionic disturbances that disturb the bidirectional coupling between the membrane potential and intracellular calcium. This has been identified as a potential precursor to the dangerous reentrant arrhythmias and SCD; yet the mechanism is not well understood [22]. 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]. Those studies used computational models to suggest specific model configuration that can produce alternans by either a calcium-dependent or voltage-dependent mechanism. They then used rapid pacing in guinea-pig ventricular myocyte and observed calcium-dependent alternans under conditions of control of action potential alternans. They suggested that in these myocytes, a voltage-driven mechanism was possible due to the steep action potential duration restitution curves.

In this study we use our local-control model of the rat ventricular myocyte to study the interplay of calcium-dependent and voltage-dependent mechanisms to give rise to rapid pacing induced alternans. Through dissection of calcium dynamics and sarcolemmal ion channel activity at the release site level we demonstrate that for alternans brought about by increased pacing rate, there is a synergy between these two mechanisms. The modeling suggests that restitution of the action potential due to incomplete recovery of the Na⁺ channels that sets the stage for alternans by a calcium dependent mechanism. Here the small action potential reduces the number of L-type calcium channels that are activated creating the substrate for stochastic activation of the ryanodine receptors that is favored by elevated myoplasmic calcium.

2. Materials and Methods

2.1. The Model

In order to study the mechanisms of alternans, we developed a novel local-control model for excitation-contraction coupling in the rat ventricular myocytes with 20,000 Ca²⁺ release sites that improves upon our previous work [26,27,28] (Figure 1). Notable improvements of this model over previous work include the following: (1) explicit buffer dynamics in the subspace, (2) an improved L-type calcium channel that incorporate calcium-bound calmodulin dependent inactivation, (3) updated formulation of Ito, and (4) updated parameters based on newer experimental data from rat ventricular myocytes [29]. In addition, by using our patented Ultra-fast Monte Carlo Method and GPU technology, it now allows us to do larger scale simulations that provide insights into calcium dynamics. The model equations are detailed in Appendix B.

2.1.1. 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].

2.1.2. 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²⁺ (${\left[{\mathrm{Ca}}^{2+}\right]}_{ds}^i$) and junctional SR Ca²⁺ (${\left[{\mathrm{Ca}}^{2+}\right]}_{jsr}^i$) and $k_{ryr}^{+}$ being a constant. The RyR2s are arranges in a cluster and display coupled gating as using our previous formulation [27].

$$C\begin{array}{c}\stackrel{k_{ryr}^{+}}{\to }\\ \underset{k_{ryr}^{-}}{\leftarrow }\end{array}O$$

(1)

2.1.3. L-Type Ca²⁺ Channel Model

The 6-state L-type Ca²⁺ 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.

The inactivation of the LCC was modelled via 2 separate pathways: VDI (voltage-dependent inactivation (O2→C5)) and CDI (Ca²⁺-dependent inactivation (O2→C4)). At each release site, the Ca²⁺ 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²⁺ (CaCaM) is the effector for Ca²⁺-dependent inactivation (CDI) of L-type Ca²⁺ 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²⁺ channels (Cav1.2), CDI is triggered by the binding of two Ca²⁺ ions to the C-lobe of CaM [33,34,35]. It means that to model CDI, two Ca²⁺ 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²⁺]_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.

2.1.4. 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² which is in the range 3–25 mS/cm² [39,48].

2.1.5. 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].

2.1.6. Sarcoplasmic Reticulum Ion Pumps

The sarco (endo)plasmic reticulum Ca²⁺-ATPase (SERCA) pump re-sequesters Ca²⁺ 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].

2.1.7. 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).

2.1.8. Membrane Potential

During AP, the dynamics of the membrane voltage V_m are governed by the ionic currents described above.

2.2. Numerical Methods

We used our patented Ultra-fast Monte Carlo Simulation Method to solve for the states of the RyR2 and L-type Ca²⁺ 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.

3. 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²⁺ visible and invisible leak by flux through the RyR2 channels in the forms of Ca²⁺ sparks and non-spark openings, respectively (Figure A1) [27,58,59]. Details of the ionic currents are shown in Figure A2.

3.1. Dynamics of Calcium during a Twitch-Relaxation Cycle

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.

3.2. 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²⁺]_myo and V_m alternating in phase. There are also alternating small and large amplitude and duration of other observable quantities such as [Ca²⁺]_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²⁺]_myo and [Ca²⁺]_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²⁺ 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²⁺]_jsr local depletion for individual CRUs as shows that the CRUs (red, blue, green, and black) activate at alternating beats.

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²⁺ channels assuming the Ca²⁺ 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²⁺ after a release event can contribute to the development of alternans. To assess the SR Ca²⁺ 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²⁺]_jSR at these release sites, right at before the next beat to start, (2) the nadir of [Ca²⁺]_SR during calcium release (Figure 8). Larger SR [Ca²⁺] comes with larger Ca²⁺ release and the beats with larger release have more activated CRUs. We can see that the ensemble initial [Ca²⁺]_SR is not an indicator of alternans formation. If we look at the nadir [Ca²⁺]_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²⁺ 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²⁺ 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²⁺] increases the stochastic activation of Ca²⁺ 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²⁺ 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%.

We then explored the hypotheses that a high level of [Na⁺]_i increases the likelihood of alternans through the modulation of Na⁺/Ca²⁺ 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²⁺]. In other words, increasing Na⁺ lead to more reverse mode NCX, leading to increases [Ca²⁺]_myo and [Ca²⁺]_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).

Simulations were performed with [Na⁺]_i held at 1 mM, 12 mM, and 20 mM. As [Na⁺]_i increase the diastolic (peak) [Ca²⁺]_nsr and the [Ca²⁺]_myo both increase (Figure 11A,B). As expected, with increasing [Na⁺]_i the reverse NCX current increases bringing more Ca²⁺ into the cell (Figure 11C). The forward NCX (Ca²⁺ extrusion) also increases due to the higher [Ca²⁺]_myo. When the NCX rate is increased intracellular [Ca²⁺] falls due to the increase extrusion (Figure 11D). Increasing [Na⁺]_i also attenuates alternans (Figure A4 and Figure 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²⁺] 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.

4. 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²⁺ 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²⁺ sparks; recruitment of a Ca²⁺ spark by neighboring Ca²⁺ sparks; and refractoriness of Ca²⁺ 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²⁺-dependent inactivation of the L-type Ca²⁺ channel varies beat to beat during alternans. Figure 6D shows more L-type Ca²⁺ during large beats than small beats during alternans. If this is caused by increased Ca²⁺ 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²⁺ is more elevated displays a faster decay than the small beat.

Refractoriness in Ca²⁺ release is demonstrated to be crucial to alternans by our statistical analysis of sparks at each release site (Ca²⁺ release unit). Our studies also show that an uncoupling of depolarization to Ca²⁺ release occurs as the L-type current is reduced both through Ca²⁺-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²⁺ during rapid pacing, in part produced by more sparks, promotes Ca²⁺ sparks at additional Ca²⁺ 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²⁺ 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²⁺ entry through reverse-NCX, it will suppress the greater role of Ca²⁺ 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²⁺ load increases, yet alternans are attenuated. There is a decrease in the action potential duration with increasing [Na⁺]_i as shown in the Figure A3, Figure A4 and Figure 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²⁺]_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²⁺ sensitivity is conferred though the binding of Ca²⁺ 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²⁺-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²⁺ 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²⁺ dependent inactivation of the L-type Ca²⁺ channel [80]. Furthermore, at high pacing rates where Ca²⁺ 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²⁺ 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²⁺ oscillations have been observed in cardiac ventricular myocytes under Ca²⁺ overload conditions. Similar to deterministic systems such as our previous model, this model is capable of Ca²⁺ oscillations due to the dynamics of the Ca²⁺ 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²⁺ release events (calcium oscillations) similar to experiment [87]. Furthermore, we have developed a spatial model of the rat ventricular myocyte that displays oscillatory Ca²⁺ 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²⁺. 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²⁺ 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.

5. 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²⁺]_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²⁺ 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²⁺ alternans [92,93,94].