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

David Conesa; Blas Echebarria; Angelina Peñaranda; Inmaculada R. Cantalapiedra; Yohannes Shiferaw; Enrique Alvarez-Lacalle

doi:10.1371/journal.pcbi.1007572

Abstract

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.

Author summary

Cardiomyocytes, upon voltage excitation, release calcium, which leads to cell contraction. However, under some pathological conditions, calcium handling is impaired. Recently, SERCA gene therapy, whose aim is to improve Ca²⁺ sequestration by the Sarcoplasmic Reticulum (SR), has failed to improve the prognosis of patients with Heart Failure. This, together with recent counterintuitive results in calcium handling, has highlighted the need for a framework to understand calcium homeostasis across species and pathologies. We show here that the proper framework is a general equilibrium approach of two independent variables. The development of this framework allows us to find a possible mechanism for the failure of SERCA gene therapy even when it manages to increase Ca SERCA expression.

Citation: Conesa D, Echebarria B, Peñaranda A, Cantalapiedra IR, Shiferaw Y, Alvarez-Lacalle E (2020) Two-variable nullcline analysis of ionic general equilibrium predicts calcium homeostasis in ventricular myocytes. PLoS Comput Biol 16(6): e1007572. https://doi.org/10.1371/journal.pcbi.1007572

Editor: Andrew D. McCulloch, University of California San Diego, UNITED STATES

Received: November 27, 2019; Accepted: May 5, 2020; Published: June 5, 2020

Copyright: © 2020 Conesa et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the manuscript and its Supporting Information files.

Funding: E.A-L, B.E, A. P, I.R.C and D. C received grants from Fundació La Marató de TV3 in 2014, under grant number 20151110, URL:http://aquas.gencat.cat/web/.content/minisite/aquas/ambits/recerca/convocatories_ex_ante/la_marato_tv3/maratotv3_2014_cor.pdf They also received grants from the Spanish Ministerio de Economía, Industria y Competitividad (MINECO) under grant number SAF2017-88019-C3-2-R and later A.P received grant. rom Ministerio de Ciencia Innovación y Universidades grant PGC2018-095456-B-I00. URL: https://sede.micinn.gob.es/stfls/eSede/Ficheros/2017/Propuesta_Resolucion_Provisional_Proyectos_Retos_2017.pdf URL:https://sede.micinn.gob.es/stfls/eSede/Ficheros/2019/Segunda_Propuesta_Resolucion_Definitiva_PGC_2018.pdf Y.S. received grant R01-119095 from National Heart, Lung and Blood Institute (NHLBI) URL: https://www.nhlbi.nih.gov/ Y.S.; B. E., and E.A-L acknowledge also that part of this work was completed in the Kavli Institute of Theoretical Physics which sponsored the summer program Integrative Cardiac Dynamics under grants NSF PHY-1748958, NIH- R25GM067110, and by the Gordon and Betty Moore Foundation (Grant No. 2919.01) https://www.kitp.ucsb.edu/activities/cardio18 The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Complex systems often present robust regulation that makes them resilient to changes in environmental conditions. Despite the presence of a large number of nonlinearly interacting elements, averaged quantities in these systems often attain predictable, and sometimes constant, values. Examples range from ecosystems and climate regulation to long-run macroeconomics [1, 2]. Our bodies present a good example of such behavior. The maintenance of homeostatic equilibrium in temperature, pH, osmolality or ionic concentration is crucial to sustaining life, and complex feedback mechanisms have been developed to achieve such goal.

Another important example is the regulation of heart function. Fundamental for survival is the ability of the heart to contract and expel more blood by increasing the heart rate [3]. In most animal species the amount of blood pumped at each beat increases rather strongly with beat rate, which stems from a positive correlation between contractile force and beating frequency [4–6]. However, the variability of the force-frequency relation is quite large across species, being quite flat for some, like mouse and rat [7, 8]. As the contraction machinery in the sarcomeres is activated by an increase in intracellular calcium, this variability is determined by differences in calcium handling [9–11]. During the action potential of the heart, calcium enters into the cell via L-type calcium channels (LCC) triggering the release of calcium stored in the sarcoplasmic reticulum (SR). Each time voltage raises, the intake of calcium due to the opening of LCC raises the probability of RyR2 to open, releasing calcium from the SR during the depolarization of voltage (systole)—the well-known Calcium-Induced Calcium-release (CICR) mechanism—which then binds to TnC, triggering contraction. Once the voltage is repolarized, calcium released from the SR is uptaken by SERCA back to the SR while the sodium-calcium exchanger (NCX) is able to extrude the calcium which entered via LCC, and contraction ceases. This process seems universal across species, although quantitatively can be very different. Not just the release of calcium is species-dependent, its reuptake into the SR by the SR Ca² ATPase (SERCA) and the extrusion of calcium from the cell by the sodium-calcium exchanger also differs. Translation of results obtained in animal models targeting different regulatory pathways into human patients is thus not trivial.

To reach homeostatic equilibrium, the amount of calcium entering the cell must be compensated by the quantity extruded [12, 13]. Equally, the amount of calcium released from the SR during the transient must be the same reuptaken by SERCA [12, 14]. Global equilibrium is thus determined by a balance of global fluxes, that encompass the complex internal machinery of calcium release, composed of thousands of interacting release units. A good understanding of calcium homeostasis is necessary to understand the effect of changes in calcium handling regulatory mechanisms. In particular, dysregulations in calcium homeostasis have important implications for arrhythmias since calcium overload of the SR often facilitates the initiation of calcium alternans, calcium waves, and arrhythmogenic EADs [13, 15–17].

Altered expression and/or functionality of proteins involved in calcium dysregulation, has often been treated with therapeutic strategies such as gene therapy [18]. In particular, SERCA gene therapy has been successfully applied to preclinical treatment of heart failure (HF) in animals, showing an increase in SERCA2a expression and improved cardiac function [19, 20]. In humans, despite early promising results [21, 22], more recent studies with wider samples of 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.

Download:

Fig 1. Structure of the computational ventricular rabbit cardiomyocyte model.

Panel (1): Scheme of the three-dimensional model where the cardiomyocyte is divided into thousands of 3d stacked Calcium Release Units (highlighted one in red). Panel (2): Depiction of the basic structure of the CaRUs showing the L-type Calcium Channels (LCC) in the vicinity of the Ryanodine Receptor (RyR), the SERCA pump and the Na-Ca Exchanger. The different states of the LCC and the RyR2 are also depicted. Panel (3): Typical cytosolic (orange) and SR (blue) free calcium transients in the rabbit ventricular model (described in Methods) at a voltage-clamped pacing rate of 2 Hz.

https://doi.org/10.1371/journal.pcbi.1007572.g001

Download:

Fig 2. Example of counterintuitive calcium homeostatic regulations.

Panel (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.

https://doi.org/10.1371/journal.pcbi.1007572.g002

Materials and methods

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³ 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 average distance between large clusters of RyR (40) or the aggregation scales of smaller ones that we take to be 500 nm [33, 34].

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²⁺ 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²⁺ 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: (1) (2) (3) (4) (5) 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³) 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²⁺ 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²⁺ exchanger depends on intracellular and extracellular Na⁺ concentrations ([Na]_i and [Na]_o respectively), which are set as constant in voltage-clamped model, as well as on extracellular Ca²⁺ 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²⁺ 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²⁺ concentrations in cytosol and network SR. We give the full expression here given the relevance for our analysis of SERCA gene therapy. (6) 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 , 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 (7) 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: (8)

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 bi-dimensional 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 , 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: (9) 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: (10) 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. (11) Notice that both ΔQ_in and ΔQ_out depend, in principle, on the value of the roughly one hundred thousand internal variables (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 at beat n follows the relation so that at steady state (ss) we have: (12) 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: (13) 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 (14)

As stated before, ventricular myocytes reproduce consistent average behavior. Therefore, fluxes in Eqs (12) and (14) depend on the relevant spatial average variables 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 . 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: (15) (16)

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: (17) (18)

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.

Results

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.

Download:

Fig 3. Schematics of the procedure to reconstruct one of the nullclines of the system.

The nullcline corresponds to the partial equilibrium where the intake of calcium into the SR equals its release, where we have multiplied all ΔQ by the volume of the cytosol, to obtain this value in femtomols (fmol). The main central graph reproduces total release ΔQ_rel and uptake ΔQ_up as a function of the initial free concentrations in the cytosol c_i(t = 0) and SR c_SR(t = 0). The nullcline is determined by the line where both surfaces cross. For each initial condition, one single transient is simulated and the release J_RyR and uptake J_SrCa computed and shown in the upper graphs. The integrated values are indicated in the bar graph below and then placed as elements of the surface. The whole surface is constructed by reproducing this procedure with multiple different initial calcium concentrations where all other initial variables are in their fast equilibrium approximation.

https://doi.org/10.1371/journal.pcbi.1007572.g003

Download:

Fig 4. General equilibrium. Prediction of steady-state from the nullclines crossing.

Panel (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.

https://doi.org/10.1371/journal.pcbi.1007572.g004

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.

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 bi-dimensional 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.

Download:

Fig 5. Shocks in general equilibrium models. Change in SERCA uptake is a two-nullcline shock.

Panel (1): 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

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

Download:

Fig 6. Homeostatic response of increasing SERCA uptake with normal and low RyR channel conductivity.

Panel (1): Dependence of pre-systolic values of free calcium concentrations in the SR and in the cytosol at steady state with maximum uptake of SERCA pump v_up in the rabbit model. On the right, the relative increase of calcium during the calcium transient as a function of the maximum SERCA uptake ν_up. We take 100% to be the transient with the standard ν_up = 0.5 μM/ms. Panel (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

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

Download:

Fig 7. Schematics of possible SERCA gene therapy outcomes.

Schematics of how SERCA gene therapy may fail in certain cells. A) During SERCA 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

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

https://doi.org/10.1371/journal.pcbi.1007572.s001

(PDF)

Acknowledgments

We want to thank L.Hove-Madsen and N.Otani for fruitful discussions.

References

1. Munasinghe M. Making economic growth more sustainable. Ecological Economics. 1995;15(2):121–124.
- View Article
- Google Scholar
2. Hicks JR. Mr. Keynes and the “classics”; a suggested interpretation. Econometrica: journal of the Econometric Society. 1937;p. 147–159.
- View Article
- Google Scholar
3. Bers D. Excitation-contraction coupling and cardiac contractile force. vol. 237. Springer Science & Business Media; 2001.
4. Endoh M. Force–frequency relationship in intact mammalian ventricular myocardium: physiological and pathophysiological relevance. European journal of pharmacology. 2004;500(1-3):73–86. pmid:15464022
- View Article
- PubMed/NCBI
- Google Scholar
5. Gwathmey J, Slawsky M, Hajjar R, Briggs G, Morgan J. Role of intracellular calcium handling in force-interval relationships of human ventricular myocardium. The Journal of clinical investigation. 1990;85(5):1599–1613. pmid:2332508
- View Article
- PubMed/NCBI
- Google Scholar
6. Janssen PM, Periasamy M. Determinants of frequency-dependent contraction and relaxation of mammalian myocardium. Journal of molecular and cellular cardiology. 2007;43(5):523–531. pmid:17919652
- View Article
- PubMed/NCBI
- Google Scholar
7. Antoons G, Mubagwa K, Nevelsteen I, Sipido KR. Mechanisms underlying the frequency dependence of contraction and [Ca2+] i transients in mouse ventricular myocytes. The Journal of physiology. 2002;543(3):889–898. pmid:12231646
- View Article
- PubMed/NCBI
- Google Scholar
8. Georgakopoulos D, Kass DA. Minimal force-frequency modulation of inotropy and relaxation of in situ murine heart. The Journal of Physiology. 2001;534(2):535–545. pmid:11454970
- View Article
- PubMed/NCBI
- Google Scholar
9. Bers DM. Calcium fluxes involved in control of cardiac myocyte contraction. Circulation research. 2000;87(4):275–281. pmid:10948060
- View Article
- PubMed/NCBI
- Google Scholar
10. Gattoni S, Røe ÅT, Frisk M, Louch WE, Niederer SA, Smith NP. The calcium–frequency response in the rat ventricular myocyte: an experimental and modelling study. The Journal of physiology. 2016;594(15):4193–4224. pmid:26916026
- View Article
- PubMed/NCBI
- Google Scholar
11. Maier LS, Bers DM, Pieske B. Differences in Ca2+-handling and sarcoplasmic reticulum Ca2+-content in isolated rat and rabbit myocardium. Journal of molecular and cellular cardiology. 2000;32(12):2249–2258. pmid:11113000
- View Article
- PubMed/NCBI
- Google Scholar
12. Eisner D, Bode E, Venetucci L, Trafford A. Calcium flux balance in the heart. Journal of molecular and cellular cardiology. 2013;58:110–117. pmid:23220128
- View Article
- PubMed/NCBI
- Google Scholar
13. Lugo CA, Cantalapiedra IR, Peñaranda A, Hove-Madsen L, Echebarria B. Are SR Ca content fluctuations or SR refractoriness the key to atrial cardiac alternans?: insights from a human atrial model. American Journal of Physiology-Heart and Circulatory Physiology. 2014;306(11):H1540–H1552. pmid:24610921
- View Article
- PubMed/NCBI
- Google Scholar
14. Eisner DA, Caldwell JL, Kistamás K, Trafford AW. Calcium and excitation-contraction coupling in the heart. Circulation research. 2017;121(2):181–195. pmid:28684623
- View Article
- PubMed/NCBI
- Google Scholar
15. Landstrom AP, Dobrev D, Wehrens XH. Calcium signaling and cardiac arrhythmias. Circulation research. 2017;120(12):1969–1993. pmid:28596175
- View Article
- PubMed/NCBI
- Google Scholar
16. Němec J, Kim JJ, Salama G. The link between abnormal calcium handling and electrical instability in acquired long QT syndrome–does calcium precipitate arrhythmic storms? Progress in biophysics and molecular biology. 2016;120(1-3):210–221. pmid:26631594
- View Article
- PubMed/NCBI
- Google Scholar
17. Alvarez-Lacalle E, Cantalapiedra IR, Peñaranda A, Cinca J, Hove-Madsen L, Echebarria B. Dependency of calcium alternans on ryanodine receptor refractoriness. PloS one. 2013;8(2):e55042. pmid:23390511
- View Article
- PubMed/NCBI
- Google Scholar
18. Chamberlain K, Riyad JM, Weber T. Cardiac gene therapy with adeno-associated virus-based vectors. Current opinion in cardiology. 2017. pmid:28169951
- View Article
- PubMed/NCBI
- Google Scholar
19. Gorski PA, Ceholski DK, Hajjar RJ. Altered myocardial calcium cycling and energetics in heart failure-a rational approach for disease treatment. Cell metabolism. 2015;21(2):183–194. pmid:25651173
- View Article
- PubMed/NCBI
- Google Scholar
20. Samuel T, Rosenberry R, Lee S, Pan Z. Correcting Calcium Dysregulation in Chronic Heart Failure Using SERCA2a Gene Therapy. International journal of molecular sciences. 2018;19(4):1086.
- View Article
- Google Scholar
21. Jaski BE, Jessup ML, Mancini DM, Cappola TP, Pauly DF, Greenberg B, et al. Calcium upregulation by percutaneous administration of gene therapy in cardiac disease (CUPID Trial), a first-in-human phase 1/2 clinical trial. Journal of cardiac failure. 2009;15(3):171–181. pmid:19327618
- View Article
- PubMed/NCBI
- Google Scholar
22. Jessup M, Greenberg B, Mancini D, Cappola T, Pauly DF, Jaski B, et al. Calcium upregulation by percutaneous administration of gene therapy in cardiac disease (CUPID) a phase 2 trial of intracoronary gene therapy of sarcoplasmic reticulum Ca2+-ATPase in patients with advanced heart failure. Circulation. 2011;124(3):304–313. pmid:21709064
- View Article
- PubMed/NCBI
- Google Scholar
23. Hulot JS, Ishikawa K, Hajjar RJ. Gene therapy for the treatment of heart failure: promise postponed. European heart journal. 2016;37(21):1651–1658. pmid:26922809
- View Article
- PubMed/NCBI
- Google Scholar
24. Hamilton S, Terentyev D. Altered Intracellular Calcium Homeostasis and Arrhythmogenesis in the Aged Heart. International journal of molecular sciences. 2019;20(10):2386.
- View Article
- Google Scholar
25. Eisner D. A., Caldwell J. L., Trafford A.W., and Hutchings D.C. The Control of Diastolic Calcium in the Heart: Basic Mechanisms and Functional Implications. Circulation Research. 2020; 126(3):395–412. pmid:31999537
- View Article
- PubMed/NCBI
- Google Scholar
26. Trafford A, Diaz M, Eisner D. Coordinated control of cell Ca2+ loading and triggered release from the sarcoplasmic reticulum underlies the rapid inotropic response to increased L-type Ca2+ current. Circulation research. 2001;88(2):195–201. pmid:11157672
- View Article
- PubMed/NCBI
- Google Scholar
27. Mahajan A, Shiferaw Y, Sato D, Baher A, Olcese R, Xie LH, et al. A rabbit ventricular action potential model replicating cardiac dynamics at rapid heart rates. Biophysical journal. 2008;94(2):392–410. pmid:18160660
- View Article
- PubMed/NCBI
- Google Scholar
28. Stern M. D., Song L. S., Cheng H., Sham J. S., Yang H. T., Boheler K. R., and Ríos E. Local control models of cardiac excitation–contraction coupling: a possible role for allosteric interactions between ryanodine receptors. The Journal of general physiology. 1999;113(3):469–489. pmid:10051521
- View Article
- PubMed/NCBI
- Google Scholar
29. Cantalapiedra I. R., Alvarez-Lacalle E., Peñaranda A., Echebarria B. Minimal model for calcium alternans due to SR release refractoriness. Chaos: An Interdisciplinary Journal of Nonlinear Science. 2017; 27(9): 093928.
- View Article
- Google Scholar
30. Alvarez-Lacalle E, Echebarria B, Spalding J, Shiferaw Y. Calcium alternans is due to an order-disorder phase transition in cardiac cells. Physical review letters. 2015;114(10):108101. pmid:25815968
- View Article
- PubMed/NCBI
- Google Scholar
31. Galice S, Xie Y, Yang Y, Sato D, Bers DM. Size matters: Ryanodine receptor cluster size affects arrhythmogenic sarcoplasmic reticulum calcium release. Journal of the American Heart Association. 2018;7(13):e008724. pmid:29929992
- View Article
- PubMed/NCBI
- Google Scholar
32. Rog-Zielinska EA, Johnston CM, O’Toole ET, Morphew M, Hoenger A, Kohl P. Electron tomography of rabbit cardiomyocyte three-dimensional ultrastructure. Progress in biophysics and molecular biology. 2016;121(2):77–84. pmid:27210305
- View Article
- PubMed/NCBI
- Google Scholar
33. Macquaide N, Tuan HTM, Hotta Ji, Sempels W, Lenaerts I, Holemans P, et al. Ryanodine receptor cluster fragmentation and redistribution in persistent atrial fibrillation enhance calcium release. Cardiovascular research. 2015;108(3):387–398. pmid:26490742
- View Article
- PubMed/NCBI
- Google Scholar
34. Soeller C. Ryanodine receptor cluster size sets the tone in cerebral smooth muscle. Proceedings of the National Academy of Sciences. 2018;115(41):10195–10197.
- View Article
- Google Scholar
35. Wagner J., and Keizer J. Effects of rapid buffers on Ca2+ diffusion and Ca2+ oscillations. Biophysical Journal. 1994;67(1):447–456. pmid:7919018
- View Article
- PubMed/NCBI
- Google Scholar
36. Smith G. D., Wagner J., and Keizer J. Validity of the rapid buffering approximation near a point source of calcium ions. Biophysical Journal. 1996;70(6):2527–2539. pmid:8744292
- View Article
- PubMed/NCBI
- Google Scholar
37. Shiferaw Y., Watanabe M. A., Garfinkel A., Weiss J. N., and Karma A. Model of intracellular calcium cycling in ventricular myocytes. Biophysical journal. 2003;85(6):3666–3686. pmid:14645059
- View Article
- PubMed/NCBI
- Google Scholar
38. Shannon TR, Wang F, Puglisi J, Weber C, Bers DM. A mathematical treatment of integrated Ca dynamics within the ventricular myocyte. Biophysical journal. 2004;87(5):3351–3371. pmid:15347581
- View Article
- PubMed/NCBI
- Google Scholar
39. Haken H. Synergetics: introduction and advanced topics. Springer Science & Business Media; 2013.
40. Qu Z., Liu M. B., and Nivala M. A unified theory of calcium alternans in ventricular myocytes. Scientific reports. 2016;6:35625. pmid:27762397
- View Article
- PubMed/NCBI
- Google Scholar
41. Huertas MA, Smith GD and Györke S. Ca2+ alternans in a cardiac myocyte model that uses moment equations to represent heterogeneous junctional SR Ca2+. Biophysical journal. 2010:99(2):377–387. pmid:20643055
- View Article
- PubMed/NCBI
- Google Scholar
42. Shannon T. R., Ginsburg K. S., and Bers DM. Reverse mode of the sarcoplasmic reticulum calcium pump and load-dependent cytosolic calcium decline in voltage-clamped cardiac ventricular myocytes. Biophysical Journal, 2000, 78.1: 322–333. pmid:10620296
- View Article
- PubMed/NCBI
- Google Scholar
43. Keynes JM. The general theory of employment, interest, and money. Springer; 2018.
44. Sato D, Shiferaw Y, Garfinkel A, Weiss JN, Qu Z, Karma A. Spatially discordant alternans in cardiac tissue: role of calcium cycling. Circulation research. 2006;99(5):520–527. pmid:16902177
- View Article
- PubMed/NCBI
- Google Scholar
45. Shiferaw Y, Aistrup GL, Wasserstrom JA. Intracellular Ca2+ waves, afterdepolarizations, and triggered arrhythmias. Oxford University Press; 2012.
46. Izu LT, Xie Y, Sato D, Bányász T, Chen-Izu Y. Ca2+ waves in the heart. Journal of molecular and cellular cardiology. 2013;58:118–124. pmid:23220129
- View Article
- PubMed/NCBI
- Google Scholar
47. Li M. and Otani N. F. Ion channel basis for alternans and memory in cardiac myocytes. Annals of biomedical engineering. 2003;31(10):1213–1230. pmid:14649495
- View Article
- PubMed/NCBI
- Google Scholar
48. Veasy J., Lai Y. M., Coombes S., and Thul R. Complex patterns of subcellular cardiac alternans. Journal of Theoretical Biology. 2019;478:102–114. pmid:31220466
- View Article
- PubMed/NCBI
- Google Scholar

[ref1] 1. Munasinghe M. Making economic growth more sustainable. Ecological Economics. 1995;15(2):121–124.
View Article
Google Scholar

[2] View Article

[3] Google Scholar

[ref2] 2. Hicks JR. Mr. Keynes and the “classics”; a suggested interpretation. Econometrica: journal of the Econometric Society. 1937;p. 147–159.
View Article
Google Scholar

[5] View Article

[6] Google Scholar

[ref3] 3. Bers D. Excitation-contraction coupling and cardiac contractile force. vol. 237. Springer Science & Business Media; 2001.

[ref4] 4. Endoh M. Force–frequency relationship in intact mammalian ventricular myocardium: physiological and pathophysiological relevance. European journal of pharmacology. 2004;500(1-3):73–86. pmid:15464022
View Article
PubMed/NCBI
Google Scholar

[9] View Article

[10] PubMed/NCBI

[11] Google Scholar

[ref5] 5. Gwathmey J, Slawsky M, Hajjar R, Briggs G, Morgan J. Role of intracellular calcium handling in force-interval relationships of human ventricular myocardium. The Journal of clinical investigation. 1990;85(5):1599–1613. pmid:2332508
View Article
PubMed/NCBI
Google Scholar

[13] View Article

[14] PubMed/NCBI

[15] Google Scholar

[ref6] 6. Janssen PM, Periasamy M. Determinants of frequency-dependent contraction and relaxation of mammalian myocardium. Journal of molecular and cellular cardiology. 2007;43(5):523–531. pmid:17919652
View Article
PubMed/NCBI
Google Scholar

[17] View Article

[18] PubMed/NCBI

[19] Google Scholar

[ref7] 7. Antoons G, Mubagwa K, Nevelsteen I, Sipido KR. Mechanisms underlying the frequency dependence of contraction and [Ca2+] i transients in mouse ventricular myocytes. The Journal of physiology. 2002;543(3):889–898. pmid:12231646
View Article
PubMed/NCBI
Google Scholar

[21] View Article

[22] PubMed/NCBI

[23] Google Scholar

[ref8] 8. Georgakopoulos D, Kass DA. Minimal force-frequency modulation of inotropy and relaxation of in situ murine heart. The Journal of Physiology. 2001;534(2):535–545. pmid:11454970
View Article
PubMed/NCBI
Google Scholar

[25] View Article

[26] PubMed/NCBI

[27] Google Scholar

[ref9] 9. Bers DM. Calcium fluxes involved in control of cardiac myocyte contraction. Circulation research. 2000;87(4):275–281. pmid:10948060
View Article
PubMed/NCBI
Google Scholar

[29] View Article

[30] PubMed/NCBI

[31] Google Scholar

[ref10] 10. Gattoni S, Røe ÅT, Frisk M, Louch WE, Niederer SA, Smith NP. The calcium–frequency response in the rat ventricular myocyte: an experimental and modelling study. The Journal of physiology. 2016;594(15):4193–4224. pmid:26916026
View Article
PubMed/NCBI
Google Scholar

[33] View Article

[34] PubMed/NCBI

[35] Google Scholar

[ref11] 11. Maier LS, Bers DM, Pieske B. Differences in Ca2+-handling and sarcoplasmic reticulum Ca2+-content in isolated rat and rabbit myocardium. Journal of molecular and cellular cardiology. 2000;32(12):2249–2258. pmid:11113000
View Article
PubMed/NCBI
Google Scholar

[37] View Article

[38] PubMed/NCBI

[39] Google Scholar

[ref12] 12. Eisner D, Bode E, Venetucci L, Trafford A. Calcium flux balance in the heart. Journal of molecular and cellular cardiology. 2013;58:110–117. pmid:23220128
View Article
PubMed/NCBI
Google Scholar

[41] View Article

[42] PubMed/NCBI

[43] Google Scholar

[ref13] 13. Lugo CA, Cantalapiedra IR, Peñaranda A, Hove-Madsen L, Echebarria B. Are SR Ca content fluctuations or SR refractoriness the key to atrial cardiac alternans?: insights from a human atrial model. American Journal of Physiology-Heart and Circulatory Physiology. 2014;306(11):H1540–H1552. pmid:24610921
View Article
PubMed/NCBI
Google Scholar

[45] View Article

[46] PubMed/NCBI

[47] Google Scholar

[ref14] 14. Eisner DA, Caldwell JL, Kistamás K, Trafford AW. Calcium and excitation-contraction coupling in the heart. Circulation research. 2017;121(2):181–195. pmid:28684623
View Article
PubMed/NCBI
Google Scholar

[49] View Article

[50] PubMed/NCBI

[51] Google Scholar

[ref15] 15. Landstrom AP, Dobrev D, Wehrens XH. Calcium signaling and cardiac arrhythmias. Circulation research. 2017;120(12):1969–1993. pmid:28596175
View Article
PubMed/NCBI
Google Scholar

[53] View Article

[54] PubMed/NCBI

[55] Google Scholar

[ref16] 16. Němec J, Kim JJ, Salama G. The link between abnormal calcium handling and electrical instability in acquired long QT syndrome–does calcium precipitate arrhythmic storms? Progress in biophysics and molecular biology. 2016;120(1-3):210–221. pmid:26631594
View Article
PubMed/NCBI
Google Scholar

[57] View Article

[58] PubMed/NCBI

[59] Google Scholar

[ref17] 17. Alvarez-Lacalle E, Cantalapiedra IR, Peñaranda A, Cinca J, Hove-Madsen L, Echebarria B. Dependency of calcium alternans on ryanodine receptor refractoriness. PloS one. 2013;8(2):e55042. pmid:23390511
View Article
PubMed/NCBI
Google Scholar

[61] View Article

[62] PubMed/NCBI

[63] Google Scholar

[ref18] 18. Chamberlain K, Riyad JM, Weber T. Cardiac gene therapy with adeno-associated virus-based vectors. Current opinion in cardiology. 2017. pmid:28169951
View Article
PubMed/NCBI
Google Scholar

[65] View Article

[66] PubMed/NCBI

[67] Google Scholar

[ref19] 19. Gorski PA, Ceholski DK, Hajjar RJ. Altered myocardial calcium cycling and energetics in heart failure-a rational approach for disease treatment. Cell metabolism. 2015;21(2):183–194. pmid:25651173
View Article
PubMed/NCBI
Google Scholar

[69] View Article

[70] PubMed/NCBI

[71] Google Scholar

[ref20] 20. Samuel T, Rosenberry R, Lee S, Pan Z. Correcting Calcium Dysregulation in Chronic Heart Failure Using SERCA2a Gene Therapy. International journal of molecular sciences. 2018;19(4):1086.
View Article
Google Scholar

[73] View Article

[74] Google Scholar

[ref21] 21. Jaski BE, Jessup ML, Mancini DM, Cappola TP, Pauly DF, Greenberg B, et al. Calcium upregulation by percutaneous administration of gene therapy in cardiac disease (CUPID Trial), a first-in-human phase 1/2 clinical trial. Journal of cardiac failure. 2009;15(3):171–181. pmid:19327618
View Article
PubMed/NCBI
Google Scholar

[76] View Article

[77] PubMed/NCBI

[78] Google Scholar

[ref22] 22. Jessup M, Greenberg B, Mancini D, Cappola T, Pauly DF, Jaski B, et al. Calcium upregulation by percutaneous administration of gene therapy in cardiac disease (CUPID) a phase 2 trial of intracoronary gene therapy of sarcoplasmic reticulum Ca2+-ATPase in patients with advanced heart failure. Circulation. 2011;124(3):304–313. pmid:21709064
View Article
PubMed/NCBI
Google Scholar

[80] View Article

[81] PubMed/NCBI

[82] Google Scholar

[ref23] 23. Hulot JS, Ishikawa K, Hajjar RJ. Gene therapy for the treatment of heart failure: promise postponed. European heart journal. 2016;37(21):1651–1658. pmid:26922809
View Article
PubMed/NCBI
Google Scholar

[84] View Article

[85] PubMed/NCBI

[86] Google Scholar

[ref24] 24. Hamilton S, Terentyev D. Altered Intracellular Calcium Homeostasis and Arrhythmogenesis in the Aged Heart. International journal of molecular sciences. 2019;20(10):2386.
View Article
Google Scholar

[88] View Article

[89] Google Scholar

[ref25] 25. Eisner D. A., Caldwell J. L., Trafford A.W., and Hutchings D.C. The Control of Diastolic Calcium in the Heart: Basic Mechanisms and Functional Implications. Circulation Research. 2020; 126(3):395–412. pmid:31999537
View Article
PubMed/NCBI
Google Scholar

[91] View Article

[92] PubMed/NCBI

[93] Google Scholar

[ref26] 26. Trafford A, Diaz M, Eisner D. Coordinated control of cell Ca2+ loading and triggered release from the sarcoplasmic reticulum underlies the rapid inotropic response to increased L-type Ca2+ current. Circulation research. 2001;88(2):195–201. pmid:11157672
View Article
PubMed/NCBI
Google Scholar

[95] View Article

[96] PubMed/NCBI

[97] Google Scholar

[ref27] 27. Mahajan A, Shiferaw Y, Sato D, Baher A, Olcese R, Xie LH, et al. A rabbit ventricular action potential model replicating cardiac dynamics at rapid heart rates. Biophysical journal. 2008;94(2):392–410. pmid:18160660
View Article
PubMed/NCBI
Google Scholar

[99] View Article

[100] PubMed/NCBI

[101] Google Scholar

[ref28] 28. Stern M. D., Song L. S., Cheng H., Sham J. S., Yang H. T., Boheler K. R., and Ríos E. Local control models of cardiac excitation–contraction coupling: a possible role for allosteric interactions between ryanodine receptors. The Journal of general physiology. 1999;113(3):469–489. pmid:10051521
View Article
PubMed/NCBI
Google Scholar

[103] View Article

[104] PubMed/NCBI

[105] Google Scholar

[ref29] 29. Cantalapiedra I. R., Alvarez-Lacalle E., Peñaranda A., Echebarria B. Minimal model for calcium alternans due to SR release refractoriness. Chaos: An Interdisciplinary Journal of Nonlinear Science. 2017; 27(9): 093928.
View Article
Google Scholar

[107] View Article

[108] Google Scholar

[ref30] 30. Alvarez-Lacalle E, Echebarria B, Spalding J, Shiferaw Y. Calcium alternans is due to an order-disorder phase transition in cardiac cells. Physical review letters. 2015;114(10):108101. pmid:25815968
View Article
PubMed/NCBI
Google Scholar

[110] View Article

[111] PubMed/NCBI

[112] Google Scholar

[ref31] 31. Galice S, Xie Y, Yang Y, Sato D, Bers DM. Size matters: Ryanodine receptor cluster size affects arrhythmogenic sarcoplasmic reticulum calcium release. Journal of the American Heart Association. 2018;7(13):e008724. pmid:29929992
View Article
PubMed/NCBI
Google Scholar

[114] View Article

[115] PubMed/NCBI

[116] Google Scholar

[ref32] 32. Rog-Zielinska EA, Johnston CM, O’Toole ET, Morphew M, Hoenger A, Kohl P. Electron tomography of rabbit cardiomyocyte three-dimensional ultrastructure. Progress in biophysics and molecular biology. 2016;121(2):77–84. pmid:27210305
View Article
PubMed/NCBI
Google Scholar

[118] View Article

[119] PubMed/NCBI

[120] Google Scholar

[ref33] 33. Macquaide N, Tuan HTM, Hotta Ji, Sempels W, Lenaerts I, Holemans P, et al. Ryanodine receptor cluster fragmentation and redistribution in persistent atrial fibrillation enhance calcium release. Cardiovascular research. 2015;108(3):387–398. pmid:26490742
View Article
PubMed/NCBI
Google Scholar

[122] View Article

[123] PubMed/NCBI

[124] Google Scholar

[ref34] 34. Soeller C. Ryanodine receptor cluster size sets the tone in cerebral smooth muscle. Proceedings of the National Academy of Sciences. 2018;115(41):10195–10197.
View Article
Google Scholar

[126] View Article

[127] Google Scholar

[ref35] 35. Wagner J., and Keizer J. Effects of rapid buffers on Ca2+ diffusion and Ca2+ oscillations. Biophysical Journal. 1994;67(1):447–456. pmid:7919018
View Article
PubMed/NCBI
Google Scholar

[129] View Article

[130] PubMed/NCBI

[131] Google Scholar

[ref36] 36. Smith G. D., Wagner J., and Keizer J. Validity of the rapid buffering approximation near a point source of calcium ions. Biophysical Journal. 1996;70(6):2527–2539. pmid:8744292
View Article
PubMed/NCBI
Google Scholar

[133] View Article

[134] PubMed/NCBI

[135] Google Scholar

[ref37] 37. Shiferaw Y., Watanabe M. A., Garfinkel A., Weiss J. N., and Karma A. Model of intracellular calcium cycling in ventricular myocytes. Biophysical journal. 2003;85(6):3666–3686. pmid:14645059
View Article
PubMed/NCBI
Google Scholar

[137] View Article

[138] PubMed/NCBI

[139] Google Scholar

[ref38] 38. Shannon TR, Wang F, Puglisi J, Weber C, Bers DM. A mathematical treatment of integrated Ca dynamics within the ventricular myocyte. Biophysical journal. 2004;87(5):3351–3371. pmid:15347581
View Article
PubMed/NCBI
Google Scholar

[141] View Article

[142] PubMed/NCBI

[143] Google Scholar

[ref39] 39. Haken H. Synergetics: introduction and advanced topics. Springer Science & Business Media; 2013.

[ref40] 40. Qu Z., Liu M. B., and Nivala M. A unified theory of calcium alternans in ventricular myocytes. Scientific reports. 2016;6:35625. pmid:27762397
View Article
PubMed/NCBI
Google Scholar

[146] View Article

[147] PubMed/NCBI

[148] Google Scholar

[ref41] 41. Huertas MA, Smith GD and Györke S. Ca2+ alternans in a cardiac myocyte model that uses moment equations to represent heterogeneous junctional SR Ca2+. Biophysical journal. 2010:99(2):377–387. pmid:20643055
View Article
PubMed/NCBI
Google Scholar

[150] View Article

[151] PubMed/NCBI

[152] Google Scholar

[ref42] 42. Shannon T. R., Ginsburg K. S., and Bers DM. Reverse mode of the sarcoplasmic reticulum calcium pump and load-dependent cytosolic calcium decline in voltage-clamped cardiac ventricular myocytes. Biophysical Journal, 2000, 78.1: 322–333. pmid:10620296
View Article
PubMed/NCBI
Google Scholar

[154] View Article

[155] PubMed/NCBI

[156] Google Scholar

[ref43] 43. Keynes JM. The general theory of employment, interest, and money. Springer; 2018.

[ref44] 44. Sato D, Shiferaw Y, Garfinkel A, Weiss JN, Qu Z, Karma A. Spatially discordant alternans in cardiac tissue: role of calcium cycling. Circulation research. 2006;99(5):520–527. pmid:16902177
View Article
PubMed/NCBI
Google Scholar

[159] View Article

[160] PubMed/NCBI

[161] Google Scholar

[ref45] 45. Shiferaw Y, Aistrup GL, Wasserstrom JA. Intracellular Ca2+ waves, afterdepolarizations, and triggered arrhythmias. Oxford University Press; 2012.

[ref46] 46. Izu LT, Xie Y, Sato D, Bányász T, Chen-Izu Y. Ca2+ waves in the heart. Journal of molecular and cellular cardiology. 2013;58:118–124. pmid:23220129
View Article
PubMed/NCBI
Google Scholar

[164] View Article

[165] PubMed/NCBI

[166] Google Scholar

[ref47] 47. Li M. and Otani N. F. Ion channel basis for alternans and memory in cardiac myocytes. Annals of biomedical engineering. 2003;31(10):1213–1230. pmid:14649495
View Article
PubMed/NCBI
Google Scholar

[168] View Article

[169] PubMed/NCBI

[170] Google Scholar

[ref48] 48. Veasy J., Lai Y. M., Coombes S., and Thul R. Complex patterns of subcellular cardiac alternans. Journal of Theoretical Biology. 2019;478:102–114. pmid:31220466
View Article
PubMed/NCBI
Google Scholar

[172] View Article

[173] PubMed/NCBI

[174] Google Scholar

Figures

Abstract

Author summary

Introduction

Materials and methods

Rabbit computational model structure

General structure.

Rate equations in the CaRU.

Pacing by external clamped potential: Calcium transient and recovery.

Currents of rabbit model.

Calcium buffering implementation.

Parameters of the rabbit model and its modified version.

General equilibrium method to study calcium homeostasis

Results

Validation

General shocks

Consequences of shocks

Change in SERCA function is a double shock.

Physiological mechanism for SERCA gene therapy failure.

Discussion and conclusion

Supporting information

S1 File. Supporting information file with appendix.

Acknowledgments

References