Instabilities of the resting state in a mathematical model of calcium handling in cardiac myocytes

Mathematical Biosciences 236 (2012) 97–107

Contents lists available at SciVerse ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Instabilities of the resting state in a mathematical model of calcium handling in cardiac myocytes Aslak Tveito a,b, Glenn Terje Lines a,c,⇑, Johan Hake a,b, Andrew G. Edwards b a

Center for Biomedical Computing, Simula Research Laboratory, and Center for Cardiological Innovation, Oslo University Hospital, Norway Department of Bioengineering, University of California, San Diego, CA, USA c Department of Informatics, University of Oslo, Norway b

a r t i c l e

i n f o

Article history: Received 16 August 2011 Received in revised form 9 February 2012 Accepted 15 February 2012 Available online 3 March 2012 Keywords: Calcium waves Cardiomyocyte Stability Eigenvalue analysis

a b s t r a c t We analyze a recently published model of calcium handling in cardiac myocytes in order to find conditions for the presence of instabilities in the resting state of the model. Such instabilities can create calcium waves which in turn may be able to initiate cardiac arrhythmias. The model was developed by Swietach, Spitzer and Vaughan-Jones [1] in order to study the effect, on calcium waves, of varying ryanodine receptor (RyR)-permeability, sarco/endoplasmic reticulum calcium ATPase (SERCA) and calcium diffusion. We study the model using the extracellular calcium concentration ce and the maximal velocity of the SERCA-pump v SERCA as control parameters. In the ðce ; v SERCA Þ-domain we derive an explicit function v  ¼ v  ðce Þ, and we claim that any resting state based on parameters that lie above the curve (i.e. any pair ðce ; v SERCA Þ such that with v SERCA > v  ðce Þ) is unstable in the sense that small perturbations will grow and can eventually turn into a calcium wave. And conversely; any pair ðce ; v SERCA Þ below the curve is stable in the sense that small perturbations to the resting state will decay to rest. This claim is supported by analyzing the stability of the system in terms of computing the eigenmodes of the linearized model. Furthermore, the claim is supported by direct simulations based on the non-linear model. Since the curve separating stable from unstable states is given as an explicit function, we can show how stability depends on other parameters of the model. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction The contraction of cardiac myocytes is controlled by a process termed excitation–contraction coupling [2]. An action potential causes influx of calcium, through L-type calcium channels (LCCs), and this ‘‘trigger’’ calcium elicits further calcium release from the sarcoplasmic reticulum (SR) via the ryanodine receptors (RyR). Calcium then binds to troponin C, which initiates contraction. Finally, calcium is removed from cytosol to facilitate relaxation. Most of the calcium is reuptaken into the SR by the sarco-endoplasmic reticulum calcium ATPase (SERCA) pump. However a fraction of all calcium is taken out of the cell through the sodium-calcium exchanger (NCX). The NCX is a protein that exchanges 1 calcium with 3 sodium ions, creating a net flow of positive charge into the cell. Under experimental conditions where the extracellular calcium level is increased, pathological calcium waves can be induced [3,4]. The elevated extracellular calcium is transported across the cell membrane and eventually pumped into SR, via the SERCA pump. Inside the SR, calcium sensitizes the RyRs, which, together with ⇑ Corresponding author at: Center for Biomedical Computing, Simula Research Laboratory, and Center for Cardiological Innovation, Oslo University Hospital, Norway. E-mail address: [email protected] (G.T. Lines). 0025-5564/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2012.02.005

the elevated cytosolic calcium, will eventually trigger calcium release from SR [5,6]. Calcium overload can also be induced by stimulation of the b-adrenergic pathway [7], which would mimic more natural conditions for the myocyte. Whatever the cause of calcium overload, if the size of this spontaneous release is large enough, a propagating wave will be generated [8]. An SR calcium threshold for waves has been identified and named store-overload-induced calcium release [9]. The released calcium is extracted through the NCX which generates a net flow of positive charge into the cell. If the magnitude of the release is sufficiently large, it can trigger an action potential; a so called ectopic beat. Such beats are believed to play a decisive role in the initiation of cardiac arrhythmias [10,11]. It has been shown that a critical balance between sensitized RyRs, and altered SERCA activity is required to trigger calcium waves [12,7]. Computational models have been used to shed light on the exact balance between these two parameters for the generation of calcium waves [13–15]. A model of calcium wave generation and propagation has recently been suggested by Swietach, Spitzer and Vaughan-Jones [1] and used by Stokke et al. [16] to probe the generation of calcium waves during pathological conditions. The model is a deterministic system of partial differential equations which is simple enough to be analyzed but also complex enough to represent the processes underlying wave generation. Importantly, features necessary for a complete view of myocyte

98

A. Tveito et al. / Mathematical Biosciences 236 (2012) 97–107

electrophysiology, such as the generation and shape of an action potential, are not included. However, we are studying generation of calcium waves, which are important for delayed after depolarisations, and such waves are triggered during the resting phase of the action potential. Our aim is to understand the model from a mathematical point of view. More specifically, we want to assess the stability of the resting state as a function of model parameters that are known to be important for the generation of calcium waves. If the resting state is unstable, a propagating calcium wave can be generated. Our prime control variables are the extracellular calcium concentration ce and the strength of the SERCA pump v SERCA . For standard parameters we theoretically derive a function v  in the ðce ; v SERCA Þ space dividing stable and unstable parameters. This curve explains how the stability of the cell depends on the main control parameters. Furthermore, we demonstrate how this curve is altered as a function of the RyR sensitivity parameter K a . The model of Swietach, Spitzer and Vaughan-Jones [1] is a system of reaction–diffusion equations. As shown in e.g. [17–19] the resting state of such systems can be analyzed by eigenmode analysis. We compute a resting state of the system and derive a linearized model around that state. The linearized model is discretized using a finite difference method leading to a linear system of ordinary differential equations. It is well known, see e.g. Chicone [20], that the stability of resting states for a system of ordinary differential equations can be analyzed by invoking the real parts of the eigenvalues of the system matrix. By using this method for the present model, we show that the curve v  ¼ v  ðce Þ separates stable from unstable parameters. Since this method is based on a linearization, we present computations based on the full non-linear model showing that the curve we have derived separates stable from unstable states also in the non-linear case. We emphasize, however, that our findings are entirely based on the mathematical model and we do not know if similar properties can be seen in experiments. The present paper is organized as follows: The mathematical model of [1,16] is presented in the next section where we also derive an analytical expression for the curve v  ¼ v  ðce Þ separating stable and unstable parameters. In Section 3 we present numerical computations based on linear analysis of the system of reaction–diffusion equations, and simulations based on stochastic initial conditions. In Section 4, we present a generalized model for which a similar analysis can be performed, and a conclusion is given in Section 5. 2. The mathematical model The mathematical model under consideration was developed by Swietach, Spitzer and Vaughan-Jones [1], and applied in a slightly modified version by Stokke et al. [16]. Here we use the form of the model presented in the online supplement of [16]. The model consists of a one dimensional system of partial differential equations of the following form:

@ccyt @ 2 ccyt 1  J cytbuf þ ðJ þ J RyR  J SERCA Þ; ¼ Dcyt V cyt sl @t @x2 2 @cSR @ cSR 1  J SRbuf þ ðJ RyR þ J SERCA Þ; ¼ DSR V SR @t @x2 2 @bcyt @ bcyt þ J cytbuf ; ¼ Dbuf @t @x2 @bSR ¼ JSRbuf ; @t C SRbuf  bSR @p ¼ kon2 ccyt ð1  pÞ  koff 2 p: @t C SRbuf  bSR0

J sl ¼ ksl ðce  ccyt Þ  v NCX

ccyt ; ccyt þ K NCX

ð6Þ

and the RyR release/leak is given by

J RyR ¼

  kRyR max 0;

  ccyt  p þ kleak ðcSR  ccyt Þ: ccyt þ K a

ð7Þ

In this system, Eqs. (1)–(4) are on standard form for models of calcium dynamics, see e.g. [21,22]. The dynamics of the RyR current are modeled in terms of a threshold; whenever the value of the fast c gate, given by ccytcyt is larger than a threshold, given by the slow þK a inactivation gate, p, there is a current. The difference between the fast and the slow gate will then give the size of the current. The fast gate is not included as a state variable in the system of equations, as it is assumed to be in quasi steady state with cytosolic calcium. The slow gate acts as an inactivation gate. From (5) we see that this gate is activated by increased cytosolic calcium ccyt together with the amount of unbound calcium buffer in SR, ðC SRbuf  bSR Þ. The former models inactivation of the RyR by cytosolic calcium and the latter models RyR activation of SR calcium and explains why a calcium overloaded cell is more prone to generate calcium waves [5,6,23]. The SERCA-pump brings calcium from cytosol into the SR and is modeled by

J SERCA ¼ v SERCA

c2cyt  ðcSR =dÞ2 c2cyt þ ðcSR =dÞ2 þ K 2SERCA

:

ð8Þ

Here K SERCA is the dissociation constant for cytosolic calcium, d is the ratio between the dissociation constants of calcium in SR and cytosol, and v SERCA is the maximum velocity of the SERCA-pump. By increasing the activity of the SERCA pump, we increase the SR calcium load, as the increased SERCA activity has to be balanced during rest by an increase in SR calcium leak. The latter can only be achieved by an increase in SR calcium load. Finally the transition rates from unbound to bound cytosolic and SR buffer are given by

J cytbuf ¼ kon ðccyt C cytbuf  ðccyt þ K cytbuf Þbcyt Þ; J SRbuf ¼ kon ðcSR C SRbuf  ðcSR þ K SRbuf ÞbSR Þ:

ð9Þ ð10Þ

Here C cytbuf and C SRbuf refer to the total amount of calcium buffer in cytosol and SR respectively. The parameters of the model are taken from [16] and are listed in the Table 1. In our analysis the extracellular calcium concentration ce , and the strength of the SERCA-pump v SERCA will be used as control parameters; their default values taken from [16] are ce =1500 lM and v SERCA ¼ 650 lM s1 . We start our analysis by considering the system given by Eqs. (1)–(5) with the fluxes and constants described above. In Section 4, we will present a similar analysis for a larger class of models that can be written in the form Eqs. (1)–(5). 2.1. Constant solutions at rest

ð1Þ ð2Þ

A constant solution in space and time of the system above satisfies the following algebraic system of equations

ð3Þ

0 ¼ J cytbuf þ

1 ðJ þ J RyR  J SERCA Þ; V cyt sl

ð11Þ

0 ¼ J SRbuf þ

1 ðJ RyR þ J SERCA Þ; V SR

ð12Þ

ð4Þ ð5Þ

Here ccyt and cSR denote the calcium concentration in cytosol and SR; similarly bcyt and bSR denote the buffered calcium concentrations, and p denotes the occupancy of the inactivating gate of the RyR. The net sarcolemmal calcium influx is given by

0 ¼ J cytbuf ;

ð13Þ

0 ¼ J SRbuf ;

ð14Þ

0 ¼ kon2

C SRbuf  bSR ccyt ð1  pÞ  koff 2 p: C SRbuf  bSR0

ð15Þ

99

A. Tveito et al. / Mathematical Biosciences 236 (2012) 97–107

where /  104 . We get the equation

Table 1 The table gives the constants used in the mathematical model. The constants are taken from Stokke et al. [16]. 1

1

kon

200 lM

Dbuf Dcyt DSR

70 lm2 s1 500 lm2 s1 20 lm2 s1 36 lM

kNCX

v NCX

s

2400 lM s1 851.1586 lM

bSR0

kon2 V cyt C cytbuf C SRbuf K SERCA Ka d

1

10 lM s 0.65 175 lM 4400 lM 0.25 lM 0.35 lM 2000

1

koff2

ksl /2 c2e þ ðv NCX þ ksl K NCX  ce ksl Þ/ce  ce ksl K NCX ¼ 0;

1

V SR K cytbuf K SRbuf

3s 0.035 2 lM 630 lM

ksl kRyR kleak

0:0066 s1 60 s1 0:8 s1

and thus, by ignoring higher order terms, we get

/

ksl K NCX

v NCX þ ksl K NCX  ce ksl ksl K NCX



ksl K NCX

v NCX þ ksl K NCX

þ 2:722  1010 ce

þ 7:487  1016 c2e

It follows from this that

ð16Þ

/¼

J sl ¼ 0:

ð22Þ

and thus we use the approximation

and

ð17Þ

ksl K NCX

ð23Þ

v NCX þ ksl K NCX

so

From (17) and (6) we get

ksl c2cyt

ð21Þ

:

From the parameters given in Table 1, we have

v NCX þ ksl K NCX  ce ksl J SERCA ¼ J RyR

ð20Þ

þ ðv NCX þ ksl K NCX  ce ksl Þccyt  ce ksl K NCX ¼ 0:

ð18Þ

ccyt ðce Þ ¼ /ce ¼

By using the parameters values given in Table 1 for ksl ; K NCX and v NCX , and then graph ccyt as a function of ce , we note that ccyt  104 ce (see Fig. 1). Based on this observation it is reasonable, in general, to seek an approximate relation of the form

ccyt ¼ /ce

ksl K NCX

v NCX þ ksl K NCX

ce :

ð24Þ

In Fig. 1 we have graphed ccyt as a function of ce given as the solution of (18) and compared it to the linear approximation given by (24) for four different sets of the parameters ðksl ; K NCX ; v NCX Þ. We observe that the linear approximation given by (24) is accurate for all four parameter sets and it will therefore be used below.

ð19Þ

ksl

1.5 1

c

cyt

(μM)

2

0.5 0 0

1000

2000

3000

4000

5000

6000

7000

8000

9000

6000

7000

8000

9000

6000

7000

8000

9000

VNCX

ccyt (μM)

2 1.5 1 0.5 0 0

1000

2000

3000

4000

5000

KNCX

ccyt (μM)

2 1.5 1 0.5 0 0

1000

2000

3000

4000

5000

c (μM) e

Fig. 1. The accuracy of the linear approximation (24) is illustrated using three sets of parameters. Upper panel: The green graphs show the exact solution (solid line) of (18) and the linear approximation (dotted line) given by (24) for standard parameters. Similar results are given by the red graphs where the value of ksl is multiplied by two, and in the blue graphs where the value of ksl is divided by two. Center panel: Same as the upper panel, but the value of v NCX is varied; standard (green), two times standard (red), and half of standard (blue). Lower panel: Similar to upper and center but the values of K NCX is varied. We conclude that for all parameter sets, the approximation given by (24) is accurate. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

100

A. Tveito et al. / Mathematical Biosciences 236 (2012) 97–107

In the system Eqs. (11)–(15) we seek a solution that is constant in space and time, and we seek a solution at a point in the statespace where the RyR-release is about to activate. From (7) it follows that the transition from inactive to active RyR-release is characterized by the equation

ccyt  p ¼ 0; ccyt þ K a

ð25Þ

/ce /ce þ K a

ð26Þ

where / is given by (23). From (13) and (9) we have

ccyt C cytbuf  ðccyt þ K cytbuf Þbcyt ¼ 0

ð27Þ

so by applying the approximation (24) we get

bcyt ðce Þ ¼

C cytbuf /ce /ce þ K cytbuf

ð28Þ

where again / is given by (23). Furthermore, from (15), we get

kon2 C SRbuf  bSR p ccyt ¼ 1p koff 2 C SRbuf  bSR0

ð29Þ

c

and since p ¼ ccytcyt (see (25)) we have þK a

kon2 C SRbuf  bSR 1 ¼ koff 2 C SRbuf  bSR0 K a

ð30Þ

We define

bSR1 ¼ C SRbuf 

koff 2 ðC SRbuf  bSR0 Þ: kon2 K a

ð31Þ

and conclude that

ð32Þ

which is a constant independent of ce . From (14) and (10) we have

cSR C SRbuf  ðcSR þ K SRbuf ÞbSR ¼ 0

ð33Þ

and then

K SRbuf bSR C SRbuf  bSR

ð34Þ

Since bSR ¼ bSR1 , we have

cSR ¼ cSR1

Above, we computed a curve v  in the ðce ; v SERCA Þ-domain and we claim that this curve divides parameters leading to stable and unstable resting states. The main purpose of this section is to provide computational evidence that this claim holds. We do that by assessing the stability of the resting state by invoking the real parts of the eigenvalues of an associated linear problem. This technique has recently been applied to assess the stability of resting state in cardiac tissue, see [17–19]. We start this section by repeating the basic steps of linear analysis of a system of reaction diffusion equations. In the rest of the section we present numerical experiments discussing the validity of our claim regarding the stability of the problem. In the numerical computations, we start by using linear theory to investigate the stability as a function of extracellular calcium concentration and the strength of the SERCA-pump. Next we consider the influence of molecular diffusion, and the magnitude of K a representing the dissociation constant of the fast activating RyR-site. We study how the magnitude of the principal eigenvalue depends on the RyR channel flux constant kRyR , and how a change in sign of an eigenvalue alters the behavior of the solution. Finally, we discuss the accuracy of linear analysis. 3.1. Linear analysis revisited

zt ¼ dzxx þ f ðzÞ

K SRbuf bSR1 C SRbuf  bSR1

ð36Þ

v SERCA ¼

kleak ðcSR 

ccyt Þðc2cyt

2

þ ðcSR =dÞ þ

c2cyt  ðcSR =dÞ2

K 2SERCA Þ

:

ð37Þ

We have found the critical strength of the SERCA-pump as a function of the extracellular calcium concentration ce at the point in state-space where the RyR-release activates. From the approximation (24 ) and the fact that cSR is a constant given by (36) we get

v  ðce Þ ¼

kleak ðcSR1  /ce Þðð/ce Þ2 þ ðcSR1 =dÞ2 þ K 2SERCA Þ ð/ce Þ2  ðcSR1 =dÞ2

ð38Þ

ð40Þ

where z is a vector containing the variables ðccyt ; cSR ; bcyt ; bSR ; pÞ, the function f carries the fluxes and functions defined in Eqs. (6)–(10) and d denotes a diagonal matrix of the form d ¼ diagðDcyt ; DSR ; Dbuf ; 0; 0Þ. We discretise the system by introducing the nodes xj ¼ jDx, where Dx ¼ L=n. Here L denotes the length of the cell (set to 100 l m in our computations) and n denotes the number of computational nodes (set to 50 in our computations). Let zj ðtÞ denote the approximate solution at x ¼ xj defined as the solution of the following system of ordinary differential equations,

z0j ðtÞ ¼ d

which is also independent of ce . It follows from (16) that J RyR ¼ J SERCA and because of (25) it follows from (7) and (8) that

where

In the ðce ; v SERCA Þ domain, we claim that any resting state based on parameters below this graph given by v SERCA ¼ v  ðce Þ is stable and any resting state based on parameters above the graph is unstable. The validity of this claim will be addressed in computations below.

ð35Þ

where

cSR1 ¼

ð39Þ

The system Eqs. (1)–(5) can be written in the form

bSR ¼ bSR1

cSR ¼

ksl K NCX

v NCX þ ksl K NCX

3. Numerical computations

and thus, from (24), we get

pðce Þ ¼

/¼

zj1  2zj þ zjþ1 ðDxÞ2

þ f ðzj Þ

ð41Þ

for j ¼ 0; . . . ; n. Suppose a stationary solution is given by z0 , i.e. f ðz0 Þ ¼ 0. Then we can introduce y ¼ z  z0 , and find that, up to linear terms, y is governed by a system of the form

y0 ¼ Ay:

ð42Þ

It is well known from the theory of ordinary differential equations (see e.g. Chicone [20]) that the stability of the resting state z0 is determined by the real parts of the eigenvalues of the system matrix A. More specifically, the system is termed unstable if the largest real part is greater than zero. Observe that if k is a real-valued eigenvalue of A and r is the associated eigenvector, then ekt r solves the system (42) when the initial condition is given by yð0Þ ¼ r. Thus an initial perturbation in terms of an eigenvector will grow or decay at exponential rate depending on the sign of the associated eigenvalue. Due to the inherit rhythm of cardiomyocytes, the magnitude of the eigenvalues are of great importance; all the cells in the tissue are reset after a few hundred milliseconds and thus a perturbation that grows slowly, will not be able to create a calcium wave. In the

101

A. Tveito et al. / Mathematical Biosciences 236 (2012) 97–107

1400

V

SERCA

(μM/s)

1200

1000

800

600

400

200 2000

3000

4000

5000

6000

7000

8000

ce (μM) Fig. 2. The solid line represents the estimated transition from stable to unstable equilibrium point given by (38). The dark grey area represents a region where the real part of all eigenvalues of the linearized system is negative, and in the light grey area, at least one eigenvalue has a positive real part. Furthermore, the transition from negative to positive values are found by using a root-finding function in Matlab and are marked by . We observe that the curve computed by (38) coincides completely with the curve separating negative and positive principal eigenvalues.

computations presented below we will see that in the present model the eigenvalues of the unstable cases are rather big and can therefore be able to initiate waves. 3.2. Stability as a function of extracellular calcium and the strength of the SERCA-pump To assess the susceptibility of the model to generating calcium waves we examine how the level of extracellular calcium ce and the strength of the SERCA-pump v SERCA impact the stability of the 900

associated resting state. For the parameters given in Table 1, we show the stability as a function of ce and v SERCA in Fig. 2. Given a point ðce ; v SERCA Þ, we compute the resting state of the system Eqs. (1)–(5) and the largest real part of the eigenvalues of the associated linearized system (42). In the figure, the light grey area (upper part) denotes an area where the principal eigenvalue has a positive real part. Thus for sufficiently large values of ce and v SERCA the resting state is unstable. Correspondingly, the dark grey area represents a stable area in the ðce ; v SERCA Þ space. These areas are found by using the linearization described above. Furthermore, the transition from negative to positive values are found by using a root-finding function in Matlab and are marked by . The solid line is given by the function (38) derived above and we observe that the solid line coincides with the transition from the stable to the unstable region, and thus provides support for the claim that the curve defined by (38) separates stability from instability. The curve in Fig. 2 illustrates that increasing either ce or v SERCA (or both) will make the model less stable. In the unstable region, a small perturbation will grow and may turn into a full calcium wave. By increasing ce we increase ccyt bringing the RyR current closer to the release threshold, see (7). Although increasing v SERCA does not alter ccyt , it lowers the threshold variable, p, by increasing steady state SR calcium load. The latter is indirectly accomplished by lowering the amount of unbound calcium buffer in SR, which regulates p. 3.3. Stability as a function of diffusion In Fig. 3 we study how the stability of a resting state depends on the diffusion coefficients of the system Eqs. (1)–(5) given by the three parameters ðDcyt ; DSR ; Dbuf Þ. More specifically, we have plotted the third largest real part of the eigenvalues (the two largest seem to be independent of diffusion) for three values of molecular diffusion; one standard given in Table 1 and denoted D in the figure, and two cases where we either multiply or divide D ¼ ðDcyt ; DSR ; Dbuf Þ by 20

λ1 λ ,D/10

800

3

λ3

15

λ ,D*10 3

600

largest real part of the eigenvalues

largest real part of the eigenvalues

700

500 400 300 200 100 00

−100 2000

10

5

0

−5

4000

ce (μM)

6000

8000

2722.55

2722.6

2722.65

ce (μM)

Fig. 3. Left panel: The figure shows the magnitude of the largest and third largest real part of the eigenvalues, and how they depend on the diffusion coefficients. The largest and second largest real part of eigenvalues are independent of diffusion, but the third largest depends on diffusion; the magnitude is reduced as diffusion is increased. Right panel: Same data as on the left, zoomed in around the critical interval.

102

A. Tveito et al. / Mathematical Biosciences 236 (2012) 97–107

10.0. We observe that the stable region seems to be rather insensitive to the magnitude of molecular diffusion. We graph the real part of the principal eigenvalue as a function of ce and we observe that as ce passes the critical point from the stable region to the unstable region, the eigenvalue changes sign and the magnitude increases dramatically. This indicates that for unstable parameters, even a small perturbation will diverge rapidly from the resting state. The sudden change in largest eigenvalue is a result of the threshold for the RyR current being passed. Although the sign of the real part of the eigenvalue seems to be almost independent of the size of the molecular diffusion, the magnitude of the real parts of the eigenvalues may still depend on it. This is also illustrated in Fig. 3 were we observe that the magnitude of the third largest eigenvalue in the unstable region decreases as the magnitude of molecular diffusion is increased. This is because a larger diffusion will lower the regenerative effect of the RyR current as calcium is faster transported away from the local release site. So more diffusion stabilizes the stationary solutions of Eqs. (1)–(5), which is generally expected to be the case, see e.g. [17]. This effect is highlighted in the right panel of Fig. 3. In summary we have observed that the magnitude of molecular diffusion can affect the magnitude of the real part of the eigenvalue, but it does not seem to affect the sign of the real part of the principal eigenvalue.

region is sensitive to changes in K a ; the size of the stable region is significantly increased by increasing the value of K a . This is because the activating gate becomes less sensitive to calcium and a larger ccyt and hence larger ce is needed to trigger the RyR current. In the right panel of Fig. 4 we observe, again, a significant jump in the magnitude of the principal eigenvalue as the extracellular concentration is changed from the stable to the unstable region. Interestingly, the largest eigenvalue gets larger as K a increases. This is because a higher K a will cause an elevated SR calcium content, which results in a larger gradient for the RyR current once the release starts. 3.5. The magnitude of the principal eigenvalue as a function of the RyR channel flux constant kRyR The RyR current flux constant kRyR determines the size of the release current once it is triggered, see (7 ). A larger release will, through positive feedback via an elevated ccyt , act to re-enforce itself. Hence we would expect the principal eigenvalue of an unstable system to increase with kRyR . In Fig. 5 we illustrate that this turns out to be the case. 3.6. Stability as a function of the strength of the Na/Ca-exchanger In the model above, the Na/Ca-exchanger (NCX) flux is given by

3.4. Stability as a function of the dissociation constant of the fast, activating RyR-site; K a

J NCX ¼ v NCX

In Fig. 4 we study how the regions of stability/instability depend on the dissociation constant of the fast, activating RyR-site called K a ; see (7). This constant represents the calcium affinity of the fast activating gate. The default value is given by K a ¼ 0:35 (green graph), and in the figure we also consider the case of K a ¼ 0:9  0:35 (blue graph) and K a ¼ 1:1  0:35 (red graph). We first note that the results obtained by the theoretical estimate (38) again coincide with the results obtained through the eigenvalue analysis described above. Second, we observe that the stable

ccyt ; ccyt þ K NCX

ð43Þ

see latter term of Jsl in Equation (6). The strength of the flux is governed by the two parameters v NCX and K NCX , and it is of interest to understand how the stability of the problem depends on these parameters. The stability curve (38) is given by

v  ðce ; K NCX ; v NCX Þ ¼

kleak ðcSR1  /ce Þðð/ce Þ2 þ ðcSR1 =dÞ2 þ K 2SERCA Þ ð/ce Þ2  ðcSR1 =dÞ2 ð44Þ

900 Kα=0.9*0.35

1200

theoretical Kα=0.35 theoretical K =1.1*0.35 α

VSERCA (μM/s)

1000

theoretical

800

600

largest real part of the eigenvalues

800 700 600 500 400 300 200

400 100 0

200 2000

4000

6000

c (μM) e

8000

−100

2000

4000

6000

8000

c (μM) e

Fig. 4. Left panel: The figure shows that the stable region depends heavily on the dissociation constant of the fast, activating RyR-site K a ; the stable region is increased when the value of K a is increased. The figure also illustrates that the theoretical estimate (solid line) provides a good approximation to the transition from stable to unstable states;  illustrates points in the ðce ; v SERCA Þ-domain where the largest real part of the eigenvalues goes from negative to positive. Right panel: The figure shows the largest real part of the eigenvalues as a function of the extracellular Calcium concentration ce . Standard parameters are used (see Table 1), and v SERCA ¼ 650 l Ms1 . Observe that there is a large jump in the magnitude of the eigenvalue as it goes from being mildly negative to strongly positive.

103

A. Tveito et al. / Mathematical Biosciences 236 (2012) 97–107

1400 A: Default B: 1.5*K

NCX

5000

1200

4000

1000

3000

2000

1000

800

600

400

0 20

40

60

k

RyR

80

3000

4000

5000

6000

7000

8000

ce (μM)

(1/s)

Fig. 5. The figure illustrates how the largest real part of the eigenvalues depends on the RyR channel flux constant kRyR ; it is strongly increasing. The  marks the default value of kRyR . We have used the standard parameters given in Table 1, and we have put v SERCA ¼ 650 l Ms1 .

where we recall that

ksl K NCX =v NCX : 1 þ ksl K NCX =v NCX

200 2000

100

ð45Þ

From the form of /, we immediately see that the stability curve is a function of the ratio K NCX =v NCX . The changes in the stable region as a function of K NCX and v NCX are analyzed using numerical computations in Fig. 6. We consider four cases; A : ðK NCX ; v NCX Þ; B : ð1:5K NCX ; v NCX Þ; C : ðK NCX ; 1:5v NCX Þ; D : ð1:5K NCX ; 1:5v NCX Þ where the default case is given by ðK NCX ; v NCX Þ ¼ ð36; 2400Þ. The numerical computations show that the stability regions are identical for the two cases A and D. This is what we expected since the value of / is equal for these cases. Furthermore, v SERCA is smaller than default for case B and larger than default for case C, so the stable region is smaller than default in case B and larger than default in case C. Again, we note that the theoretical estimate of the boundary between stable and unstable regions coincides with the curve where the largest real part of the eigenvalues changes sign.

Fig. 6. The figure shows how changes in the parameters K NCX and v NCX (see Eq. (43)) governing the strength of the NCX-current affect the stable region. We observe that case A (default parameters) and D (both K NCX and v NCX multiplied by 1.5) are identical. If we increase K NCX by a factor of 1.5 and keep v NCX at the default value (case B), we observe that the stable region becomes smaller, and thus the stability decreases when K NCX is increased. On the other hand, if we increase v NCX by a factor of 1.5 and keep K NCX constant (case C), we observe that the stable region becomes larger, and hence the stability increases when v NCX is increased. Similarly to Fig. 2, the solid lines show the theoretical estimate, while transitions from negative to positive eigenvalues are marked with .

2 1.8

deviation from equilibrium

0

/¼

C: 1.5*VNCX D: 1.5*Default

VSERCA (μM/s)

largest real part of the eigenvalues

6000

λ100 =2.1

1.6

λ

101

=−0.59

1.4 1.2 1 0.8 0.6 0.4

3.7. An illustration of the effect of change of sign of an eigenvalue 0.2

In Fig. 7 we consider two possible perturbations of the resting state. In order to do the comparison, we introduce the L2  norm given by

kzk ¼

X Z i

L

1=2 ðzi ðxÞÞ2 dx

0

where i denotes the components of the vector z ¼ ðccyt ; cSR ; bcyt ; bSR ; pÞ. The blue graph illustrates the L2  norm of the deviation between the resting state and the solution of the system (1)–(5) when the initial condition is perturbed using the eigenvector associated with eigenvalue number 100 (arranged according to the size of real part). The eigenvalue is positive, and the solution deviates from equilibrium at exponential rate. This computation is based on the solution of the fully non-linear model (1)–(5) and illustrates that the linear analysis provides reasonable prediction about the solution of the problem. In the green line we illustrate the non-linear solution of a problem where the initial condition is defined by a perturbation of the resting state using the eigenvector associated with eigenvalue number 101. We observe that the eigenvalue is negative and a rather big perturbation decays in agreement with linear theory.

0

0

0.5

1

1.5

2

time (s) Fig. 7. The figure illustrates how a perturbation develops as a function of time measured in the L2  norm. In the computations we have used ðce ; v SERCA Þ ¼ ð3000; 650Þ. If we perturb the constant state using an eigenvector associated with an eigenvalue with a positive real part, the magnitude of the perturbation grows. And similarly, if we perturb the resting state using an eigenvector associated with an eigenvalue with a negative real part, the magnitude of the perturbation decays. This behavior is what we expect based on linear theory.

3.8. Accuracy of linear analysis In general terms, linear analysis is believed to provide accurate predictions for small perturbations and over small time-intervals. This was illustrated in the example mentioned above. However, in the present case it is extremely hard to provide accurate analytical assessments of how small the perturbation must be and for how small time-intervals the analysis holds. To shed some light on this we compare the development of the deviation from equilibrium, measured in the L2  norm, for the solution of the non-linear

104

A. Tveito et al. / Mathematical Biosciences 236 (2012) 97–107

K = 0.35

deviation from equilibrium

α

1 0.95 non−linear

0.9

linear

0.85 0.8 0.75 0.7

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

3

3.5

4

4.5

1

deviation from equilibrium

Kα = 0.35*0.9

10 8

non−linear linear

6 4 2 0

0

0.5

1

1.5

2

2.5

5 −3

time (s)

x 10

Fig. 8. The figure illustrates the accuracy of linear analysis. In the upper panel we compare the development of the L2  norm of the solution of the non-linear system Eqs. (1)– (5) and the L2  norm as predicted by linear theory. The initial condition is defined by the eigenvector, scaled to unit length, associated with the principal eigenvalue. In the upper case, we assess a stable problem and the two graphs are almost identical. In the unstable case below, we note that the estimate based on linear theory is good initially and then deviates more and more from the non-linear solution.

system (1)–(5) and the deviation as predicted by linear analysis. We use an eigenvector r as the initial perturbation and then linear theory predicts that the deviation will evolve like ekt r where k is the associated eigenvalue. In Fig. 8 we compare the deviation from

Default

V

=350

Serca

equilibrium measured in the L2  norm for the linear and non-linear cases. We observe that linear analysis is extremely accurate in the stable case (upper panel), and in the unstable case it gives accurate prediction for the first few milliseconds.

K

*0.9

alpha

Diffusion/10

Fig. 9. The simulations are based on ce ¼ 1; 2; 3; 4; 5 mM from top to bottom; time is horizontal and space is the vertical axis. Each simulation uses initial conditions equal to the equilibrium solution for the given parameter choice (indicated in the titles, relative default values). Plot shows deviation from the initial calcium distribution where blue is zero and red is 1 lM. Each 10 ms, ccyt is perturbed by adding a uniformly distributed variable in the interval plus/minus 0.05 lM. The simulations go on for 50 ms, the vertical axis is 100 lM, and Dx ¼ 2 lM. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

105

A. Tveito et al. / Mathematical Biosciences 236 (2012) 97–107

3.9. Assessment of stability in terms of non-linear analysis Finally, we use direct simulations based on the system Eqs. (1)–(5) to evaluate whether the stability predictions given above hold in the fully nonlinear case. In Fig. 9 we consider simulations for four different data-sets using five different values of the extracellular calcium concentration ce . The initial condition is provided by the associated equilibrium solutions. For each 10 ms the numerical approximation of ccyt is perturbed, at each computational node, by adding a uniformly distributed random variable in the interval plus/minus 0.05 lM. In the left column we have used default data (see Table 1) and we have used v SERCA ¼ 650 lM s1 . By solving the algebraic equation v  ðce Þ ¼ 650 we find that the problem should become unstable for ce > c  2743 lM, or  2724 lM if we use the fully-non-linear version of ccyt given as the solution (18). We observe from the figure that for small values of ce the resting state is stable in the sense that no wave is created from the perturbations. This means that each perturbation simply decays and the solution converges to equilibrium, which now defines a resting state. But when ce P 3 mM, the equilibrium state becomes unstable and we observe that perturbations lead to calcium waves that do not return to the equilibrium. In the second column we have reduced the value of v SERCA and the theory predicts that we should have unstable solutions for ce > c  4500 lM. From the figure, we observe that perturbations do not create waves for ce smaller than 4 mM, and that we do get calcium waves for ce ¼ 5 mM. In the two latter columns we change other parameters and we observe that these changes in the nonlinear case behave as predicted by linear theory. The computations presented in Fig. 9 support the claim that the curve given in (38) separates stable and unstable resting states of the system Eqs. (1)–(5).

where kRyR ¼ kRyR ðzÞ can be any strictly positive function of any of the variables; z denotes the vector ðccyt ; cSR ; bcyt ; bSR ; pÞ, and where the coefficient m can be any positive real number. Furthermore, the function kleak ¼ kleak ðzÞ can be any positive function of all the variables. 3. There is a function F 1 such that FðbÞ ¼ b0 implies b ¼ F 1 ðb0 Þ. 4. The SERCA pump is given by

J SERCA ¼ v SERCA Sðccyt ; cSR Þ;

ð52Þ

where S can be any strictly positive function. Note that in the original paper [1] it is pointed out that in some cases kleak is given by:

 c 1:7 SR kleak ðcSR Þ ¼ 0:8 : 500

ð53Þ

We therefore only use kleak ¼ kleak ðcSR Þ, but the more general case of kleak ¼ kleak ðzÞ can be handled in the same manner. Furthermore, we can also have S ¼ SðzÞ as long as S is strictly positive. 4.1. Constant solutions at rest Again we start by seeking a constant and stationary solution of the system

1 ðJ þ J RyR  J SERCA Þ; V cyt sl 1 þ ðJ RyR þ J SERCA Þ; V SR

0 ¼ J cytbuf þ

ð54Þ

0 ¼ J SRbuf

ð55Þ

0 ¼ J cytbuf ; 0 ¼ J SRbuf ; 0 ¼ kon2 FðbSR Þccyt ð1  pÞ  koff 2 p:

ð56Þ ð57Þ ð58Þ

As above, it follows from this that 4. A generalized mathematical model

J SERCA ¼ J RyR

We have derived a sharp estimate of a stable region for the model stated above. In this section, we will consider the problem of generalizing the result to other models that can be written on the same form. We start with the system of the form

and

2

@ccyt @ ccyt 1  J cytbuf þ ðJ þ J RyR  J SERCA Þ; ¼ Dcyt V cyt sl @t @x2

ð46Þ

@cSR @ 2 cSR 1  J SRbuf þ ðJ RyR þ J SERCA Þ; ¼ DSR V SR @t @x2 2 @bcyt @ bcyt þ J cytbuf ; ¼ Dbuf @t @x2 @bSR ¼ J SRbuf ; @t @p ¼ kon2 FðbSR Þccyt ð1  pÞ  koff 2 p: @t

ð47Þ

J sl ¼ 0:

ð60Þ

Therefore, by Assumption 1, we have

ccyt ðce Þ ¼ J 1 sl ðc e Þ:

ð49Þ ð50Þ

Compared to the system Eqs. (1)–(5), we allow more general functions and fluxes to be involved. We keep the form of J cytbuf and JSRbuf unchanged (see (9), (10)) and for the remaining fluxes and functions we make the following assumptions:

m cm cyt þ K a

 p ¼ 0;

ð62Þ

and thus

pðce Þ ¼ ¼

cm cyt ðc e Þ m ccyt ðce Þ þ K m a J m sl ðc e Þ m J m sl ðc e Þ þ K a

ð63Þ ð64Þ

1 where we have used the convention that J m ðcÞÞm . sl ðcÞ ¼ ðJ From (56) we have

J cytbuf ¼ 0

1. There is a function J 1 sl such that

J sl ðce ; ccyt Þ ¼ 0

ð65Þ

and then from (9) we get

implies that

ccyt C cytbuf  ðccyt þ K cytbuf Þbcyt ¼ 0

ccyt ðce Þ ¼ J 1 sl ðc e Þ:

ð66Þ

so

2. The RyR-term is given by

J RyR ðzÞ ¼ kRyR ðzÞmax 0;

ð61Þ

Since we are seeking stationary and constant in space solutions where the RyR-release is about to activate, it follows from (51) that

cm cyt ð48Þ

ð59Þ

cm cyt

!

!

bcyt ðce Þ ¼

m  p þ kleak ðzÞ ðc SR  c cyt Þ;

cm cyt þ K a

ð51Þ

¼

C cytbuf ccyt ccyt þ K cytbuf C cytbuf J 1 sl ðc e Þ 1 J sl ðce Þ þ K cytbuf

ð67Þ :

ð68Þ

106

A. Tveito et al. / Mathematical Biosciences 236 (2012) 97–107

Furthermore, from (58), we get

kon2 FðbSR Þccyt ð1  pÞ ¼ koff 2 p

ð69Þ

Then, by using (61) and (62) we get

kon2 p J m ðce Þ ¼ sl m FðbSR ÞJ m sl ðc e Þ ¼ 1p koff 2 Ka

ð70Þ

so, using Assumption 3, we arrive at

bSR ¼ F 1



 koff 2 : kon2 K m a

ð71Þ

We observe that also in this quite general case, bSR is again a constant not depending on ce . From (57) we have

cSR C SRbuf  ðcSR þ K SRbuf ÞbSR ¼ 0

ð72Þ

so

cSR ¼

K SRbuf bSR C SRbuf  bSR

ð73Þ

and thus

cSR ¼

K SRbuf F

1



C SRbuf  F 1

koff 2 kon2 K m a





koff 2 kon2 K m a

:

ð74Þ

We note that also cSR is a constant independent of ce . It follows from (59) that J RyR ¼ J SERCA and because of (62) it follows from (51) that

v SERCA Sðccyt ; cSR Þ ¼ kleak ðzÞðcSR  ccyt Þ; so

v SERCA ðce Þ ¼

kleak ðzÞðcSR  ccyt Þ Sðccyt ; cSR Þ

kleak @

v  ðce Þ ¼

K SRbuf F 1



C SRbuf F 1

koff 2 kon2 K m a



0

 10

koff 2 kon2 K m a

A@

S@J1 sl ðc e Þ;

5.2. The Ryanodine Receptor Mutation of the ryanodine receptor is capable of promoting calcium waves and ventricular arrhythmia. The best-described pathology of this type is the autosomal dominant form of catecholaminergic polymorphic ventricular tachycardia (CPVT). In this syndrome the RyR pool exhibits heightened sensitivity to cytosolic calcium, as reflected by increased open probability [26]. This reduces the SR calcium load required to initiate a calcium wave, and causes these cardiac cells to be unstable during b-adrenergic stimulation. Our modifications to K a compare well with experimental data describing CPVT. This can be seen in Fig. 4, where increasing K a (decreasing RyR calcium sensitivity) increases the stable range of cellular and SR calcium load. 5.3. SERCA

Again, we have found the strength of the SERCA-pump as a function of the extracellular calcium concentration ce at the point in statespace where the RyR-release activates,

0

of myocyte calcium-handling, and it is useful to frame the following discussions with a comment on those differences and how they may impact stability. First, rodents exhibit a greater reliance upon SR calcium release for contraction than do higher order species, including rabbits and humans. This increased SR release is fueled largely by a higher SR calcium load, which itself results from enhanced SERCA expression, and reduced NCX function [2]. Both of these differences shift the competitive balance between NCX-mediated calcium extrusion and SR reuptake strongly in favor of reuptake. Based upon our model results, this shift in balance would be expected to promote a less stable diastolic state in the rodent, but it has been clearly shown that these events also occur in higher order species, particularly the rabbit [24,25]. So our overall expectation is that the model analysed here may become unstable with less severe changes to v SERCA and ce , than in a rabbit- or human-specific model. However, rabbit and human models would also be expected to exhibit a greater arrhythmogenic current response to those instabilities due to greater forward mode Na/Ca exchange.

K SRbuf F 1



C SRbuf F 1

K SRbuf F 1



C SRbuf F 1

koff 2 kon2 K m



a

koff 2 kon2 K m

koff 2 kon2 K m a



 

 a1

koff 2 kon2 K m a

1 A J 1 sl ðc e Þ :

ð75Þ

A

This demonstrates that we can compute the critical curve for the general class of models written in the form given by Eqs. (54)– (58) provided that Assumptions 1–4 hold. We have, however, not performed numerical investigations covering cases other than those reported above.

Experimental manipulations that enhance SR calcium reuptake and SERCA function promote calcium waves, but generally require simultaneous calcium overload [27,28]. Both of these alterations occur during b-adrenergic stimulation, and for this reason it is not surprising that b-adrenergic agonists produce calcium waves and arrhythmia in CPVT and heart failure [26,24]. The results of the model suggest that when extracellular calcium concentration is increased  30% increased from baseline, maximal SERCA pump velocity must still more than double before calcium handling becomes unstable. This finding may explain why recent experimental investigations involving SERCA overexpression [29,30] did not notice increased wave propensity. These studies achieved approximately 37% and 150% increases in SERCA protein, and may therefore not have achieved an SR calcium load sufficient to elicit waves. 5.4. Sodium-Calcium Exchange

5. Discussion The physiological relevance of our findings are limited to the physiological relevance of the model we have analyzed, which as mentioned above is relatively simple. However, this model was constructed specifically to interrogate the determinants of calcium waves, and at least qualitatively, our analysis regarding the sensitivity of stability transition to specific parameter modifications are in good agreement with experimental data. 5.1. Species-Specific Calcium Stability Many of the parameter alterations we have imposed are analogous to the differences underlying species-specific characteristics

The manipulation of K NCX applied in Fig. 6B is roughly analogous to diminishing the thermodynamic gradient for forward-mode Na/ Ca exchange (as would occur under conditions of sodium overload). It is well known that intracellular sodium overload promotes calcium retention and elevates SR calcium load, particularly in rodent and failing cardiac myocytes [31,32]. Both of these effects would be expected to reduce diastolic calcium-handling stability, in line with our analytic observations. The prescribed manipulation of v NCX in Fig. 6C is more directly analogous to increased NCX protein expression, which is also a well-described manifestation of heart failure [33,34] and clearly contributes to arrhythmia associated with calcium waves [35,24,25]. Interestingly, our simple analysis lends quantitative support to the idea that elevated NCX expression is

A. Tveito et al. / Mathematical Biosciences 236 (2012) 97–107

an effective early compensation in heart failure. By inhibiting calcium wave initiation (i.e. extending the stable region in ce ; V SERCA ), increased NCX expression may improve diastolic stability in the short term. However, under acute b-adrenergic challenge, or with further pathological development, diastolic stability may still be exceeded, and in this circumstance increased v NCX will clearly magnify the resulting arrhythmogenic current development. 6. Conclusion We have analyzed a mathematical model of calcium handling in cardiomyocytes and found explicit formulas separating parameters leading to stable and unstable resting states of the model. The suggested stability curve is verified through a series of computations based on linear analysis where the eigenvalues of a linearized model are analyzed. Furthermore, the result is confirmed by using fully non-linear computations. The model we have analyzed is intentionally simple and does not include the action potential. By choosing this simple model we have been able to describe the transition from an unstable to a stable resting state for calcium handling with an analytic expression. Acknowledgements This project was supported by the National Center for Research Resources (5P41RR008605–19), the National Institute of General Medical Sciences (8 P41 GM103426–19) from the National Institutes of Health, and a fellowship from the Heart Rhythm Society (AGE). References [1] P. Swietach, K.W. Spitzer, R.D. Vaughan-Jones, Modeling calcium waves in cardiac myocytes: importance of calcium diffusion, Front. Biosci. 15 (2010) 661. [2] D.M. Bers, Excitation-Contraction Coupling and Cardiac Contractile Force, second ed., Kluwert Academic, Dordrecht, The Netherlands, 2001. [3] W.G. Wier, M.B. Cannel, J.R. Berlin, E. Marban, W.J. Lederer, Cellular and subcellular heterogeneity of ½Ca2þ i in single heart cells revealed by fura-2, Science 235 (1987) 325. [4] L.A. Venetucci, A.W. Trafford, S.C. O’Neill, D.A. Eisner, The sarcoplasmic reticulum and arrhythmogenic calcium release, Cardiovasc. Res. 77 (2) (2008) 285. [5] R. Sitsapesan, A.J. Williams, Regulation of current flow through ryanodine receptors by luminal Ca2þ , J. Membr. Biol. 159 (3) (1997) 179. [6] V. Lukyanenko, S. Subramanian, I. Györke, T.F. Wiesner, S. Györke, The role of luminal Ca2þ in the generation of Ca2þ waves in rat ventricular myocytes, J. Physiol 518 (Pt 1) (1999) 173. ´ az, S.C. O’Neill, D.A. Eisner, Reducing [7] L.A. Venetucci, A.W. Trafford, M.E. Dı ryanodine receptor open probability as a means to abolish spontaneous Ca2þ 2þ release and increase Ca transient amplitude in adult ventricular myocytes, Circ. Res. 98 (10) (2006) 1299. ´ az, A.W. Trafford, S.C. O’Neill, D.A. Eisner, Measurement of sarcoplasmic [8] M.E. Dı reticulum Ca2þ content and sarcolemmal Ca2þ fluxes in isolated rat ventricular myocytes during spontaneous Ca2þ release, J. Physiol. 501 (Pt 1) (May 1997) 3. [9] D. Jiang, B. Xiao, D. Yang, R. Wang, P. Choi, L. Zhang, H. Cheng, S.R.W. Chen, RyR2 mutations linked to ventricular tachycardia and sudden death reduce the threshold for store-overload-induced Ca2þ release (soicr), Proc. Natl. Acad. Sci. USA 101 (35) (2004) 13062–13067. [10] M. Rubart, D.P. Zipes, Mechanisms of sudden cardiac death, J. Clin. Invest. 115 (9) (2005) 2305. [11] J.N. Weiss, A. Karma, Y. Shiferaw, P-S. Chen, A. Garfinkel, Z. Qu, From pulsus to pulseless: the saga of cardiac alternans, Circ. Res. 98 (2006) 1244. [12] S.C. O’Neill, L. Miller, R. Hinch, D.A. Eisner, Interplay between SERCA and sarcolemmal Ca2þ efflux pathways controls spontaneous release of Ca2þ from

[13]

[14]

[15]

[16]

[17] [18]

[19]

[20] [21] [22]

[23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31] [32] [33] [34]

[35]

107

the sarcoplasmic reticulum in rat ventricular myocytes, J. Physiol. 559 (Pt 1) (2004) 121. J. Okada, S. Sugiura, S. Nishimura, T. Hisada, Three-dimensional simulation of calcium waves and contraction in cardiomyocytes using the finite element method, Am. J. Physiol. Cell Physiol. 288 (3) (2005) C510. P. Swietach, K.W. Spitzer, R.D. Vaughan-Jones, Ca2þ -mobility in the sarcoplasmic reticulum of ventricular myocytes is low, Biophys. J. 95 (3) (2008) 1412. H.R. Ramay, M.S. Jafri, W.J. Lederer, E.A. Sobie, Predicting local sr Ca2þ dynamics during Ca2þ wave propagation in ventricular myocytes, Biophys. J. 98 (11) (2010) 2515. M.K. Stokke, S.J. Briston, G.F. Jlle, I. Manzoor, W.E. Louch, L. yehaug, G. Christensen, D.A. Eisner, A.W. Trafford, O.M. Sejersted, I. Sjaastad, Ca2þ wave probability is determined by the balance between SERCA2-dependent Ca2þ reuptake and threshold sr Ca2þ content, Cardiovasc. Res. 90 (3) (2011) 503. A. Tveito, G.T. Lines, A condition for setting off ectopic waves in computational models of excitable cells, Math. Biosci. 213 (2008) 141. A. Tveito, O. Skavhaug, G.T. Lines, R. Artebrant, Computing the stability of steady-state solutions of mathematical models of the electrical activity in the heart, Comput. Biol. Med. 41 (8) (2011) 611. A. Tveito, G.T. Lines, R. Artebrant, O. Skavhaug, M.M. Maleckar, Existence of excitation waves for a collection of cardiomyocytes electrically coupled to fibroblasts, Math. Biosci. 230 (2) (2011) 79. C. Chicone, Ordinary Differential Equations with Applications, volume 34 of Texts in Applied Mathematics. Springer, 1999. M.D. Stern, Buffering of calcium in the vicinity of a channel pore, Cell Calcium 13 (3) (1992) 183. G.D. Smith, J.E. Keizer, M.D. Stern, W.J. Lederer, H. Cheng, A simple numerical model of calcium spark formation and detection in cardiac myocytes, Biophys. J. 75 (1) (1998) 15. S. Györke, I. Györke, V. Lukyanenko, D. Terentyev, S. Viatchenko-Karpinski, T.F. Wiesner, Regulation of sarcoplasmic reticulum calcium release by luminal calcium in cardiac muscle, Front. Biosci. 7 (2002) d1454. S.M. Pogwizd, K. Schlotthauer, L. Li, W. Yuan, D.M. Bers, Arrhythmogenesis and contractile dysfunction in heart failure: roles of sodium-calcium exchange inward rectifier potassium current and residual beta-adrenergic responsiveness, Circ. Res. 88 (11) (2001) 1159. K. Schlotthauer, D.M. Bers, Sarcoplasmic reticulum Ca2þ release causes myocyte depolarization. underlying mechanism and threshold for triggered action potentials, Circ. Res. 87 (9) (2000) 774. M. Fernández-Velasco, A. Rueda, N. Rizzi, J-P. Benitah, B. Colombi, C. Napolitano, S.G. Priori, S. Richard, A.M. Gómez, Increased Ca2þ sensitivity of the ryanodine receptor mutant ryr2r4496c underlies catecholaminergic polymorphic ventricular tachycardia. Circ. Res. 104 (2) (2009) 201–209, 12p following 209. J.A. Wasserstrom, Y. Shiferaw, W. Chen, S. Ramakrishna, H. Patel, J.E. Kelly, M.J. O’Toole, A. Pappas, N. Chirayil, N. Bassi, L. Akintilo, M. Wu, R. Arora, G.L. Aistrup, Variability in timing of spontaneous calcium release in the intact rat heart is determined by the time course of sarcoplasmic reticulum calcium load. Circ. Res. 107 (9) (2010) 1117–1126. Katsuji Fujiwara, Hideo Tanaka, Hiroki Mani, Takuo Nakagami, Tetsuro Takamatsu, Burst emergence of intracellular Ca2þ waves evokes arrhythmogenic oscillatory depolarization via the Naþ -Ca2þ exchanger: simultaneous confocal recording of membrane potential and intracellular Ca2þ in the heart, Circ. Res. 103 (5) (2008) 509. D.L. Baker, K. Hashimoto, I.L. Grupp, Y. Ji, T. Reed, E. Loukianov, G. Grupp, A. Bhagwhat, B. Hoit, R. Walsh, E. Marban, M. Periasamy, Targeted overexpression of the sarcoplasmic reticulum Ca2þ -atpase increases cardiac contractility in transgenic mouse hearts, Circ. Res. 83 (12) (1998) 1205. R.J. Hajjar, J.X. Kang, J.K. Gwathmey, A. Rosenzweig, Physiological effects of adenoviral gene transfer of sarcoplasmic reticulum calcium atpase in isolated rat myocytes, Circulation 95 (2) (1997) 423. D.M. Bers, Altered cardiac myocyte Ca2þ regulation in heart failure, Physiology (Bethesda) 21 (2006) 380–387. D.M. Bers, S. Despa, J. Bossuyt, Regulation of Ca2þ and Naþ in normal and failing cardiac myocytes, Ann. N. Y. Acad. Sci. 1080 (2006) 165–177. T.R. Shannon, D.M. Bers, Integrated Ca2þ management in cardiac myocytes, Ann. N. Y. Acad. Sci. 1015 (2004) 28–38. S.M. Pogwizd, M. Qi, W. Yuan, A.M. Samarel, D.M. Bers, Upregulation of Naþ / Ca2þ exchanger expression and function in an arrhythmogenic rabbit model of heart failure, Circ. Res. 85 (11) (1999) 1009. M.K. Stokke, K. Hougen, I. Sjaastad, W.E. Louch, S.J. Briston, U.H. Enger, K.B. Andersson, G. Christensen, D.A. Eisner, O.M. Sejersted, A.W. Trafford, Reduced SERCA 2 abundance decreases the propensity for Ca2+ wave development in ventricular myocytes, Cardiovasc. Res. 86 (2010) 63.