A multiscale computational model of spatially resolved calcium cycling in cardiac myocytes: from detailed cleft dynamics to the whole cell concentration profiles

Mathematical modeling of excitation-contraction coupling (ECC) in ventricular cardiac myocytes is a multiscale problem, and it is therefore difficult to develop spatially detailed simulation tools. ECC involves gradients on the length scale of 100 nm in dyadic spaces and concentration profiles along the 100 μm of the whole cell, as well as the sub-millisecond time scale of local concentration changes and the change of lumenal Ca²⁺ content within tens of seconds. Our concept for a multiscale mathematical model of Ca²⁺ -induced Ca²⁺ release (CICR) and whole cardiomyocyte electrophysiology incorporates stochastic simulation of individual LC- and RyR-channels, spatially detailed concentration dynamics in dyadic clefts, rabbit membrane potential dynamics, and a system of partial differential equations for myoplasmic and lumenal free Ca²⁺ and Ca²⁺-binding molecules in the bulk of the cell. We developed a novel computational approach to resolve the concentration gradients from dyadic space to cell level by using a quasistatic approximation within the dyad and finite element methods for integrating the partial differential equations. We show whole cell Ca²⁺-concentration profiles using three previously published RyR-channel Markov schemes.

1. Introduction

Cardiomyocyte muscle filament shortening and lengthening is a Ca²⁺ dependent process. The timing of contraction is controlled through electrical excitation via a process known as excitation-contraction-coupling (ECC). ECC is mediated through Ca²⁺, and is facilitated through an amplification process known as Ca²⁺-induced Ca²⁺ release (CICR). In ventricular myocytes, CICR is controlled locally by the colocalization of L-type Ca²⁺-channels (LCCs) in the T-tubule membrane on the one side of a dyadic cleft [aka Ca²⁺ release unit (CRU)] and ryanodine receptor channels (RyRs) in the junctional sarcoplasmic reticulum (jSR) membrane on the other side. Depolarization of the plasma membrane leads to the activation of LCCs, which causes Ca²⁺ entry from the extracellular space into the dyadic space. The influx of Ca²⁺ activates RyRs, which release Ca²⁺ from the sarcoplasmic reticulum (SR) (Fabiato and Fabiato, 1975). Within Ca²⁺ -release-units, Ca²⁺ dynamics are distinct from the myoplasm, with steep [Ca²⁺] gradients and several-fold higher [Ca²⁺] (Stern, 1992). The local release events (sparks) of the ~20,000 individual CRUs sum to produce a macroscopic whole-cell Ca²⁺ transient. The CRUs are coupled through myoplasmic Ca²⁺ diffusion, through SR Ca²⁺ diffusion, and through the spatially homogeneous membrane voltage.

Fundamental properties of ECC/CICR are: (i) that through a mechanism of local control of CICR, a graded LCC current produces a graded RyR Ca²⁺ release (Barcenas-Ruiz and Wier, 1987; Cannell et al., 1987; Stern, 1992); (ii) high whole cell CICR-gain at negative membrane potential V and lower whole cell CICR-gain at positive V (Altamirano and Bers, 2007); (iii) CICR is achieved by local Ca²⁺ release in highly localized cell-subcompartments, with dyadic space [Ca²⁺] rising to ~1000 × myoplasmic [Ca²⁺] (Koh et al., 2006) and strong gradients inside the dyads (Tanskanen et al., 2007; Hake and Lines, 2008; Hake et al., 2012); (iv) stochastic transition rates of LCCs depending on membrane potential (Cleemann et al., 1998) and small fluctuations in the number of LCCs can result in variability of AP duration and in early after depolarization formation (Tanskanen et al., 2005).

These properties illustrate that multiple length scales (tens of nanometers in the dyadic space to 100 μm cell size) and time scales (sub-millisecond for [Ca²⁺] changes in the dyad to tens of seconds for SR dynamics) are involved. To account for that multi-scale character, recently several models with spatially distributed Ca²⁺ release sites have been developed (Restrepo et al., 2008; Restrepo and Karma, 2009; Hatano et al., 2011; Williams et al., 2011; Nivala et al., 2012; Walker et al., 2014). These realistic models of CICR release aim at reproducing the fundamental properties of Ca²⁺ dynamics and at gaining independence from model-simplifying assumptions as far as possible with reasonable effort. Many of these models represent the dyadic space by a single compartment not resolving concentration gradients. Detailed spatially resolved models of the CRU have been developed which represent the steep local [Ca²⁺] gradients in the cleft (Koh et al., 2006; Schendel and Falcke, 2010; Hake et al., 2012; Schendel et al., 2012; Cannell et al., 2013; Stern et al., 2013 and others reviewed in Williams et al., 2010). However, a detailed CRU model has not yet been coupled to whole cell calcium dynamics and to whole cell electrophysiology, with a detailed representation of the spatial distributions of CRUs, challenges we address in this study.

We use mathematical multiscale techniques (Green Function, quasistatic approximation) to simulate a computationally efficient mathematical model of CRUs with spatially resolved [Ca²⁺] and stochastic state dynamics of all individual LC- and RyR-channels. The dynamics of up to 5120 CRUs is then embedded into simulations of the cellular concentration fields for Ca²⁺ and Ca²⁺-binding molecules as well as membrane potential time course.

2. Mathematical Modeling

The mathematical model comprises a system of partial differential equations for the cytosolic and sarcoplasmic concentration dynamics, N_c models for the individual CRUs and a system of ordinary differential equations for the electrophysiology (see Figure 1). We present the individual modules first and then describe their coupling to a whole cell model. All parameter values are listed in Tables 1–5.

FIGURE 1

Figure 1. The modules of the model and their interaction. The mathematical model comprises a set of partial differential equations (PDEs) for the bulk concentrations of cytosolic and sarcoplasmic free Ca²⁺, cytosolic and sarcoplasmic mobile buffers and a cytosolic stationary buffer. The N_c Ca²⁺ release units (CRUs) are simulated all individually and are source terms in the bulk concentration dynamics PDEs. The state dynamics of each of their LC- or RyR-channels is a continuous time Markov chain. The concentration profile in the dyadic space is modeled in spatial detail with a quasistationary approximation, the dynamics of the concentrations of free Ca²⁺ and buffer in the jSR are determined by release into the cleft and refilling from the network SR (nSR). The electrophysiology model has been developed by Mahajan et al. (2008a). The LCC current in the CRUs and the Na⁺/Ca²⁺-exchanger flux couple the membrane potential dynamics directly to the concentration dynamics.

TABLE 1

Table 1. Buffering and diffusion parameters.

TABLE 2

Table 2. Exchanger and uptake parameters.

TABLE 3

Table 3. Physical constants and ionic concentrations.

TABLE 4

Table 4. LCC and RyR parameters for (Stern et al., 1999; Cannell et al., 2013; Walker et al., 2014) model.

TABLE 5

Table 5. Ionic current conductances.

2.1. PDE Model

The dynamics of the cytosolic Ca²⁺ concentration, c, comprises of plasma membrane transport, release and uptake by the SR and binding to buffers. Plasma membrane transport is controlled by ion channels and the Na⁺/Ca²⁺-exchanger. The T-tubule network is an interface to extracellular fluid in the bulk of the cytosol due to which membrane fluxes like the Na⁺/Ca²⁺-exchanger contribute to bulk concentration dynamics (J_NaCa). The Na⁺/Ca²⁺-exchanger flux through the plasma membrane ($J_{NaCa}^{pm}$) enters the partial differential equations as flux boundary condition. The flux J_pump describes the pumping of Ca²⁺ by SERCAs into the SR. A leak flux (often denoted by J_leak) is not included as stochastic RyR openings during diastole account for SR leak. The Ca²⁺-binding molecules (b_j, j = s, m, f) in the cytosol include stationary and mobile buffers and Fluo-4. The total concentration $b_j^{tot}$ is spatially homogeneous for all of the buffers. The partial differential equations for the cytosolic concentration fields and their boundary conditions are:

$\begin{array}{cc}\frac{\partial c}{\partial t}=\nabla ·(D\nabla c)+J_{cru}+J_{NaCa}-J_{pump}-J_{buf}\hfill & \left(1\right)\hfill \end{array}$

$\begin{array}{cc}\begin{array}{l}\frac{\partial b_j}{\partial t}=\nabla ·(D_b^j\nabla b_j)+R_j\left(c,b_j\right),j=s,m,f\hfill \\ {\overrightarrow{n}·D\nabla c|}_{\Gamma }=J_{NaCa}^{pm}\hfill \\ {\overrightarrow{n}·D_b^j\nabla b_j|}_{\Gamma }=0,j=s,m,f,\hfill \end{array}& \left(2\right)\hfill \end{array}$

where D and ${\mathbf{\text{D}}}_b^j$ are diagonal diffusion matrices. We denote the plasma membrane by Γ, and the cell volume by ${\nu }_{cell}\subset {ℝ}^3$. We simulate between 1 and 16 z-discs in a subvolume Ω < ν_cell. The expressions for the fluxes are:

$\begin{array}{c}{J}_{buf}={\displaystyle \sum _{j=1}^{N_B}{R}_j}\left(c,b_j\right),R_j\left(c,b_j\right)=k_j^{+}\left(b_j^{tot}-b_j\right)c-k_j^{-}{b}_j,\\ j=s,m,f& \left(3\right)\hfill \end{array}$

$\begin{array}{cc}{J}_{cru}={\displaystyle \sum _{i=1}^{N_c}\frac{{\nu }_{cell}}{{\nu }_{cyt}}}\frac{\Theta (R_{cru}^i-\left|\overrightarrow{r}-\overrightarrow{r_i}\right|)}{\frac{4}{3}\pi {\left(R_{cru}^i\right)}^3}\left[{\displaystyle \sum _{j=1}^{N_{LCC}^i}{I}_{LCC}^{i,j}}(t)+{\displaystyle \sum _{j=1}^{N_{RyR}^i}{I}_{RyR}^{i,j}}(t)\right]& \left(4\right)\hfill \end{array}$

$\begin{array}{cc}{J}_{pump}=\frac{V_P^{max}{c}^2}{K_P^2+c^2},& \left(5\right)\hfill \end{array}$

where Θ(x) denotes the step function

$\begin{array}{cc}\Theta (x)=\{\begin{array}{cc}0& x<0\\ 1& x\ge 0\end{array}& \left(6\right)\hfill \end{array}$

Ca²⁺ influx and release through LCC ($I_{LCC}^{i,j}(t)$) and RyR ($I_{RyR}^{i,j}(t)$) channels occurs in dyadic clefts. They are represented in the bulk dynamics by spherical source volumes with a radius determined by the numbers of LCCs ($N_{LCC}^i$) and RyRs ($N_{RyR}^i$) in the cleft. Their dependence on time t is caused by their stochastic behavior described in detail in Section 2.3. If an individual channel is open, its current obeys Equations (18 or 19), and is 0 otherwise.

We use the bidomain concept (Keener and Sneyd, 1998) for modeling cytosol and SR processes. It perceives both compartments to fill the same volume continuously with volume ratio ν_cyt/ν_cell and ν_sr/ν_cell, resp. We include one buffer B_sr in the SR lumen. The partial differential equation for the SR Ca²⁺ concentration S and Ca²⁺ buffer B_sr are

$\begin{array}{cc}\begin{array}{l}\frac{\partial S}{\partial t}=\nabla ·(D_S\nabla S)-k_{sr}^{+}\left(B_{sr}^{tot}-B_{sr}\right)S+k_{sr}^{-}{B}_{sr}-J_{jsr}\hfill \\ +\frac{{\nu }_{cyt}}{{\nu }_{sr}}(J_{pump})\hfill \end{array}& \left(7\right)\hfill \end{array}$

$\begin{array}{cc}\begin{array}{l}\frac{\partial B_{sr}}{\partial t}=\nabla ·(D_B\nabla B_{sr})+k_{sr}^{+}\left(B_{sr}^{tot}-B_{sr}\right)S-k_{sr}^{-}{B}_{sr}\hfill \\ {\overrightarrow{n}·D_S\nabla S|}_{\Gamma }={\overrightarrow{n}·D_B\nabla B_{sr}|}_{\Gamma }=0.\hfill \end{array}& \left(8\right)\end{array}$

D_S and D_B are diffusion matrices. Release appears as the flux J_jsr in the S-dynamics, which refills the individual junctional sarcoplasmic reticulae. Their concentration ($c_{jsr}^i$) dynamics are described in the context of the CRU models (Equation 21). The volume of the ith junctional SR is ${\nu }_{jsr}^i$. The refill flux is localized within a spherical volume with radius $R_{jsr}^i$ in the network SR:

$\begin{array}{cc}{J}_{jsr}={\displaystyle \sum _{i=1}^{N_c}\frac{{\nu }_{cell}{\nu }_{jsr}^i}{{\nu }_{sr}}}\frac{\Theta (R_{jsr}^i-|\overrightarrow{r}-{\overrightarrow{r}}_i|)}{\frac{4}{3}\pi {\left(R_{jsr}^i\right)}^3}\frac{S({\overrightarrow{r}}_i)-c_{jsr}^i}{{\tau }_{refill}}.& \left(9\right)\hfill \end{array}$

2.2. Electrophysiology Model

In this subsection we briefly introduce the membrane potential model of the rabbit ventricular myocyte by Mahajan et al. (2008a). We assume that the membrane potential is uniform over the computational domain. The dynamics of the membrane potential V is described by:

$\begin{array}{cc}\frac{dV}{dt}=-\left(I_{ion}+I_{stim}\right).& \left(10\right)\hfill \end{array}$

I_stim is the stimulus current to depolarize the cell. I_ion comprises (Mahajan et al., 2008a):

$\begin{array}{cc}\begin{array}{l}{I}_{ion}=I_{Na}+I_{to,f}+I_{to,s}+I_{Kr}+I_{Ks}+I_{K1}+I_{NaK}\\ +I_{CaL}+I_{NaCa}.\end{array}& \left(11\right)\hfill \end{array}$

I_Na is the fast Na⁺ current, I_K1 is the inward rectifier current, I_{to, f} is the fast component of the rapid outward K⁺ current, I_{to, s} is the slow component of the rapid outward K⁺ current, I_Kr is the rapid component of the delayed rectifier current, I_Ks is the slow component of the delayed rectifier current and I_NaK is Na⁺/K⁺ pump current. These currents and the equations defining them are described in detail in Mahajan et al. (2008a). The Na⁺/Ca²⁺-exchanger current I_NaCa is determined by the spatial integral of the corresponding fluxes in the cytosolic concentration dynamics

$\begin{array}{cc}{I}_{NaCa}=\frac{{\nu }_{cyt}F}{{\nu }_{cell}{C}_m}\left(\frac{{\nu }_{cell}}{\Omega }{\displaystyle {\int }_{\Omega }{J}_{NaCa}}d\Omega +\frac{A_{cell}}{\Gamma }{\displaystyle {\int }_{\Gamma }{J}_{NaCa}^{pm}} d\Gamma \right).& \left(12\right)\hfill \end{array}$

A prefactor f_NaCa for each finite element has been introduced to simulate the distribution of the Na⁺/Ca²⁺-exchanger according to Jayasinghe et al. (2009). The LCC-current I_CaL is the sum of all individual single channel currents (Equation 19) from the CRU models

$\begin{array}{cc}{I}_{CaL}=-2\alpha {\displaystyle \sum _{i=1}^{N_c}{\displaystyle \sum _{j=1}^{N_{LCC}^i}{I}_{LCC}^{i,j}}}.& \left(13\right)\hfill \end{array}$

The factor α = Fν_cell/C_mΩ converts the current in terms of ions/s into pA/pF and at the same time scales the currents up from values obtained with the simulation volume Ω, which is set by the number of simulated z-discs, to values corresponding to a realistic cell volume ν_cell. Similarly, the factor A_cell/Γ in Equation (12) scales the plasma membrane component from simulated (Γ) to whole cell area (A_cell). C_m is the cell membrane capacitance and F the Faraday constant.

2.3. The CRU Model

The CRU model was based on a model previously developed by Schendel et al. (2012). We adapted it in a few key points to properly interact with the PDE model described above. We omit the index numbering the CRUs in this section for simpler notation. The dyadic space is modeled as a cylinder of variable radius and a height of 15 nm containing LCCs and RyRs, which are placed in regular arrays at the bottom and top of the cylinder (see Figure 2). The number of RyRs in each CRU was randomly determined using an exponential distribution with mean 50. For each four RyRs there was one LCC. The diameter of the cylinder was determined for each CRU such that there was a margin of 60 nm between the outermost channel and the boundary of the cleft.

FIGURE 2

Figure 2. Two examples for channel placement in the CRU model. For 30 RyRs and 8 LCCs (left) and 16 RyRs and 5 LCCs (right). The number of RyR channels per CRU obeys an exponential distribution with an average of 50 across all CRUs.

2.3.1. Stochastic Channel Gating

The state dynamics of the LC- and RyR-channels are simulated with Markov models. For the RyRs we explored three models, the first is a model originally developed by Stern et al. (1999) in it's modified version by Schendel et al. (2012). It is a four state model [see Figure 3A(i)], including one open, one resting and two inactivated states. It features RyR inhibition in case of high cytosolic or low jSR Ca²⁺concentrations. For an in-depth description of this model see (Schendel et al., 2012). The second model is the two-state model [see Figure 3A(ii)] developed by Cannell et al. (2013), in which termination of CICR is mediated through the steep Ca²⁺dependence of the RyR closed time. With a small decline in jSR [Ca²⁺] the Ca²⁺flux via open RyRs declines, causing a decline in local dyadic [Ca²⁺], which in turn causes a decrease in the open probability of neighboring RyRs, a process known as induction decay. This RyR model does not rely on experimentally un-substantiated biophysical mechanisms for CICR termination, such as dyadic/cytoplasmic RyR Ca²⁺-dependent inactivation, or RyR-lumenal Ca²⁺-dependent inactivation. The third model is the two-state model [see Figure 3A(iii)] developed by Walker et al. (2014) (adapted from Williams et al., 2011) which incorporates modulation of the RyR-opening rate by junctional RyR-lumenal [Ca²⁺] ($c_{jsr}^i$). This model has a fixed closing rate, and, in accordance with the experimental data of Cannell et al. (2013), there is only weak regulation of the RyR opening-rate when ($c_{jsr}^i$) is < 1 mM. It is of note that the opening rate of this last RyR-Model is theoretically unbound, however in our simulations opening rates larger than 0.7 ms⁻¹ were not encountered.

FIGURE 3

Figure 3. State scheme for the RyR and LCC Markov Models. A(i) 4-state RyR-channel model of Stern et al. (1999), O denotes the open, I the inactivated, and R the resting state. The activation rate k_o is a fourth order Hill function of dyadic Ca²⁺, the inactivation rate k_i is a first order Hill function of Ca²⁺. The Ca²⁺ dependent rates are influenced by the Ca²⁺ concentration in the jSR. The model was originally developed by Stern et al. (1999) and was modified by Schendel et al. (2012). A(ii) 2-state RyR-channel model of Cannell et al. (2013) where O denotes the open and C the closed state and k_open and k_close are polynomial functions of [Ca²⁺]. A(iii) 2-state RyR-channel model of Walker et al. (2014) where O denotes the open and C the closed state, k_open is a polynomial function of both [Ca²⁺] and [Ca²⁺] _jsr (i.e., k_open is influenced by the Ca²⁺concentration in the jSR) and k_close is a constant. (B) 7-state scheme for the LCC according to Mahajan et al. (2008a) with the states O open, I1_Ca, I2_Ca Ca²⁺ dependent inactivated states, C1, C2 closed states, I1_Ba, I2_Ba Ca²⁺ independent inactivated states (for details see Mahajan et al., 2008a).

For the LCCs we used the 7-state model developed by Mahajan et al. (2008a) (see Figure 3B) which exhibits both Ca²⁺ dependent and Ca²⁺ independent inactivation.

Each channel's state was chosen randomly according to the steady state distribution for the initial values. The transition times between different states were determined using the Gillespie Algorithm for time dependent transition rates. That algorithm determines the time of an event as the time of the crossing of a random threshold by an integrated propensity for each individual Markov chain. With the given rate of changes in membrane potential and Ca²⁺ concentration, propensities are strongly time dependent. Hence, before we finally accept a calculated update of V and c, we verify that we have taken into account all events triggered by propensity crossings by the update during the iteration step. The precise definition of the transition rates are given in Schendel et al. (2012), Cannell et al. (2013), and Walker et al. (2014) for the three RyR-Models and in Mahajan et al. (2008a) and the accompanying CellML data (Mahajan et al., 2008b) for the LCC Model (there are slight discrepancies in the equations between the CellML data and the original paper; the equations we use are from the CellML data). In some cases we had to adjust the constants involved in the equations, and these are given in Table 4.

The frequency of transitions becomes very large with a large number of CRUs and with many channels in each CRU. While this in itself does not have a significant impact on the speed of the stochastic algorithm, it can force very small timesteps on the PDE model. In order to alleviate this problem and allow for longer timesteps we take advantage of the properties of the stochastic process. Those channel state transitions which do not represent either openings or closings are (within a given iteration step) independent stochastic events. We can therefore execute an arbitrary number of such transitions in a single iteration. This approach scales almost independently of the total number of states of the channel model. It is only dependent on the frequency of transitions from open to closed state or vice versa.

2.3.2. Ca²⁺profile in the CRU

The Ca²⁺concentration in the cleft (c_di) is modeled by the following equation in cylindrical coordinates:

$\begin{array}{cc}\begin{array}{l}\frac{\partial c_{di}}{\partial t}={\displaystyle \sum _{i=1}^{N_{LCC}}{J}_{\text{LCC}}^i}+{\displaystyle \sum _{j=1}^{N_{RyR}}{J}_{\text{RyR}}^j}+D_c{\Delta }_{r,\phi }{c}_{di}(\overrightarrow{r})-\frac{\partial }{\partial z}{J}_z.\hfill \end{array}& \left(14\right)\hfill \end{array}$

Because the time to reach the stationary concentration profile upon opening or closing of a channel is about one order of magnitude shorter than the timescale of channel state transitions, we assume steady state for c_di. We make a separation ansatz (yielding only a negligible error as discussed in Schendel and Falcke, 2010):

$c_{di}(r,\phi ,z)=f(r,\phi )·Z(z).$

The large aspect ratio radius/height of the dyadic space renders gradients in the z-direction of the cylinder due to channel fluxes negligible (Schendel et al., 2012). However, electro-diffusion due to differential charge densities on T-tubule and jSR membrane causes a substantial gradient. Soeller and Cannell (1997) determined the charge density profile shaping Z(z) (Soeller and Cannell, 1997). With their result we obtain

$Z(z)=e^{-2{\varphi }_0{e}^{-\kappa z}}.$

The values of ϕ₀ and κ are listed in Table 3. To find f we consider (Equation 14) for a given configuration of open channels. The flux through the channels is constant due to the steady state assumption. We describe the channel fluxes as Dirac-δ-functions in space with the channel currents as pre-factors, integrate the equation over z and obtain

$\begin{array}{cc}\begin{array}{l}{\Delta }_{r,\phi }f(r,\phi )=-\frac{1}{D_c{h}^{\ast }}[{\displaystyle \sum _{i=1}^{N_{LCC}}{I}_{LCC}^i}\delta (r-r_i,\phi -{\phi }_i)\\ +{\displaystyle \sum _{j=1}^{N_{RyR}}{I}_{RyR}^j}\delta (r-r_j,\phi -{\phi }_j)]\text{},\end{array}& \left(15\right)\hfill \end{array}$

where $h^{*}={\int }_0^{h}Z(z)dz$¹. Solving this using the Green function G((r, φ), (r_i, φ_i)) we find

$\begin{array}{cc}\begin{array}{l}f(r,\phi )-c_{bulk}=\frac{1}{h^{\ast }{D}_c}{\displaystyle \sum _{i=1}^{N_{LCC}}{I}_{LCC}^i}G((r,\phi ),(r_i,{\phi }_i))\\ +{\displaystyle \sum _{j=1}^{N_{RyR}}{I}_{RyR}^j}G((r,\phi ),(r_j,{\phi }_j))\\ G(r,r_i)=\frac{1}{2\pi }\log \left(\frac{\parallel \widehat{r_i}-r\parallel }{q\parallel r_i-r\parallel }\right)\text{ for }{r}_i\ne 0\\ G(r,0)=\frac{1}{2\pi }\log \left(\frac{R}{\parallel r\parallel }\right)\end{array}& \left(16\right)\hfill \end{array}$

where R is the radius of the dyadic space, $q=\frac{R}{\parallel r_i\parallel }$, and $\widehat{{\mathbf{\text{r}}}_i}=q^2{\mathbf{\text{r}}}_i$. c_bulk is the Ca²⁺concentration at the boundary of the cleft (cylinder barrel) as computed by the PDE model.

With

$\eta (\overrightarrow{r_i},\overrightarrow{r})=\frac{Z(z)}{h^{\ast }{D}_c}G((r,\phi ),(r_i,{\phi }_i))$

c_di can be written as

$\begin{array}{cc}{c}_{di}(\overrightarrow{r})=Z(z_i)c_{bulk}+{\displaystyle \sum _{i=1}^{N_{LCC}}{I}_{LCC}^i}\eta (\overrightarrow{r_i},\overrightarrow{r})+{\displaystyle \sum _{i=1}^{N_{RyR}}{I}_{RyR}^i}\eta (\overrightarrow{r_i},\overrightarrow{r}).& \left(17\right)\hfill \end{array}$

The currents $I_{LCC}^i$ and $I_{RyR}^i$ depend themselves on the local Ca²⁺ concentration at the channel mouth:

$\begin{array}{cc}{I}_{RyR}^i=g(c_{jsr}-c_{di}({\overrightarrow{r}}_i))\text{ for }i\in N_{\text{RyR}}& \left(18\right)\hfill \end{array}$

$\begin{array}{cc}{I}_{LCC}^i=J_L\delta V\frac{c_{\text{ext}}{e}^{-\delta V}-c_{di}({\overrightarrow{r}}_i)}{1-e^{-\delta V}}\text{ for }i\in N_{LCC}& \left(19\right)\hfill \end{array}$

with δ = 2F/(RT). Here, ${\overrightarrow{r}}_i$ denotes the location of the mouth of the ith channel. Inserting Equations (18, 19) into Equation (17) and evaluating it at each channel mouth defines a system of linear equations for the concentration values setting the currents for given values of c_jsr and c_bulk. The coefficients of that system can be calculated in advance of a simulation from the cleft geometry, which renders the simulation very efficient.

2.3.3. Buffers in the Dyadic Space

In order to employ the two state RyR-models by Cannell et al. (2013) and Walker et al. (2014) (see the paragraph on Stochastic Channel Gating below) we found it necessary to introduce a form of Ca²⁺buffering in the dyadic space. This is to be expected since in the original papers the models were fitted to data from experiments conducted in the presence of buffers such as Calmodulin, Fluo-4, and ATP. Cannell et al. (2013) used also in their simulations high buffer concentrations within the dyadic space, and in order to adapt our model to this degree of buffering we introduce a linear buffer factor β. The coupling factor η now takes the form:

$\eta (\overrightarrow{r_i},\overrightarrow{r})=\frac{1}{\beta }\frac{Z(z)}{h^{\ast }{D}_c}G((r,\phi ),(r_i,{\phi }_i))$

We adjusted β so that the Ca²⁺ concentration at the channel mouths and the average dyadic Ca²⁺ concentrations matched the values given in the paper by Cannell et al. (2013). An approximation of [Ca²⁺]_i that would be measured by a single-wavelength Fluo-4 experimental recording (denoted ${[C{a}^{2+}]}_i^{exp}$) was calculated from the Ca²⁺-bound Fluo-4. The in vitro calibration approach described in Takahashi et al. (1999) was used:

$\begin{array}{cc}{[C{a}^{2+}]}_i^{exp}=K_d\times (F-F_{min})/(F_{max}-F),& \left(20\right)\hfill \end{array}$

where K_d is the dissociation constant of Fluo-4, F is the experimentally measured fluorescence intensity (here the [Ca²⁺-bound Fluo-4]), F_max is the measured fluorescence intensity in Ca²⁺-saturated dye (here this is set as the maximum [Ca²⁺-bound Fluo-4], i.e., $b_{Fluo-4}^{tot}$), and F_min is the measured fluorescence intensity in the absence of Ca²⁺ (here set to zero) (see Table 1).

2.3.4. Modeling jSR Dynamics

Each dyadic space is paired with its own junctional SR (jSR) compartment. We assume spatially uniform Ca²⁺concentration in the jSR. Ca²⁺dynamics in the jSR depends on the release flux through the RyRs and a refill flux from the network SR (nSR), which we assume to depend simply on the concentration difference S − c_jsr. The buffering by Calsequestrin is modeled using the fast buffer approximation.

$\begin{array}{cc}\begin{array}{l}\frac{d{c}_{jsr}}{dt}={\beta }_{jsr}\left(\frac{S(\overrightarrow{r})-c_{jsr}}{{\tau }_{refill}}-\frac{1}{{\nu }_{jsr}}{\displaystyle \sum _{j=1}^{N_{RyR}}{I}_{RyR}^j}\right)\hfill \\ {\beta }_{jsr}={\left(1+\frac{n{K}_{csqn}{B}_{csqn}}{{(K_{csqn}+c_{jsr})}^2}\right)}^{-1}\hfill \end{array}& \left(21\right)\hfill \end{array}$

$S(\overrightarrow{r})$ denotes the free nSR Ca²⁺concentration at the location $\overrightarrow{r}$ of the spherical volume representing refilling in the S-dynamics (Equation 7). This equation is solved numerically by linearization with very small time steps (< 10⁻⁴ ms).

2.4. Module Interactions

The modules interact on a time scale longer than a single iteration step by the dependencies of the dynamics on Ca²⁺, V, and other state variables. Instantaneous interactions are important for the algorithmic realization of the time integration. We use for the completely coupled system an Euler forward method. The instantaneous interactions are between the PDE system and the CRU models mediated by Ca²⁺ currents and concentrations, between the CRU models and the electrophysiological model by the LCC current, and between the PDE system and the electrophysiological model via the bulk and plasma membrane components of the Na+/Ca²⁺-exchanger flux (see also Figure 1).

2.4.1. Bulk Concentrations—Electrophysiology

We use the Ca²⁺ concentration spatially resolved at the plasma membrane for calculating local values of $J_{NaCa}^{pm}$ and analogously in the bulk for J_NaCa, and then average (Equation 12) to obtain the current entering the membrane potential dynamics. Vice versa, $J_{NaCa}^{pm}$ serves as boundary condition and J_NaCa as bulk source term for the PDEs.

2.4.2. Bulk Concentrations—CRUs

The value of the Ca²⁺ concentration c determined by the PDEs is averaged along the rim of a cleft volume and then serves as boundary condition c_bulk in Equation (17) to determine the RyR and LCC currents ($I_{LCC}^{i,j}$ and $I_{RyR}^{i,j}$) in each CRU. Vice versa, these currents determine the source terms for the bulk concentration dynamics in Equation (1).

2.4.3. Electrophysiology—CRUs

The membrane potential affects the LCC currents (Equation 19) and vice versa the sum of all individual LCC currents the membrane potential (Equation 13).

2.5. Numerical Approach

We used a piecewise bi-linear finite element method for the solution of the spatially three-dimensional reaction-diffusion model including the complex distribution of CRUs at multiple z-discs. The first challenge was the fine scale resolution of the computational grid to resolve the strong concentration gradients at the boundary of the CRUs. We take the equidistant tetrahedral elements with the size of 0.05 μM in our computations. The next challenge is to deal with the adaptive time stepping schemes for solving the reaction-diffusion systems. Due to the fast transitions of the channel opening/closings in a CRU, the time scales vary from tens of microseconds to milliseconds. To resolve such rapid changes adaptive and higher order time steppings are inevitable to treat the very smooth diffusion effects as efficiently as possible. To this end, we use higher order linearly implicit Runge-Kutta methods for time discretization of reaction-diffusion systems, see (Lang, 2001; Chamakuri et al., 2012). These belong to a large class of methods which try to avoid the nonlinear system and replace it by a sequence of linear ones. Also, this allows the use of adaptive timescales in the simulations. Specifically, we employed a second order Rosenbrock method called ROWDA (Lang, 2001), and we avoid the time discretization implications of this problem as described in Chamakuri and Rüdiger (2012, Section 4.2).

Our parallel implementation of the discretization routines are based on the public domain package DUNE (Bastian et al., 2008), especially the dune-pdelab discretization module. The parallel linear solvers depend on the dune-istl module. Based on this, we developed a finite element simulator to solve the whole cell calcium cycling model. The internal parallel Cartesian (called Yasp) grid in DUNE is used for the parallel grid constructions. In our domain decomposition approach, the original domain is partitioned into subdomains and each subdomain is assigned to a single processor. In our computations, we used a non-overlapping domain decomposition approach to solve the discretized PDE model. In this regard, we parallelized our code by using minimum global communications for solving the stochastic part of the problem. The time discretization results in a system of linear equations which can be solved using efficient iterative solvers. The BiCGSTAB (van der Vorst, 1994) with Jacobi preconditioner is used as the linear solver and the relative tolerance of 10⁻⁶ is used as the stopping criteria for the linear solver at each step of the ODE time integrator.

Here we propose a novel technique to determine the new timestep during the stochastic opening of many CRUs. Due to the presence of the large numbers of CRUs, the stochastic algorithm that governs the timestep for the next channel transition plays an important role for the computations. As mentioned before, channel state transitions which do not represent openings or closings are, within a given iteration step, independent stochastic events. We can therefore execute several of them within a single iteration. Additionally, typical time steps during an AP are in the range of 0.01 ms, i.e., they are shorter than the diffusion time between neighboring CRUs. Consequently, conductance changing events in different CRUs are statistically independent on the time scale of a single iteration and we can allow for several of them in (distinct) CRUs within one time step.

We introduce two time steps: first the deterministic timestep τ_det (which is allowed by the numerical integration of the PDEs) and second the stochastic timestep τ_stoc. We propose the following algorithm.

The bulk calcium cycling PDE model and the electrophysiology model are integrated from t to t + τ_det, where τ_det is the accepted deterministic timestep of the PDE solver. Then, the stochastic channel transitions are predicted from t to t + τ_det. Suppose there were N_s conductance changing stochastic events at times t + τ_i where i = 1…N_s, τ_i ≤ τ_det and 0 ≤ N_s ≤ N_c. Here, the time of the stochastic event is τ_i for the i^th CRU. In case that there is no stochastic event for a CRU, τ_i is set to τ_i = τ_det. The stochastic timestep τ_stoc is determined from the τ_i as the time by which a maximum number of acceptable transitions is reached. The maximum number has been determined empirically to be sufficiently small with 0.1 N_s to cause no essential difference to simulations with τ_stoc sufficiently small to guarantee N_s = 1. Now all the occurring events in the CRUs up to t + τ_stoc are set to take place at time t + τ_stoc. By doing so, we avoid time steps which are too small for acceptable simulation time. A schematic illustrating a single iteration and how the time steps are determined can be found in Figure 4.

FIGURE 4

Figure 4. A schematic illustrating the basic steps for a single iteration. Orange marks all values produced and all work done by the CRU-Model. The work and results of the PDE-Model are marked in red, everything involved with the Electrophysiology ODE-Model (Electrophys. in the diagram) is marked in brown. A superscript t marks the value of a given quantity at time t, while a superscript t + τ denotes the predicted value of that quantity at time t + τ, e.g., the prediction step for the PDE-Model uses the flux values from the CRU and ODE-Model at time t. The CRU Transition Prediction step uses both the current values at time t and the predicted values at time t + τ from the PDE- and ODE-model.

3. Results

Our main result is the fully coupled simulation tool. Our motivation was to be able to take stochastic channel state dynamics for each channel in each CRU and the concentration profile within CRUs into account while executing the simulation of (partial) differential equations for other state variables. Figure 5 shows such a concentration profile inside the dyadic space calculated with Equation (17) for a typical RyR current. Gradients are substantial such that distant RyRs and LCCs experience much smaller concentrations and Ca²⁺-dependent transition rates than channels close to an open one.

FIGURE 5

Figure 5. Simulation of a typical Ca²⁺ gradient in a CRU. The stationary dyadic [Ca²⁺] profile of a single open RyR is shown.

We simulated 16 z-discs, each with 320 CRUs, and the Ca²⁺ dynamics for this sub-cellular region (~30% of the cardiomyocyte) were coupled to the whole-cell electrophysiology ODE model. It took 64.2 h to solve a single action potential on 848 Intel Xeon E5-2650 v2 2.60 GHz CPUs (central processing units). The Ca²⁺ concentration profile 70.0 ms after stimulus at a single z-disc is shown in Figure 6. The local Ca²⁺ dynamics and corresponding whole cell electrophysiology are shown in Figure 7A and in Supplementary Movie 1. The SR free [Ca²⁺] is visualized in Figure 7B and in Supplementary Movie 2. Currents are shown in Figure 8.

FIGURE 6

Figure 6. Free [Ca²⁺]_i and SR free [Ca²⁺] at 70.0 ms after activation, using the Walker et al. (2014)-RyR-model. (A) Myoplasm [Ca²⁺]_i. (B) nSR [Ca²⁺] for the 8th z-disc numbered from the bottom in Figure 7. The concentration is color-coded according to the color scale shown. There are 320 CRU per z-disc, with an average of 50 RyR and 12.5 LCC per CRU. See Supplementary Movies 1, 2 for the evolution of [Ca²⁺]_i and nSR [Ca²⁺] through an AP.

FIGURE 7

Figure 7. Three dimensional visualization of spatially resolved [Ca²⁺] at 70.0 ms after activation, using the (Walker et al., 2014)-RyR-model. (A) Isosurfaces show [Ca²⁺]_i in green for [Ca²⁺]_i = 0.6 μM and red for [Ca²⁺]_i = 2.4 μM. (B) The yellow isosurface shows free SR-[Ca²⁺] = 430 μM. There are 320 CRU per z-disc, amounting to 5160 CRUs in total, with an average of 50 RyR and 12.5 LCC per CRU. See Supplementary Movies 1, 2 for the evolution of [Ca²⁺]_i and nSR [Ca²⁺] through an AP.

FIGURE 8

Figure 8. Membrane potential, ionic currents, and concentrations after a stimulus of 40 mV for 1 ms with the (Walker et al., 2014)-RyR model. (A) Membrane potential. (B) Currents I_K1, I_Kr, I_Ks, I_NaK, and I_{to, s}. (C) LCC current (I_CaL), Na⁺/Ca²⁺-exchanger current (I_NaCa), I_{to, s}, and I_Na (truncated). (D) Average cytosolic [Ca²⁺]_i and ${[C{a}^{2+}]}_i^{exp}$ as defined in Equation (20). (E) Ca²⁺ -fluxes: J_rel, J_Ca, J_NaCa, J_up. (F) Buffer-bound [Ca²⁺] ([b_m], [b_s], [B_sr]), and nSR lumenal [Ca²⁺] ([Ca_sr]) are plotted. [b_m] and [b_s] are conventional concentrations in units of μM whereas [Ca_sr] and [B_sr] are expressed for simplicity in units of μMν_sr/ν_cyt (μmol/l cytosol). Plots are of the first AP from a simulation with 16 z-discs.

The simulation shows that we can reproduce the whole sequence of Ca²⁺ transients from the initial concentration profile in the dyadic space, via sparks to the whole cell transient. The cytoplasmic [Ca²⁺] is rather heterogeneous, due to randomness of release events as well as variations in CRU size. The concentration of 0.6 μM is reached in almost the whole volume, but 2.4 μM only in the proximity of the the z-discs. The currents reproduce basic features of the Mahajan-model (Mahajan et al., 2008a), and likewise reproduce experimental plots of rabbit Ca²⁺ currents well (Weber et al., 2002), but show some discrepancies mostly related to Na⁺/Ca²⁺-exchanger current.

We explored the use of three RyR models which all produce realistic action potentials (Figure 9). When using the (Stern et al., 1999) model there were slow kinetics of Ca²⁺ release and a low RyR maximum open percentage (2.7%). Gain (the ratio of J_rel to J_Ca) was ~4 during the first 50 ms of the AP and ~3 for the remainder of the AP (AP duration = 265 ms, basic cycle length BCL = 350 ms). With the (Cannell et al., 2013) model many RyRs open and close (spark-like) within 20 ms of the start of the AP, and gain was high in this early phase of the AP (~10) and moderate in the later AP (~6) and there was a corresponding fast rise in the [Ca²⁺]i with an early peak at ~10 ms, followed by a reduction of [Ca²⁺]_i and a second [Ca²⁺]_i rise and fall through the AP. With the (Walker et al., 2014) model gain was ~6 in the early AP, then ~2 in the remainder of the AP, and there was a relatively high RyR maximum open percentage (32%). The (Cannell et al., 2013; Walker et al., 2014) 2-state RyR-schemes have markedly faster (and more physiological) kinetics of Ca²⁺ release than the (Stern et al., 1999) 4-state scheme. With the (Walker et al., 2014)-RyR scheme our model displays expected restitution properties with shorter APD and [Ca²⁺]_i transient with shorter BCL (Figure 10).

FIGURE 9

Figure 9. Comparison of 4-state RyR model (Stern et al., 1999, top panels), 2-state RyR model (Cannell et al., 2013, middle panels), and 2-state RyR model (Walker et al., 2014, bottom panels). (A) Stern et al. (1999) membrane potential. (B) Stern et al. (1999) average cytosolic [Ca²⁺]_i and ${[C{a}^{2+}]}_i^{exp}$ (as defined in Equation 20). (C) Cannell et al. (2013) membrane potential. (D) Cannell et al. (2013) average cytosolic [Ca²⁺]_i and ${[C{a}^{2+}]}_i^{exp}$. (E) Walker et al. (2014) membrane potential. (F) Walker et al. (2014) average cytosolic [Ca²⁺]_i and ${[C{a}^{2+}]}_i^{exp}$. Plots are of the first AP from a simulation of a single z-disc.

FIGURE 10

Figure 10. Restitution properties with the (Walker et al., 2014)-RyR model. (A) Membrane potential. (B) Average cytosolic [Ca²⁺]_i]. (C) Average cytosolic ${[C{a}^{2+}]}_i^{exp}$ as defined in Equation (20). Solid line: basic cycle length (BCL) = 350 ms; dashed line: BCL = 300 ms; dotted line: BCL = 250 ms. At cycle length 350 ms the DI was 177 ms. At cycle length 300 ms the DI was 132 ms. At cycle length 225 ms the DI was 87 ms. (D) Constant BCL Restitution Curve (Otani and Gilmour, 1997). For clarity, the AP with the shortest DI in (D) is not shown in (A–C). These simulations have been carried out with a single z-disc.

4. Discussion

This proof of concept study demonstrates a new multi-scale model of CICR linking three spatial scales: (1) detailed molecular stochastic modeling of the CRU at a continuous spatial scale; (2) a whole-cell ODE electrophysiology model, which describes the potassium channels, the voltage gated sodium channels, the sodium-potassium pump, and integrates all membrane fluxes to derive the total membrane current and voltage; (3) a PDE FEM calcium diffusion model representing myoplasmic and nSR Ca²⁺ diffusion between CRUs and between z-discs. One rationale for this approach is that it allows the removal of the artificial compartment for the sub-membrane space, and hence provides a more quantitative modeling of diffusion processes and gradients.

Taking gradients in the dyadic cleft into account in a simulation with thousands of CRUs would be impossible without the use of the stationary Green function. An estimate of the advantage of the quasi-static approximation with respect to computational speed can be obtained from estimating the number of operations per cleft. Our approach requires to solve a linear system of equations whenever channels open or close or the boundary condition or the jSR concentration has changed significantly. That means ${(N_{RyR}+N_{LCC})}^2$ operations each. Simulating the concentrations inside the cleft would require about 1000 grid points and a time step in the microsecond range (see Thul and Falcke, 2004). Consequently, it would be at least a few hundred times slower than the use of the stationary Green function. Applying the dynamic Green function would entail similar computational efforts, since it converges extremely slowly (Bentele and Falcke, 2007).

To our knowledge, our model is the only whole cell AP model which includes CRUs clustered around the z-discs (a proxy for explicit T-tubules) and spatial resolution of intra-CRU [Ca²⁺] with realistic steep [Ca²⁺] gradients. Other approaches include the multiscale model of Restrepo and Karma (2009) which simulates 20,000 CRUs with Markov models for individual LC- and RyR-channels in a three-dimensional grid interacting with a membrane voltage ODE-model (Mahajan et al., 2008a). This was developed further by Nivala et al. (2012) and was modified to a bidomain approach with finite elements belonging to cytoplasm or SR. A related model of Williams et al. (2011) incorporates energetic coupling of the RyRs, and non-junctional RyRs, but is used only for the analysis of the diastolic RyR dynamics without an ODE electrophysiology model. Walker et al. (2014), based on Williams et al. (2011), developed a finite volume model of a single CRU which includes modeling of intra-CRU calcium concentration differences, diffusion of Ca²⁺ and mobile buffers, T-tubules, Markov chain models for the channels and excitation-contraction coupling gain. Cannell et al. (2013) investigated spark dynamics and termination using a two-state RyR model with modeling of intra-CRU calcium concentration differences, but until our current study the Cannell et al. (2013) model had not been coupled to a whole cell geometry or to an action potential model. Torres et al. (2014) modeled a modified local control mechanism to simulate confocal image data and to explore the relationship between Ca²⁺ transients and the activation of non-junctional RyR clusters of myocytes with sparse T-systems. A finite-element model has been implemented by Hatano et al. (2011) which incorporates subcellular structures including the T-tubule system, SR, myofibrils and mitchondria, where electrophysiology and local calcium dynamics are coupled to a local myofibril contraction model to simulate local and global contraction. They model the spatial detail of CICR release and local contraction for three myofibrils of one sarcomere length but do not model the CRU [Ca²⁺] in spatially resolved detail.

We explored the use of three RyR models: the 4-state model of Stern et al. (1999), the 2-State induction decay model of Cannell et al. (2013), and the 2-State model of Walker et al. (2014) [similar to the Cannell et al., 2013-model except with the addition of modulation of the RyR open probability by junctional RyR-lumenal [Ca²⁺] ($c_{jsr}^i$)]. With the (Stern et al., 1999) RyR-model peak [Ca²⁺]_i equals 0.8 μM with time to peak of 200 ms. With the (Cannell et al., 2013) model, dual peaks in [Ca²⁺]_i were observed with an early peak of [Ca²⁺]_i = 0.6 μM at time to peak ~10 ms, with the principal peak [Ca²⁺]_i transient of ~0.8 μM and time to peak 100 ms. With the (Walker et al., 2014) peak [Ca²⁺]_i was ~1.0 μM and time to peak 130 ms. The time to peak [Ca²⁺]_i for both the 2-state models is comparable to that reported for the rabbit by Weber et al. (2002). The early [Ca²⁺]_i local maximum with the (Cannell et al., 2013) model may be explained from the origins of the (Cannell et al., 2013) model, which is from fits to experimental data from Ca²⁺ sparks and blinks, and has never been previously used in a model of prolonged CICR resulting from an AP. Indeed the authors of this model state that the model equations are not designed to capture behavior at resting [Ca²⁺]_i with fully loaded SR. The RyR-dynamics of this model with its current parameters are therefore reproducing spark like behavior in an AP. We had to introduce strong buffering inside the dyadic space and jSR-depletion in order to reach closing rates sufficiently fast for termination of release at the end of an action potential. However, this proof of concept investigation does not provide evidence that induction decay alone is an insufficient mechanism to cause the termination of CICR in an AP. Rather, the model provides a platform for investigating AP CICR. Careful experimentally led tuning of the 2-state model parameters alongside Ca²⁺ buffering parameters and the jSR refill flux will allow assessment of the feasibility of induction decay as the sole CICR termination mechanism in an AP. A second possible explanation for the early local maximum in [Ca²⁺]_i with the (Cannell et al., 2013)-model is that in our model we do not explicitly represent the [Ca²⁺]_i uptake and release by the mitochondria. Although mitochondrial [Ca²⁺]_i uptake and release is controversial on such fast timescales (Boyman et al., 2014), there is evidence that fast early mitochondrial [Ca²⁺]_i uptake can act like a fast stationary buffer, and that inhibition of fast uptake through a specific inhibitor of the mitochondrial Ca²⁺ uniporter (MCU) can result in an early local maximum of the [Ca²⁺]_i transient (Maack et al., 2006). Furthermore, the data in Figure 9 show that the [Ca²⁺]_i recorded in experiments via a fluorescent Ca²⁺ probe such as Fluo-4 [approximated as (${[C{a}^{2+}]}_i^{exp}$), as defined in Equation 20] would smooth-out any early [Ca²⁺]_i peak, if this were a genuine experimental feature of cardiomyocyte CICR. In AP simulations with the 2-state RyR models, some local high calcium transients at CRU-sites remained into the diastolic period, associated with CRUs where some RyRs remained open. A similar phenomenon has been described in the setting of spontaneous Ca²⁺-sparks by Stern et al. (2013), who used an RyR-scheme similar to Walker et al. (2014) [a 2-state model, with modulation of the open probability by junctional RyR-lumenal [Ca²⁺] ($c_{jsr}^i$)].

The discrepancies between our simulations and the ODE-Mahajan model with respect to Na⁺/Ca²⁺-exchanger current illustrate the value of spatially resolved modeling in exploring detailed properties of Ca²⁺ dynamics. Using the same Na⁺/Ca²⁺-exchanger model as Mahajan et al. (2008a) we were not able to reproduce the same time course of the total Na⁺/Ca²⁺-exchanger current. In both models the current has the same general profile, with a short period in reverse mode (positive current) in the first phase of an action potential, followed by forward mode (negative current) during the V-plateau and decline. However, in our approach the I_NaCa current turns negative much later during the plateau than in the ODE model. The Na⁺/Ca²⁺-exchanger depends on the [Ca²⁺] in sub-membrane space in the ODE model and on the cytoplasmic concentration c in our simulations, including of course the regions with high c-values at and close to CRUs. The main reason for the differences is that the sub-membrane [Ca²⁺] in the ODE model exhibits a sharp rise and decline while the average myoplasmic concentration shows slower dynamics. We have experimented with localizing the Na⁺/Ca²⁺-exchanger molecule density in proximity to the CRUs to reflect the observed higher abundance of Na⁺/Ca²⁺-exchanger near the CRU (Scriven et al., 2000, 2002; Jayasinghe et al., 2009; Scriven and Moore, 2013). However, that has not fully resolved the discrepancies with respect to the ODE results yet. We expect to gain detailed insights into the factors shaping the Na⁺/Ca²⁺-exchanger current through investigation of these problems.

5. Conclusions and Limitations

We showed on a level of proof-of-concept that multiscale modeling of cardiomyocyte ECC from sub-dyadic scales to many z-discs using full partial differential equations is possible. We demonstrate that the model produces realistic physiology on these scales and has the potential to provide new insight into subcellular mechanisms and structures. Maybe the most severe limitation of this modeling approach is the requirement for high performance computing (HPC) to run it. We expect some improvement of simulation efficiency by more specifically tailored numerical methods, however the requirement for HPC will remain. The use of the Green function inside the dyadic cleft requires linearity of the reaction diffusion equations there. That excludes non-linear buffering terms and allows for linear buffering only. Some of the limitations of our model arise from its early state of development and will be removed with inclusion of more detail like e.g., anisotropic diffusion, mitochondria, more detailed buffering, and non-junctional RyR. We used modular programming as far as possible to ease the exploration of a variety of ion channel models and other species specific membrane potential dynamics.

Author Contributions

NC implemented the PDE model and coupled this with the CRU model and cell electrophysiology ODE model. JV co-implemented the cell electrophysiology ODE and PDE model. WN and SG implemented the RyR Markov-models. WN implemented the CRU model and co-implemented the PDE model. NC, WN, JV, SG carried out simulations, processed data from simulations and produced results figures. SG, NC, WN, JV, MF drafted the manuscript. MF conceived the CRU model, the PDE model and the coupling of the model scales. All authors approved and assisted in redrafting the final manuscript.

Funding

This study has been supported by DFG grant FA350/9-1, DFG GRK 1772, BMBF grant eMed:SMART and a collaborative grant from the German Centre for Cardiovascular Research to MF and S. Luther, MPI DS Göttingen.

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

We acknowledge the work of Thomas Schendel who developed the CRU model which was adapted to the CRU model used in this study (Schendel et al., 2012). NC gratefully acknowledges support from the Austrian Academy of Sciences (OÄW).

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys.2015.00255

Supplementary Movie 1. Simulation of the first AP with 16 z-discs using the (Walker et al., 2014)-RyR model. (A) Shows the three-dimensional 16 z-disc myoplasmic [Ca²⁺] by means of 0.6 μM and 2.4 μM free myoplasmic [Ca²⁺]_i isosurfaces; (B) shows the membrane-potential and the mean myoplasmic [Ca²⁺]_i transient; (C) is a section through the PDE solution at the 8th z-disc from the top with the image colored by free myoplasmic [Ca²⁺].

Supplementary Movie 2. Simulation of the first AP with 16 z-discs using the (Walker et al., 2014)-RyR model. (A) Shows the three-dimensional 16 z-disc [Ca_sr] by means of a 430 μM isosurface; (B) shows the membrane-potential and the mean myoplasmic [Ca²⁺]_i transient; (C) is a section through the PDE solution at the 8th z-disc from the top with the image colored by [Ca_sr].

Footnotes

1. ^A factor of h on the rhs of Equation (15) cancels with the same factor in the denominator of Equations (18, 19) turning point sources into line densities.

References

Altamirano, J., and Bers, D. M. (2007). Voltage dependence of cardiac excitation-contraction coupling: unitary Ca²⁺ current amplitude and open channel probability. Circ. Res. 101, 590–597. doi: 10.1161/CIRCRESAHA.107.152322

PubMed Abstract | CrossRef Full Text | Google Scholar

Barcenas-Ruiz, L., and Wier, W. G. (1987). Voltage dependence of intracellular [Ca²⁺]_i transients in guinea pig ventricular myocytes. Circ. Res. 61, 148–154.

PubMed Abstract | Google Scholar

Bastian, P., Blatt, M., Dedner, A., Engwer, C., Klöfkorn, R., Kornhuber, R., et al. (2008). A generic grid interface for parallel and adaptive scientific computing. Part II: Implementation and tests in DUNE. Computing 82, 121–138. doi: 10.1007/s00607-008-0004-9

CrossRef Full Text | Google Scholar

Bentele, K., and Falcke, M. (2007). Quasi-steady approximation for ion channel currents. Biophys. J. 93, 2597–2608. doi: 10.1529/biophysj.107.104299

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyman, L., Chikando, A. C., Williams, G. S., Khairallah, R. J., Kettlewell, S., Ward, C. W., et al. (2014). Calcium movement in cardiac mitochondria. Biophys. J. 107, 1289–1301. doi: 10.1016/j.bpj.2014.07.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Cannell, M. B., Berlin, J. R., and Lederer, W. J. (1987). Effect of membrane potential changes on the calcium transient in single rat cardiac muscle cells. Science 238, 1419–1423.

PubMed Abstract | Google Scholar

Cannell, M. B., Kong, C. H., Imtiaz, M. S., and Laver, D. R. (2013). Control of sarcoplasmic reticulum Ca2+ release by stochastic RyR gating within a 3D model of the cardiac dyad and importance of induction decay for CICR termination. Biophys. J. 104, 2149–2159. doi: 10.1016/j.bpj.2013.03.058

PubMed Abstract | CrossRef Full Text | Google Scholar

Chamakuri, N., and Rüdiger, S. (2012). Whole-cell simulations of hybrid stochastic and deterministic calcium dynamics in 3d geometry. J. Comp. Interdiscipl. Sci. 3, 3–18. doi: 10.6062/jcis.2012.03.01.0047

CrossRef Full Text | Google Scholar

Chamakuri, N., Rüdiger, S., Warnecke, G., and Falcke, M. (2012). Adaptive space and time numerical simulation of reactiondiffusion models for intracellular calcium dynamics. Appl. Math. Comput. 218, 10194–10210. doi: 10.1016/j.amc.2012.03.100

CrossRef Full Text | Google Scholar

Cleemann, L., Wang, W., and Morad, M. (1998). Two-dimensional confocal images of organization, density, and gating of focal Ca²⁺ release sites in rat cardiac myocytes. Proc. Natl. Acad. Sci. U.S.A. 95, 10984–10989.

PubMed Abstract | Google Scholar

Fabiato, A., and Fabiato, F. (1975). Contractions induced by a calcium-triggered release of calcium from the sarcoplasmic reticulum of single skinned cardiac cells. J. Physiol. 249, 469–495.

PubMed Abstract | Google Scholar

Hake, J., Edwards, A. G., Yu, Z., Kekenes-Huskey, P. M., Michailova, A. P., McCammon, J. A., et al. (2012). Modelling cardiac calcium sparks in a three-dimensional reconstruction of a calcium release unit. J. Physiol. 590(Pt 18), 4403–4422. doi: 10.1113/jphysiol.2012.227926

PubMed Abstract | CrossRef Full Text | Google Scholar

Hake, J., and Lines, G. T. (2008). Stochastic binding of Ca²⁺ ions in the dyadic cleft; continuous versus random walk description of diffusion. Biophys. J. 94, 4184–4201. doi: 10.1529/biophysj.106.103523

PubMed Abstract | CrossRef Full Text | Google Scholar

Hatano, A., Okada, J., Washio, T., Hisada, T., and Sugiura, S. (2011). A three-dimensional simulation model of cardiomyocyte integrating excitation-contraction coupling and metabolism. Biophys. J. 101, 2601–2610. doi: 10.1016/j.bpj.2011.10.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Jayasinghe, I. D., Cannell, M. B., and Soeller, C. (2009). Organization of ryanodine receptors, transverse tubules, and sodium-calcium exchanger in rat myocytes. Biophys. J. 97, 2664–2673. doi: 10.1016/j.bpj.2009.08.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Keener, J., and Sneyd, J. (1998). Mathematical Physiology. New York, NY: Springer.

Google Scholar

Koh, X., Srinivasan, B., Ching, H. S., and Levchenko, A. (2006). A 3D Monte Carlo analysis of the role of dyadic space geometry in spark generation. Biophys. J. 90, 1999–2014. doi: 10.1529/biophysj.105.065466

PubMed Abstract | CrossRef Full Text | Google Scholar

Lang, J. (2001). Adaptive Multilevel Solution of Nonlinear Parabolic PDE Systems, Vol. 16 of Lecture Notes in Computational Science and Engineering. Berlin: Springer-Verlag.

Maack, C., Cortassa, S., Aon, M. A., Ganesan, A. N., Liu, T., and O'Rourke, B. (2006). Elevated cytosolic Na+ decreases mitochondrial Ca2+ uptake during excitation-contraction coupling and impairs energetic adaptation in cardiac myocytes. Circ. Res. 99, 172–182. doi: 10.1161/01.RES.0000232546.92777.05

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahajan, A., Shiferaw, Y., Sato, D., Baher, A., Olcese, R., Xie, L. H., et al. (2008a). A rabbit ventricular action potential model replicating cardiac dynamics at rapid heart rates. Biophys. J. 94, 392–410. doi: 10.1529/biophysj.106.98160

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahajan, A., Shiferaw, Y., Sato, D., Baher, A., Olcese, R., Xie, L. H., et al. (2008b). CellML of the Mahjan 2008 Action Potential Model. Available online at: https://models.cellml.org/workspace/mahajan_shiferaw_sato_baher_olcese_xie_yang_chen_restrepo_karma_garfinkel_qu_weiss_2008. (Accessed July 13, 2015).

Nivala, M., de Lange, E., Rovetti, R., and Qu, Z. (2012). Computational modeling and numerical methods for spatiotemporal calcium cycling in ventricular myocytes. Front. Physiol. 3:114. doi: 10.3389/fphys.2012.00114

PubMed Abstract | CrossRef Full Text | Google Scholar

Otani, N. F., and Gilmour, R. F. (1997). Memory models for the electrical properties of local cardiac systems. J. Theor. Biol. 187, 409–436.

PubMed Abstract | Google Scholar

Restrepo, J. G., and Karma, A. (2009). Spatiotemporal intracellular calcium dynamics during cardiac alternans. Chaos 19, 037115. doi: 10.1063/1.3207835

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Schendel, T., and Falcke, M. (2010). Efficient and detailed model of the local Ca²⁺ release unit in the ventricular cardiac myocyte. Genome Inform. 22, 142–155. doi: 10.1142/9781848165786_0012

PubMed Abstract | CrossRef Full Text

Schendel, T., Thul, R., Sneyd, J., and Falcke, M. (2012). How does the ryanodine receptor in the ventricular myocyte wake up: by a single or by multiple open L-type Ca²⁺ channels? Eur. Biophys. J. 41, 27–39. doi: 10.1007/s00249-011-0755-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Scriven, D. R., and Moore, E. D. (2013). Ca²⁺ channel and Na⁺/Ca²⁺ exchange localization in cardiac myocytes. J. Mol. Cell. Cardiol. 58, 22–31. doi: 10.1016/j.yjmcc.2012.11.022

CrossRef Full Text | Google Scholar

Scriven, D. R. L., Dan, P., and Moore, E. D. W. (2000). Distribution of proteins implicated in excitation-contraction coupling in rat ventricular myocytes. Biophys. J. 79, 2682–2691. doi: 10.1016/S0006-3495(00)76506-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Scriven, D. R. L., Klimek, A., Lee, K. L., and Moore, E. D. W. (2002). The molecular architecture of calcium microdomains in rat cardiomyocytes. Ann. N. Y. Acad. Sci. 976, 488–499. doi: 10.1111/j.1749-6632.2002.tb04783.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Soeller, C., and Cannell, M. B. (1997). Numerical simulation of local calcium movements during L-type calcium channel gating in the cardiac diad. Biophys. J. 73, 97–111.

PubMed Abstract | Google Scholar

Stern, M. D. (1992). Theory of excitation-contraction coupling in cardiac muscle. Biophys. J. 63, 497–517.

PubMed Abstract | Google Scholar

Stern, M. D., Rios, E., and Maltsev, V. A. (2013). Life and death of a cardiac calcium spark. J. Gen. Physiol. 142, 257–274. doi: 10.1085/jgp.201311034

PubMed Abstract | CrossRef Full Text | Google Scholar

Stern, M. D., Song, L. S., Cheng, H., Sham, J. S., Yang, H. T., Boheler, K. R., et al. (1999). Local control models of cardiac excitation-contraction coupling. A possible role for allosteric interactions between ryanodine receptors. J. Gen. Physiol. 113, 469–489.

PubMed Abstract | Google Scholar

Takahashi, A., Camacho, P., Lechleiter, J. D., and Herman, B. (1999). Measurement of intracellular calcium. Physiol. Rev. 79, 1089–1125.

PubMed Abstract | Google Scholar

Tanskanen, A. J., Greenstein, J. L., Chen, A., Sun, S. X., and Winslow, R. L. (2007). Protein geometry and placement in the cardiac dyad influence macroscopic properties of calcium-induced calcium release. Biophys. J. 92, 3379–3396. doi: 10.1529/biophysj.106.089425

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanskanen, A. J., Greenstein, J. L., O'Rourke, B., and Winslow, R. L. (2005). The role of stochastic and modal gating of cardiac L-type Ca²⁺ channels on early after-depolarizations. Biophys. J. 88, 85–95. doi: 10.1529/biophysj.104.051508

PubMed Abstract | CrossRef Full Text | Google Scholar

Thul, R., and Falcke, M. (2004). Release currents of IP₃ receptor channel clusters and concentration profiles. Biophys. J. 86, 2660–2673. doi: 10.1016/S0006-3495(04)74322-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Torres, N. S., Sachse, F. B., Izu, L. T., Goldhaber, J. I., Spitzer, K. W., and Bridge, J. H. (2014). A modified local control model for Ca2+ transients in cardiomyocytes: junctional flux is accompanied by release from adjacent non-junctional RyRs. J. Mol. Cell. Cardiol. 68, 1–11. doi: 10.1016/j.yjmcc.2013.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Vorst, H. A. (1994). Bi-CGSTAB: a fast and smoothly converging variant of bi-cg for the solution of nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 13, 631–644.

Google Scholar

Walker, M. A., Williams, G. S., Kohl, T., Lehnart, S. E., Jafri, M. S., Greenstein, J. L., et al. (2014). Superresolution modeling of calcium release in the heart. Biophys. J. 107, 3018–3029. doi: 10.1016/j.bpj.2014.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Weber, C. R., Piacentino, V., Ginsburg, K. S., Houser, S. R., and Bers, D. M. (2002). Na(+)-Ca(2+) exchange current and submembrane [Ca(2+)] during the cardiac action potential. Circ. Res. 90, 182–189. doi: 10.1161/hh0202.103940

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, G. S., Smith, G. D., Sobie, E. A., and Jafri, M. S. (2010). Models of cardiac excitation-contraction coupling in ventricular myocytes. Math. Biosci. 226, 1–15. doi: 10.1016/j.mbs.2010.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar