Two-variable nullcline analysis of ionic general equilibrium predicts calcium homeostasis in ventricular myocytes

Ventricular contraction is roughly proportional to the amount of calcium released from the Sarcoplasmic Reticulum (SR) during systole. While it is rather straightforward to measure calcium levels and contractibility under different physiological conditions, the complexity of calcium handling during systole and diastole has made the prediction of its release at steady state impossible. Here we approach the problem analyzing the evolution of intracellular and extracellular calcium fluxes during a single beat which is away from homeostatic balance. Using an in-silico subcellular model of rabbit ventricular myocyte, we show that the high dimensional nonlinear problem of finding the steady state can be reduced to a two-variable general equilibrium condition where pre-systolic calcium level in the cytosol and in the SR must fulfill simultaneously two different equalities. This renders calcium homeostasis as a problem that can be studied in terms of its equilibrium structure, leading to precise predictions of steady state from single-beat measurements. We show how changes in ion channels modify the general equilibrium, as shocks would do in general equilibrium macroeconomic models. This allows us to predict when an enhanced entrance of calcium in the cell reduces its contractibility and explain why SERCA gene therapy, a change in calcium handling to treat heart failure, might fail to improve contraction even when it successfully increases SERCA expression.


Competing interests:
The authors have declared that no competing interests exist.
patients have found no statistical effects of the treatment [23]. These unsuccessful results have been related to reduced viral gene transfer in humans and the effect of neutralizing antibodies present in patients. However, one can not disregard the possibility that SERCA therapy may fail even when SERCA function is improved since its failure or success also depends on how the underlying pathology affects the general equilibrium properties of calcium homeostasis. Highlighting also the relevance of understanding calcium homeostasis in order to predict calcium levels when some channel properties are changed, recent reviews [12,24,25] have shown clear counterintuitive calcium behavior. For example, an increase in entrance of calcium via LCC is not always followed by an increase in total calcium levels at the steady state.
To clarify why the calcium behavior seems often counterintuitive and show its very complex nature we can use that precise example. Consider that, in a cell under constant pacing, the conductance of the L-type calcium channels is suddenly increased. Naively, one would expect this to lead to an increase in SR calcium load, as more calcium enters the cell, and gets accumulated into the SR. This is indeed often the case, as our simulations show for a detailed model of calcium handling in rabbit myocytes (see Fig 1 and details of the model in the next section). It is not too difficult, however, to modify the model with changes in the levels of buffers and SERCA, and obtain counterintuitive results. In this modified model, as the strength of LCC increases, the SR Ca concentration decreases as seen in Fig 2. In this case, a larger entrance of calcium into the cell produces a larger release from the SR, that then is taken more efficiently out of the cell by NCX, resulting in a net efflux of calcium from the cell. This counterintuitive result (described experimentally in [26]) is not at all unique. It seems clear that understanding and predicting calcium homeostatic levels is key to understand the contractibility behavior of cardiac ventricles.

Rabbit computational model structure
General structure. We consider the cell to be a 3-dimensional array composed of 25 × 25 × 50 Calcium Release Units (CaRUs) of volume 0.5 × 0.5 × 2μm 3 following the basic structure described by Mahajan et al [27] (see first panel in Fig 1) where the cell is divided into volumes associated with RyR clusters. Ryanodine Receptors channels (RyR) are introduced following the description by Stern et al [28], as adapted in [29] with four possible states. We do not consider, for the L-type Calcium Channels (LCC), the states sensitive to Barium in [27]. We also introduce important modifications in the way we write the equations of the model to ensure calcium mass balance, so the difference between the calcium entering and leaving the cell matches exactly the change of calcium inside the cell. As in previous models [30], each CaRU is divided into 5 compartments or subvolumes, as shown in the second panel of Fig 1: cytosol (i), subsarcolemma (s), network sarcoplasmic reticulum (nsr), junctional sarcoplasmic reticulum (jsr) and dyadic junction (d), which is the space between a group of LCC channels, in the cellular membrane, and a cluster of RyR channels in the membrane of the sarcoplasmic reticulum. The basic structure of the CaRU is a rectangular cuboid with dz being the distance between Z-planes, which we take to be roughly 2 μm [31,32], while the distance in X-Y is the  (1): Evolution of the SR free calcium concentration when, at a time marked by a black vertical line, the conductivity of the LCC is increased (left) or decreased (right). An increase in intake via the LCC leads to a cardiomyocyte with higher calcium average concentrations in the SR. Panel (2): Response of the cardiomyocyte to an increase in LCC conductivity, as in the first panel, when its buffering and SERCA uptake is different. Now, the increased intake via LCC is homeostatically regulated in a completely different way. The higher level of calcium intake leads to an even larger extrusion for the SR leading to a lower calcium load.
Rate equations in the CaRU. This model considers the diffusion of calcium ions between compartments of the same CaRU, as well as diffusion between neighboring CaRUs; moreover, there are other currents related to the dynamics of calcium buffers, thus differentiating between free and bound to buffers calcium ions in the different compartments. Finally, there are currents due to channels, exchangers or pumps. L-type calcium channels introduce calcium, j LCC , in the dyadic space. Flux across RyR also goes into the dyadic from the junctional SR, j RyR . SERCA pumps calcium j SrCa from the cytosol to the network SR and the Na + /Ca 2+ exchanger is considered to be in the membrane and t-tubules, as LCC, but not close to it. We, therefore, consider that the exchanger flux, j NCX , is sensitive to the calcium gradient between extracellular calcium and calcium concentrations in the subsarcolemma. The dynamics are deterministic except for the RyRs and LCCs which are defined by internal markovian states and whose transitions are considered to be stochastic.
The set of differential equations for Ca 2+ concentration in the different compartments on each CaRU (where we have omitted the i, j, k-th superscripts indexing the position of the CaRUs for simplicity) reads as follows, where we take no-flux boundary conditions: where j ds , j si , j tr represent internal diffusion currents and v d , v s , v i , v jsr and v nsr represent local volumes for dyadic space, subsarcolemma, cytosol, junctional SR and network SR respectively. We take v jsr to be equal to v nsr and reserve the notation c sr for the average concentration of c nsr and c jsr , which can also be defined as the free SR calcium concentration. While in the model we track both, in our analysis, given that at the time scale of one beat both c nsr and c jsr equilibrate and become roughly one and the same with c sr , we will often use the latter expressed in terms of liters of cytosol. We take the cytosolic volume to be half of the available space (0.25 μm 3 ) given the presence of mitochondria, which we neglect in this model. D s , D i and D sr are the inter-CaRU diffusion constants for subsarcolemma, cytosol, and SR respectively. In the case of subsarcolemmal diffusion, we consider that these volumes do not exist across z-planes so (D s ) z = 0. Regarding D i and D sr , they are different from each other, slower in the SR, but the same in all directions. Notice, however, that given that the characteristic length is different across the z-plane direction, the characteristic time scales of diffusion in the Z direction is different than in the X-Y direction.
The model spatial structure considers three different diffusive currents within a CaRU, which depend only on the difference of Ca 2+ concentration among the affected compartments with a given characteristic time of diffusion. These are: diffusion between dyadic junction and subsarcolemma, j ds = (c d − c s )/τ d ; diffusion between subsarcolemma and cytosol, j si = (c s − c i )/ τ s ; diffusion between network and junctional sarcoplasmic reticulum, j tr = (c nsr − c jsr )/τ tr .
Pacing by external clamped potential: Calcium transient and recovery. The above equations have to be complemented with an external clamped potential which fixes an external pacing. In this model, the pacing of the model is set by a periodic clamped action potential with a constant pacing period T and an action potential duration APD described in the Supporting information file S1 File. Each time voltage raises, the model reproduces the Calcium-Induced Calcium-release (CICR) mechanism at the global level from the stochastic behaviour of the RyR clusters, where the LCC opening raises the probability of RyR2 to open, releasing calcium from the SR. This results in an averaged behavior as the one depicted in the third panel of Fig 1. Average free calcium concentration values in the cytosol raise and drop during one beat (calcium transient), while in the SR a transient drop in the concentration of free calcium is followed by its recovery.
Currents of rabbit model. The four basic currents that allow for the normal calcium transient, namely intake of calcium by the LCC, extrusion by NCX, release by the RyR and uptake by SERCA, are fully detailed in the SM. Here we present their basic dependencies.
The current associated with the Na + /Ca 2+ exchanger depends on intracellular and extracellular Na + concentrations ([Na] i and [Na] o respectively), which are set as constant in voltageclamped model, as well as on extracellular Ca 2+ concentration ([Ca] o ), also constant, and subsarcolemmal calcium concentration c s .
The current through the LCC channels, carrying calcium from the extracellular medium towards the intracellular space, depends on the transmembrane voltage, the extracellular Ca 2+ concentration and the calcium concentration in the dyadic volume close to both the RyR cluster and the LCC channels. The current is proportional to the number of LCC channels in the open state (O LCC ). We use the expression given in [27] for the properties of LCC in the rabbit, not taking into account those states related to barium inactivation. This leaves us with five possible states: O for the open state, C1 and C2 for closed states, and I1 and I2 for the inactivated states. In each CaRU, we consider a group of 5 LCC channels, as shown in the second panel of Fig 1. We track the state of each one of the channels using the length rule to establish a transition between states with the uniformly distributed random number generator RANMAR.
The calcium release current, through the RyRs, is diffusive, but depends not only on the difference among c jsr and c d but also on the conductivity of the RyR, which is set by the number of these proteins in the open state (O RyR ): j RyR = g rel O RyR (c jsr − c d ). In each CaRU we consider a cluster of N RyR = 40 RyR, where each one of the receptors can be in one of the 4 possible states shown in the second panel of Fig 1: O, for the open state; C is the closed state, while I1 and I2 are two inactivated/terminated states. Their dynamics are treated stochastically as a function of transition rates which depend, both for opening and inactivation, on calcium concentrations in both the dyadic and the jSR spaces (a detailed description of the transition rates in the RyR can be found in the Supporting information S1 File).
Finally, the current given by the SERCA pump is thermodynamically limited. It is considered to depend only on Ca 2+ concentrations in cytosol and network SR. We give the full expression here given the relevance for our analysis of SERCA gene therapy.
where K i and K sr are half occupation binding constants, and ν up is the maximum uptake strength whose effects on the model will be analyzed in detail.

Calcium buffering implementation.
Buffers are located by construction in the cytosol and in the junctional SR where calcium can bind to them. Therefore, the above equations have to be complemented with the equivalent dynamics equations for the calcium binding to the buffers. In this respect our model presents a very important difference with previous models. We do not consider the fast-buffering approximations as described in [35,36] and used, for instance, in [37], since this leads to a loss in the mass balance in the dynamical equations and it does not particularly speed up the code. In our analysis, it is critical to have an algorithm that strictly balances mass to all orders. This leads us to consider jSR to be the volume around the RyR where calsequestrin is present and to compute the evolution of c jsr TOT , which includes both free calcium in jSR (c jsr ) and calcium bound to calsequestrin in that compartment. The reason is that calsequestrin is the only buffer that is fast enough to slow down the performance of the code if the binding of calcium to the buffer is computed dynamically. Once the total amount of calcium is computed in the jSR, this calcium is split between the free and buffered parts using the fast-buffering approximation where B CSQ is the concentration of calsequestrin in the jsr and K CSQ its dissociation constant. The free calcium concentration c jsr can be obtained solving the quadratic equation: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ffiffi All the other buffers have their own dynamical equation where the calcium taken out from its free level goes to the calcium bound to the buffer. Time scales and affinities are mostly taken from Shannon et al. [38]. A full description is incorporated in the Supporting information S1 File.
Parameters of the rabbit model and its modified version. A full list of parameters of the rabbit model including buffers is provided in the appendix of the Supporting information S1 File. Furthermore, this appendix also includes the changes in the parameters of the modified rabbit model used in Fig 2. This modification of the rabbit model allows us to show in the introduction that the effects of increasing the conductivity of LCC can depend strongly on the animal model. We must stress that we change the parameters of our rabbit model but not its structure. More specifically, our modified rabbit model has different buffer levels in order to change the effect on the level reached by calcium in the transient and a slightly reduced maximum SERCA uptake to mimic a slight reduction in its expression.

General equilibrium method to study calcium homeostasis
We develop here a novel method to study calcium homeostasis. We proceed to show that a bidimensional general equilibrium approach provides the proper framework to understand and predict intracellular calcium levels at steady state. The actual homeostatic calcium level observed in an in-silico rabbit ventricular model can be predicted using a dramatic reduction of dimensionality, where the general equilibrium of thousands of variables is effectively reduced into a general equilibrium of two variables, namely, the level of free calcium in the cytosol and the level of free calcium in the SR. The reduction of the dimensionality allows us to show that calcium homeostasis can be understood in terms of shocks as in basic macroeconomic models and unveil that changes in SERCA function are complex double shocks. This opens the possibility of understanding calcium homeostasis from this new perspective.
As described in the previous section, we consider a subcellular compartmental model of calcium handling, dividing the cell into CaRUs, each incorporating a cluster of RyRs and of LCC, whose gating follows a Markov chain, solved stochastically. Let us denote the current of calcium that goes into the cell through the LCC channels associated with the ith CaRU as j i LCC , where now this current is defined by liter of cytosol (μmol/(Lcyt ms)). The average flux of calcium entering the cell via LCC J LCC (μmol/(Lcyt ms)) is: where N is the number of CaRUs in the cell. Integrating during the period T we obtain the total amount of calcium that enters/intrudes the cell during one beat: We can define similarly the average flux of calcium that leaves/extrudes the cell via the Na-Ca exchanger J NCX and the total amount of calcium that leaves the cell ΔQ out .
Notice that both ΔQ in and ΔQ out depend, in principle, on the value of the roughly one hundred thousand internal variables � i j (where j stands for the number of open LCCs, amount of Ca bound to buffers, etc, at each of the i different CaRUs) that define the state of the cell at the beginning of the beat. We assume that, despite the stochastic nature of LCC and RyR2 channels, average values are well-behaved. This being the case, we can establish that the total amount of calcium Q n T at beat n follows the relation Q nþ1 T ¼ Q n T þ DQ n in À DQ n out so that at steady state (ss) we have: The total amount of calcium released via RyR2 from the SR to the cytosol, ΔQ rel , and the total calcium uptaken by SERCA, ΔQ up , can be equally defined from the average fluxes of all individual RyR clusters, J RyR , and SERCA pumps, J SrCa : This gives a beat-to beat change in the SR calcium concentration given by: In steady state, the amount of calcium in the SR must be constant, leading to As stated before, ventricular myocytes reproduce consistent average behavior. Therefore, fluxes in Eqs (12) and (14) depend on the relevant spatial average variables � � j and not on their spatial distribution.
Except for average calcium concentrations, all the average quantities that define the state of the cell under voltage-clamped conditions have relaxation time scales which are shorter than the fastest pacing period reached. This typically goes from T min = 300 − 400 ms in humans to T min = 150 − 200 ms in rats. With this in mind, we notice that LCC, exchanger, SERCA and buffer characteristic time scales are typically faster than T min . Recovery times of the RyR2 are probably around 100 ms, safely below T min for most species.
Thus, we claim that it is possible to do a separation of time scales, such that the dynamics of the homeostatic process occurs in a slow manifold, where the dynamics of the fast variables slaves to that of the slow variables (or order parameters, in synergetic terminology [39]). We, thus, divide the variables into fast, ϕ fast , and slow, ϕ slow . The latter equilibrate in a time scale of several beats. At each individual beat, however, the former attain their equilibrium values, that will depend on the state of the slow variables, so � n fast ¼ hð� n slow Þ. This leads us to the conclusion that the surface which determines homeostasis suffers a huge effective reduction of its dimensionality since memory is mainly eliminated by the external pacing. Flows in one beat are then determined only by the diastolic free SR concentration c sr and cytosolic concentrations c i , that act as slow variables, such that one can compute effective maps for the total amount of calcium in the cell and in the SR, as: Maps as the ones indicated in Eq (16) for the concentration of calcium in the SR have been already constructed and used, for instance, to study the stability of the dynamics to alternans [37,40,41]. However, in previous work, this map does not include the dependence on cytosolic calcium concentrations since they do not consider homeostatic equilibrium. The models have one global variable and one map to equilibrate. In order to treat the homeostatic problem, we need to supplement it with the map for the global calcium level in Eq (15) and clarify the double dependence on cytosolic and SR calcium levels.
Then, at steady state, we have: These are two equilibrium conditions for two independent global average variables. If a steady state exists, it corresponds to the simultaneous fulfillment of the two global equilibrium conditions. Alternatively, the solution to these equations may either not exist or, if it does, be unstable, corresponding to, for instance, alternans or chaotic behavior. In this paper, we address the first scenario with a stable fixed point.
We provide further information regarding these maps in the Supporting information S1 File where we show that the total calcium concentrations in the cell and the SR could also be treated as the fundamental analytical variables instead of free calcium concentrations. The reason is that end-diastolic free calcium concentrations can be related to total calcium concentrations assuming the equilibrium concentration of calcium bound to buffers. So, even though in the mathematical model we only use the fast buffering approximation for calsequestrin, the other buffers equilibrate so fast that they attain the equilibrium values by the end of each period.

Validation
From Eqs (17) and (18) we can predict the average diastolic calcium level at steady-state from measurements away from steady-state. Figs 3 and 4 sketch the basic procedure to predict the steady state from single beat measurements. We simulate the evolution of calcium concentrations, c i and c sr , during one single beat, taking different initial calcium diastolic values (i.e. different initial conditions), with all the local variables unrelated to calcium concentrations set initially at their fast-equilibrium approximation. From this evolution we can compute how much calcium leaves the SR, ΔQ rel , or enters the SR, ΔQ up , during this beat, and repeat the calculation for a different initial condition. This gives new values of ΔQ rel and ΔQ up . If we repeat this procedure for multiple combinations of initial values of c i and c sr we can reconstruct the dependence of the currents with the diastolic calcium concentrations. We obtain the functions ΔQ rel (c sr , c i ) and ΔQ up (c sr , c i ), which are nothing else than surfaces in a diastolic calcium concentration space as shown in Fig 3. It is rather intuitive that release and uptake depend on cytosolic and SR calcium concentration at the beginning of one beat. Less intuitive is that intake and extrusion also depends on c sr even though the L-type Calcium current and exchanger expressions do not depend instantaneously on it. The reason for this strong dependence on the SR calcium concentration at the beginning of the beat is because a larger or smaller initial calcium SR load leads to larger or smaller cytosolic calcium transients during the beat. A different calcium transient results, in general, in a different inactivation in the LCC and NCX extrusion during that beat.
The line where both surfaces cross (see ΔQ rel = ΔQ up in Fig 3), or g-nullcline, fulfills the partial equilibrium indicated by Eq (17). To obtain the partial equilibrium given by Eq (18) we construct the surfaces for ΔQ in (c sr , c i ) and ΔQ out (c sr , c i ) using the same procedure indicated in Fig 3. The line where they cross is the f-nullcline, which represents the equilibrium condition that at steady state the same amount of calcium that enters the cell must leave it. The point where both nullclines intersect gives us the concentration of diastolic calcium that fulfills both equilibrium conditions. This sets the steady state, i.e., it predicts the diastolic homeostatic calcium level and its internal distribution between the SR and cytosol of the cell.  (1): The first two surfaces indicate the total amount of calcium entering ΔQ in and leaving ΔQ out the cell, computed following the same procedure explained in Fig 3. Below we compute the surfaces corresponding to release and uptake of calcium from the SR. On the right, we plot the two nullclines that correspond to the crossing of the previous surfaces. The f-nullcline is the set of initial c i and c sr where ΔQ in = ΔQ out while the g-nullcline corresponds to the crossing between release ΔQ rel and uptake ΔQ up . The point where both nullclines cross gives the steady state of the system (c i , c SR ) where concentrations return to the same pre-systolic values after a stimulation. Notice that the parameters of the model fix this crossing precisely at the expected value of roughly 150nM in the cytosol and 40 μmol/Lcyt in the SR, which correspond to a local concentration of free calcium in the SR c sr at roughly 500 μM, given the volume ratio between SR and cytosol. We notice that this free calcium concentration corresponds to 100 μmol/Lcyt of total calcium in the SR, free and bound to CSQ, as reported in [42]. Panel (2): Steady state from single beat measurements at 2 Hz described in the top panel with modification of the rabbit model where the conductivity of LCC is increased. The prediction of thesteady state obtained from the global equilibrium as a function of the LCC conductivity fits very well the computed steady-state obtained after letting the system evolve for 100 beats. We have checked that the prediction given by the crossing of the nullclines agrees with the calcium concentration values obtained letting the system evolve to steady-state (bottom panel of Fig 4). This agreement confirms the validity of the reduction in dimensionality. This reduction is not at all trivial. It comes from the fact that average values are not very sensitive to the inherent noise of the system in ventricle myocytes and that the processes involved in calcium handling have time scales that are faster than the pacing period. Regarding calcium homeostasis, each period of the system provides a global reset of the fast variables leaving only the slow concentration variables in play.

General shocks
Our key result has deep implications for the analysis of cardiomyocyte contractibility. We can use this framework to infer how the cell behavior will be modified upon a change in channel properties or the increase of buffers levels. More specficially, we can use the insights from bidimensional general equilibrium problems from other fields and understand the counterintuitive response of cells to changes at the channel level, as described in the introduction. We will focus our discussion on changes in SERCA function given its possible relevance to understand gene therapy failure. We will show that, in terms of shocks, as they are normally called in the literature, a change of SERCA strength is actually a double shock. With this realization in mind we proceed to unveil and discuss a physiological mechanism for SERCA gene therapy failure where, even in the case of SERCA improved function, calcium transient and contractibility are not improved.
Analysis of stability based on partial equilibrium conditions is common in other fields, such as macroeconomics. In fact, this framework can be mapped one-to-one to the Investment-Saving, Liquidity preference-Money supply (IS-LM) model of macroeconomics developed by J. Hicks [2] as an explanation of Keynes general theory of Employment, Interest and Money [43]. In the IS-LM, the steady-state of a large complex closed economy is defined by its ouput (average GDP per capita) and by the average interest rate, just like in our case we define our steady state by two variables, the average calcium concentration in the cytosol and the SR. Equally, the general equilibrium conditions that must be met are two. First, Investment has to balance Saving, leading to the IS curve, equal to our f-nullcline. The equivalent of our g-nullcline involves the interest rates and output for which the money market is in equilibrium.
There, a shock is a sudden change in the curves that determine equilibrium. In our case, a change in the conductivity or the property of a particular receptor or the strength of SERCA/ NCX is a shock that changes the four surfaces corresponding to the four fluxes. If the shock were only to affect strongly the surfaces related to the release and uptake of calcium from the SR (as sketched in the first panel of Fig 5), then just the g-nullcline would shift position. From this shift and the slope of the f-nullcline we can predict the effect of the shock. For example, if the uptake of calcium is raised (increasing the ΔQ up surface) then the g-nullcline will shift down. Given that, in this example, the slope of the f-nullcline is negative, the downward shift of the g-nullcline will move the crossing of the two nullclines, and the steady state, to higher c sr and lower c i .
However, as we proceed to describe in the next sections, changes in key functioning elements of calcium homeostasis normally affect two or all four of the surfaces. To be specific, we are going to consider a change in SERCA function.

Consequences of shocks
We aim to understand how an increase in SERCA function affects the homeostatic state and whether it systematically leads to improved calcium transients and contraction. We will model We show a hypothetical change in the cell that only affects one of the four surfaces. This produces a change in one of the nullclines. As an example, we pick an increase of the ΔQ up surface. This leads to a shift of the g-nullcline to the right compared with the previous case, resulting in a lower diastolic level in the cytosol but larger in the SR. Panel (2): Effect of changing the maximum uptake of the SERCA pump from 0.15 μ M/ ms to 0.3 μ M/ms in our rabbit model. Increasing the function of SERCA does not only affect the g-nullcline, expected since the uptake is larger, but also affects how well the exchanger works shifting also the f-nullcline to the right. The end result is similar to the previous case.
https://doi.org/10.1371/journal.pcbi.1007572.g005 this increase in SERCA function with an increase in its maximum uptake rate v up . This change mimics increases/decreases in the number of pumps present in each of the CaRUs. This is related to the number of pumps expressed in the SR membrane and should increase if SERCA therapy works properly. We proceed to answer whether a change of SERCA function always leads to larger/broader calcium transients.
Change in SERCA function is a double shock. We observe that changes in v up affect not only the surface directly related to the uptake ΔQ up , but also ΔQ out . Thus, a change in SERCA function is a double shock, since it changes two nullclines, as shown in the bottom panel of Fig  5. The reason is rather intuitive. A change in the uptake does not have a strong effect in the behavior of the LCC channels and hardly changes the release, but a more efficient SERCA reduces quickly cytosolic calcium levels, which reduces the efficiency of the exchanger. The effect is that the shock associated with larger uptake moves both nullclines to the right. Given the structure of the nullclines, shown in Fig 5, where the f-nullcline has a negative slope while that of the g-nullcline is positive, it is thus straightforward that they will cross at a higher level of SR. In fact, we observe that the level of calcium in the SR increases monotonically, while the level in the cytosol decreases also monotonically (left of the first panel of Fig 6). In the Supporting information S1 File, we show that the total amount of calcium in the cell, this is, the sum of calcium in the cytosol and SR, increases slightly as we increase the SERCA uptake given that the increase in calcium in the SR compensates a reduction of calcium levels in the cytosol. While the increase in SR load depends on the slope of the nullclines, whether pre-systolic calcium in the cytosol increases or decreases depends on how far both nullclines move with the shock. Thus, the fact that the nullclines cross at lower cytosolic calcium is not a structural result but it depends on the particulars of the system.
The shift of the nullclines due to the increase in SERCA function is a quantification of how SERCA interacts with the exchanger and, indirectly, via changes in the release, with the LCC. A faster SERCA uptake reduces the ability of the exchanger to extrude calcium since SERCA pushes down calcium in the cytosol faster during diastole, making the gradient of calcium between the cytosol and extracellular space larger (see Fig 7). Given that the exchanger does not work as well with a large gradient in calcium, one could think that it will extrude less calcium. However, this cannot be the case. Since ΔQ in is a rather flat surface, calcium entering via LCC, as SERCA increases, remains roughly unchanged. Thus, to reach equilibrium, calcium extruded by the exchanger has to remain almost constant too. The only possibility is that the system reaches a different homeostatic equilibrium in which the release and calcium transient adapt so that the exchanger can extrude the same as it did before the shock (Fig 7).
The resulting changes in the calcium transient when we increase the maximum SERCA uptake at steady-state are summarized on the right graph of the first panel in Fig 6 and explained schematically in Fig 7. The transient amplitude must change if diastolic cytosolic calcium decreases. We observe that the peak of calcium transient barely changes, but given that the minimum of the transient decreases, the amplitude raises appreciably. The end result is that in this animal model, a higher v up increases the contractibility of the cell, given the larger transient.
Physiological mechanism for SERCA gene therapy failure. Now we aim to answer a simple question. Will a change of SERCA function always lead to larger/broader calcium transients? For that, we will reproduce a possible malfunction. Given that there are calcium-binding proteins that affect both SERCA and RyR [19], we consider a rabbit cell where the conductivity of the RyR is reduced without changes in the exchanger or the LCC. The details of the rabbit cell with reduced RyR activity are explained in S1 File where we show that the slopes of both nullclines are still the same as in the original rabbit model. However, as the second panel of Fig 6 shows, there are crucial differences with the wild-type scenario, mainly that the calcium transient is actually reduced leading to lower contraction when SERCA maximum uptake is increased.
Given that the structure of the nullcline is the same, increasing the SERCA uptake in the RyR2 modified rabbit cell leads to higher diastolic SR load (see Fig 6). However, as shown in the left graph of the second panel in Fig 6, higher v up leads to an increase of the diastolic  (2): We present the same graphs as in the first panel but for a cell where the single channel conductivity of the RyR has been reduced. In this situation, the calcium transient is always lower than before and it does not improve as the uptake of SERCA is larger.
https://doi.org/10.1371/journal.pcbi.1007572.g006 cytosolic calcium. The reason is that, in this case, the shift in the f-nullcline is larger than in the g-nullcline (see S1 File for details). This increased level of calcium makes the exchanger work more efficiently. On the other hand, the surface ΔQ in is very flat, so the calcium influx remains almost constant. This means that the exchanger has to extrude roughly the same amount of calcium no matter what the SERCA uptake is. This directly leads to a decrease of the transient, as shown at the graph on the right of the second panel in Fig 6 and explained schematically in Fig 7. Actually, both the minimum is increased and the maximum is reduced at steady-state, leading to a strongly reduced transient, exactly the opposite of the wild-type case.
The sketch of the physiological mechanism for SERCA gene therapy failure is described in detail in Fig 7. Upon an increase in SERCA uptake, initially calcium intake via LCC and release via RyR remain unaffected. However, the NCX will be affected by a larger gradient in calcium leading to a reduction in its function. The system is no longer in equilibrium and must rebalance loading calcium (Fig 7A). This change in calcium levels affects all elements of calcium handling, including those that were not initially affected: intake and release. The specifics of how this increase in calcium is distributed depends on the particulars of the general equilibrium state of the cell. The slope of the nullclines will determine if SR will be more loaded as it is normally the case. The effect of the shock on the nullclines will determine whether cytosolic calcium is increased, because the cell will load until NCX function is restored thanks to a reduced gradient, or decreased, in case the NCX function is restored thanks to a larger gene therapy a sudden increase in maximum SERCA uptake does not change initially the LCC or the release but it decreases NCX function as it decreases cytosolic calcium. This necessarily leads to a state of calcium unbalance that must be rebalanced again. B) Homeostasis will always fix and rebalance the system. However, the particulars of how this is done are not universal. They depend on the nullcline structure and reaction to shocks. For instance, in a healthy rabbit, the cell increases total calcium content by increasing considerably the calcium in the SR while diminishing cytosolic calcium levels. This increase in SR calcium results in a higher transient that allows a higher extraction of calcium through the NCX, even if diastolic calcium is reduced. However, it can happen that the RyR has a weak sensitivity to SR calcium content, so release remains almost constant. Since uptake is increased, the diastolic cytosolic level has to be increased to allow calcium extrusion through the NCX, then the transient is decreased, which results in decreased contractibility despite a stronger SERCA.
https://doi.org/10.1371/journal.pcbi.1007572.g007 transient provided by a more loaded SR (Fig 7B). Both are perfectly possible outcomes of the shifts in the nullclines given the shock given to the cell.

Discussion and conclusion
Detailed models of calcium handling have been able to reproduce spatiotemporally complex behavior such as discordant alternans [44] or calcium waves [45,46], and have been used, for instance, to understand the onset of alternans as a synchronization transition [30]. Surprisingly given the complexity of the dynamics, the equilibrium conditions can be obtained using, as in very simple models, average variables and balance of fluxes in the cytosol and the SR. This result is not trivial. It requires that all internal variables, except cytosolic and SR calcium, equilibrate at the time scale of the pacing period. Also, that all CaRUs behave on average as the averaged variable, to prevent the appearance of intracellular inhomogeneities with different behavior.
As we have discussed earlier in this paper, this dimensionality reduction is not uncommon in other complex systems, and has been used previously to study the appearance of alternans in calcium cycling [37,40,41]. We understand the equilibrium point as satisfying two partial equilibrium conditions: calcium in and out of the cell and in and out of the SR must balance. We can then use this description to understand the effects of shocks, i.e., changes in parameters that result in shifts in the curves denoting the partial equilibrium conditions. We observe that the effect of the shocks can be very different depending on the form of these curves. This is, the same recipe may have very different outcomes depending on the state of the system.
In this paper, we have considered a clamped AP because the APD generally adapts to changes in currents and clamping does not affect this main insight. If we were to take the proper shape of the APD for a given frequency as our clamped AP in the model, the nullclines depicted here would not be affected. However, for future work, it would be interesting to study the interplay of calcium homeostasis with changes in the AP at different frequencies in order to address other cardiac properties such as the structure of force-frequency relations in different animal models. Since most of the adaptation of the AP is fast, we expect it to be generally slaved to the dynamics of calcium, at least, at the time scale of seconds. In this respect, a full linear stability analysis of the periodic transient, with an explicit calculation of the most unstable (or less stable) eigenvalues and eigenmodes [47,48] would help to validate this point. However, there are slower times scales, associated with the long term accumulation of ions, that will become another dynamical variable, increasing the complexity of the general equilibrium problem. More specifically, one should expect the slow change of potassium and sodium concentrations to affect the LCC and NCX, which in turn will affect the nullclines which will feedback into the APD and back again into the ionic balance. The general structure of our approach will hold, but a full analysis of the relevant slow variables in the full model will be needed.
The analysis of a ventricular model in terms of shocks is thus robust, and we expect it to be useful to understand how these shocks may produce counterintuitive results. As, for instance, a decrease in SR calcium concentration as the strength of the LCC channels is increased. This can be well understood with the shifts in the partial equilibrium curves. Clearly, a further investigation of what determines the slope of the curves and how they depend on species or underlying pathologies should be undertaken.
Another relevant example is the effect of a change in SERCA. Upregulation of SERCA gene expression, or SERCA gene therapy, has been used to treat patients with heart failure. Even if this therapy is often successful, in some occasions it does not seem to work [22,23]. From our results, one could think of a possible explanation as sketched in Fig 7. Under some modifications on the parameters, an increase of SERCA strength does not lead to an increase in the cytosolic calcium transient, but rather the opposite. This is an in-silico failure of SERCA therapy due to how the cell reacts to the double shock that the changes in SERCA represent. We have thus unveiled a possible physiological mechanism for SERCA therapy failure and how its success or failure depends on the basic structure and movement of the nullclines determining homeostatic steady-state upon changes in SERCA function.
Supporting information S1 File. Supporting information file with appendix. In this supplemental information file, we give a full description of the rabbit model, develop the nullcline analysis, mentioned in the results section, using total calcium concentrations in the cell instead of free calcium concentrations. We also provide additional details on how SERCA uptake changes affect cell behavior concerning the physiological mechanism for SERCA gene therapy failure, as indicated in the discussion section. This supplemental file also includes an appendix with tables of the different parameters used in the models. (PDF)