Calcium buffers and L-type calcium channels as modulators of cardiac subcellular alternans

In cardiac myocytes, calcium cycling links the dynamics of the membrane potential to the activation of the contractile filaments. Perturbations of the calcium signalling toolkit have been demonstrated to disrupt this connection and lead to numerous pathologies including cardiac alternans. This rhythm disturbance is characterised by alternations in the membrane potential and the intracellular calcium concentration, which in turn can lead to sudden cardiac death. In the present computational study, we make further inroads into understanding this severe condition by investigating the impact of calcium buffers and L-type calcium channels on the formation of subcellular calcium alternans when calcium diffusion in the sarcoplasmic reticulum is strong. Through numerical simulations of a two dimensional network of calcium release units, we show that increasing calcium entry is proarrhythmogenic and that this is modulated by the calcium-dependent inactivation of the L-type calcium channel. We also find that while calcium buffers can exert a stabilising force and abolish subcellular Ca2+ alternans, they can significantly shape the spatial patterning of subcellular calcium alternans. Taken together, our results demonstrate that subcellular calcium alternans can emerge via various routes and that calcium diffusion in the sarcoplasmic reticulum critically determines their spatial patterns.


Intoduction
Cardiac arrhythmias constitute a leading public health problem and cause most cases of sudden cardiac death. In the US alone, sudden cardiac death accounts for approximately 300, 000 − 450, 000 lives every year [1]. Among the many forms of cardiac arrhythmias, cardiac alternans feature prominently. This 5 rhythm disturbance at the level of a single cardiac myocyte is characterised by alternating patterns of the membrane potential and the intracellular calcium (Ca 2+ ) concentration on successive beats. For instance, at one beat, a long action potential duration (APD) is accompanied by a large intracellular Ca 2+ transient, while on the next beat, the APD is shortened concomitant with a 10 small amplitude Ca 2+ transient. As a consequence, contractile efficiency is impaired, which in turn can cause a detrimental reduction in blood flow. In early experimental studies, the intracellular Ca 2+ concentration was averaged across a cardiac myocyte. The advent of high-resolution microscopy revealed that alternating Ca 2+ dynamics were already present at individual Ca 2+ release units 15 (CRU). While one CRU follows a pattern of large-small-large Ca 2+ transients, neighbouring CRUs exhibit small-large-small Ca 2+ transients. Crucially, both CRUs experience the same membrane potential. These findings gave rise to the concept of subcellular Ca 2+ alternans [2][3][4][5][6][7][8][9][10] and illustrated that nonlinear processes govern cardiac dynamics across multiple scales: the cell wide membrane 20 potential and the Ca 2+ fluxes restricted to single dyadic clefts.
The existence of subcellular Ca 2+ alternans reinforces the notion of cardiac myocytes as a network of networks. Each CRU can be conceptualised as a network of interacting components such as L-type Ca 2+ channels, sodium-calcium exchangers (NCXs) and ryanodine receptors (RyRs). These local networks are 25 then coupled via Ca 2+ diffusion through both the cytosol and the sarcoplasmic reticulum (SR). This interconnectedness offers multiple explanations for the origin of subcellular Ca 2+ alternans. On the one hand, we have previously shown that Ca 2+ alternans can emerge purely through coupling [11]. In other words, an isolated CRU displays a regular period-1 orbit, but upon increasing the cou- 30 pling strength between CRUs, period-2 orbits characteristic of Ca 2+ alternans emerge. Crucially, the shape of the Ca 2+ alternans may depend on the form of coupling. In a recent model of Ca 2+ cycling, Ca 2+ alternans emerge for dominant cytosolic coupling via the traditional period-doubling bifurcation, where an eigenvalue of the associated linearised map exits the unit disk through −1 35 along the real axis [12,13]. Here, each node exhibits alternating Ca 2+ dynamics, and neighbouring nodes (or nodes in different parts of a cell) oscillate out-ofphase with each other. For dominant luminal coupling, there is a saddle-node bifurcation, where the leading eigenvalue leaves the unit disk at +1 along the real axis. In this case, each node follows a period-1 orbit, but the amplitudes 40 of neighbouring CRUs varies. On the other hand, changes to the molecular components of a CRU can induce Ca 2+ alternans, exemplified by weakening sarco-endoplasmic Ca 2+ ATP (SERCA) pumps or increasing Ca 2+ flux through L-type Ca 2+ channels.
To date, investigations on how Ca 2+ alternans emerge due to modifications 45 at the CRU level have almost exclusively focussed on dominant cytosolic coupling [3,[14][15][16][17][18][19][20]. However, the question as to whether Ca 2+ diffusion in the SR is slow or fast -and hence weak or strong -is still unanswered [21][22][23]. Here, we focus on strong SR Ca 2+ diffusion and explore the impact of two modifiers of the local Ca 2+ dynamics on the genesis of subcellular Ca 2+ alternans: L-type 50 Ca 2+ channels and Ca 2+ buffers.
The L-type Ca 2+ channel has received significant attention due to its central role in excitation-contraction coupling [24][25][26][27]. Its contribution to the formation of Ca 2+ alternans is ambiguous though [28]. On the one hand, several studies have provided compelling evidence that altering the dynamics of L-type Ca 2+ 55 channel through e.g. cooperative gating or reducing the current can either promote or inhibit Ca 2+ alternans [29,30]. On the other hand, Ca 2+ alternans have been observed with clamped membrane voltage, thus limiting the degree of control that L-type Ca 2+ channels can exert on the genesis of Ca 2+ alternans [31]. Here, we investigate the role of Ca 2+ -dependent inactivation of the L-type 60 Ca 2+ channel on the dynamics of a CRU, which occurs in addition to voltagedependent activation and inactivation. [32,33]. We find that Ca 2+ -dependent inactivation affects the formation of subcellular Ca 2+ alternans in a nontrivial manner that depends on the unitary current of the L-type Ca 2+ channel.

Materials and methods
We consider a two-dimensional network of 15×10 CRUs, where the dynamics of a node with label µ is governed by the Shiferaw-Karma model [31] dc µ The Ca 2+ concentrations in the subsarcolemmal space and in the cytosolic bulk are denoted by c µ s and c µ i , respectively, while the total Ca 2+ concentration in the SR and the Ca 2+ concentration in the unrecruited SR are given by c µ j and c µ u , respectively. The Ca 2+ release current from the unrecruited SR into the subsarcolemmal space is I µ r , and we refer to the L-type Ca 2+ current, the NCX current and the SERCA uptake current by I µ CaL , I µ NCX and I µ up , respectively. The model contains four diffusive currents with timescales τ s , τ c , τ sr and τ a , describing coupling between the subsarcolemmal space and the cytosolic bulk, through the cytosolic bulk between neighbouring CRUs (indexed by I n ), between the total and unrecruited SR, and through the SR between neighbouring CRUs (indexed by I n ), respectively. In some instances, we report the network coupling strengths as inverse of the timescales, i.e. σ Of particular interest for the present study is where c e sets the EC 50 value, i.e. the value of the subsarcolemmal Ca 2+ concentration c s at which q ∞ equals 0.5, and γ controls the sensitivity of Ca 2+dependent inactivation. Essentially, the larger γ the more step-like the inactivation around a Ca 2+ concentration of c e . Buffering is modelled based on the fast-buffer approximation [50,51] yielding where B SR denotes the total buffer concentration in the SR and K SR the asso-

Results
To establish a baseline for our findings, we first investigate the dynamics  Table 1 and σc = 15s −1 , σsr = 3s −1 .
row from the bottom. When we decrease T p , we observe the emergence of subcellular Ca 2+ alternans as depicted in Fig For dominant luminal coupling, where τ c τ sr (σ c σ sr ), we again find stable synchrony at low pacing frequencies (see Fig. 2A). Indeed, the space-  Table 1 and σc = 2s −1 , σsr = 15s −1 .
time plot of the subsarcolemmal Ca 2+ concentration is identical to the one in we find subcellular Ca 2+ alternans that emerge via a saddle-node bifurcation at the network level, in contrast to a period doubling bifurcation for the former.
As Fig. 2B highlights, each CRU follows a period-1 orbit, but this orbit differs 115 amongst the CRUs in the network. Figures 2C and 2D illustrate that here, CRUs on the left form large Ca 2+ amplitude transients, while the transients are smaller towards the right. In the following we will use Figures 2B to 2D as a reference case and contrast them with the network behaviour when we alter the behaviour of the L-type Ca 2+ channel and that of Ca 2+ buffers.  Table 1.
when γ is small. On the other hand, as γ is increased, subcellular Ca 2+ alternans emerge via a saddle-node bifurcation as shown in Figs. 3B-3D. Figure 3B Table 1.
where c µ,i s is the maximum of c µ s on the ith beat. Then, we determine the maximum of all θ µ across the CRU network, θ = max µ θ µ . When i Ca is small, θ 4000 6000 8000 10000  Table 1.
larger values of γ, we observe a sharp transition from synchrony to subcellular Ca 2+ alternans as the L-type Ca 2+ channel becomes stronger. There appears to be an L-shape stability boundary in that for a large range of γ, subcellular Ca 2+ 165 alternans appear for approximately the same value of i Ca , while for a large range of i Ca , alternans set in for approximately the same small value of γ. We also note that the transition from stable synchrony to subcellular Ca 2+ alternans is quite abrupt, as indicated by the sharp transition from blue to yellow. Taken together, our findings provide strong evidence that the L-type Ca 2+ channel can 170 initiate subcellular Ca 2+ alternans, either via its Ca 2+ -dependent inactivation or the strength of its unitary Ca 2+ current.

Buffers
All results so far were obtained for constant buffer contributions. In other words, we set β(c µ i ) and β(c µ s ) to constants β i and β s , respectively, consistent 175 with earlier work [11]. In this way, we eliminate any time-dependent modulation of the Ca 2+ dynamics through binding and unbinding to Ca 2+ buffers. Under more general conditions, however, Eq. (4) entails that β(c µ i ) and β(c µ s ) oscillate with the same frequency as c µ i and c µ s , respectively. Figure 6A illustrates that in this case, subcellular Ca 2+ alternans can be abolished and synchrony is stable.

180
This behaviour needs to be contrasted with that depicted in Fig. 4B, which we would obtain with the parameter values used in Fig. 6A upon replacing the dynamic buffers with the constant buffers used in Fig. 4B. In other words, while the dynamics of the L-type Ca 2+ channel can induce subcellular Ca 2+ alternans (as demonstrated in Fig. 4), dynamic Ca 2+ buffers can rescue this pathologi-  Table 1 and K SR = 6.0µM, K T = 600.0µM, K Cd = 7.0µM, B SR = 250.0µ mol/1 cytosol, B T = 12000.0µ mol/1 cytosol, B Cd = 1.0µ mol/1 cytosol. prompted us to explore another form of non-responsive buffers. The sensitivity of buffers is usually determined by their dissociation constants, which in the present study are the three constants K SR , K Cd and K T in Eq. (4), as well as the corresponding concentration of binding sites B SR , B Cd and B T . By choos- 190 ing appropriate values, we can effectively "desensitise" the Ca 2+ buffers. As Figs. 6B-6D illustrate for the desensitised dynamics, subcellular Ca 2+ alternans re-emerge consistent with a saddle-node bifurcation. Figure 6B Table 1.
Ca 2+ alternans, we next explore the impact of the buffer time course on the 220 formation of subcellular Ca 2+ alternans. To do this in a controlled fashion, we extract the time course of both β(c µ s ) and β(c µ i ) from the full nonlinear model and then clamp the buffer time courses at each node to these profiles. In other words, each node experiences nonlinear buffer dynamics, but the buffers do not alternate from node to node. As Fig. 8A reveals, we obtain subcellular 225 Ca 2+ alternans that differ from those reported so far in this study. Here, every node in the network follows the same period-2 orbit characteristic of subcellular Ca 2+ alternans that emerge via a period-doubling bifurcation. Figures 8B and   8C illustrate the uniform behaviour across the network that alternates between successive beats. This spatial pattern is known as spatially concordant Ca 2+  Table 1.
alternans. When we change the coupling strengths, but keep all other parameter values unaltered, we observe spatially discordant Ca 2+ alternans as plotted in Figs. 8D-8F. Every node follows a period-2 orbit, but different parts of the network oscillate out-of-phase with each other.
Our results so far strongly suggest that the time course and amplitude of 235 Ca 2+ buffers significantly impacts on the genesis of subcellular Ca 2+ alternans. Figure 9 shows results from an in silico experiment in which we tune the Ca 2+ buffer dynamics from constant (ε = 0) to fully nonlinear (ε = 1). As a measure for the strength of subcellular Ca 2+ alternans, we report the maximal beat-tobeat variation θ as defined after Eq. (5). As ε increases, we find a monotonic  Table 1.

Discussion
Subcellular Ca 2+ alternans have been firmly linked to the genesis of cardiac arrhythmias. Despite this crucial connection, we still lack a complete picture of 245 how the dynamics of the intracellular Ca 2+ concentration transitions from its healthy period-1 orbit to its various pathological forms.
Our focus has been on understanding subcellular Ca 2+ alternans in tubulated myocytes, such as ventricular myocytes. The presence of t-tubules in these cells gives rise to well-defined CRUs, which form a network where nearest neigh- to subcellular Ca 2+ alternans as is manifest from the abrupt colour change from blue to yellow. It remains to be seen whether this behaviour can be understood more formally in terms of a phase transition.
All results for the L-type Ca 2+ channel were obtained with constant buffer contributions. However, since the concentration of Ca 2+ bound buffers directly 300 depends on the intracellular Ca 2+ concentration, the buffer function β in Eq. (4) should evolve over time. Using the full nonlinear buffers, we find that subcellular Ca 2+ alternans are extinguished (Fig. 6A) Fig. 8, we broke this connection and clamped the Ca 2+ buffer dynamics in such a way that each node exhibits the same nonlinear orbit. In other words, Ca 2+ buffers alternate at each node, but there is no spatial variation of the buffer dynamics. In this regime, the patterns of the subcellular Ca 2+ alternans vary drastically from the 315 ones we observed so far. We found spatially concordant alternans, which can transition into spatially discordant alternans upon altering the coupling strength of cytosolic and SR diffusion. While we employed buffers to induce this pattern change, it is conceivable that such dynamics could originate from other dynam-ical variables of cardiac Ca 2+ cycling. In this case, our results point to more subtle dependencies in that the nonlinear dynamics of cardiac Ca 2+ cycling can be easily disturbed into new dynamic regimes, potentially inducing a plethora of cardiac arrhythmias. It is therefore astonishing that cardiac Ca 2+ dynamics more often than not behaves completely regularly; a fact that certainly deserves more attention. 325 We first reported the emergence of subcellular Ca 2+ alternans via a saddlenode bifurcation in a PWL caricature of an established Ca 2+ cycling model [12]. One might wonder if this novel form of subcellular Ca 2+ alternans is a consequence of the approximations used in the derivation of the PWL model.
The results presented here show that this is not the case. The fully nonlinear 330 model exhibits the same instabilities. This provides further evidence that PWL models are valuable in exploring the behaviour of complex nonlinear systems and thus adds to earlier success stories such as the McKean model, which represents a PWL version of the Fitzugh-Nagumo model for the propagation of neural action potentials [53][54][55]. The advantage of PWL models is that the 335 majority of the analysis can be performed semi-analytically, which greatly facilitates the exploration of the associated parameter space. In turn, this allows for a more comprehensive classification of the possible dynamics. In contrast, fully nonlinear systems can often only be dissected via direct numerical simulations, which is often only done for a small subset of parameter values. In this respect, 340 PWL models can provide guidance for the analysis of the nonlinear systems and where to explore in parameter space for interesting behaviour.
The last point becomes especially pertinent for the exploration of the different spatial patterns that emerge via a saddle-node bifurcation. As Figs. 1, 2, 6 and 8 illustrate, the Ca 2+ profiles across the network exhibit significant vari-345 ability. In a PWL model, these patterns can be classified and understood from a linear stability analysis, which can be performed in closed form [12,13]. On the other hand, the nonlinear model requires direct simulations, which are computationally more expensive and limited in scope as to what parameter values to sample.

350
As stated above, our focus here is on tubulated myocytes. However, Ca 2+ alternans have also been observed in non-tubulated cells such as atrial myocytes and failing ventricular myocytes [56][57][58][59][60][61][62]. In these cells, L-type Ca 2+ channels are only located at the cell periphery, where they trigger Ca 2+ release from the SR through the RyR. A Ca 2+ wave then propagates centripetally from the pe-355 riphery via diffusion and Ca 2+ induced Ca 2+ release [63,64]. Conceptually, it therefore makes sense to distinguish junctional CRUs (that contain L-type Ca 2+ channels) and non-junctional CRUs (that lack L-type Ca 2+ channels). Due to the stronger reliance on Ca 2+ diffusion, it will be interesting to explore how differences in the diffusive coupling between CRUs and the fact there are two 360 classes of CRUs shape subcellular Ca 2+ alternans and whether the bifurcation structure observed for tubulated myocytes carries over to non-tubulated ones.
Answering this question will not only unravel further similarities or differences between tubulated and non-tubulated myocytes, it will also advance our understanding of atrial fibrillation, which is projected to become epidemic with an 365 ageing population [65].
We here provide the parameter values used in the study unless otherwise stated.