## Abstract

Intracellular calcium (Ca) alternans in cardiac myocytes have been shown in many experimental studies, and the mechanisms remain incompletely understood. We recently developed a “3R theory” that links Ca sparks to whole cell Ca alternans through three critical properties: randomness of Ca sparks; recruitment of a Ca spark by neighboring Ca sparks; and refractoriness of Ca release units. In this study, we used computer simulation of a physiologically detailed mathematical model of a ventricular myocyte couplon network to study how sarcoplasmic reticulum (SR) Ca load and other physiological parameters, such as ryanodine receptor sensitivity, SR uptake rate, Na-Ca exchange strength, and Ca buffer levels affect Ca alternans in the context of 3R theory. We developed a method to calculate the parameters used in the 3R theory (i.e., the primary spark rate and the recruitment rate) from the physiologically detailed Ca cycling model and paced the model periodically to elicit Ca alternans. We show that alternans only occurs for an intermediate range of the SR Ca load, and the underlying mechanism can be explained via its effects on the 3Rs. Furthermore, we show that altering the physiological parameters not only directly changes the 3Rs but also alters the SR Ca load, having an indirect effect on the 3Rs as well. Therefore, our present study links the SR Ca load and other physiological parameters to whole cell Ca alternans through the framework of the 3R theory, providing a general mechanistic understanding of Ca alternans in ventricular myocytes.

- calcium cycling
- calcium spark
- modeling

t-wave alternans have been closely associated with cardiac arrhythmias and sudden cardiac death (3, 50, 51, 66, 69). Intracellular calcium (Ca) alternans, as a potential cause of T-wave alternans (15, 46), have been shown in many experimental studies (2, 5, 7, 16–19, 25, 35, 41, 68, 71). However, the underlying mechanism of Ca alternans remains incompletely understood. Adler et al. (1) proposed a theoretical mechanism linking a nonlinear sarcoplasmic reticulum (SR) Ca release function to mechanical alternans. Evidence from experimental studies (17–19, 71) shows that Ca alternans could be explained by a steep nonlinear dependence of SR Ca release upon the diastolic SR Ca load (Ca_{DSRL}) immediately preceding the release (a steep fractional release-load relationship). This mechanism requires that Ca_{DSRL} alternate concomitantly with SR Ca release, but the refractoriness of SR Ca release is not explicitly required. Subsequent mathematical analyses and simulation studies either potentiated this theory (30, 45, 57, 64, 69, 71) or used it to study subcellular Ca alternans dynamics (21, 56). On the other hand, experimental studies in rabbit ventricular myocytes by Picht et al. (41) and in cat atrial myocytes by Huser et al. (25) showed that Ca alternans can occur without Ca_{DSRL} alternans. In other words, the SR Ca concentration ([Ca]_{SR}) reaches the same level every beat even though the myoplasmic Ca transient amplitude (and thus the SR Ca depletion amplitude) undergoes alternans. Since Ca_{DSRL} is the same in each beat, the refractoriness of SR Ca release is considered as the primary mechanism of alternans (41), and the importance of refractoriness in Ca alternans has been demonstrated in later studies (32, 58). Therefore, the above experimental studies raise the issue of the roles of SR Ca and refractoriness of SR Ca release in the genesis of Ca alternans. In addition, it is well known that the elementary Ca signaling events are Ca sparks (12, 13), which are Ca release events of discrete Ca release units (CRUs) or couplons. A typical ventricular myocyte contains ∼20,000 CRUs (20). Since the openings of the L-type Ca channels (LCCs) and the ryanodine receptors (RyRs) are random, Ca sparks exhibit randomness, even when elicited by action potentials (8, 26). However, most of the previous theories of Ca alternans used the whole cell Ca transient and SR Ca content as the order parameters (or variables) without taking into account the properties of Ca sparks. In other words, the link between Ca sparks and Ca alternans has not been made in these theories.

Ca cycling models taking into account Ca sparks have been developed recently to study Ca alternans in computer simulations (14, 24, 48, 49, 52). In our studies (14, 52), we proposed a theory that links the properties of Ca sparks to whole cell Ca alternans. We called this theory the “3R theory” in which alternans arises as a result of the interactions of three critical properties: randomness of Ca sparks; recruitment of a Ca spark by its neighboring CRUs; and refractoriness of the CRU. A simple mathematical formulism of this theory was derived using a mean-field approach in which we assumed that the system is well mixed so that the firing probability of a CRU is uniform in space. The formulism links the number of Ca sparks in the present beat (*N*_{k+1}), to the number of sparks in the previous beat (*N*_{k}) as:
*f* is a nonlinear function describing the recruitment rate. A specific formulism was derived in our previous studies as (14, 52):

In these equations, *N*_{0} is the total number of CRUs; α is the probability of a primary Ca spark, triggered by the opening of the LCCs or occurring spontaneously in a CRU; γ is the probability that a primary spark recruits a neighboring CRU to spark; and β is the probability that a CRU activated on the previous beat remains refractory during the present beat. *M* is the number of neighbors of a CRU. In *Eq. 1*, *N*_{0} − β*N*_{k} is the number of CRUs recovered from the previous beat, α(*N*_{0} − β*N*_{k}) is the number of primary sparks, and (1 − α)(*N*_{0} − β*N*_{k})*f* is the number of recruited sparks. Ca alternans occurs when β and γ are large and α is in the intermediate range (Fig. 1). Note that in this theory, Ca_{DSRL} is not an explicit parameter required for Ca alternans. We also demonstrated in computer simulation that although in the case of free running SR, Ca alternans is accompanied by Ca_{DSRL} alternans, when [Ca]_{SR} was clamped at a constant value Ca alternans could still occur (52). This indicates that Ca_{DSRL} alternans is not required for Ca alternans, providing a possible explanation for the observations by Picht et al. (41) and Huser et al. (25).

Even though Ca_{DSRL} alternans is not required for Ca alternans, this does not necessarily mean that Ca_{DSRL} is not important for the genesis of Ca alternans. While Ca_{DSRL} is not an explicit parameter in the 3R theory, it affects α, β, and γ by regulating the RyR open probability, refractoriness, and the amount of Ca released when an RyR channel opens. Therefore, there is no doubt that SR Ca load is important for Ca alternans; the question remaining is how it affects Ca alternans. Furthermore, there are many detailed physiological parameters [e.g., the LCC open probability, the Na-Ca exchange current (*I*_{NCX}), the RyR open probability, the extracellular Ca concentration ([Ca]_{o}), the sarcoplasmic/endoplasmic reticulum Ca ATPase (SERCA) pump, the CRU spacing, Ca buffers, etc.] that are altered in diseases and experiments. These parameters not only directly affect α, β, and γ but also affect Ca_{DSRL} and thus indirectly affect α, β, and γ. For example, a stronger SERCA pump rate causes a higher Ca_{DSRL}, but this also removes Ca from the myoplasm faster so that less Ca can reach the neighboring CRUs, which reduces the recruitment rate γ. Therefore, the effects of these parameters on Ca alternans may not be straightforward but rather complex, which may be responsible for the seemingly contradictory observations in many experimental studies. For example, partially reducing LCCs promoted alternans in some experiments (35, 38) but suppressed them in other experiments (37). In mouse models of catecholaminergic polymorphic ventricular tachycardia, alternans were both suppressed and potentiated by leaky RyRs (32, 34, 53). Since in cardiac diseases, such as heart failure, many parameters are remodeled, a general theory that links these parameter changes to alternans is important for understanding the underlying mechanisms. The 3R theory provides such a theoretical framework in which one can link a physiological parameter and SR Ca load to α, β, and γ in the 3R theory and thus understand how the physiological parameters affect Ca alternans.

The purpose of the present study is twofold: *1*) to investigate the roles of SR Ca load on Ca alternans and understand the mechanisms using the 3R theory; and *2*) to link the detailed physiological parameters to the parameters (α, β, and γ) in the 3R theory and thus understand the underlying mechanisms by which they affect Ca alternans. We used a two-CRU system to calculate how SR Ca load and other physiological parameters affect α and γ. We then carried out simulations of the full model to study how SR Ca load and the detailed physiological parameters affect Ca alternans. By comparing their effects on α and γ (and thus the effects on alternans in the 3R theory) to their effects on Ca alternans in the detailed model, we link the SR Ca load and other physiological parameters to whole cell Ca alternans through the framework of the 3R theory, providing a general mechanistic understanding of Ca alternans in ventricular myocytes.

## METHODS

### Spatially Distributed Ca Cycling Model

We improved our previous two-dimensional CRU network model (52) into a three-dimensional (3D) model that includes 100 × 20 × 10 CRUs (Fig. 2), representing a ventricular myocyte (40). Each CRU contains (Fig. 2*A*) a junctional SR (jSR) that is diffusively connected to the network SR (SR) and a dyadic space (DS) that is diffusively connected to the myoplasm (Myo). Extracellular Ca enters the DS through the voltage-gated LCCs, which open stochastically and are simulated by a Markov model (Fig. 2*B*) (39). Ca is released from the jSR through its associated cluster of RyRs to the DS. The RyRs also open stochastically and are simulated using a Markov model (Fig. 2*C*) (61) in which activation and inactivation of RyRs are regulated by Ca in the DS. Ca is either extruded from the cell via the Na-Ca exchanger or taken back up into the SR via the SERCA pump. The CRUs are locally coupled by Ca diffusion throughout the Myo and SR domains, and a full myocyte model is simulated as a 3D grid of 100 × 20 × 10 CRUs (Fig. 2*D*). Each CRU contains 100 RyRs and 10 LCCs.

The differential equations describing the Ca concentrations in different spaces are as follows:
*c*_{m}(*x*,*y*,*z*,*t*) and *c*_{s}(*x*,*y*,*z*,*t*) are the local Ca concentrations in the Myo and the SR space, respectively, and *c*_{d}^{(i)} and *c*_{j}^{(i)} are the Ca concentrations in the *i*^{th} DS and jSR space, respectively. The β_{m}, β_{s}, β_{d}, and β_{j} are Ca buffering functions, which are instantaneous functions of the Ca concentration in each space, following Wagner and Keizer (67). The *J*_{m}, *J*_{s}, *J*_{d}^{(i)}, and *J*_{j}^{(i)} are the net Ca flux for each space. *D*_{m} and *D*_{s} are the Ca diffusion constants in the Myo space and the SR space, respectively, and we set them as *D*_{m} = *D*_{s} = 0.3 μm^{2}/ms. The detailed mathematical formulations of the buffering constants and the Ca fluxes are based on the two-dimensional model by Rovetti et al. (52) and were described in more detail for the 3D model in our recent publication (40).

### Measuring α and γ

Since the parameter β is largely controlled by the pacing cycle length (PCL), we set out to measure the primary spark probability (α) and the spark recruitment probability (γ) directly during the simulations of the detailed model. To measure α and γ for a given set of physiological parameters, we used the following protocol. We coupled two CRUs, say *CRU #1* and *CRU #2* (Fig. 2*E*). We removed all the LCCs from *CRU #2* so that the sparks occurring in *CRU #1* are primary ones only and the sparks occurring in *CRU #2* are recruited ones only. For a given set of parameter values, we paced *CRU #1* at PCL = 1,000 ms, which is much larger than the time required for the RyRs to recover from the previous beat. This allows us to isolate the effect of recruitment without any impact from refractoriness. After pacing the system to the steady state, data were collected for 100 beats. We then measured α and γ directly as follows: α = number of sparks in *CRU #1*/total number of stimuli, and γ = number of sparks in *CRU #2*/number of sparks in *CRU #1*.

### Pacing Protocol

The model is paced periodically with a clamped voltage waveform (Fig. 2*F*) in which the voltage jumps from −80 to 0 mV, maintains at 0 mV for 150 ms, and then decreases linearly to −80 mV in 50 ms. This waveform is more similar to an action potential than a rectangular wave form used in voltage clamp studies. The PCL is the time interval from the beginning of the previous voltage jump to the beginning of the following voltage jump.

### Simulation Methods

The SR and Myo domains of a CRU are discretized into 3D spatial grids of 5 × 5 × 5 grid points. The distance between two grid points in the Myo is 0.2 μm with the one for SR space scaling with the corresponding volume ratio (40). A ventricular myocyte model of 100 × 20 × 10 CRUs contains two sets of 500 × 100 × 50 grid points. The two domains are coupled through the collection of jSR and DS domains at a spacing of 1.0 μm. The equations are simulated using an operator splitting method (44) by advancing first the diffusion step and then the flux steps using a first-order forward Euler method with a time step of 0.01 ms. The stochastic transitions of LCCs and RyRs are simulated using the Gillespie's stochastic simulation algorithm (22), modified to handle time-dependent rates (40). All computations were performed on an Intel Xeon 2.53 GHz processor using graphical processing unit parallel computing with an NVIDIA Tesla C2050. Typical runtime is 10 min for a 1-s simulation in the full model of 100 × 20 × 10 CRUs. However, it is still computationally nontrivial for collecting large data sets, such as the bifurcation diagrams used in this study. We found that reducing the model to 20 × 20 × 5 CRUs still produced nearly identical results to the full 100 × 20 × 10 CRU model (Fig. 3; also see Fig. 5*A* for the comparison of the onset of bifurcation). Therefore, we used the reduced model to produce results that required lengthy simulations in this study.

## RESULTS

### Cycle Length Dependence of Ca Alternans

To determine how PCL affects Ca alternans in our spatially distributed Ca cycling model, we paced the model with different PCLs (Fig. 4). In Fig. 4, *A* and *B*, we plotted bifurcation diagrams with respect to the PCL for the peak average Ca concentration ([Ca]_{i}) in the myoplasmic space and Ca_{DSRL}, respectively. For each PCL, the total time simulated was 10 s and the last 10 beats were plotted in the graphs. Ca alternans occurs at around PCL = 300 ms and persists until PCL = 200 ms, which is the shortest PCL that can be used in this model (since the clamped duration is 200 ms). In Fig. 4, *C* and *D*, we plotted the [Ca]_{i} (solid line) and [Ca]_{SR} (dashed line) for two different PCLs, PCL = 400 ms (Fig. 4*C*), and PCL = 250 ms (Fig. 4*D*), along with snapshots of the Ca concentration distribution in the myoplasmic space at the time points marked. We observed that during alternans, local Ca waves occurred during the large Ca transient but rarely during the small Ca transient (see Supplemental Online Movies; Supplemental Material for this article is available online at the *Am J Physiol Heart Circ Physiol* website), in agreement with the experimental observations of Diaz et al. (17, 18).

Although it is well known that Ca alternans occurs as PCL decreases, the exact reasons are still not well understood. Based on the 3R theory, alternans can only occur when β is large. Decreasing PCL increases β, since at faster pacing more CRUs will be in their refractory period following their previous firing beat. In addition, as PCL decreases, the SR Ca load increases (Fig. 4*B*), which may increase α and γ to potentiate Ca alternans (see results below).

### Effects of SR Load on Ca Alternans

To study the effects of SR Ca load on Ca alternans, we set the extracellular Ca, [Ca]_{o}, at various levels and paced the model at PCL = 250 ms (Fig. 5). We found that Ca alternans occurred at [Ca]_{o} = 2 mM and terminated at [Ca]_{o} = 8 mM (Fig. 5*A*), corresponding to SR Ca loads roughly equal to 800 and 950 μM, respectively (Fig. 5*B*). This shows that at both low and high Ca load, Ca alternans cannot occur. However, the reason why alternans only occurs for intermediate values of SR Ca load is not clear and we explore this below.

Changing [Ca]_{o} does not isolate the effects of SR Ca load because ([Ca]_{o}) also affects other fluxes in the model such as the NCX and LCC fluxes. To isolate the role of SR Ca load, we clamped [Ca]_{SR} at fixed values in the model. We first clamped [Ca]_{SR} at various fixed values and measured α and γ in the two-CRU system (Fig. 2*E*). The parameter γ increases steeply from 0 to 1 over a narrow range of [Ca]_{SR} values, and α increases from 0 to 1 in a similar but wider range of [Ca]_{SR} values (Fig. 5*C*). Next, we clamped [Ca]_{SR} at various fixed values in the whole model and paced it at PCL = 250 ms (Fig. 5*D*). Ca alternans begins between 700 and 800 μM and terminates between 1,300 and 1,400 μM. Based on the 3R theory, alternans occurs when γ is large and α is at the intermediate values (Fig. 1). At low SR Ca load, both α and γ are small; thus our simulation results agree with the prediction of the 3R theory that no Ca alternans can occur. At a higher SR Ca load, both α and γ reach one; thus our simulation results also agree with the 3R theory prediction that no alternans can occur since α is too large. Therefore, alternans only occurs in the intermediate values of SR Ca load.

Note that the SR Ca load levels for which α is in the intermediate range (as measured from the two-CRU system) does not match the SR Ca load levels for which alternans occurs in the whole CRU network system. There are many possible reasons for this difference. First, the CRUs in the two-CRU system face a different environment from the ones in the full model (e.g., both CRUs in the two-CRU system are at the boundary and one of the CRUs has no LCCs). Second, α and γ in Fig. 5*C* were measured at PCL = 1,000 ms to avoid the confluent effects of refractoriness, while the bifurcation in Fig. 5*D* was obtained at PCL = 250 ms. If we also measure α and γ at PCL = 250 ms, they do shift to higher SR Ca values (Fig. 6).

### Effects of RyR Sensitivity on Ca Alternans

To increase RyR sensitivity, we increased the rate constant (k_{ap}) for the transition from the close state “C” to the open state “O” in the Markov RyR model (Fig. 2*C*) by a factor σ_{RyR}:k_{ap}*c*_{d}^{2}→σ_{RyR}k_{ap}*c*_{d}^{2} (we also increased the rate constant for the transition from the refractory state “R” to the inactive state “I” by the same factor to preserve detailed balance). As in Fig. 5*C*, Fig. 7*A* shows α and γ measured from the two-CRU system at different clamped [Ca]_{SR} values, but with k_{ap} increased by 40% (σ_{RyR} = 1.4). Both curves shifted to the left compared with the ones shown in Fig. 5*C*. Due to this shift, one would expect Ca alternans to begin and terminate at lower SR Ca loads than for the control value of RyR sensitivity. To test this hypothesis in the whole model, we clamped the SR load at various fixed values and paced the cell at PCL = 250 ms (Fig. 7*B*). Ca alternans begins ∼600 μM and terminates between 900 and 1,000 μM, which are lower than the values obtained with the control RyR sensitivity (Fig. 5*D*). Therefore, increased RyR sensitivity may cause Ca alternans to occur at a lower SR Ca load and thus promote Ca alternans. Note that in this study, by promoting or suppressing alternans we generally mean that changing a parameter causes an earlier onset (or a lower threshold) or termination (or a higher threshold) of alternans and are not referring to a change in the amplitude of alternans.

To study how RyR sensitivity affects Ca alternans under free SR, we set [Ca]_{o} = 2.5 mM and paced the couplon network model at PCL = 250 ms. Figure 7, *C* and *D*, shows peak [Ca]_{i} and Ca_{DSRL} vs. σ_{RyR} when the SR Ca was set free in our simulations. As the RyR sensitivity was increased from control, the magnitude of Ca alternans reduces until the RyR sensitivity was increased by >50%, at which point Ca alternans disappears. Increased RyR sensitivity caused alternans at lower Ca_{DSRL}, (at control RyR sensitivity, alternans occurs when Ca_{DSRL} is above 800 μM; Fig. 5*B*), agreeing with the prediction of 3R theory that increasing RyR sensitivity promotes alternans by causing alternans to occur at lower SR Ca loads. However, as RyR sensitivity is increased, α eventually reaches one, at which point Ca alternans terminates. On the other hand, if one reduces RyR sensitivity from its control value, the magnitude of Ca alternans increases, but the total Ca release becomes more and more irregular (Fig. 7, *E* and *F*). The cause of this irregularity can be understood as follows. As RyR sensitivity reduces, the incoming Ca flux through the LCCs is no longer able to trigger the opening of RyRs with high probability (i.e., low α). Therefore, only a few sparks occur on a paced beat causing the SR to load to a higher level. At some point the SR loads to a level high enough such that the incoming Ca flux through the LCCs is able to trigger the opening of RyRs with high probability due to the effect of SR Ca load on α. Therefore, a large Ca release occurs, depleting the SR and starting the process anew. This irregularity may agree with a recent experimental observation that partial inhibition of SR Ca release elevated SR Ca, favoring beat-to-beat irregular responses (17, 37).

### Effects of SERCA Pump on Ca Alternans

Overexpressing SERCA2a, which enhances Ca reuptake through the SERCA pump, has been shown to suppress Ca alternans in cardiac myocytes by causing alternans to occur at faster pacing rates (16, 71). To understand how the accelerated SERCA pump suppresses Ca alternans in the framework of 3R theory, we first measured α and γ in the two-CRU system with fixed [Ca]_{SR} for various strengths of the SERCA pump multiplied by a factor σ_{SERCA}: *J*_{SERCA}→σ_{SERCA}*J*_{SERCA}. We found a steep decrease in γ from 1 to 0 as σ_{SERCA} was increased from half its control value (σ_{SERCA} = 1) to double its control value, while α remained nearly unchanged (Fig. 8*A*). The reason that increasing SERCA pump strength reduces γ is that a stronger SERCA pump enhances local Ca reuptake and thus reduces the amount of Ca diffusing to neighboring CRUs, therefore, weakening the coupling of the CRUs. Based on the 3R theory, increasing SERCA pump strength will suppress alternans. To demonstrate this, we performed simulations in the couplon network model with clamped [Ca]_{SR} at 900 μM. Ca alternans were terminated when the strength of the SERCA pump was increased to 10 times its control value (Fig. 8*B*). Note that γ decreased to almost 0 when σ_{SERCA} was doubled, which indicates that Ca alternans should be terminated at a much lower SERCA pump strength than the one in the simulation in Fig. 8*B*. As we noted above, this may be due to that the conditions in two-CRU system are different from the full system.

On the other hand, enhancing the SERCA pump also increases SR Ca load, which increases α and γ as shown in Fig. 5. Therefore, enhancing SERCA pump results in complex competing effects on Ca alternans: *1*) it reduces γ to suppress alternans; *2*) it increases SR Ca load, which increases γ to potentiate alternans; and *3*) it can also suppress alternans by increasing SR Ca load to increase α to large values. To study the effects of the SERCA pump on Ca alternans in the couplon network model with an unclamped SR load, we performed simulations with different values of σ_{SERCA} for [Ca]_{o} = 5 mM and PCL = 250 ms (Fig. 8, *C* and *D*). At the control SERCA pump strength (σ_{SERCA} = 1), Ca alternans occurs. Ca alternans disappears at both large and small values of σ_{SERCA}. However, compared with the control case in Fig. 5*B* for which Ca alternans began at a Ca_{DSRL} of 800 μM, alternans were maintained at much lower Ca_{DSRL} values when the SERCA pump strength was lower than control (Fig. 8*D*). This can be understood as follows. Although reducing SERCA lowers Ca_{DSRL}, which tends to reduce γ, reducing SERCA pump also directly increases γ. The two competing effects may result in a net effect of increasing γ to potentiate alternans. As the SERCA pump strength is reduced further, both α and γ may be reduced due to further lowering of Ca_{DSRL}, which eventually eliminates alternans. On the other hand, when SERCA pump strength is increased, it may reduce γ to suppress alternans, and it may also increase α (due to increase of SR load) to suppress alternans. Therefore, based on our results shown in Fig. 8, *C* and *D*, enhancing the SERCA pump can either suppress or potentiate alternans depending on the specific state of the system.

### Effects of Na-Ca Exchange on Ca Alternans

The effects of Na-Ca exchange on Ca alternans are more straightforward. When we fixed SR Ca load and measured α and γ in the two-CRU system with different strengths of *J*_{NCX}: *J*_{NCX}→σ_{NCX}*J*_{NCX}, there was almost no change in α, while γ decreased moderately as σ_{NCX} increased (Fig. 9*A*). This suggests that the main role of NCX in Ca alternans is through its regulation of the average SR Ca load. To test this hypothesis, we unclamped the SR Ca load. As seen in Fig. 9*A*, increasing σ_{NCX} with an unclamped SR load caused a sharp decrease in γ without much effect on α. In the couplon network model with free SR, Ca alternans occurred when we set [Ca]_{o} = 4 mM and paced the system at PCL = 250 ms. Setting σ_{NCX} = 2 eliminated these alternans (Fig. 9*B*). Therefore, increasing the strength of Na-Ca exchange suppresses Ca alternans through its synergistic effects on SR Ca load and γ.

### Effects of Ca Buffers on Ca Alternans

We also studied the effects of Ca buffers on Ca alternans. We found that with a fixed SR load, increasing the concentration of the SR-bound buffer in the myoplasm, B_{SR}→σ_{Buff}B_{SR}, caused a large decrease in γ and only a small decrease in α (Fig. 9*C*). This is due to the fact that Ca released from a primary spark is buffered in the myoplasmic space before it can cause its neighbors to fire, causing γ to decrease. However, the myoplasmic buffering of Ca only has a small effect on Ca diffusion out of the dyadic space and thus has little effect on α. We were able to induce alternans at [Ca]_{o} = 1.8 mM and PCL = 250 ms by setting σ_{Buff} = 0.5 in the couplon network model (Fig. 9*D*), even though the SR load was reduced. In other words, reducing myoplasmic Ca buffer promotes Ca alternans mainly via its direct effect on increasing γ.

## DISCUSSION

In previous studies (14, 52), we developed a theoretical framework for Ca alternans, called the 3R theory, which links the Ca spark properties to whole cell Ca alternans in cardiac myocytes. The present study links the detailed cellular and subcellular physiological parameters to whole cell Ca alternans via the 3R theory by comparing their effects on the parameters in the 3R theory to their effects on Ca alternans in the spatially distributed Ca cycling model. Our observations and conclusions can be summarized as follows (illustrated in Fig. 10).

### Effects of SR Ca Load

The SR Ca load affects Ca alternans by directly affecting α, β, and γ. Specifically, higher SR Ca load facilitates Ca-induced Ca release (CICR), making it easier for LCC to trigger Ca sparks and thus increases the primary spark probability α. For the same reason, it also increases the recruitment parameter γ. Our observation that Ca alternans only occurs in an intermediate range of SR Ca loads for a fixed PCL can be well explained by the effects of SR Ca load on α and γ in the 3R theory (Fig. 1), i.e., for low SR Ca load, both α and γ are too small for alternans to occur, while for high SR Ca load, α is too large. Note that in our RyR model, the RyR open probability is not directly regulated by SR Ca. The reason that SR Ca affects α and γ is that once one or more of the RyRs open, the Ca flux through these channels is larger with a higher SR Ca concentration, facilitating CICR and thus increasing α and γ. However, many experimental studies (23, 29, 43, 65) have shown that the RyR open probability is also affected by luminal Ca, i.e., for a higher SR Ca load, the RyR open probability is higher. This should further enhance the effects of SR Ca load on α and γ.

### Effects of Heart Rate

The most direct effect of heart rate on Ca alternans is its effect on the refractory parameter β, since shortening PCL increases the probability (β) of a CRU remaining in the refractory state after its previous firing beat. Another effect of increased heart rate is that the LCC open probability is reduced due to the incomplete recovery of LCCs. This reduces α to affect Ca alternans. Heart rate can also affect Ca alternans by altering SR Ca load. It has been shown (6) that in large animals, such as rabbit, SR Ca load increases as heart rate increases, while in small animals, such as mouse and rat, SR Ca load decreases as heart rate increases. Therefore, heart rate affects Ca alternans through all the three parameters of the 3R theory, but its main effect is through β.

### Effects of the Cellular and Subcellular Physiological Parameters

The detailed cellular and subcellular physiological parameters, such as the RyR properties (open probability, refractoriness, etc.), the SERCA pump rate, the RyR cluster distribution, and the LCC-RyR interaction, etc., directly affect α, β, and γ, as well as indirectly by affecting the SR Ca load (Fig. 10). The direct effects may work synergistically with the indirect effects or competitively to promote or suppress Ca alternans. Next, we discuss the effects of some specific parameters.

#### RyR sensitivity to myoplasmic Ca.

Increased RyR sensitivity by redox modification has been shown to increase susceptibility to alternans (5). Our study shows that increasing the RyR sensitivity to myoplasmic Ca increases both α and γ when SR Ca is fixed. On the other hand, an increase in RyR sensitivity causes the SR to be more leaky through small numbers of RyR channels opening randomly at a higher rate, as well as more spontaneous Ca sparks. Thus the average SR load is decreased, which causes a decrease in α and γ. Therefore, the net effect of increasing RyR sensitivity is complex. As we show in this study, alternans can be suppressed either by increasing or decreasing RyR sensitivity to myoplasmic Ca, depending on the particular state of the system. This may explain the contradictory observations in mouse models of catecholaminergic polymorphic ventricular tachycardia that alternans was either suppressed or potentiated by leaky RyRs in these animal models (32, 34, 53).

#### SERCA pump rate.

Altering SERCA pump rate mainly affects γ (when the SR Ca is fixed) because as the strength of SERCA is increased, Ca released by a primary spark is removed from the myoplasm more quickly; thus, less Ca is diffused to the neighboring CRUs, decreasing the recruitment rate γ. On the other hand, increased SERCA also causes a higher SR Ca load, which then increases both α and γ. These competing effects may either promote or suppress alternans depending on the state of the system. We observe that increased SERCA from the control value in our model suppressed alternans, agreeing with experimental observation that overexpressing SERCA2a suppresses Ca alternans (16, 71) and reducing SERCA2a enhances Ca alternans (68). However, we also show that severely attenuating the SERCA pump may also abolish Ca alternans due to reduced SR Ca load. A recent study (62) showed that Ca waves were suppressed in myocytes of SERCA2a knockout mice, indicating that the recruitment (and thus γ) is reduced by the diminished SERCA2a abundance.

#### Ca buffering.

We show that increased Ca buffering in the myoplasm mainly decreased γ when SR Ca is fixed. This is due to the fact that more of the Ca released from a primary spark is buffered and thus less Ca diffuses to its neighbors, causing γ to decrease. Buffering has little effect on α since Ca dynamics in the dyadic space are extremely fast and the Ca concentration is much higher than that in the myoplasmic space, increasing myoplasmic Ca buffering only slightly affects Ca diffusion from the dyadic space to myoplasmic space, and thus has only a small effect on CICR in the dyadic space. However, increased Ca buffering makes a spontaneous spark less likely to cause its neighbors to fire, reducing SR Ca leak, and thus causes a higher SR Ca load that increases γ. Our present study shows that a decrease in the buffering concentration promotes alternans, even though the average SR load is lowered, indicating the direct effect of buffering on γ dominates. It should be kept in mind that there are many buffers in the myoplasmic space with different binding and unbinding kinetics (55), which may affect Ca alternans differently from the buffers in our model, which are instantaneous. In addition, SR Ca buffers may also affect Ca alternans in potentially different ways than the effects of Ca buffers in the myoplasmic space. We expect that these different effects can still be explained via the framework of the 3R theory.

#### NCX and LCCs.

*I*_{NCX} and *I*_{Ca,L} are the two sarcolemmal ionic currents that directly affect Ca cycling and thus alternans. The effect of *I*_{NCX} on Ca alternans is straightforward, i.e., increasing NCX reduced the Ca load of the cell, therefore, mainly affecting alternans via SR Ca load. The effect of LCCs is more complex. In this study, although we did not explore the parameter space of LCCs, we needed to set the LCC conductance in a proper range in order for alternans to occur in our model. In the previous study (52), we showed in simulation that Ca alternans disappeared when the LCC conductance is either high or low. We attributed this mainly to the effect of the LCC open probability or conductance on the primary spark probability α. However, increasing (or decreasing) the LCC open probability or conductance increases (or decreases) Ca load of the cell, which will also affect γ. Therefore, how the LCC open probability affects alternans is not straightforward, i.e., increasing the open probability can either suppress or promote alternans. These theoretical considerations may explain why alternans was induced in some experiments (35, 38) but suppressed in other experiments (37) by partially reducing LCCs.

#### Other parameters.

Many other physiological parameters or processes at the subcellular level affect Ca alternans, such as RyR cluster distribution, the LCC-RyR interaction (or the volume of DS), the myofilament sensitivity to Ca, the T-tubule distribution, etc. Although we did not carry out simulations to study their effects in this study, we believe that their effects can still be explained via the 3R theory: their direct effects on α, β, and γ, and their indirect effects through SR Ca load. For example, reduced CRU spacing, which occurs in failing myocytes (10), increases γ, which may contribute to the promotion of Ca alternans and Ca waves (27) in failing hearts. Increasing the Ca diffusion rate also increases γ, which has been shown in computer simulations to promote Ca alternans (49). These physiological parameters, similar to the parameters discussed above (e.g., RyR sensitivity and SERCA pump rate), may affect alternans in a complex manner, but the underlying mechanisms can be understood through the framework of the 3R theory.

In the present study, we only directly measured α and γ, not the refractoriness parameter β. In our RyR model (Fig. 2*C*), the refractoriness is mainly determined by the transition rate from state “R” to state “C,” which was fixed in our study. However, as shown in experimental studies, the refractoriness is also affected by many subcellular physiological parameters. For example, increasing the sensitivity of RyR to myoplasmic Ca appears to speed the recovery of Ca spark triggering (47). Increased SERCA pump function may accelerate the recovery from refractoriness, as shown experimentally at both the cellular (63) and the spark (47) levels. In addition, myoplasmic Ca buffers may slow refilling of SR Ca and thereby slow recovery from refractoriness (54). The effects of altering the physiological parameters on the refractoriness β need to be considered together with their effects on α and γ to gain a holistic mechanistic understanding of Ca alternans through the 3R theory.

### Summary and Conclusions

In addition to SR Ca load and heart rate, all Ca-related physiological parameters and regulating processes may affect Ca cycling dynamics and Ca alternans. As illustrated in Fig. 10, even though the effects of these parameters in Ca alternans may be very complex, they have direct effects on α, β, and γ and indirect effects by altering SR Ca load. In other words, by analyzing their direct and indirect effects on α, β, and γ, one can understand how these parameters affect Ca alternans via the 3R theory. Since changing any of the physiological parameters may alter SR Ca load, SR Ca load plays an extremely important role in the genesis of Ca alternans. As shown in this study, varying a specific physiological parameter may either promote Ca alternans or suppress Ca alternans, and the 3R theory may serve as a unifying theory for explaining seemingly different experimental observations that perturb the same physiological parameters to study Ca alternans. In conclusion, our present study shows that the complex effects of the SR Ca load and other physiological parameters on Ca alternans can be understood through the framework of the 3R theory, which provides a general mechanistic understanding of Ca alternans in ventricular myocytes.

### Limitations

There are several limitations in this study, which are discussed below.

First, since the 3R theory was derived under a mean-field approximation, we do not expect it to be quantitatively accurate for predicting Ca alternans in the physiologically detailed Ca cycling model. In addition, the calculated α and γ from the two-CRU system may not be the same as the true values of α and γ in the full system, another limitation that prevents us from accurate predictions. However, we believe that the two-CRU system still gives qualitatively correct trends in a parameter's effects on α and γ. The 3R theory can only take into account Ca alternans, and is not able to predict the irregular beat-to-beat responses seen in experimental studies (33, 37) and in computer simulations of the detailed model for severely reduced RyR sensitivity (Fig. 7 in the present study) or LCC conductance [Fig. 7*A* in Rovetti et al. (52)]. It should be noted that the 3R theory is not likely the only mechanism of intracellular Ca alternans; other mechanisms are expected to exist and need to be revealed in future studies.

Second, we used a RyR model (61) that includes myoplasmic Ca-dependent inactivation, which may not be what is happening in the real cell (36). In addition, no luminal Ca regulation of RyR is present in this model, the existence of which has been demonstrated in experimental studies (23, 29, 43, 65). However, the regulation of RyRs is still an issue in debate, which is summarized in a recent review article by Sobie and Lederer (60). Nevertheless, many RyR models (11, 28, 31, 49, 59, 61, 70, 72) have been developed and used to study Ca cycling dynamics in cardiac myocytes. We expect that a specific RyR model may give specific predictions of the effects of a physiological parameter on Ca alternans, but we also believe that these specific effects can be mechanistically unified using the 3R theory.

Third, our detailed Ca cycling model is a uniform CRU network with the CRU spacing in both longitudinal and transverse direction being the same (1 μm). However, we also assumed that the Ca diffusion constants in both directions are identical. Due to the complex structure of cardiac myocytes, the diffusion rate in the longitudinal direction may be faster than in the transverse direction (42), and anisotropic Ca diffusion in myoplasm has been used in previous models (27, 49). Taking into account the faster diffusion in the longitudinal direction, our model is equivalent to a longer CRU spacing in the longitudinal direction. In addition, the RyR cluster distribution inside a cell is heterogeneous (4, 9), and new mechanisms of Ca alternans may emerge due to the heterogeneous firing properties of the CRUs, which need to be revealed in future studies.

Fourth, our detailed Ca cycling model is a generic model, which does not take into account species-specific properties. To our knowledge, no species-specific spatially distributed Ca cycling model has been developed so far for ventricular myocytes, and the experimental information available for distinguishing the species-specific differences is still limited. In a recent publication (40), we coupled the rabbit ventricular action potential model (39) with the spatially distributed Ca cycling model, which, with further improvements, may be used to study Ca cycling dynamics in rabbit ventricular myocytes. Nevertheless, the goal of using animal models is to study human diseases, such as alternans and arrhythmias, and reveal the underlying common mechanisms, which may eventually lead to novel therapeutic strategies. In that light, a generic model is an advantage in terms of understanding common mechanisms.

Despite the limitations, our present study links the SR Ca load and other cellular and subcellular physiological parameters to whole cell Ca alternans through the framework of the 3R theory, providing a general mechanistic understanding of Ca alternans in ventricular myocytes.

## GRANTS

This work is supported by National Heart, Lung, and Blood Institute Grants P01-HL-078931 and R01-HL-093205.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## AUTHOR CONTRIBUTIONS

Author contributions: M.N. and Z.Q. conception and design of research; M.N. performed experiments; M.N. and Z.Q. analyzed data; M.N. and Z.Q. interpreted results of experiments; M.N. prepared figures; M.N. and Z.Q. drafted manuscript; M.N. and Z.Q. edited and revised manuscript; M.N. and Z.Q. approved final version of manuscript.

- Copyright © 2012 the American Physiological Society

## REFERENCES

- 1.↵
- 2.↵
- 3.↵
- 4.↵
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 64.↵
- 65.↵
- 66.↵
- 67.↵
- 68.↵
- 69.↵
- 70.↵
- 71.↵
- 72.↵