Buffering and total calcium levels determine the presence of oscillatory regimes in cardiac cells

Calcium oscillations and waves are often behind instances of extra depolarization in cardiac cells, eventually giving rise to life-threatening arrhythmias. In this work, we study the conditions for the appearance of calcium oscillations in both a detailed subcellular model of calcium dynamics and a minimal model that takes into account just the minimal ingredients of the calcium toolkit. To avoid the effects of homeostatic changes and the interaction with the action potential we consider the somewhat artificial condition of a cell without pacing and with no calcium exchange with the extracellular medium. This permits us to isolate the main reasons responsible for the oscillations by controlling externally the total calcium content of the cell. We find that as the calcium content is increased, the system transitions between two stationary states, corresponding to one with closed ryanodine receptors (RyR) and most calcium in the cell stored in the sarcoplasmic reticulum (SR), and another, with open RyRs and a depleted SR. In between these states, calcium oscillations may appear. This transition depends very sensitively in the amount of buffering in the cell. We find, for instance, that at high values of calsequestrin (CSQ) oscillations disappear, while they are present for a broad range of parameters at low values of CSQ. Using the minimal model, we can relate the stability of the oscillating state to the nullcline structure of the system, and find that its range of existence is bounded by a homoclinic and a Hopf bifurcation. Author summary In cardiac cells, calcium plays a very important role. An increase in calcium levels is the trigger used by the cell to initiate contraction. Besides, calcium modulates several transmembrane currents, affecting the cell transmembrane potential. Thus, dysregulations in calcium handling have been associated with the appearance of arrhythmias. Often, this dysregulation results in the appearance of periodic calcium waves or global oscillations, providing a pro-arrhythmic substrate. In this paper, we study the onset of calcium oscillations in cardiac cells using both a detailed subcellular model of calcium dynamics and a minimal model that takes into account just the minimal ingredients of the calcium toolkit. Both reproduce the main experimental results and link this behavior with the presence of different steady-state solutions and bifurcations that depend on the total amount of calcium in the cell and in the level of buffering present. We expect that this work will help to clarify the conditions under which calcium oscillations appear in cardiac myocytes and, therefore, will represent a step further in the understanding of the origin of cardiac arrhythmias.

deterministic model of calcium in a cardiac cell (or in a CaRU). Within this model, we 48 study the existence and stability of different solutions. We show that oscillations 49 typically appear at high global calcium concentration and/or high RyR open probability. 50 Their appearances depend on a delicate balance between the total calcium level in the 51 cell and the level of buffering of calcium available. For instance, at high values of 52 calsequestrin (CSQ), the system presents a transition from a low concentration, 53 excitable state, to a high concentration state. Such a transition has been proposed to be 54 the basis of complex states, such as long-lasting sparks [28]. At low concentrations of 55 CSQ, in between these two stable states, oscillations appear. We study this transition 56 using a minimal model, that includes the concentration of dyadic and SR calcium and 57 the open probability of the RyR and show that it suffices to explain the appearance of 58 oscillations. A further reduction to a minimal two-dimensional model allows us to 59 explain the transition to the oscillatory regime in terms of the nullcline structure of the 60 system.  Detailed subcellular calcium model 72 We model the spatial structure of the cell as in a previous model of a cardiomyocyte 73 presented in Marchena and Echebarria [27], which has been modified to add the effects 74 of calsequestrin. The equations of the model read: dc bi (r, t) dt = J bi (r, t) where c i is the calcium concentration in the cytosol, c tot sr the total calcium concentration 76 in the SR, and c bi represents the concentration of a given buffer in the cytosol. Besides, 77 J rel and J up are the release flux from the SR and the uptake by SERCA, respectively, 78 and J bi represents the binding of free calcium to the different buffers in the cytosol 79 (TnC, SR binding buffer and CaM). These currents are given by: J up = g up (c i /K i ) 2 − (c sr /K sr ) 2 1 + (c i /K i ) 2 + (c sr /K sr ) 2 (5) The spatial structure of the model includes cytoplasmic and SR spaces, with a spatial 81 discretization of 100 nm.
Assuming fast binding, the stationary condition for c bSQ (ċ bSQ = 0) is: where K SQ = k of f SQ /k onSQ is the dissociation constant. From this, the concentration 106 of free calcium can be obtained solving Eq. (2) for the total amount of calcium in the 107 SR, c tot Solving this equation we obtain the value of free calcium in the SR, The advantage of using this formulation of the rapid buffer approximation over the 110 more usual, for example in [32], is that it conserves mass exactly.

111
Under physiological conditions, the total amount of calcium in the cell at steady For a better comparison with the results from a reduced calcium model, described later, 121 we will consider as our control parameter the average calcium content of the cell 122 c T = Q T /(v i + v sr ). Thus, in our simulations,c T is a constant value that is determined 123 by the initial conditions for cytosolic and luminal calcium (free and bound to buffers). 124 Reduced calcium model 125 The minimal model for the local dynamics of calcium is based on the schematics shown 126 in Fig. 2. We consider a simplified description of the system, with dynamics of the total 127 dc tot with the currents given by A detailed derivation of these equations and their range of validity can be found in 131 Appendix A. Notice that, as in the full model, we consider a situation where no external 132 pacing is imposed. In this sense, neither external intake from the LCC is considered, nor 133 any extrusion via the sodium-calcium exchanger.

134
For simplicity, we consider a SERCA pump without an equilibrium condition, that 135 always pumps calcium from the cytosol to the SR. This gives a basal solution at 136 c i = c d = 0, instead of the physiological value of ∼ 100nM. However, given that, at basal 137 conditions, c sr ∼ 1mM, this is a reasonable simplification. As in the detailed subcellular 138 model, we assume the approximation of rapid CSQ buffer, so we can compute the 139 amount of free luminal calcium c sr from the total luminal calcium c tot sr from Eq. (9).

140
To close the system we should introduce an extra equation for calcium concentration 141 in the cytosol c i . However, as we assume that the total calcium content in the cell is 142 constant, then we have a conservation equation. Therefore, we can compute c i solving 143 the following quadratic equation for the conservation ofc T where v cell is the unit volume defined as v cell ≡ v i + v sr + v d , and B b is the 145 concentration of a generic buffer in the cytosol.

146
To simplify the analysis, we proceed to work with the assumption that the dynamics 147 of the RyRs is faster than that of calcium concentration (Ṗ o 0), obtaining then a 148 minimal two-variable model. This will be our base-line minimal model. However, we  In this case, we assume that the open and close dynamics of the RyR are fast, so we can 155 assume that it is in a quasi-steady state (Ṗ o 0). Then, from Eq. (14), we obtain: where the parameter K 2 o = k m /k p is the ratio of the open and close rates of the RyR. Then, with these assumptions, the simplified model becomes Fast dyadic calcium dynamics 159 In order to test the robustness of the analysis, we also consider a simplified model given 160 by Eqs. (12)- (14), in the limit of fast dynamics in the dyadic space and takeċ d 0.

161
Then, from Eq. (12): Substituting this expression in Eqs. (13) and (14), we obtain another minimal model, 163 given by where again, c i must be computed solving the quadratic equation for the conservation of 165 massc T [Eq. (16)]. For simplicity we will consider the case when no calsequestrin is

168
We first present the results of the numerical simulations of both the full detailed model 169 and the minimal model of calcium cycling. Both produce the same basic scenarios for 170 intracellular calcium dynamics, with three different dynamical behaviors, which we then 171 proceed to analyze. The goal of the development of the minimal model is, precisely, to 172 be able to perform this analytical treatment and check how the behavior depends on 173 total calcium and buffering levels.  When the calcium load increases, the system starts to spontaneously show calcium 187 waves. These waves persist in time with different shapes and durations, giving rise to a 188 nearly periodic oscillation in the global calcium signal. Roughly, we observe one calcium 189 wave per second (Fig. 4). Waves are normally initiated at different sites each time but 190 they appear systematically indicating a strong oscillation at the substrate level that we 191 will address in the discussion.  We must point out that a similar trend has been observed experimentally by Stevens 204 et al [23], even if in the experimental preparation the control parameter was the amount 205 of cytosolic calcium, and not total calcium, as in our simulations. Oscillations appear as 206 the amount of calcium in the cell increases, giving rise to a state with depleted SR 207 calcium (and RyRs in the open state), at high calcium concentrations. Furthermore, 208 experimentally it has been shown that changes in buffering levels can have also 209 important effects on this transition. More specifically, Stevens et. al [23] have shown 210 that the reduction of CSQ in the SR bulk enhances the appearance of oscillations. We 211 have checked if this situation is also present in our simulation and found this to be the 212 case. As shown in Fig. 4, when we reduce the CSQ concentration, the oscillations 213 appear at lower values ofc T , they have a higher frequency and the range of oscillations 214 in terms ofc T becomes broader.

215
Minimal model without calsequestrin 216 We proceed to explain the results obtained in the minimal model where we find the 217 same qualitative behavior as in the results obtained with the full subcellular stochastic 218 model. We have performed simulations in the approximation of fast RyR dynamics at 219 different values of the cell average calcium concentrationc T . We consider first the case 220 when no calsequestrin is present, B SQ = 0. As we observed in the full subcellular model, 221 at low values ofc T the system remains in a low concentration steady state (Fig. 5). In 222 this state the system is excitable, so the fixed point is locally stable, but a large enough 223 perturbation produces an increase in calcium concentration, resulting in a calcium 224 transient typical from CICR. On the other hand, at high calcium loads, there is a new 225 fixed point with high cytosolic calcium concentration, that has been related to the 226 appearance of long-lasting sparks [28]. As in the subcellular model, in between the low 227 concentration fixed point and the permanently open state, the system presents . After a transient, the system ends up in either a steady state which is excitatory at low levels of total calcium in the cell with observed low levels of calcium in the cytosol, in an oscillatory state with intermediate levels of total calcium in the cell, or in a state of high total levels of calcium in the cell with observed high cytosolic calcium levels.  Indeed, the number of stationary solutions changes with the calcium concentration 230 c T . The fixed points of the system can be found from: together with: Eqs. (23)-(25) represent three algebraic equations that give the concentrations c i , c d and c sr as a function of total calcium concentration in the cellc T . When no calsequestrin is present, B SQ = 0, it is easy to obtain that Introducing these expressions into Eq. (23)  To calculate the stability of the stationary solutions, we have computed the value of 239 the eigenvalues of the Jacobian matrix, corresponding to Eqs. (18)- (19). We find that, 240 while the lower branch is always stable, the other branch of solutions is unstable for a 241 large range of parameters (Fig. 7), due to the appearance of oscillations. The stability of 242 the corresponding periodic orbit has been calculated using XPPAUT [33] (Fig. 8). We 243 obtain that, asc T is increased, a limit cycle appears in a global homoclinic bifurcation, 244 with zero frequency (Fig. 8c). Increasingc T , this limit cycle finally disappears in a Hopf 245 bifurcation, at which the upper fixed point becomes stable (Fig. 8b).   In addition, at this point, the system presents five fixed points where only the lowest 259 and uppermost are stable. It is important to note that this shows that CSQ enhances 260 the elimination of oscillations. Finally, for high values ofc T , the system has again three 261 fixed points.

263
The mathematical tractability of the minimal model allows us to get a better 264 understanding of the transition to the upper state via oscillations. A first insight can be 265 obtained by plotting the nullclines of the system (Fig. 10), which can help us understand 266 the main mechanisms behind the transition to the oscillatory state. In particular, we 267 obtain the critical average calcium concentration for the onset of oscillations, which 268 depends on buffering levels, and the conditions for the appearance of the upper state.   (Fig. 11c). Due to the emergence of the pinch-off, the system 281 dynamics follows the lower nullcline up to the tip, at the largest value of c sr , where, due 282 Furthermore, from Eq. (25), we can write c i in terms of c sr (assuming B SQ = 0, Assuming that the concentration of bound cytosolic calcium is much larger than that of 296 free cytosolic calcium, we obtain The same equation can be written as Once the pinch-off has been produced (Fig. 10c) respectively. Just at the pinch-off these two points merge. This allows to establish the 304 critical value at which the oscillation appears as the one that makes zero the 305 discriminant of the second order polynomial of c * sr . After some algebra, the condition 306 for the criticalc * T becomes 12. Dependency of the onset of oscillations with buffer parameters. The filled dots represent the control values B b = 80µM, K b = 0.5µM, given in Table 1. which can be written as This gives Using the parameters in Table 1, this expression gives a critical value for the onset of which gives a value of c * sr = 0.86 mM. At the onset of oscillations, there is thus a sudden 316 decrease in basal SR calcium concentration, to less than half its previous value before 317 the oscillations.

318
An increase in the quantity of cytosolic buffers (higher B b ) results in a delay in the 319 onset of oscillations, that would occur at higher calcium load (Fig. 12a). A higher 320 calcium affinity (lower K b ) on the contrary, would result in oscillations at lower loads 321 (Fig. 12b). It is interesting to notice, too, that the strength of SERCA does not affect 322 the onset of oscillations.
A fixed point will be stable provided f 1 + g 2 < 0. When f 1 g 2 − f 2 g 1 > 0 the stable fixed 328 point corresponds to a node and if f 1 g 2 − f 2 g 1 < 0 to a stable spiral. We can use this to 329 relate the slope of the nullclines to the stability of the upper fixed point. Let us denote 330 α ≡ċ d and β ≡ċ sr the time derivatives of the independent variables.

331
At large values of the total concentration (see Fig. 13a, withc T = 75 µM), the slope 332 of the black nullcline (α = 0) at the fixed point is positive, while the slope of the orange 333 nullcline (β = 0) is negative. Then, increasing c d at constant c sr , α goes from being 334 positive to negative. This means that f 1 ≡ ∂α/∂c d < 0. Using the same argument, it is 335 easy to check that also g 2 ≡ ∂β/∂c sr < 0, and therefore the fixed point is stable. At slope of the black nullcline becomes negative. Thus, in this case, while g 2 is still 338 negative, f 1 becomes positive, and it is not possible to determine the stability of the 339 fixed point. It will depend on the speed of rate of c d and c sr close to the fixed point. If 340 the dynamics of c d is faster, then one expects this point to be unstable, if it is c sr that 341 varies fast, then stable. One would then expect that buffers that change the dynamics 342 either in the cytosol or in the SR would effect the stability of the fixed point and, 343 therefore, the range of existence of the limit cycle.

344
Robustness of the results. Fast dyadic calcium dynamics 345 Fig 14. Solutions as a function of total calcium concentrationc T , with calsequestrin concentration set to zero when the fast dyadic approximation is used. We obtain the same type of structure as expected.The system can be in a monostable state, which is excitatory (low load), in an oscillatory state (intermediate load), or in a bistable state (high loads), where it usually ends in a state of open RyR and depleted SR calcium concentration.
It is useful to check that the basic points of our discussion hold when different  (Fig. 14) and test that we find the same basic behavior: at low values 352 ofc T the system remains in a low concentration steady state (Fig. 15), but oscillations 353 appear for a range ofc T , up to a limit where an upper state becomes stable.

354
More importantly, the basic structure of the nullclines determines the possible 355 solutions and again, oscillations appear when pinch off is produced (Fig. 14). The fixed 356 points in this case are the same as in Fig. 7, since they correspond to fixed points of  [34], to the control of the pacemaker rhythm 365 in both early embryonic heart cells [35,36] and in sinoatrial nodal pacemaker cells 366 (SANCs) [37][38][39]. Spontaneous calcium releases or oscillations are also related to the 367 appearance of early or delayed-afterdepolarizations [40], which may give rise to open state [23,41]. Calcium overload can be obtained, for instance, by inhibition of the 373 Na + -K + pump current I N aK , that results in [Na + ] i overload. The consequent build-up 374 of [Na + ] i reduces the effectiveness of the Na + -Ca 2+ exchanger at removing calcium 375 from the cell and intracellular calcium concentrations become elevated. A similar effect 376 is observed in models of hypercalcemia [42]. The effect of elevated [Na + ] i has been studied in computational models, finding calcium oscillations [43], that, depending on 378 the model, appear via a supercritical Hopf [24] or a homoclinic bifurcation [25,26]. 379 We have shown, using a full subcellular model, that under global calcium overload, 380 SR oscillations appear, leading finally to a state with permanently open RyR and 381 depleted SR. In these simulations, oscillations give rise to periodic calcium waves, that 382 propagate along the myocyte. To obtain a better understanding of the origin and 383 parameter dependence of these transitions, we have studied them in a simplified model 384 of the calcium dynamics. Despite not including all the physiological details (or because 385 of that), minimal models are often useful to gain a better understanding of the origin of 386 complex calcium rhythms [44], as oscillations [45]. We have thus analyzed the different 387 transitions within a minimal model, that takes into account the RyR dynamics, as well 388 as fluxes among different compartments. This allows us to explain the origin of 389 oscillations using a nullcline analysis, as well as to give analytic expressions for the 390 different transitions. One interesting conclusion is that buffers affect heavily the 391 dynamics. The effect of buffers on oscillations has been studied previously [32], and 392 experimentally a change in CSQ levels has been observed to alter the range of values of 393 cytosolic calcium at which oscillations are observed [23]. Here we find that an increase 394 in the levels of CSQ prevents altogether the oscillations, obtaining a direct transition to 395 an open state, that also occurs at larger values of total calcium content. We have shown 396 that this picture is robust and independent on which is faster, whether the opening of 397 the RyR2 or the diffusion of calcium near it.

398
The fact that oscillations appear both in a fully detailed stochastic model and in a 399 simple deterministic model seems another proof of the robustness of the behavior. This 400 robustness is not unexpected since it has already been observed in other calcium 401 signaling processes, such as calcium oscillations in hepatocytes [46,47]. Nevertheless, 402 while the minimal deterministic model presents a sudden drop in SR calcium content at 403 the onset of oscillations, the subcellular stochastic model presents a gradual decrease in 404 the level of basal SR, in full agreement with the experimental results obtained by 405 Stevens et al [23]. Thus, the comparative analysis of our results seems to suggest that 406 the origin of oscillations in experiments arises from the smearing out of a sharp 407 transition to the oscillatory state, due to the stochasticity of the RyR in an extended In this appendix, we show how to derive the simplified model in Eqs. (12)- (14) from a detailed deterministic model of calcium handling. We consider cytosolic and luminal spaces, each separated into different compartments, i.e, dyadic and cytosolic by the one side, and junctional and network SR, by the other. We also consider the effect of two buffers in the cytosol, TnC and SR, and Calsequestrin (CSQ) in the SR. The effect of the buffers is actually very relevant, as they change the structure of possible solutions of the system.
With this, the dynamics can be described by the following set of deterministic equations for the calcium concentration at the different compartments (40) where c d , c j , c sr , c i stand for the concentration in dyadic space, free luminal, SR network and cytosol, P o is the fraction of RyRs in the open state, τ i and τ sr are the diffusion time constants out of the dyadic space and SR network, v i , v d , v jsr and v sr are the volumes associated with each compartment, g up is the strength of SERCA pump, K s is the concentration at which SERCA closes and g is the strength of the release current. The dynamics of the buffers are given by linear reactions with the following set of ODEs Since the total amount of calcium has just a variation of about by a 5% or 10% over a calcium cycle, we assume that the total calcium concentration,c T , is fixed. In this way, one does not have to solve a differential equation for the calcium concentration in the cytosol. Rather, it is derived from the algebraic equation: Homeostatic behavior of the cell will eventually load the system more or less, increasing or decreasingc T .
Gating of the RyR can been described by Markov models that describe the transitions among different conformations of the channel. Thus, for the dynamics of the RyR we consider a phenomenological four state model [48].
with the last equation given by the condition that the sum of probabilities is equal to 1 P R and P o are the ratios of local RyR in the recovered and open states. P I A and P I B stand for the terminated states. To make the model treatable we assume several hypotheses.
1. We consider rapid equilibrium of c j (ċ j = 0). Thuṡ 2. At first order, we approximate c j as c sr in the right hand side of Eq. (50). This approximation is valid when τ sr (51) where we have eliminated the dependence with c j .
3. We apply the rapid buffer approximation in the SR because CSQ is very fast [29,30]. The derivation of c sr in terms of c tot sr have been already shown in the main text (see Eqs. 7-10). 4. We assume that all the buffers in the cytosol are in equilibrium. Then 5. Besides, we combine both TnC and SR buffers in only one buffer with a total concentration of B b and an affinity of K b , which will be 6. Finally, we consider that the inactivated states of the RyRs do not play a role in the mechanisms that produce oscillations. For that reason, we reduce the four state model to a two state model The parameters of these equations are taken from the literature where, except for those of the RyR, are well documented. The ratio of SR to cytosol in the cell is roughly 1 to 20-30. The order of magnitude of the volume of the cleft where calcium is released is around 10 −3 µm 3 . SERCA is roughly activated at around 2-5 µM and closed at around 15-20 µM with the number of buffers relevant to absorption around 50-100 µM , taking an average affinity of around 0.5 µM . The whole set of parameters is given in Table 1.