Computational Modeling and Numerical Methods for Spatiotemporal Calcium Cycling in Ventricular Myocytes

Intracellular calcium (Ca) cycling dynamics in cardiac myocytes is regulated by a complex network of spatially distributed organelles, such as sarcoplasmic reticulum (SR), mitochondria, and myofibrils. In this study, we present a mathematical model of intracellular Ca cycling and numerical and computational methods for computer simulations. The model consists of a coupled Ca release unit (CRU) network, which includes a SR domain and a myoplasm domain. Each CRU contains 10 L-type Ca channels and 100 ryanodine receptor channels, with individual channels simulated stochastically using a variant of Gillespie’s method, modified here to handle time-dependent transition rates. Both the SR domain and the myoplasm domain in each CRU are modeled by 5 × 5 × 5 voxels to maintain proper Ca diffusion. Advanced numerical algorithms implemented on graphical processing units were used for fast computational simulations. For a myocyte containing 100 × 20 × 10 CRUs, a 1-s heart time simulation takes about 10 min of machine time on a single NVIDIA Tesla C2050. Examples of simulated Ca cycling dynamics, such as Ca sparks, Ca waves, and Ca alternans, are shown.

1. Introduction

Calcium (Ca) signaling is fundamental to many biological functions (Berridge et al., 2000). In cardiac myocytes, Ca cycling is regulated by a diffusively coupled network of Ca release units (CRUs), which plays a central role in cardiac excitation-contraction coupling and heart rhythms (Bers, 2002; Lakatta et al., 2010). The elementary Ca cycling event is a Ca spark (Cheng et al., 1993; Cheng and Lederer, 2008). Ca sparks are discretized Ca release events due to random and collective openings of the ryanodine receptor (RyR) channels clustered in a CRU. A typical cardiac myocyte includes about 10,000–20,000 CRUs, and the spatial arrangement of CRUs varies widely across different myocyte types and changes in diseased conditions (Franzini-Armstrong et al., 1999; Soeller and Cannell, 1999; Chen-Izu et al., 2007; Soeller et al., 2007; Wei et al., 2010). The mathematical modeling of subcellular Ca cycling presents many numerical challenges that are unique to the Ca subsystem.

One challenge is the fact that Ca concentration gradients between the myoplasmic space and other intra and intercellular subspaces are extremely large. Ca concentrations outside the cell in the extracellular space are roughly 1.8 mM, but prolonged exposure to high Ca concentrations in the myoplasm are toxic to the cell. Thus, intracellular Ca is highly buffered (Bers, 2002) and the cell expends energy through membrane pumps to maintain a steady state free Ca level of about 0.1 μM in the myoplasm. However, high Ca levels are necessary (and need to be accessed quickly) in order to interact with the myofibrils and cause contraction. The cell’s solution is an intracellular storage space called the sarcoplasmic reticulum (SR). Steady state levels of free Ca in the SR are in the range of 1 mM, about the same order of magnitude of Ca concentration that exists outside the cell in the extracellular space, and around four orders of magnitude larger than in the myoplasm. The large concentration gradient results is an extremely high Ca flux whenever a Ca channel opens to connect the myoplasm to the SR or extracellular spaces. This is optimal for cellular function since a large amount of Ca enters the myoplasm quickly, as is necessary for contraction. But this causes problems for numerical simulation as the high fluxes cause the differential equations to be numerically stiff, requiring a very small time step. In addition, there exist multiple time scales in the various fluxes, and Ca is released into subspaces of the myoplasm which have extremely small volumes, presenting further numerical challenges.

Another challenge is that the dynamics of subcellular Ca are spatially dependent. Unlike other ionic concentrations, such as sodium (Na) and potassium (K), for which whole-cell averages provide reasonable approximations, subcellular Ca release is locally controlled through Ca-induced Ca release (CICR). This becomes especially important when trying to understand arrhythmogenic phenomena such as delayed afterdepolarizations which are triggered by Ca waves. To gain insight into the initiation of Ca waves through local Ca releases and to understand the complex wave dynamics that ensue, any realistic mathematical model of Ca dynamics must include the underlying spatial structure of the Ca signaling network. This implies that instead of ordinary differential equations (ODEs) one must resort to partial differential equations (PDEs) in order to capture the spatial diffusion of Ca in both the myoplasm and the SR spaces.

A third major challenge is that subcellular Ca release is inherently stochastic. A CRU consists of 50–200 RyR channels and 5–10 L-type Ca channels (LCCs) (Soeller et al., 2007; Bers, 2008). The small number of channels involved means that averages do not well approximate the dynamics governed by the law of mass action. This suggests that instead of deterministic differential equations, one should use continuous-time discrete-state Markov models, which are notorious for their cost to solve numerically. Furthermore, the transition rates are typically dependent on state variables such as Ca concentrations in the myoplasm and SR, and are thus implicitly time-dependent, further adding to the difficulty in numerical simulation.

Despite the challenges in modeling cardiac Ca cycling and excitation-contraction coupling, mathematical, and computational models have been developed at many scales using different mathematical and computational means (Luo and Rudy, 1994; Jafri et al., 1998; Greenstein and Winslow, 2002; Sobie et al., 2002; Shiferaw et al., 2003; Wang et al., 2005; Izu et al., 2006; Groff and Smith, 2008; Restrepo et al., 2008; Williams et al., 2008, 2011; Huertas et al., 2010; Ramay et al., 2010; Rovetti et al., 2010; Wasserstrom et al., 2010), including: (1) low-dimensional deterministic models described by ODEs; (2) single CRU models with stochastically simulated individual RyRs; (3) ODE models with stochastically simulated individual RyRs coupled to a common myoplasmic pool; (4) high-dimensional deterministic models described by PDEs and stochastically simulated CRU firing (fire-diffusion-fire type models), and (5) coupled CRU network models with stochastically simulated individual RyRs. As the complexity of the model increases, the computational demands also increase tremendously due to the aforementioned issues. However, most of the above models have been implemented on central processing units (CPUs) and have used direct numerical methods for the simulation of stochastic channel openings. These numerical methods limit the use of highly resolved spatial discretizations and scale of simulations (for example, to simulate 20,000 CRUs in a whole-cell with accurate spatial diffusion of Ca). Therefore, the development of computationally efficient models and advanced numerical algorithms is important.

Here we present a recently developed three-dimensional (3D) spatially distributed Ca cycling model in which CRUs are locally coupled by Ca diffusion throughout the myoplasmic (Myo) and sarcoplasmic reticulum (SR) domains, improved from our previous CRU network model (Rovetti et al., 2010). The Ca cycling model has been coupled to the action potential (AP) model of the rabbit ventricular myocyte (Mahajan et al., 2008), allowing for the study of excitation-contraction coupling. The major improvements to the Ca cycling model from Rovetti et al. (2010) include: (1) 100 × 20 × 10 CRUs which corresponds to the number of CRUs in a typical myocyte and the Myo and SR spaces are modeled as true 3D spaces. (2) The phenomenological LCC model (Rovetti et al., 2010) has been replaced by a 7-state model based on experimental patch-clamp data obtained in isolated rabbit ventricular myocytes (Mahajan et al., 2008). In addition, we used a time-dependent version of Gillespie’s method for the fast computation of the stochastic opening of the RyRs and LCCs. All numerical algorithms were implemented using graphical processing units (GPUs). These advanced numerical and computational methods allow us to use highly resolved Myo and SR spaces for accurate Ca diffusion. We first present the mathematical details of the model, and then the advanced numerical algorithms we developed and/or employed for their solution. Finally we show how these algorithms are implemented using GPUs and the computational efficacy is discussed. Simulation results of Ca cycling and excitation-contraction coupling are shown as examples.

2. Materials and Methods

2.1. Mathematical Model

A ventricular myocyte is composed of a spatially distributed complex network of consisting of the SR, mitochondria, and myofibrils, among other organelles, as in Figure 1A. A t-tubular system (Figure 1B) facilitates effective communication of this network with the extracellular space. Ca diffuses in both the Myo and the SR, the latter of which is an interconnected network inside the cell distinguished as junctional SR (jSR) and network SR (Figure 1A). A unifying cardiac excitation-contraction picture is illustrated in Figure 1C (Bers, 2002). In a normal action potential (see Figure 1D), voltage-dependent opening of the LCCs brings Ca into the dyadic space (DS), a very small space between the LCC cluster and jSR (shaded area in Figure 1C). Elevated Ca concentrations in the vicinity of the LCCs causes their inactivation. The RyR channels open stochastically and their open probability is sensitive to Ca in the DS, a process called CICR. Therefore, the RyR channels can be triggered by Ca entry from the LCCs, high Myo and SR Ca, and Ca diffusing from neighboring CRUs. Ca entered from the LCCs and released from the SR diffuses to the myofibrils to signal contraction and participates in many other signaling processes in the myocytes. Ca is pumped back into the SR by the sarcoplasmic-endoplasmic reticulum Ca ATPase (SERCA) pump, and extruded by Na-Ca exchange (NCX). Ca is also uptaken by mitochondria through the mitochondrial uniporter and released from mitochondria via NCX in the mitochondrial membrane and opening of other channels, such as the mitochondrial permeation transition pore. LCC and NCX couple Ca and voltage bi-directionally, but all other currents also affect this coupling either indirectly via their effects on voltage or directly via Ca regulation of the ion channels.

FIGURE 1

Figure 1. (A) Detailed intracellular structure in a cardiac myocyte (Katz, 2011). (B) T-tubule network of a rat ventricular myocyte (Soeller and Cannell, 1999). (C) Illustration of cardiac excitation-contraction system. (D) A typical action potential and Ca transient.

Modeling the complete detailed structure and Ca signaling of a myocyte is much too complicated and challenging, even using the most advanced computational technologies. Instead of modeling the complex detailed structure of the cell, we used a simplified approach in which we model the cell using a two-domain structure, the Myo and SR domains (see Figure 2). We assume Ca freely diffuses throughout Myo and SR domains, which is mathematically modeled by the diffusion equation. In computer simulation, the Myo and SR domains are discretized depending upon the spatial accuracy desired (see Section 2). The Myo and SR domains are coupled via SR release and uptake. Each CRU contains a Myo space, a SR space, a jSR which is diffusively connected to the SR, and a DS which is diffusively connected to the Myo. Extracellular Ca enters the DS through voltage-gated LCCs, which open stochastically and are simulated by a Markov model (see Figure 3). Ca is released from the jSR through its associated cluster of RyRs to the DS. The RyRs also open stochastically and are simulated using a Markov model in which activation and inactivation of RyRs are regulated by Ca in the DS (see Figure 3). Ca is either extruded from the cell via the Na-Ca exchanger, or taken back up into the SR via the SERCA pump.

FIGURE 2

Figure 2. The spatially distributed Ca cycling model. (A) Plot of a 3D coupled CRU network. (B) Detailed illustration of a CRU.

FIGURE 3

Figure 3. The Markov models for LCCs and RyRs. (A) 7-state LCC model. (B) 4-state RyR model.

To couple the Ca cycling to voltage dynamics, we assume that the voltage is uniform across the cell membrane. The ionic currents are taken from the Mahajan model (Mahajan et al., 2008), except the Ca related currents which are computed through summing the local Ca currents from the spatial Ca cycling model. The Ca cycling parameters are mainly based on Rovetti et al. (2010) with modifications to account for the 3D geometry (see Tables 1–5). The parameters for the ionic currents used in the AP model and transition rates used in the LCC model are from Mahajan et al. (2008).

TABLE 1

Table 1. Parameter values for diffusion and buffering equations.

TABLE 2

Table 2. Parameter values for jSR and DS fluxes.

TABLE 3

Table 3. Parameter values for SERCA (uptake), NCX, and background leak fluxes.

TABLE 4

Table 4. Parameter values for LCCs.

TABLE 5

Table 5. Parameter values for RyR channels.

2.1.1. Ca cycling model

The time evolution of the concentration of Ca in the 3D Myo and SR domains is modeled using a system of reaction-diffusion equations. Reaction terms come in two types: continuous flux terms couple the two domains and extrude Ca from the Myo domain. Flux terms at discrete locations couple the two domains via the CRUs through a jSR compartment and a DS compartment.

The system of equations for the evolution of Ca in the cell is

$\begin{array}{rcll}{\beta }_m\left(c_m\right)\frac{\partial c_m}{\partial t}& =D_m{\nabla }^2{c}_m+J_m,& & \text{(1)}\\ {\beta }_s\left(c_s\right)\frac{\partial c_s}{\partial t}& =D_s{\nabla }^2{c}_s+J_s,& & \text{(2)}\\ {\beta }_d\left(c_d\right)\frac{d{c}_d^{\left(i\right)}}{dt}& =J_d^{\left(i\right)},& & \text{(3)}\\ {\beta }_j\left(c_j\right)\frac{d{c}_j^{\left(i\right)}}{dt}& =J_j^{\left(i\right)},& & \text{(4)}\end{array}$

where c_m(x, y, z, t) and c_s(x, y, z, t) are the local Ca concentrations in the Myo and SR, respectively, and $c_d^{\left(i\right)}$ and $c_j^{\left(i\right)}$ are the Ca concentration in the ith DS and ith jSR, respectively. The Myo and SR domains have diffusion coefficients D_m and D_s. We assume Ca is buffered in the Myo (β_m), SR (β_s), DS (β_d), and jSR (β_j) by both calmodulin and SR-bound proteins in the Myo and DS, and calsequestrin in the SR and jSR (Shannon et al., 2004). We also assume that such buffering occurs rapidly enough to be modeled with instantaneous capacity functions (Wagner and Keizer, 1994)

$\begin{array}{rcll}{\beta }_m\left(c\right)& =1+\frac{B_{sr}{K}_{sr}}{{\left(c+K_{sr}\right)}^2}+\frac{B_{cd}{K}_{cd}}{{\left(c+K_{cd}\right)}^2},& & \text{(5)}\\ {\beta }_d\left(c\right)& =1+\frac{B_{sr}^{\prime }{K}_{sr}^{\prime }}{{\left(c+K_{sr}^{\prime }\right)}^2}+\frac{B_{cd}^{\prime }{K}_{cd}^{\prime }}{{\left(c+K_{cd}^{\prime }\right)}^2},& & \text{(6)}\\ {\beta }_s\left(c\right)& =1+\frac{B_{cal}{K}_{cal}}{{\left(c+K_{cal}\right)}^2},& & \text{(7)}\\ {\beta }_j\left(c\right)& =1+\frac{B_{cal}^{\prime }{K}_{cal}^{\prime }}{{\left(c+K_{cal}^{\prime }\right)}^2}.& & \text{(8)}\end{array}$

The terms J_m, J_s, $J_d^{\left(i\right)},$ and $J_j^{\left(i\right)}$ represent the net current for each Ca space and are specified below. In order to account for the different volumes of the Myo, DS, SR, and jSR, the magnitudes of various fluxes between spaces are rescaled by the ratio of the appropriate volume elements v_m, v_d, v_s, and v_j, the local volumes of the Myo, DS, SR, and jSR respectively.

Myoplasm flux. Ca enters and leaves the Myo due to uptake, exchange, and background leak fluxes, with net flux

$\begin{array}{rcll}{J}_m& =-J_{up}-J_{NCX}+J_{bg}+J_{Srleak}& & \\ & +\frac{v_d}{v_m}\sum _i\delta \left(x-x^{\left(i\right)},y-y^{\left(i\right)},z-z^{\left(i\right)}\right)J_{ds}^{\left(i\right)}& & \text{(9)}\end{array}$

The proteins mediating the SR uptake, NCX, and cell background and SR leak fluxes are generally found distributed evenly throughout their respective regions, thus we model the corresponding fluxes as spatially continuous functions of position.

The SERCA uptake pump is modeled to be driven solely by the Myo Ca with the simple Hill function

$J_{up}=v_{up}\frac{c_m^2}{c_m^2+k_{up}^2}$(10)

with a maximum velocity v_up and a half-maximal concentration k_up.

The Ca flux through the NCX pump, J_NCX, known to be a reversible flux sensitive to both Myo Ca and membrane voltage, is modeled following the physiological model of Weber et al. (2001):

$\begin{array}{rcll}{J}_{NCX}& =v_2\frac{\left(e^{\frac{\eta VF}{RT}}{\left[Na\right]}_i^2{\left[Ca\right]}_o-e^{\frac{\left(\eta -1\right)VF}{RT}}{\left[Na\right]}_o^3{c}_m\right)}{Deno{m}_{1}Deno{m}_2},& & \text{(11)}\\ Deno{m}_1& =\left(1+k_{sat}{e}^{\frac{\left(\eta -1\right)VF}{RT}}\right)\left(1+{\left(\frac{K_{mC{a}_{act}}}{c_m}\right)}^3\right),& & \text{(12)}\\ Deno{m}_2& =K_{mC{a}_o}{\left[Na\right]}_i^3+K_{mN{a}_o}^3{c}_m+K_{mN{a}_i}^3{\left[Ca\right]}_o& & \\ & \times \left(1+\frac{c_m}{K_{mC{a}_i}}\right)+K_{mC{a}_i}{\left[Na\right]}_o^3\left(1+\frac{{\left[Na\right]}_i^3}{K_{mN{a}_i}^3}\right)& & \\ & +{\left[Na\right]}_i^3{\left[Ca\right]}_o+{\left[Na\right]}_o^3{c}_m,& & \text{(13)}\end{array}$

where V is the net membrane potential and [Na]_i is the intracellular Na concentration, both of which are dynamic variables (see Section 2). The symbols [Ca]_o and [Na]_o represent the extracellular Ca and Na concentrations, respectively, which are assumed to be constant.

A background Ca leak flux J_bg is driven by the membrane potential and Myo and external Ca concentrations in the form of the Nernst equation taken from Luo and Rudy (1994)

$J_{bg}=-g_{bg}\frac{F}{2RT}\left(V-log\left(\frac{{\left[Ca\right]}_o}{c_m}\right)\right).$(14)

It has been shown experimentally that there is little to no voltage drop across the SR membrane and the Myo, therefore a passive Ca leak flux J_SRleak out of the SR and into the Myo is modeled by a simple diffusive current dependent upon the Ca gradient

$J_{SRleak}=g_{SRleak}\left(c_s-c_m\right).$(15)

The ionic flux $J_{ds}^{\left(i\right)}$ between the local Myo and the ith DS, defined below, is a spatially discrete flux, present only at position (x⁽ⁱ⁾, y⁽ⁱ⁾, Z⁽ⁱ⁾).

SR flux. Ca enters and leaves the network SR due to uptake, leak, and diffusion, with net value

$J_s=\frac{v_m}{v_s}\left(J_{up}-J_{SRleak}\right)-\frac{v_j}{v_s}\sum _i\delta \left(x-x^i,y-y^i,z-z^i\right)J_{jsr}^{\left(i\right)}.$(16)

The equation for the spatially discrete jSR flux $J_{jsr}^{\left(i\right)}$ is given below.

Junctional SR flux. The net flux of Ca for the ith jSR is

$J_j^{\left(i\right)}=J_{jsr}^{\left(i\right)}-\frac{v_d}{v_j}{J}_{ryr}^{\left(i\right)},$(17)

in which $J_{ryr}^{\left(i\right)}$ is described below, and

$J_{jsr}^{\left(i\right)}=g_{jsr}\left(c_s\left(x^{\left(i\right)},y^{\left(i\right)},z^{\left(i\right)}\right)-c_j^{\left(i\right)}\right)$(18)

gives a first-order refilling of Ca form the network SR to the jSR with refilling rate g_jsr.

Dyadic space flux. The net flux of Ca for the ith DS, including the fluxes through the RyRs and LCCs, is

$J_d^{\left(i\right)}=-J_{ds}^{\left(i\right)}+J_{lcc}^{\left(i\right)}+J_{ryr}^{\left(i\right)}.$(19)

The flux from the DS to the local Myo is dependent on the local Ca gradient with flux rate g_ds:

$J_{ds}^{\left(i\right)}=g_{ds}\left(c_d^{\left(i\right)}-c_m\left(x^{\left(i\right)},y^{\left(i\right)},z^{\left(i\right)}\right)\right).$(20)

Ca flux into the DS, via either the LCC or RyR channel, is the major route of entry of Ca into the Myo. Due to the fact that this entry is governed by a small number of interacting ion channels in each DS, we represent the opening of both types of channels as simple Markov models, as described below.

L-type Ca channel model. Each DS contains a cluster of LCCs. The kinetics of each channel are represented with a 7-state Markov model, developed in Mahajan et al. (2008) based on experimental data from rabbit ventricular myocytes (see Figure 3). It includes 4 inactive states (I2_Ca, I1_Ca, I2_Ba, I1_Ba), 2 closed states (C2, C1), and one open state (O). The transition rates are dependent upon the membrane potential V and the Ca concentration in the local DS, as opposed to an averaged Ca concentration in the submembrane space as in the original model (Mahajan et al., 2008). There are n_L LCCs in a single DS, and we assume that each channel operates independently of the others. Let $N_L^{\left(i\right)}\in \left\{0,1,\dots ,n_L\right\}$ be a stochastic variable indicating the number of LCCs in the open state (O) at any given time in the ith DS. Then the total current through all open LCCs in the ith DS is given by the voltage-dependent Goldman-Hodgkin-Katz equation

$i_{lcc}^{\left(i\right)}=N_L^{\left(i\right)}{P}_{Ca}\frac{4F^2}{RT}V\frac{{ĉ}_d^{\left(i\right)}exp\left(\frac{2FV}{RT}\right)-{\beta }_{Ca}{\left(Ca\right)}_o}{exp\left(\frac{2FV}{RT}\right)-1},$(21)

where ${ĉ}_d^{\left(i\right)}$ is the Ca concentration in the ith dyadic space (in units of nM). The current $i_{lcc}^{\left(i\right)}$ is converted to the flux $J_{lcc}^{\left(i\right)},$ with units of μM/ms, through

$J_{lcc}^{\left(i\right)}=\frac{-i_{lcc}^{\left(i\right)}}{2F{v}_d}.$(22)

The dynamics of the stochastic variable $N_L^{\left(i\right)}$ are described in Section 4.

RyR channel model. Each jSR contains a cluster of RyR channels, facing the DS. The kinetics of each channel are represented with a 4-state Markov model, adapted from Stern et al. (1999; see Figure 3). The model describes the channel as having closed (C), open (O), inactivated (I), and refractory (R) states. The transition rates for the two horizontal components (C to O, and R to I) are equivalent, with forward rate $k_a^{+}{\left(c_d^{\left(i\right)}\right)}^2$ dependent on the square of the DS Ca, and constant reverse rate $k_a^{-}.$ Similarly, the transition rates for the two vertical transitions (O to I, and C to R) are equivalent, with forward rate $k_b^{+}{c}_d^{\left(i\right)}$ and reverse rate $k_b^{-}.$ In this model the RyR open probability is regulated explicitly by the Ca concentration in the DS.

There are n_R RyR channels in a single jSR, and we assume that each channel operates independently of the others. Let $N_R^{\left(i\right)}\in \left\{0,1,\dots ,n_R\right\}$ be a stochastic variable indicating the number of RyR channels in the open state (O) at any given time in the ith DS. We model the total flux flowing through the RyR cluster as dependent on the Ca gradient between the local jSR and DS given by

$J_{ryr}^{\left(i\right)}=N_R^{\left(i\right)}\frac{g_{ryr}}{v_d}\left(c_j^{\left(i\right)}-c_d^{\left(i\right)}\right),$(23)

where g_ryr is the single RyR channel flux rate. The dynamics of the stochastic variable $N_R^{\left(i\right)}$ are described in Section 4.

2.1.2. Ca and voltage coupling

We couple the above Ca cycling model with voltage by using the ionic current formulation of the rabbit ventricular myocyte model by Mahajan et al. (2008). The voltage, assuming it is uniform across the cell membrane, is described by the differential equation

$\begin{array}{rcll}{C}_m\frac{dV}{dt}& =I_{Na}+I_{to,f}+I_{to,s}+I_{Kr}+I_{Ks}+I_{K1}+I_{NaK}+I_{stim}& & \\ & +I_{lcc}+I_{NCX}+I_{bg},& & \text{(24)}\end{array}$

where C_m is the cell membrane capacitance. The currents I_Na, I_to,f, I_to,s, I_Kr, I_Ks, I_K1, I_NaK, and I_stim are assumed to be uniform across the cell membrane and are taken directly from their formulation in Mahajan et al. (2008). The current I_Ks depends on Ca and, for simplicity, we use the average whole-cell Ca which we compute directly from the spatial Ca cycling model. The currents I_lcc, I_NCX, and I_bg represent the average LCC, NCX, and Ca background currents, respectively. We compute these from the local LCC, NCX, and background fluxes in the Ca cycling model through $I_{lcc}=\alpha ∕N_{cru}{\sum }_i{J}_{lcc}^{\left(i\right)},$$I_{NCX}=\alpha ∕N_{cru}{v}_{cru}\int \int \int J_{NCX}dxdydz,$ and $I_{bg}=\alpha ∕N_{cru}{v}_{cru}\int \int \int J_{bg}dxdydz,$ where v_cru is the cytosolic volume (in units of μm³) of a single CRU and N_cru is the total number of CRUs. Here the fluxes in equations (21) and (11) are converted to currents by the conversion factor α = −aFv_cell/C_m, where a is the ionic charge of the current carrier, v_cell is the volume of the myocyte in μl, and F is Faraday’s constant.

2.2. Numerical Algorithms

To integrate the reaction-diffusion system described by equation (1), we use an operator splitting scheme, described in Qu and Garfinkel (1999), applying first the flux terms and then the diffusion terms. The remaining state variables in the model are the stochastic ion channel states which are updated using Gillespie’s method (Gillespie, 2007), a new formulation of which is introduced here to handle time-dependent rates. We use the first-order forward Euler method to integrate the fluxes and diffusion (adapted to the multiple time scales involved). For the AP model, we first employ the standard Rush-Larson method (Rush and Larsen, 1978) to the quasi-linear ODEs before applying the Euler step. We choose this brute force approach because of its ease of implementation on GPUs. This is due to the fact that forward Euler is explicit and thus easily parallelized in an efficient manner. More complex implicit methods (see Strang, 1968; Qu and Garfinkel, 1999), for instance) are typically used in a non-GPU setting, however, for most problems the speed gained through implementation on a GPU typically outweighs the speed gained from a larger time step allowed by a higher-order method (Sato et al., 2009). For example, we found that the implicit second-order accurate Crank-Nicholson method (Haberman, 2004) allowed for a much larger time step than that imposed by the stability condition of the forward Euler discretization, and was thus much faster on a normal CPU. However, this method involves the inversion of a tridiagonal matrix which does not lend itself to an efficient parallelization, thus, the explicit method of forward Euler remains faster in the GPU setting.

2.2.1. Flux and diffusion

For a time step of size Δt, the Ca fluxes are updated according to the rule

$\begin{array}{rcll}& c_m\to c_m+\Delta t{\beta }_m{\left(c_m\right)}^{-1}{J}_m,& & \text{(25)}\\ & c_s\to c_s+\Delta t{\beta }_s{\left(c_s\right)}^{-1}{J}_s,& & \text{(26)}\\ & c_d^{\left(i\right)}\to c_d+\Delta t{\beta }_d{\left(c_d\right)}^{-1}{J}_d,& & \text{(27)}\\ & c_j^{\left(i\right)}\to c_j+\Delta t{\beta }_j{\left(c_j\right)}^{-1}{J}_j.& & \text{(28)}\end{array}$

Assuming a uniform spatial discretization (x, y, z) = (iΔx, jΔx, kΔx), and defining

$\begin{array}{rcll}{F}_{-}^x=c^{i-1,j,k}-c^{i,j,k},F_{+}^x=c^{i+1,j,k}-c^{i,j,k},& & & \text{(29)}\\ F_{-}^y=c^{i,j-1,k}-c^{i,j,k},F_{+}^y=c^{i,j+1,k}-c^{i,j,k},& & & \text{(30)}\\ F_{-}^z=c^{i,j,k-1}-c^{i,j,k},F_{+}^z=c^{i,j,k}-c^{i,j,k+1},& & & \text{(31)}\end{array}$

the spatial diffusion step is applied to c_m and c_s according to

$\begin{array}{rcll}{c}^{i,j,k}& \to c^{i,j,k}+\frac{D_m\Delta t}{\Delta x^2}{\beta }_m{\left(c^{i,j,k}\right)}^{-1}\left(F_{-}^x+F_{+}^x+F_{-}^y+F_{+}^y+F_{-}^z+F_{+}^z\right),& & \text{(32)}\end{array}$

with c = c_m, c_s, respectively. We found a mesh spacing of Δx = 0.2 μm to provide a good realization of Ca diffusion.

There exist multiple time scales in the reaction-diffusion system described by equation (1): time scales associated with the spatially dependent Ca fluxes, time scales associated with Ca diffusion, time scales associated with the spatially independent ionic concentrations and membrane potential in the AP model, and time scales associated with the stochastic opening and closing of the RyR and L-type Ca channels. Thus, we would like our numerical integration scheme to take advantage of the multiple time scales involved. We choose a global time step Δt_g equal to the largest time scale involved, and introduce the smaller time steps (which are evenly divisible into the global time step Δt_g):

$\begin{array}{rcll}& \text{Ca}\text{Fluxes}\to \Delta t_J=\frac{\Delta t_g}{m_J},& & \text{(33)}\\ & \text{Ca}\text{Diffusion}\to \Delta t_D=\frac{\Delta t_g}{m_D},& & \text{(34)}\\ & \text{AP}\to \Delta t_{AP}=\frac{\Delta t_g}{m_{AP}},& & \text{(35)}\end{array}$

for some integers m_J, m_D, and m_AP. During each global time step, we update the AP variables m_AP times, apply the Ca flux m_J times according to equation (25), and apply the Ca diffusion m_D times according to equation (32), with each smaller time step nested within the next larger time step. The time step associated with the stochastic Ca channels is more complicated and is discussed in Section 4. The fastest time scale involved in the Ca fluxes is the flux into the DS, which we are able to eliminate using a quasi-steady state approximation (see Section 3 below). We find that taking a global time step of Δt_g= 0.01 ms and m_AP= 1, m_D= 1, and m_J= 1 provide stable and fairly accurate results, while m_AP= 1, m_D= 2, and m_J= 10 provide very accurate results while maintaining reasonable simulation times (see Section 3 for performance details). Note: not only does the size of the Ca diffusion time step depend on Δx, but so does the size of the Ca flux time step since the local volume elements for the Myo and SR, v_m and v_s, depend upon the spatial discretization.

2.2.2. Boundary conditions

We assume that Ca cannot diffuse past the cell membrane, and thus impose first-order no flux boundary conditions. For a rectangular cell 0 ≤ x ≤ N_x, 0 ≤ y ≤ N_y, 0 ≤ z ≤ N_z this amounts to the conditions $F_{-}^x=0$ if x = 0, $F_{+}^x=0$ if x = N_x, $F_{-}^y=0$ if y = 0, $F_{+}^y=0$ if y = N_y, $F_{-}^z=0$ if z = 0, and $F_{+}^z=0$ if z = N_z.

2.2.3. Quasi-steady state approximation

Due to the small volume of the DS, fluxes into this space are typically large and thus numerically costly to integrate. This can be seen directly by letting V_d → 0 in the flux equation for the DS. A standard way to overcome this issue (Murray, 2002) is to assume that the DS equilibrates on a much faster time scale than that of the other fluxes, known as a quasi-steady state approximation. Thus, instead of numerically integrating the equation for the DS, one sets the DS flux equal to zero

$J_d^{\left(i\right)}=0,$(36)

which can be analytically solved for the DS concentration $c_d^{\left(i\right)}.$

2.2.4. Ion channel kinetics

The state variables N_R and N_L, representing the number of RyRs and LCCs in the open state for a single CRU, are modeled as continuous-time discrete-state Markov process. In what follows we drop the index i (which denoted the ith CRU) with the understanding that the algorithm is independently executed for each CRU.

Each dyadic space is treated as a single-pool element, thus, each ion channel (of a given type) in the same CRU experiences the same dyadic space environment. Therefore, each channel in a cluster has identical transition rates and we can treat each channel type as a group, tracking only the number of channels in each possible state (Gillespie, 2007). Let N_s represent the number of states in the underlying Markov model with N_t possible transitions between states. Furthermore, let λ_jk with l ≤ j, k ≤ N_s, j ≠ k, denote the transition rate between state j and state k, and let N_j be the number of channels in the jth state.

Direct stochastic simulation method. The most direct method for the simulation of the above Markov process is to assume a sufficiently small time step dt such that at most one transition will take place during time dt, i.e., there exists exactly one j and one k such that N_j → N_j− 1 and N_k → N_k+ 1 or no change takes place. One then partitions the unit interval into N_t− 1 partitions of length N_jλ_jkdt and one partition of length 1 − ΣN_jλ_jkdt, where N_jλ_jkdt represents the transition probability from state j to state k during time dt. Finally, one chooses a number uniformly distributed on the unit interval to determine which transition takes place.

There are two problems with the above method. First, it requires the choice of a small time step dt, which must be small enough to handle the fastest time scale involved. For time-dependent rates this is inefficient since it does not allow one to take advantage of possible time periods for which the transition rates are low. For example, the local Ca concentrations during a Ca spark can be as much as 200 times higher than during periods of inactivity, thus, the Ca dependent transition rates (some of which depend on the square of the Ca concentration in the dyadic space) can vary by 4 orders of magnitude. In Rovetti et al. (2010), a modified version of the direct method was introduced which allows for a variable time step during a single global time step, however, the method still requires one to chose a small time scale, and thus still suffers from the second problem which affects the efficiency of every direct method: it only converges to the appropriate statistics in the limit dt → 0. Therefore, we would like to use a better approach known as Gillespie’s method. Gillespie’s method is exact in the sense that it does not require the choice of a small time step dt, yet it reproduces the appropriate statistics exactly. However, the classic version of Gillespie’s method only applies to time-independent processes.

Below we introduce a modification of Gillespie’s method to handle time-dependent transition rates in an efficient manner based on an integral formulation. A similar method based on the integral approach was first used to simulate single-channel kinetics in the Hodgkin and Huxley model of nerve membrane ion conductances (Clay and DeFelice, 1983). Here we provide an explicit algorithm for the integral approach along with a new formulation which makes use of the total transition rate to calculate the wait time to the first transition, as opposed to keeping track of the state of each channel separately. In the setting of subcellular Ca modeling, other (non-direct) methods have been implemented to account for time-dependent transition rates (see, for instance, Restrepo et al., 2008; Williams et al., 2008).

Time-adaptive Gillespie’s method. Suppose the above Markov chain with time-dependent total transition rate λ(t) = ΣN_iλ_ij is in a given state at time t = 0, and let T be the random variable representing the wait time from t = 0 to the first transition out of the current state. Then it is straight forward to show that the wait time T is exponentially distributed with probability density function

$p\left(t\right)=\lambda exp\left(-{\int }_0^t\lambda \left(\tau \right)d\tau \right).$(37)

It follows from the standard inversion generating method of Monte Carlo theory (Gillespie, 2007) that one can generate a realization of T satisfying the above density function by selecting a random number r that is uniformly distributed on the unit interval, and solving the equation

${\int }_0^T\lambda \left(\tau \right)d\tau =-ln\left(r\right)$(38)

for T. For time-independent transition rates the above can be explicitly solved for T

$\frac{d\lambda }{dt}=0\Rightarrow T=-\frac{ln\left(r\right)}{\lambda },$(39)

the classic wait time in Gillespie’s method for time-independent rates.

The above method is exact in the sense that no approximations have been made in the derivation of equation (38), thus, T exactly follows the statistics given by the probability density function in equation (37). In practice, however, the total transition rate λ is an implicit function of time through its dependence on other time-dependent state variables, and we must numerically approximate the integral in equation (38) in order to solve for T. Therefore, assume that we are using a global time step of size Δt to numerically integrate the Ca fluxes in our operator splitting scheme. Over a time step of length Δt we can treat λ(t) as a constant and use the normal Gillespie algorithm to calculate the wait time to the next transition. Let t_loc represent the local time during a single global time step Δt (see Figure 4). For t_loc< Δt, we continue to update the local time step according to Gillespie’s algorithm using equation (39). Once we the local time exceeds Δt we must then numerically integrate equation (38).

FIGURE 4

Figure 4. Adaptive time step for Gillespie’s Method with time-dependent rates.

The details of the algorithm are as follows:

1. Draw a random number r uniformly distributed on the unit interval and let L_r= −ln(r).

2. Set t_loc= 0. Numerically solve equation (38), which is equivalent to evaluating the quantity

$m=L_r-\lambda \left(t\right)\Delta t.$(40)

If m ≤ 0 then a transition occurs during the current global time step after a local time step of length

$T=\frac{L_r}{\lambda \left(t\right)}.$(41)

Let t_loc= T and proceed to Step 4. Otherwise proceed to Step 3.

3. Since m > 0, the local wait time is greater than the global time step Δt. On the next global time step continue to numerically solve equation (38), which is equivalent to letting

$m\to m-\lambda \left(t\right)\Delta t.$(42)

If m ≤ 0 then a transition occurs during the current global time step after a local time step of length

$T=\frac{m+\lambda \left(t\right)\Delta t}{\lambda \left(t\right)}.$(43)

Let t_loc= T and proceed to Step 4. Otherwise repeat the current step.

4. Calculate the individual transition probabilities N_iλ_ij and the total transition probability λ = ΣN_iλ_ij. Partition the interval [0,λ] into N_t partitions of length N_iλ_ij. Draw a random number s uniformly distributed on the interval [0, λ]. Make the transition corresponding to where s falls. Go to the next step.

5. Draw a random number r uniformly distributed on the unit interval and calculate

$T=-\frac{ln\left(r\right)}{\lambda \left(t\right)}.$(44)

If t_loc+ T < Δt then let t_loc → t_loc+ T and go back to Step 4. Otherwise, integrate what remains in the interval by letting

$L_r=-ln\left(r\right)-\left(\Delta t-t_{loc}\right)\lambda \left(t\right),$(45)

and start at Step 2 on the next global time step.

There are many advantages to the above method. For example, there may be times for which a transition does not take place for multiple global time steps, thus this method is numerically efficient.

2.3. GPU computation

2.3.1. Hardware

The integration of the complete mathematical model was done on GPUs. We performed simulations on single NVIDIA Tesla C2050 high-performance Fermi-based GPGPU supporting ECC and double precision. Four Tesla cards were installed in a system with two quad-core 2.53 GHz Intel Xeon processors and 16 GB of RAM. The programs were written in C++ using the CUDA API and we used the GNU C++ compiler version 4.4.3 and NVIDIA CUDA version 4.0.

For the implementation of the model we started from the basic approach outlined in our previous study on the use of GPU to simulate cardiac tissue.

2.3.2. Single precision versus double precision

One potential problem with previous generations of GPU cards was that they could only handle single precision problems. The single precision did not pose any problem in our previous study on cardiac tissue simulation, apart from some measures that had to be introduced to avoid the singularity in the Goldman-Hodgkin-Katz equation, whenever this expression was used to coumpute the equilibrium potential of an ion channel.

In the subcellular Ca cycling model however it can happen in some special cases, e.g., when very accurate results are required, that the time step becomes too small for single precision computation. In these cases we used the double precision capability of the Fermi-based GPUs. We found the use of double precision to increase simulation time by approximately 20%.

2.3.3. Coupling to the action potential model

To couple the subcellular Ca cycling model to the action potential model in Mahajan et al. (2008), we have to compute the average of the various Ca fluxes and Ca concentrations (see Sec. 2.2) for each time step of the AP model. When this averaging is performed on the CPU the simulation time is completely dominated by this computation (>86% of the simulation time is spent in computing the averages). We therefore used the parallel reduction algorithm provided in NVIDIA’s CUDA SDK (Harris, 2007) to compute the averages on the GPU. The ODEs of the AP model are evaluated on the CPU. Taking advantage of the fact that the CPU can compute parallel to GPU computations, the actual solution of the AP model does not add to the simulation time, therefore the simulation time taken by the AP model (see below; Sec. 4) is equal to the execution time of the parallel reduction algorithm.

2.3.4. Performance

One second of simulated subcellular Ca cycling with uniform time step equal to 0.01 ms, including computation of the whole-cell action potential model, took 466.8 s to run on a single Tesla C2050. In Figure 5 we tested the accuracy and speed of the simulations with respect to different time steps. We took a uniform time step equal to 0.0001 ms as the benchmark (smaller time steps produced virtually no change in the results) to which we compared the accuracy and simulation speed with respect to other time step schemes. We found that for all time steps smaller than 0.01 ms the steady state Ca concentration was extremely accurate, to within 0.001% of the benchmark steady state. However, depending upon the time step used, the accuracy in the peak Ca concentration during depolarization varied. We found that taking a uniform time step of 0.001 gave highly accurate results to within 3% of the benchmark peak Ca concentration while the speed increased by a factor of 10. Taking a time step of 0.01 ms produced results that were within 18% of the benchmark peak Ca concentration while the speed increased by a factor of 94. Taking advantage of the multiple time scales that exist within the system, we found that taking an AP time step of 0.01 ms, a diffusion time step of 0.005 ms, and a Ca flux time step of 0.001 ms gave accurate results (within 5% of the benchmark peak Ca concentration) while maintaining a reasonable simulation speed (17 times speed up from benchmark). For all time steps tested we found the voltage to be nearly identical to the benchmark except for sleight variations in the AP duration due to the differences in peak Ca transients (see Figure 5B). It should be noted that for time steps on the order of 0.0001 ms the use of double precision is required, which causes a 20% reduction in speed. For all time steps greater than or equal to 0.001 we found single precision to be adequate, though we used double precision for the time comparisons above.

FIGURE 5

Figure 5. Comparison of various time discretization schemes. (A) Average whole-cell Ca concentrations for four different time discretization schemes [see (C) for details of schemes]. (B) Corresponding voltage traces. (C) Errors and simulation speeds for each of the time discretization schemes. Note: double precision is necessary when using a time step equal to 0.0001 ms, which slows the simulation time by 20% when compared to single precision.

Contrary to our 2D and 3D cardiac tissue simulations (Sato et al., 2009), the diffusion part of the problem was no longer the absolute bottleneck, with the simulation time spread over diffusion, reaction, and the whole-cell current computation. In our Ca cycling simulations, 50% of the simulation time was spent in the (intracellular) reaction part, of which 70% computing the Ca flux and 30% updating the Ca channels; 31% of the total time was spent in the diffusion computation and 15% in the whole-cell current computation (due to the parallel reduction algorithm), including the computation of whole-cell Ca concentration by reduction. The remaining 4% was spent in data transfer. Figure 6 shows a comparison between the 2D cardiac tissue simulation reported in Sato et al. (2009) and our new subcellular Ca dynamics simulations.

FIGURE 6

Figure 6. Ratio of the execution time of different parts of the code of the subcellular Ca cycling model, compared to our implementation of cardiac tissue simulation on a GPU (Sato et al., 2009).

For the cardiac tissue simulations reported in Sato et al. (2009) we were able to speed up the simulations considerably by using texture memory for large two-dimensional problems (by up to 27% for 1000 × 1000 cells). For our subcellular calcium simulations however the use of texture memory only sped up the simulations by about 2.5%. The explanation is that the individual dimensions (500 × 100 × 50) of the subcellular calcium simulations are too small for texture memory to be efficient.

3. Results

We implemented the 3D Ca cycling model with 100 × 20 × 10 = 20,000 identical CRUs (500 × 100 × 50 grid points), simulating the CRU network corresponding to a complete cardiac myocyte with dimensions of 100 μm × 20 μm × 10 μm. Using this model, we were able to generate well-known features of Ca and AP dynamics of cardiac myocytes. For example, we could replicate the well-known Ca signaling hierarchy (Cheng and Lederer, 2008; Weiss et al., 2011): Ca quarks, Ca sparks, macrosparks (clusters of Ca sparks), and Ca waves (see Figure 7). Under pacing, we were able to reproduce the transition from regular Ca transients and AP durations at slow pacing cycle lengths to Ca and AP duration alternans during rapid pacing, see Figure 8 (parameter values can be found in the Tables 1–5).

FIGURE 7

Figure 7. Ca signaling hierarchy reproduced by the Ca cycling model. At t = 0 ms there are quarks (q), sparks (s), and clusters of sparks (c), commonly referred to as macrosparks. The cluster of sparks (c) at t = 0 ms eventually evolves into a Ca wave (w) by t = 120 ms. Note that it is a random process, as not all clusters of sparks evolve into Ca waves.

FIGURE 8

Figure 8. Average whole-cell Ca concentration and voltage superimposed on linescans of subcellular Ca. (A) Slow pacing resulting in spatially uniform Ca release and consistent Ca release from beat to beat. (B) Rapid pacing resulting in spatially non-uniform Ca release and alternans.

4. Discussion

In this manuscript, we presented a mathematical model of Ca cycling and its coupling to membrane voltage, as well as numerical algorithms for effective computer simulation. Our model is spatially detailed, simulating a CRU network that represents a ventricular myocyte. Using this model, we are able to recapitulate the well-known Ca signaling hierarchy and excitation-contraction coupling dynamics. With advanced numerical algorithms and computational technologies, we can perform simulations at the cellular scale, while maintaining a fine spatial resolution at the intracellular scale, with reasonable computational speed.

However, we would like to note that there are several limitations. (1) Our model is a CRU network arranged in a regular 3D grid, which is much simpler than the real structure of a ventricular myocyte. (2) Mitochondria and myofibrils, which are also spatially distributed in the cell, play important roles in Ca cycling dynamics (O’Rourke and Blatter, 2009), but are not included in the model. (3) We used a RyR model that includes dyadic space Ca dependent inactivation, which may not be physiologically correct (Liu et al., 2012). Different RyR models (Sobie et al., 2002; Restrepo et al., 2008; Chen et al., 2009) have been developed for simulating Ca sparks and whole-cell Ca dynamics. Although our model can reproduce Ca alternans and the well-known Ca signaling hierarchy in general, different RyR models may affect specific predictions of Ca dynamics. (4) The formulation of certain Ca fluxes may need to be improved. For example, the explicit SR leak term, such as equation (15), may not be needed since there is evidence that random RyR openings are sufficient to account for SR leak in real myocytes (Williams et al., 2011). (5) The RyR distribution inside a cell is heterogeneous (Baddeley et al., 2009), and new Ca dynamics may emerge due to the heterogeneous firing properties of the CRUs. To overcome these limitations, our model needs to be improved through corroboration with experimental data and the incorporation of more physiological details.

Nonetheless, our present study has developed advanced mathematical and computational methods, which can be used to effectively simulate spatially detailed Ca cycling and its coupling to membrane voltage over long time scales. We hope to further improve the model and the computational algorithms to allow for the study of intracellular Ca cycling in a tissue-scale environment, and eventually contribute to the development of multi-scale modeling approaches to cardiac excitation-contraction coupling and arrhythmias (Hunter and Nielsen, 2005; Bassingthwaighte and Chizeck, 2008; Qu et al., 2011), a grand challenge in biological modeling.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work is supported by NIH/NHLBI P01 HL078931 and a fellowship award for advanced researchers from the Swiss Foundation for Grants in Biology and Medicine (Enno de Lange).

References

Baddeley, D., Jayasinghe, I. D., Lam, L., Rossberger, S., Cannell, M. B., and Soeller, C. (2009). Optical single-channel resolution imaging of the ryanodine receptor distribution in rat cardiac myocytes. Proc. Natl. Acad. Sci. U.S.A. 106, 22275–22280.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bassingthwaighte, J., and Chizeck, H. (2008). The physiome projects and multiscale modeling. IEEE Signal Process. Mag. 25, 121–144.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Berridge, M., Lipp, P., and Bootman, M. (2000). The versatility and universality of calcium signalling. Nat. Rev. Mol. Cell Biol. 1, 11–21.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bers, D. (2002). Cardiac excitation-contraction coupling. Nature 415, 198–205.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bers, D. (2008). Calcium cycling and signaling in cardiac myocytes. Annu. Rev. Physiol. 70, 23–49.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Chen, W., Wasserstrom, J. A., and Shiferaw, Y. (2009). Role of coupled gating between cardiac ryanodine receptors in the genesis of triggered arrhythmias. Am. J. Physiol. Heart Circ. Physiol. 297, H171–H180.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cheng, H., Lederer, W., and Cannell, M. (1993). Calcium sparks: elementary events underlying excitation-contraction coupling in heart muscle. Science 262, 740–744.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cheng, H., and Lederer, W. J. (2008). Calcium sparks. Physiol. Rev. 88, 1491–1545.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Chen-Izu, Y., Ward, C. W., Stark, W., Banyasz, T., Sumandea, M. P., Balke, C. W., Izu, L. T., and Wehrens, X. H. T. (2007). Phosphorylation of RyR2 and shortening of RyR2 cluster spacing in spontaneously hypertensive rat with heart failure. Am. J. Physiol. Heart Circ. Physiol. 293, H2409–H2417.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Clay, J., and DeFelice, L. (1983). Relationship between membrane excitability and single channel open-close kinetics. Biophys. J. 42, 151–157.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Franzini-Armstrong, C., Protasi, F., and Ramesh, V. (1999). Shape, size, and distribution of Ca²⁺ release units and couplons in skeletal and cardiac muscles. Biophys. J. 77, 1528–1539.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gillespie, D. (2007). Stochastic simulation of chemical kinetics. Annu. Rev. Phys. Chem. 58, 35–55.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Greenstein, J. L., and Winslow, R. L. (2002). An integrative model of the cardiac ventricular myocyte incorporating local control of Ca²⁺ release. Biophys. J. 83, 2918–2945.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Groff, J. R., and Smith, G. D. (2008). Ryanodine receptor allosteric coupling and the dynamics of calcium sparks. Biophys. J. 95, 135–154.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Haberman, R. (2004). Applied Partial Differential Equations with Fourier Series and Boundary Value Problems. Upper Saddle River: Pearson Prentice Hall.

Harris, M. (2007). Optimizing Parallel Reduction in Cuda. NVIDIA Corporation.

Huertas, M. A., Smith, G. D., and Györke, S. (2010). Ca²⁺ alternans in a cardiac myocyte model that uses moment equations to represent heterogeneous junctional SR Ca²⁺. Biophys. J. 99, 377–387.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hunter, P., and Nielsen, P. (2005). A strategy for integrative computational physiology. Physiology 20, 316–325.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Izu, L. T., Means, S. A., Shadid, J. N., Chen-Izu, Y., and Balke, C. W. (2006). Interplay of ryanodine receptor distribution and calcium dynamics. Biophys. J. 91, 95–112.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jafri, M. S., Rice, J. J., and Winslow, R. L. (1998). Cardiac Ca²⁺ dynamics: the roles of ryanodine receptor adaptation and sarcoplasmic reticulum load. Biophys. J. 74, 1149–1168.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Katz, A. M. (2011). Physiology of the Heart. Philadelphia: Lippincott Williams and Wilkins.

Lakatta, E. G., Maltsev, V. A., and Vinogradova, T. M. (2010). A coupled system of intracellular Ca²⁺ clocks and surface membrane voltage clocks controls the timekeeping mechanism of the heart’s pacemaker. Circ. Res. 106, 659–673.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liu, O. Z., Lederer, W., and Sobie, E. A. (2012). Does the goldilocks principle apply to calcium release restitution in heart cells? J. Mol. Cell. Cardiol. 52, 3–6.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Luo, C., and Rudy, Y. (1994). A dynamical model of the cardiac ventricular action potential: I. Simulations of ionic currents and concentration changes. Circ. Res. 74, 1071–1096.

Pubmed Abstract | Pubmed Full Text

Mahajan, A., Shiferaw, Y., Sato, D., Baher, A., Olcese, R., Xie, L.-H., Yang, M.-J., Chen, P.-S., Restrepo, J. G., Karma, A., Garfinkel, A., Qu, Z., and Weiss, J. N. (2008). A rabbit ventricular action potential model replicating cardiac dynamics at rapid heart rates. Biophys. J. 94, 392–410.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Murray, J. D. (2002). Mathematical Biology I: An Introduction. Heidelberg: Springer.

O’Rourke, B., and Blatter, L. A. (2009). Mitochondrial Ca²⁺ uptake: tortoise or hare? J. Mol. Cell. Cardiol. 46, 767–774.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Qu, Z., and Garfinkel, A. (1999). An advanced algorithm for solving partial differential equations in cardiac conduction. IEEE Trans. Biomed. Eng. 46, 1166–1168.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Qu, Z., Garfinkel, A., Weiss, J. N., and Nivala, M. (2011). Multi-scale modeling in biology: how to bridge the gaps between scales? Prog. Biophys. Mol. Biol. 107, 21–31.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ramay, H. R., Jafri, M. S., Lederer, W. J., and Sobie, E. A. (2010). Predicting local SR Ca²⁺ dynamics during Ca²⁺ wave propagation in ventricular myocytes. Biophys. J. 98, 2515–2523.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Restrepo, J. G., Weiss, J. N., and Karma, A. (2008). Calsequestrin-mediated mechanism for cellular calcium transient alternans. Biophys. J. 95, 3767–3789.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rovetti, R., Cui, X., Garfinkel, A., Weiss, J. N., and Qu, Z. (2010). Spark-induced sparks as a mechanism of intracellular calcium alternans in cardiac myocytes. Circ. Res. 106, 1582–1591.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rush, S., and Larsen, H. (1978). A practical algorithm for solving dynamic membrane equations. IEEE Trans. Biomed. Eng. 25, 389–392.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sato, D., Xie, Y., Weiss, J., Qu, Z., Garfinkel, A., and Sanderson, A. (2009). Acceleration of cardiac tissue simulation with graphic processing units. Med. Biol. Eng. Comput. 47, 1011–1015.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Shannon, T. R., Wang, F., Puglisi, J., Weber, C., and Bers, D. M. (2004). A mathematical treatment of integrated ca dynamics within the ventricular myocyte. Biophys. J. 87, 3351–3371.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Shiferaw, Y., Watanabe, M., Garfinkel, A., Weiss, J., and Karma, A. (2003). Model of intracellular calcium cycling in ventricular myocytes. Biophys. J. 85, 3666–3686.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sobie, E. A., Dilly, K. W., dos Santos Cruz, J., Lederer, W. J., and Jafri, M. S. (2002). Termination of cardiac Ca²⁺ sparks: an investigative mathematical model of calcium-induced calcium release. Biophys. J. 83, 59–78.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Soeller, C., and Cannell, M. B. (1999). Examination of the transverse tubular system in living cardiac rat myocytes by 2-photon microscopy and digital image processing techniques. Circ. Res. 84, 266–275.

Pubmed Abstract | Pubmed Full Text

Soeller, C., Crossman, D., Gilbert, R., and Cannell, M. B. (2007). Analysis of ryanodine receptor clusters in rat and human cardiac myocytes. Proc. Natl. Acad. Sci. U.S.A. 104, 14958–14963.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stern, M. D., Song, L.-S., Cheng, H., Sham, J. S., Yang, H. T., Boheler, K. R., and Ros, E. (1999). Local control models of cardiac excitation contraction coupling. J. Gen. Physiol. 113, 469–489.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Strang, G. (1968). On the construction and comparison of difference schemes. SIAM J. Numer. Anal. 5, 506–517.

CrossRef Full Text

Wagner, J., and Keizer, J. (1994). Effects of rapid buffers on Ca²⁺ diffusion and Ca²⁺ oscillations. Biophys. J. 67, 447–456.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wang, K., Tu, Y., Rappel, W.-J., and Levine, H. (2005). Excitation-contraction coupling gain and cooperativity of the cardiac ryanodine receptor: a modeling approach. Biophys. J. 89, 3017–3025.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wasserstrom, J. A., Shiferaw, Y., Chen, W., Ramakrishna, S., Patel, H., Kelly, J. E., O’Toole, M. J., Pappas, A., Chirayil, N., Bassi, N., Akintilo, L., Wu, M., Arora, R., and Aistrup, G. L. (2010). 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, 1117–1126.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Weber, C. R., Ginsburg, K. S., Philipson, K. D., Shannon, T. R., and Bers, D. M. (2001). Allosteric regulation of Na/Ca exchange current by cytosolic Ca in intact cardiac myocytes. J. Gen. Physiol. 117, 119–132.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wei, S., Guo, A., Chen, B., Kutschke, W., Xie, Y.-P., Zimmerman, K., Weiss, R. M., Anderson, M. E., Cheng, H., and Song, L.-S. (2010). T-tubule remodeling during transition from hypertrophy to heart failure. Circ. Res. 107, 520–531.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Weiss, J. N., Nivala, M., Garfinkel, A., and Qu, Z. (2011). Alternans and arrhythmias. Circ. Res. 108, 98–112.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Williams, G. S., Chikando, A. C., Tuan, H.-T. M., Sobie, E., Lederer, W., and Jafri, M. (2011). Dynamics of calcium sparks and calcium leak in the heart. Biophys. J. 101, 1287–1296.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Williams, G. S., Huertas, M. A., Sobie, E. A., Jafri, M. S., and Smith, G. D. (2008). Moment closure for local control models of calcium-induced calcium release in cardiac myocytes. Biophys. J. 95, 1689–1703.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text