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 Ca2+ content within tens of seconds. Our concept for a multiscale mathematical model of Ca2+ -induced Ca2+ 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 Ca2+ and Ca2+-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 Ca2+-concentration profiles using three previously published RyR-channel Markov schemes.


Introduction
Cardiomyocyte muscle filament shortening and lengthening is a Ca 2+ dependent process. The timing of contraction is controlled through electrical excitation via a process known as excitationcontraction-coupling (ECC). ECC is mediated through Ca 2+ , and is facilitated through an amplification process known as Ca 2+ -induced Ca 2+ release (CICR). In ventricular myocytes, CICR is controlled locally by the colocalization of L-type Ca 2+ -channels (LCCs) in the T-tubule membrane on the one side of a dyadic cleft [aka Ca 2+ 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 2+ entry from the extracellular space into the dyadic space. The influx of Ca 2+ activates RyRs, which release Ca 2+ from the sarcoplasmic reticulum (SR) (Fabiato and Fabiato, 1975). Within Ca 2+ -release-units, Ca 2+ dynamics are distinct from the myoplasm, with steep [Ca 2+ ] gradients and several-fold higher [Ca 2+ ] (Stern, 1992). The local release events (sparks) of the ∼20,000 individual CRUs sum to produce a macroscopic whole-cell Ca 2+ transient. The CRUs are coupled through myoplasmic Ca 2+ diffusion, through SR Ca 2+ 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 2+ 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 2+ release in highly localized cellsubcompartments, with dyadic space [Ca 2+ ] rising to ∼1000× myoplasmic [Ca 2+ ] (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 2+ ] changes in the dyad to tens of seconds for SR dynamics) are involved. To account for that multiscale character, recently several models with spatially distributed Ca 2+ 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 2+ 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 2+ ] 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 2+ ] and stochastic state dynamics of all individual LCand RyR-channels. The dynamics of up to 5120 CRUs is then embedded into simulations of the cellular concentration fields for Ca 2+ and Ca 2+ -binding molecules as well as membrane potential time course.

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.

PDE Model
The dynamics of the cytosolic Ca 2+ 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 2+ -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 2+exchanger contribute to bulk concentration dynamics (J NaCa ). The Na + /Ca 2+ -exchanger flux through the plasma membrane (J pm NaCa ) enters the partial differential equations as flux boundary condition. The flux J pump describes the pumping of Ca 2+ 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 2+ -binding molecules (b j , j = s, m, f ) in the cytosol include stationary and mobile buffers and Fluo-4. The total concentration b tot j is spatially homogeneous for all of the buffers. The partial differential equations for the cytosolic concentration fields and their boundary conditions are: where D and D j b are diagonal diffusion matrices. We denote the plasma membrane by Ŵ, and the cell volume by ν cell ⊂ R 3 . We simulate between 1 and 16 z-discs in a subvolume < ν cell . The expressions for the fluxes are: where (x) denotes the step function 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 2+ , cytosolic and sarcoplasmic mobile buffers and a cytosolic stationary buffer. The N c Ca 2+ 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 2+ 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 2+ -exchanger flux couple the membrane potential dynamics directly to the concentration dynamics.
Ca 2+ influx and release through LCC (I i,j LCC (t)) and RyR (I i,j RyR (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 i LCC ) and RyRs (N i RyR ) 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 2+ concentration S and Ca 2+ buffer B sr are 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 i jsr ) dynamics are described in the context of the CRU models (Equation 21). The volume of the ith junctional SR is ν i jsr . The refill flux is localized within a spherical volume with radius R i jsr in the network SR:

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: I stim is the stimulus current to depolarize the cell. I ion comprises (Mahajan et al., 2008a): I ion = I Na + I to,f + I to,s + I Kr + I Ks + I K1 + I NaK + I CaL + I NaCa .
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 2+exchanger current I NaCa is determined by the spatial integral of the corresponding fluxes in the cytosolic concentration dynamics A prefactor f NaCa for each finite element has been introduced to simulate the distribution of the Na + /Ca 2+ -exchanger according to Jayasinghe et al. (2009). The LCC-current I CaL is the sum of all 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.

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.

Stochastic Channel Gating
The state dynamics of the LC-and RyR-channels are simulated with Markov models. For the RyRs we explored three models,  (Stern et al., 1999;Cannell et al., 2013;Walker et al., 2014) For the equations of the 4-state RyR model see (Stern et al., 1999;Schendel and Falcke, 2010).
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 2+ 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) ). 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 openingrate when (c i jsr ) 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 −1 were not encountered.
For the LCCs we used the 7-state model developed by Mahajan et al. (2008a) (see Figure 3B) which exhibits both Ca 2+ dependent and Ca 2+ independent inactivation.
Each channel's state was chosen randomly according to the steady state distribution for the initial values. The transition times  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 2+ , the inactivation rate k i is a first order Hill function of Ca 2+ . The Ca 2+ dependent rates are influenced by the Ca 2+ 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)  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 2+ 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.

Ca 2+ profile in the CRU
The Ca 2+ concentration in the cleft (c di ) is modeled by the following equation in cylindrical coordinates: 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): 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 The values of φ 0 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 where h * = h 0 Z(z)dz 1 . Solving this using the Green function G((r, ϕ), (r i , ϕ i )) we find where R is the radius of the dyadic space, q = R r i , andr i = q 2 r i . c bulk is the Ca 2+ concentration at the boundary of the cleft (cylinder barrel) as computed by the PDE model. With c di can be written as r). (17) The currents I i LCC and I i RyR depend themselves on the local Ca 2+ concentration at the channel mouth: with δ = 2F/(RT). Here, 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.

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 2+ 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: We adjusted β so that the Ca 2+ concentration at the channel mouths and the average dyadic Ca 2+ concentrations matched the values given in the paper by Cannell et al. (2013). An approximation of [Ca 2+ ] i that would be measured by a singlewavelength Fluo-4 experimental recording (denoted [Ca 2+ ] exp i ) was calculated from the Ca 2+ -bound Fluo-4. The in vitro calibration approach described in Takahashi et al. (1999) was used: where K d is the dissociation constant of Fluo-4, F is the experimentally measured fluorescence intensity (here the [Ca 2+bound Fluo-4]), F max is the measured fluorescence intensity in Ca 2+ -saturated dye (here this is set as the maximum [Ca 2+ -bound Fluo-4], i.e., b tot Fluo−4 ), and F min is the measured fluorescence intensity in the absence of Ca 2+ (here set to zero) (see Table 1).

Modeling jSR dynamics
Each dyadic space is paired with its own junctional SR (jSR) compartment. We assume spatially uniform Ca 2+ concentration in the jSR. Ca 2+ 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.
β jsr = 1 + nK csqn B csqn (K csqn + c jsr ) 2 −1 S( r) denotes the free nSR Ca 2+ concentration at the location 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 −4 ms).

Module Interactions
The modules interact on a time scale longer than a single iteration step by the dependencies of the dynamics on Ca 2+ , 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 2+ 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 2+ -exchanger flux (see also Figure 1).

Bulk Concentrations-Electrophysiology
We use the Ca 2+ concentration spatially resolved at the plasma membrane for calculating local values of J pm NaCa 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 pm NaCa serves as boundary condition and J NaCa as bulk source term for the PDEs.

Bulk Concentrations-CRUs
The value of the Ca 2+ 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)

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

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;. 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 −6 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

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 2+ -dependent transition rates than channels close to an open one.
We simulated 16 z-discs, each with 320 CRUs, and the Ca 2+ 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 2+ concentration profile 70.0 ms after stimulus at a single z-disc is shown in Figure 6. The local Ca 2+ dynamics and corresponding whole cell electrophysiology are shown in Figure 7A and in Supplementary Movie 1. The SR free [Ca 2+ ] is visualized in Figure 7B and in Supplementary Movie 2. Currents are shown in Figure 8.
The simulation shows that we can reproduce the whole sequence of Ca 2+ transients from the initial concentration profile in the dyadic space, via sparks to the whole cell transient. The cytoplasmic [Ca 2+ ] 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 2+ currents well (Weber et al., 2002), but show some discrepancies mostly related to Na + /Ca 2+ -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 2+ 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 2+ ] i with an early peak at ∼10 ms, followed by a reduction of [Ca 2+ ] i and a second [Ca 2+ ] 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 2+ 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 2+ ] i transient with shorter BCL (Figure 10).

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 2+ 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 2+ ] with realistic steep [Ca 2+ ] 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 2+ 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 2+ 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  (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.
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 2+ ] 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 Weber et al. (2002). The early [Ca 2+ ] 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 2+ 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 2+ ] 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 2+ 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 2+ ] i with the (Cannell et al., 2013)-model is that in our model we do not explicitly represent the [Ca 2+ ] i uptake and release by the mitochondria. Although mitochondrial [Ca 2+ ] i uptake and release is controversial on such fast timescales (Boyman et al., 2014), there is evidence that fast early mitochondrial [Ca 2+ ] i uptake can act like a fast stationary buffer, and that inhibition of fast uptake through a specific inhibitor of the mitochondrial Ca 2+ uniporter (MCU) can result in an early local maximum of the [Ca 2+ ] i transient (Maack et al., 2006). Furthermore, the data in Figure 9 show that the [Ca 2+ ] i recorded in experiments via a fluorescent Ca 2+ probe such as Fluo-4 [approximated as ([Ca 2+ ] exp i ), as defined in Equation 20] would smooth-out any early [Ca 2+ ] 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 2+ -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 2+ ] (c i jsr )]. The discrepancies between our simulations and the ODE-Mahajan model with respect to Na + /Ca 2+ -exchanger current illustrate the value of spatially resolved modeling in exploring detailed properties of Ca 2+ dynamics. Using the same Na + /Ca 2+ -exchanger model as Mahajan et al. (2008a) we were not able to reproduce the same time course of the total Na + /Ca 2+ -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 Vplateau 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 2+ -exchanger depends on the [Ca 2+ ] 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 2+ ] 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 2+ -exchanger molecule density in proximity to the CRUs to reflect the observed higher abundance of Na + /Ca 2+ -exchanger near the CRU (Scriven et al., 2000(Scriven et al., , 2002Jayasinghe 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 2+ -exchanger current through investigation of these problems.

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