## Abstract

In cardiac myocytes, calcium (Ca) can be released from the sarcoplasmic reticulum independently of Ca influx from voltage-dependent membrane channels. This efflux of Ca, referred to as spontaneous Ca release (SCR), is due to Ryanodine receptor fluctuations, which can induce spontaneous Ca sparks, which propagate to form Ca waves. This release of Ca can then induce delayed after-depolarizations (DADs), which can lead to arrhythmogenic-triggered activity in the heart. However, despite its importance, to date there is no mathematical model of SCR that accounts for experimentally observed features of subcellular Ca. In this article, we present an experimentally based model of SCR that reproduces the timing distribution of spontaneous Ca sparks and key features of the propagation of Ca waves emanating from these spontaneous sparks. We have coupled this model to an ionic model for the rabbit ventricular action potential to simulate SCR within several thousand cells in cardiac tissue. We implement this model to study the formation of an ectopic beat on a cable of cells that exhibit SCR-induced DADs.

- arrhythmia

during excitation-contraction coupling Ca is released from the sarcoplasmic reticulum (SR) within several thousand submicron-scale junctions distributed throughout the cell (1, 3). These junctions, also referred to as Ca-release units (CRUs), contain membrane-bound l-type Ca channels (LCC) along with a cluster of Ca-sensitive Ryanodine receptor (RyR) channels that control the flow of Ca across the local junctional SR (JSR). Under normal conditions, Ca release occurs via calcium-induced calcium release (CICR), where an LCC channel opening “triggers” Ca release from the neighboring cluster of RyR channels. Once CICR is initiated, the large gradient of Ca across the SR membrane ensures that a substantial amount of Ca is released into the cytosol. The localized release of Ca during this processes can be visualized using fluorescence imaging and is referred to as a Ca spark (1). However, under certain conditions, Ca release from the SR can occur in the absence of an LCC trigger (3, 25, 40, 44). In this case, Ca release is induced by one or a few RyR-channel random openings, which can allow sufficient Ca flux from the SR to initiate CICR, leading to a Ca spark. These local release events are referred to as spontaneous Ca sparks (1, 6). Under conditions such as Ca overload, the Ca released during a spark can also diffuse and induce Ca sparks in neighboring CRUs. This chain reaction of events can lead to Ca waves that propagate across the cell, leading to a substantial depletion of Ca from the SR (38, 40, 44). This process, referred to as spontaneous Ca release (SCR), has been studied extensively (28, 29, 35, 37) and has been implicated as potential trigger of various cardiac arrhythmias.

SCR is arrhythmogenic because the Ca released into the cytosol will stimulate Ca-sensitive membrane currents, such as the sodium-calcium exchanger, which can induce sufficient inward currents to form delayed after-depolarizations (DADs) (9, 16, 38). In cardiac tissue these events can summate to trigger ill-timed ectopic beats, leading to arrhythmogenic-triggered activity (32, 38). The role of SCR in arrhythmogenesis is highlighted by recent gene-based studies showing that mutations of Ca-cycling proteins are associated with cardiac arrhythmias (22, 24, 42, 43). In particular, Priori and coworkers have presented elegant studies (21, 22) showing that catecholaminergic polymorphic ventricular tachycardia (CPVT), a lethal familial disease, is associated with distinct mutations of the RyR molecular complex. These studies clearly demonstrated a direct connection between a specific RyR mutation and an increased propensity for DADs. Following this work, Fernandez-Velasco et al. (10) presented a detailed study of subcellular Ca in myocytes from mutant mice. They found that, upon rapid pacing, followed by a pause, myocytes from wild-type mice remained quiescent, whereas mutant mice cells exhibited spontaneous Ca waves. A comprehensive study of these mice showed a much higher incidence of spontaneous Ca sparks and Ca waves, which was enhanced by β-adrenergic stimulation and rapid pacing. This abnormal Ca activity was attributed to an increased Ca sensitivity of the mutant RyR channels in these mice.

In the clinical setting, SCR has been implicated in the genesis of a wide variety of cardiac arrhythmias. For instance cardiac-mapping studies of a rabbit model of heart failure (HF) showed that VT was clearly induced, and maintained, by focal excitations (30, 31). Electrophysiological studies of myocytes from these hearts revealed an enhanced expression of the NaCa exchanger (30), making them more prone to DADs. Thus, in this model of HF, focal excitations are due to the larger inward currents induced by Ca release during SCR. Interestingly, SCR has also been implicated in the genesis of atrial fibrillation (AF). In an elegant study, Hove-Madsen et al. (14) extracted myocytes from human patients diagnosed with AF and found that the frequency of SCR is substantially increased in these myocytes compared with control. This work is consistent with several studies showing that AF can be perpetuated by focal excitations in the pulmonary veins, which are critically sensitive to agents that interfere with the Ca-cycling machinery (7, 13). Again, studies of isolated myocytes from these regions reveal a higher frequency and incidence of SCR. These studies, along with work from other laboratories (see Ref. 23 for a recent review), have demonstrated that SCR is an important factor in the genesis of several distinct cardiac arrhythmias.

Experimentally based computational models of voltage and calcium have led to several important insights into cardiac arrhythmias at the cell and tissue level (5, 11, 18, 26, 27, 34). Also, numerical studies of subcellular Ca have provided insights on the underlying dynamics of SCR (12, 15, 33). These studies involve the detailed analysis of the statistics of spontaneous Ca sparks (12), as well as the saltatory propagation of Ca waves (15, 19). However, to date, there is no mathematical model of the voltage and Ca dynamics of a cardiac myocyte that accounts for the spatial distribution of Ca during SCR and that can also be implemented in cardiac tissue simulations. Subsequently, it is not well understood how SCR in a population of cells can lead triggered arrhythmias in the whole heart. In this study, we present a phenomenological model of SCR that accounts for experimentally observed features (41), such as the probabilistic timing of spontaneous Ca sparks and the propagation properties of Ca waves. We have coupled this model to an existing action potential (AP) model of the rabbit ventricular myocyte (27) to model the formation of triggered activity in cardiac tissue.

## MATERIALS AND METHODS

### A Spatially Distributed Model of Subcellular Ca Cycling

To model the spatiotemporal distribution of calcium in cardiac cells, we have implemented a recently developed mathematical model attributable to Restrepo et al. (33). In this model, the cardiac cell is represented as a collection of spatially distributed subcellular compartments, which contain key elements of the Ca-cycling machinery (Fig. 1). Each subcellular unit, which we will refer to as the CRU, consists of *1*) the dyadic junction, where a few LCC (1–5 in number) on the cell membrane are in close proximity to a cluster of 50–150 RyR channels attached to the JSR; *2*) the submembrane space, which represents a volume in the vicinity of the sarcolemma. It is necessary to include this space because Ca-sensitive membrane currents are regulated by Ca near the cell membrane, which varies more rapidly than the average bulk Ca concentration; *3*) the bulk myoplasm, which represents the volume of space into which Ca diffuses; *4*) the local JSR, which is the portion of the SR that is close to the cell membrane; and *5*) the network SR, which represents the bulk SR network. The Ca concentrations in these compartments are denoted as *c*_{p}^{k}, *c*_{s}^{k}, *c*_{i}^{k}, *c*_{jsr}^{k}, and *c*_{sr}^{k}, respectively, where the superscript labels the *k*^{th} CRU. To model the dynamics of each CRU we simulate exactly the stochastic time evolution of all RyR channels within each CRU. In this study, we model an array of coupled cells, where each cell is represented as a linear array of 200 CRUs. Nearest-neighbor CRUs are coupled via Ca diffusion except at cell boundaries, where we assume Ca diffusion is negligible. Membrane voltage of the *i*^{th} cell *V*_{i} is modelled using
*N* CRUs in the *i*^{th} cell. Here, *C*_{m} is the cell-membrane capacitance, *I*_{Na} is the sodium current, *I*_{K} is the total potassium current, *I*_{stim} denotes the stimulus current, and the superscript *i* denotes the cell number. Also, the local LCC current within the *k*^{th} CRU of the *i*^{th} cell is denoted as *I*_{Ca}^{i,k}, and *I*_{NaCa}^{i,k} represents the Ca-sensitive Na-Ca-exchange current. Coupling between cells is modeled using a gap-junction conductance *G*_{g}, which is chosen such that the effective voltage diffusion in cardiac tissue (*D* = *G*_{g}*L*^{2}/*C*_{m}) is *D* = 1.0 × 10^{-3} cm^{2}/ms, where the cell length is *L* = 100 μm. Here, all ion currents are formulated according to the rabbit ventricular myocyte model of Mahajan et al. (27). Model parameters for the Ca-cycling system used are the same as in Resterepo et al. (33), with minor modifications, which are described in Chen et al. (5).

## RESULTS

### Qualitative Features of SCR in Cardiac Myocytes

#### Overview of experimental findings.

Numerous studies have investigated the spatiotemporal properties of subcellular Ca during SCR. In particular, Wasserstrom et al. (41) applied confocal line-scan imaging of subcellular Ca within a single or two or three subepicardial myocytes of intact rat heart. Using this approach, it was possible to observe the spatial distribution of Ca during SCR at the level of individual sarcomeres. On the basis of this study, we can summarize several key features of SCR. First, it is evident that SCR occurs as Ca waves, which originate at multiple distinct sources or initiation sites in the cell. Following Ref. 41, we will refer to these wave-initiation sites as SCR sources (SCRS). Second, once a Ca wave is initiated, it propagates until it collides either with the cell boundary or with waves originating from other SCRS. Therefore, SCR can be divided into distinct episodes of wave initiation and propagation. Finally, as the SR load is increased we observe that *1*) the number of SCRS for each episode of SCR increases; *2*) the time interval to the first episode of SCR decreases; *3*) more episodes of SCR occur, i.e., the time interval between episodes of SCR decreases; and *4*) the speed of Ca waves increases.

The fluorescence-imaging data presented above reveal the rich spatiotemporal dynamics of subcellular Ca during SCR. Several key temporal and spatial features of SCR can be noted. First, the timing of SCR is dictated by the timing distribution of a few SCRS. More precisely, the timing of SCR is determined by the first (∼1–5) SCRS that occur in the cell. However, the location and timing of these SCRS appear stochastic in nature, i.e., SCRS occur in seemingly random locations and random times. Thus, to understand the timing of SCR, it is essential to characterize the stochastic properties (in space and time) of these discrete SCRS events. Second, the spatial distribution of Ca during SCR is determined by the number of SCRS per episode and the speed of Ca-wave propagation. The experimental images suggest that, when a SCRS occurs, the Ca waves propagate until they extinguish at the boundaries, or collide with waves from a different source. Therefore, the amount of Ca release during an episode of SCR is determined by the speed of propagation of Ca waves and the number of wave sources. These features are essential to a description of SCR because they determine the timing and the rate at which Ca is released in the cell during SCR. Thus a mathematical model of SCR used to study the occurrence of Ca-mediated triggered activity must correctly account for these features.

#### Qualitative properties of SCR using a spatially distributed stochastic model.

To model the qualitative features of SCR we will first simulate the dynamics of a linear array of coupled CRUs. Although this model does not capture the full complexity of a real cell, it can provide some key insights into the mechanisms that underlie SCR. To mimic the experimental protocol in Ref. 41, we pace our array of coupled CRUs for 10 beats followed by a 2-s pause and repeat this protocol for different values of the initial SR Ca load in the cell. In Fig. 2, we plot the space-time distribution of subcellular Ca for an array of 1,000 CRUs. Here, each CRU is diffusively coupled to its nearest neighbor. However, to model cell boundaries, we set diffusion to be zero every 200 CRUs, which is consistent with experimental results showing that Ca waves in intact heart rarely propagate between cells (20, 41). As shown in Fig. 2, we find that at a low SR load ([Ca]_{sr} = 1,400 μM) subcellular Ca release occurs via triggered Ca sparks (Fig. 2*A*) during the AP followed by a few spontaneous Ca sparks, which begin to occur after roughly 500 ms after the AP upstroke of the last paced beat. Increasing the SR Ca load ([Ca]_{sr} = 1,600 μM), we find that a larger number of spontaneous Ca sparks occur after roughly ∼500 ms following the last paced beat (Fig. 2*B*). Furthermore, a population of these spontaneous events transition to Ca waves, which can propagate across several (∼10) CRUs before they extinguish. At an even higher SR Ca load ([Ca]_{sr} = 1,800 μM), we find that spontaneous Ca sparks tend to occur at shorter coupling intervals ∼400 ms and that most spontaneous Ca sparks transition to Ca waves, which propagate across the cell (Fig. 2*C*). These waves propagate until they extinguish either by colliding with the cell boundaries or with other waves from different sources.

Our 1D stochastic simulation provides several key insights into the underlying mechanisms for SCR in cardiac myocytes. First, we note that SCR in our simulation is initiated by a spontaneous Ca spark, which propagates as a wave. These spontaneous sparks occur due to random fluctuations of RyR channels within a cluster. Now, if the SR load is sufficiently high, these random fluctuations can allow enough SR Ca to diffuse into the dyadic space to induce a Ca spark via CICR. Furthermore, as the SR load increases, the probability that spontaneous Ca sparks occur increases rapidly. For even larger SR loads the amount of Ca released is sufficient to induce Ca release in neighboring junctions. When this occurs, a fire-diffuse-fire wave can occur, which can propagate in the cell. As this wave propagates, a substantial amount of Ca can be released from the SR, leading to a rise in intracellular Ca concentration. As in the experiments (41), here, the timing of SCR is essentially dictated by the first few SCRS that occur. Furthermore, the amount of Ca released into the cell is critically dependent on the number of wave sources and the propagation speed of the ensuing Ca waves. These results highlight the spatiotemporal complexity of subcellular Ca during SCR and the necessity to account for these features to accurately model SCR.

### A Phenomenological Model of SCR

#### Modeling subcellular Ca release by computing the rate of Ca spark recruitment.

A major limitation of the spatially distributed stochastic model is that it is too computationally demanding to be applied in cardiac tissue. This is because a simulation of a single CRU involves the exact stochastic simulation of roughly 100 independent RyR channels. However, a typical cardiac myocyte consists of about 10^{4} CRUs, and, therefore, to model a small tissue with ∼10^{3} cells, will require a simulation of roughly 10^{9} channels, which is computationally prohibitive, if not intractable. Thus any model that attempts to describe SCR on the tissue scale must account for key subcellular features without a detailed accounting of the individual ion channels involved. Furthermore, SCR is due to Ca waves propagating in the complex three-dimensional geometry of a cell. Modeling the detailed propagation of these waves is computationally demanding even at the single cell level, let alone in cardiac tissue. To overcome this limitation we have developed a simplified phenomenological model of subcellular Ca dynamics which can *1*) reproduce the correct timing of SCR without explicit simulation of stochastic ion channel transitions, and *2*) account for the release of Ca attributable to propagating Ca waves, which are initiated at multiple sites in the cell. Our approach is to follow Shiferaw et al. (34) and write the current flux draining the SR at time *t* as a summation of discrete fluxes attributable to Ca sparks that have been recruited in the cell. The time dependence of the flux attributable to each spark is assumed to have a simple exponential form
*J*_{s}(*t*,*t*′) represents the flux at time *t* attributable to a spark that is initiated at time *t*′. Following Ref. 27, we assume that the current flux draining the SR is proportional to the JSR Ca concentration [*c*_{jsr}(*t*)], and where τ_{d} is the average lifetime of a Ca spark. At the whole cell level, the total flux draining the SR is simply a summation of discrete fluxes attributable to Ca sparks distributed throughout the cell. Therefore, as shown in Ref. 34, we can approximate the total flux as an integral over discrete fluxes using
*R*(*t*′) is the rate at which new sparks are recruited in the cell at time *t*′. This integral can be shown (27) to satisfy the ordinary differential equation
*T* = τ_{d}/(1 −
*dc*_{jsr}/*dt*. Therefore, the total flux of Ca draining the SR at time *t* can be computed directly from the rate of spark recruitment *R*(*t*). The factors that control the rate of spark recruitment are *1*) the number of LCC openings that deliver Ca into the CRU, which then induces a Ca spark via CICR, and *2*) the rate at which Ca sparks occur independently of LCC channel openings via RyR channel fluctuations and/or propagating Ca waves. Thus the rate of spark recruitment in the cell can be written as a sum
*R*_{ica}(*t*) is the rate at which sparks are recruited by LCC channel openings, and where *R*_{scr}(*t*) is the rate of recruitment via spontaneous Ca sparks and waves. In the absence of SCR *R*(*t*) is determined by *R*_{ica}(*t*), where the recruitment rate is dictated directly by the number of LCC openings and the probability that an LCC opening will trigger a Ca spark. In the model of Mahajan et al. (27), this rate was modeled using a phenomenological relationship that reproduced the main features of *I*_{Ca}-mediated excitation-contraction coupling. In the following section, we will present a method to compute *R*_{scr}(*t*) that accounts for Ca release attributable to spontaneous Ca sparks and waves.

#### Modeling the spark-recruitment rate using a simplified lattice model.

To model SCR it is necessary to compute the rate of spark recruitment attributable to a Ca wave. To construct such a model we rely on confocal line-scan imaging that shows that Ca waves originate from localized point sources in the cell and then propagate as an approximately spherical wave that collides first with the transverse walls and then proceeds along the longitudinal direction (2, 36, 37, 41). Because the longitudinal direction is roughly three to five times longer than the transverse, the majority of Ca sparks are recruited by the approximately planar wave front moving in this direction (37). Now, it is well known from ultrastructure studies that a cardiac myocyte can be divided into an array of pancake-like sections of width ∼1–2 μm and diameter 20–50 μm, referred to as sarcomeres. In a typical cell there are roughly 100–150 sarcomeres, and it is generally believed that each sarcomere contains roughly the same amount and distribution of proteins involved in Ca cycling. In this study, we will take advantage of this geometry to represent a cardiac myocyte as a linear array of coupled sarcomeres. Using this reduced description of the cell, we will proceed to model SCR by specifying the rate of spontaneous Ca sparks for each sarcomere and the velocity of a Ca wave front propagating along the long axis, as it moves from one sarcomere to the next.

To implement our coarse graining we will describe a cardiac myocyte as a 1D lattice of 150 sites, where each site corresponds to a sarcomere. To model the subcellular dynamics of Ca during SCR, we let each lattice site be in one of two states, which we label 1 or 0, corresponding to an active or inactive sarcomere, respectively. This reduction is motivated by the observation that a planar Ca wave propagating in the longitudinal direction will excite all CRUs within a sarcomere at roughly the same time. Thus we can treat each sarcomere as being active or inactive depending on whether all CRUs in that sarcomere have been triggered or not. To proceed we model spontaneous spark formation at the *i*^{th} sarcomere by assigning the probability per unit time that a spontaneous Ca spark originates at that sarcomere. In this study, we will treat each sarcomere as equivalent and assign a spontaneous rate Γ_{s}, so that Γ_{s}*dt* is the probability that a Ca spark occurs in a given sarcomere in the time interval *dt*. This rate will be dependent on the local Ca concentration variables such as the JSR Ca concentration *c*_{jsr}(*t*), and the average dyadic Ca concentration *c*_{p}(*t*). In a later section we will describe in detail a method to determine this transition rate directly from the available experimental data (41).

To model dynamics on the lattice we will require that, once a lattice site acquires the value 1, then it will remain active (retain the value 1) for a fixed duration τ_{d}, after which it will be transitioned back to 0. We will take τ_{d} to be roughly the average spark life time, because all CRUs in a sarcomere are activated at roughly the same time, and release Ca for a duration comparable to the average spark duration. To model Ca-wave propagation on our lattice we simply assign a rule that, if a site transitions to 1, then a neighboring site, with value 0, will also be transition to 1 after a time τ_{w} has elapsed, with a transmission probability denoted as *p*. The transmission time τ_{w} will determine the propagation velocity on the lattice and will be fixed by fitting to the longitudinal Ca-wave speed observed experimentally. Also, the nearest neighbor transmission probability *p* will be taken to be SR load dependent because the excitability of CRUs is regulated primarily by the SR Ca load. Again, details of this choice will be given in the following section. To proceed, we then determine *R*_{scr}(*t*) by counting the rate at which lattice sites are recruited on our 1D lattice. We count this rate at each time step by directly computing the spark rate
*N*(*t*,*t* − Δ) is the number of 0 to 1 transitions on our lattice in the time interval [*t* − Δ, *t*], *n*_{s} is a constant that corresponds to the number of CRUs within each sarcomere, and where Δ = 3 ms is a fixed time interval. In this way we can effectively model the rate of spark recruitment via the stochastic formation of SCRS and the propagation of Ca waves.

#### Determination of model parameters.

The important feature of the lattice model is that it can correctly account for two essential features of SCR. First, by assigning the correct spontaneous Ca spark rate Γ_{s} we can account for the timing distribution of SCRS. Second, once these waves are initiated, our model can account for the rate of spark recruitment attributable to propagating Ca waves. Furthermore, because we consider an array of coupled sarcomeres, our model can also describe the occurrence of multiple SCRS occurring in the cell. This is an important feature because it is known experimentally that the amount of Ca released into the cell is sensitive to the number of distinct wave sources (41). Also, the model is computationally fast because we model Ca dynamics only on the scale of sarcomeres, rather than at the scale of individual ion channels. The parameters of the lattice model are found by fitting directly to experimental measurements and, in the absence of data, to our 1D array of coupled CRUs.

##### MODELING THE RATE OF SPONTANEOUS CA SPARKS.

Here we describe a systematic method to determine the rate of spontaneous Ca sparks that occur in a given sarcomere. As a starting point we will first determine the statistics of spontaneous sparks for a single CRU. Using this result, we can then infer the statistics of an ensemble of CRUs to determine the statistics of spontaneous sparks within a sarcomere. To determine the spontaneous spark statistics for an isolated CRU, we will compute *p*_{i}(*t*), defined so that that *p*_{i}(*t*)*dt* is the probability that the first spontaneous Ca spark occurs in the time interval *t*,*t* + *dt* in the *i*^{th} CRU. To compute *p*_{i}(*t*), we simulate an isolated CRU by clamping all Ca concentration variables except the dyadic junction Ca concentration (*c*_{p}^{i}) and the local JSR concentration (*c*_{jsr}^{i}) (see Fig. 1). We then measured the time until a spontaneous Ca spark is induced in that CRU via a random opening of several RyR channels in the cluster. Here, our criterion for a spark is that greater than 80% of the RyR channels in the cluster transition from the closed to the open state. In Fig. 3 we plot *p*_{i}(*t*) for three simulations runs, where the initial conditions were *c*_{p}^{i} = 0.3 μM, and with *c*_{jsr}^{i} at three different Ca concentrations. From Fig. 3 we see that *p*_{i}(*t*) can be well approximated by an exponential distribution. Thus, under these conditions, we can model the first-latency distribution of a single CRU as a Poisson process with distribution
_{i} is the rate at which a spontaneous Ca spark forms at an isolated CRU. Thus the probability that a spark occurs at the *i*^{th} CRU in small time interval *dt* is simply given by γ_{i}*dt*. Now, given an ensemble of *n* CRUs, it is straightforward to show that the probability of spontaneous Ca sparks is also exponentially distributed with a rate that is Γ = *n*γ, where
_{s} = *m*γ, where *m* is the number of CRUs within a sarcomere, and where γ is the average over all CRUs in that sarcomere. In this study we will treat all sarcomeres as equivalent although this restriction can be relaxed by simply assigning a different spontaneous rate for each sarcomere.

The above analysis reveals that we can model the timing statistics of SCR if we know the rate (Γ_{s}) of spontaneous Ca sparks for each sarcomere. To determine this rate, we rely on recent experiments (41), where the first-latency distribution of spontaneous Ca sparks was measured at the level of the whole cell. In these experiments, the distribution of the first spontaneous Ca spark, which induced a Ca wave, was measured for a whole myocyte in the intact rat heart. In this case, the SR load is not constant, as a substantial amount of Ca is released from the SR during the last AP, followed by uptake into the SR, which occurs over a time scale of a few hundred milliseconds. Thus, to reproduce the experimental whole cell data, it is necessary to incorporate the dependence of Γ_{s} on the Ca concentration variables. Now, because the SR load is the most important variable that regulates the probability of spontaneous Ca sparks, we will only consider the functional dependence of Γ_{s} on the SR concentration (*c*_{sr}). Thus our strategy is to choose a functional form for the rate Γ_{s}(*c*_{sr}), which reproduces the correct experimentally observed first-latency distribution for the whole cell. In particular, we will ensure that the average and standard deviation (SD) of the experimentally measured whole cell first-latency distribution is reproduced by our lattice of 150 sarcomeres. To fit the experimental data, we use a spontaneous spark rate of the form
_{s} is given in units of sparks/ms, *x* is the SR concentration in units of mM, and where we have taken the spontaneous rate to be zero for SR loads less than 1 mM, because SCR is rare at these concentrations and will not play a role in Ca cycling. Also, for simplicity, at large SR loads greater than 2.5 mM, we set the transition rate to a constant. Here, the functional form of Γ_{s} is chosen to fit the experimental data in Ref. 41, where the first-latency distribution is measured during a pause following rapid pacing at a cycle length of 200 ms. To compare with experiments we note that the experiments were performed with different external Ca to vary SR load. We follow Ref. 41 and adjust our model parameters so that external Ca concentrations of [Ca]_{e} = 4, 5, and 6 mM correspond in our model to *c*_{sr} = 1,400, 1,600, and 1,800 μM, respectively. In Fig. 4 we plot the whole cell, first-latency distribution for the three SR Ca loads considered. In Table 1 we show the average latency and SD at the indicated SR loads and the corresponding experimental measurements reproduced from Ref. 41. As shown, at the given SR loads, the average and SD of the first-latency distribution is comparable to the experimental measurements.

##### MODELING WAVE PROPAGATION ON THE LATTICE.

To model Ca-wave propagation in the longitudinal direction, our lattice model must reproduce the experimentally observed Ca-wave speed. Here, we again rely on recent experimental data measuring the speed (in the longitudinal direction) of Ca waves as a function of SR load (41). To incorporate this data into the model, we implement a rule that, once a lattice site is excited (i.e., a 0–1 transition), then a neighboring site in *state 0* can spark (transition to value 1) only after a waiting time τ_{w}. This waiting time τ_{w} will dictate the velocity of Ca-wave propagation since τ_{w} ≈ *d*_{s}/*v*_{Ca}, where *d*_{s} is roughly the distance between the centers of adjacent sarcomeres, and where *v*_{Ca} is the Ca-wave speed. To estimate τ_{w}, we find a best-fit function to the experimental measurements of *v*_{Ca} vs. SR load (see Fig. 5 caption) and compute the waiting time for a fixed sarcomere spacing (*d*_{s} = 1 μm). In this way we ensure that the propagation speed on our 1D lattice matches exactly with the longitudinal Ca-wave speed measured experimentally.

##### NEAREST-NEIGHBOR TRANSMISSION PROBABILITY.

Observation of the data presented in Ref. 41 shows that, once a Ca wave is initiated, it propagates until it reaches the ends of the cell or collides with waves originating from a different source. This result indicates that, for the experimental conditions considered, the longitudinal Ca waves induced during SCR are robust and do not abort. This is in contrast to our detailed stochastic simulations shown in Fig. 2, where, at moderate SR Ca loads, most Ca waves failed before colliding with other waves or with the cell boundary. This result is due to the fact that, in our simulations, the probability that a Ca spark at one CRU will trigger Ca release from a neighboring CRU is strongly dependent on the SR load. In fact we find that, indeed, in our system of coupled CRUs, the transmission probability *p* is a sigmoid function of the SR load, so that at low loads Ca waves do not propagate, whereas at high loads Ca waves can propagate a substantial distance in the cell. This strong dependence of Ca-wave propagation on SR load is likely to apply in the complex 3D geometry of the cell although the precise functional dependence on the SR load is not known. Here, we will model the experimental data by taking the transmission probability to be a sigmoid function of the SR load, where the threshold for propagation is taken to be less than the smallest SR Ca load considered in the experiment in Ref. 41, i.e., with [Ca]_{e} = 4 mM, which in our model corresponds to an SR concentration of *c*_{sr} = 1.4 mM. In this way, at the SR loads that correspond to the experimental conditions, Ca waves propagate robustly in the longitudinal direction. Thus we will pick the onset of wave propagation to occur at an SR load of 1.0 mM. To proceed we chose a function of the form *p*(*x*) = {tanh [λ(x − x_{o})] + 1}/ 2, where x_{o} = 1.0 mM gives the SR threshold for the spark to wave transition, and where λ = 5.6 mM^{-1} is the width of this threshold, which we have taken from fits to our detailed stochastic simulation given the absence of experimental data.

### Incorporating SCR in a Mathematical Model of the Rabbit Ventricular Myocyte

We have coupled our lattice model of SCR to an AP model of the rabbit ventricular myocyte attributable to Mahajan et al. (27). In this model Ca release is modeled by computing the total rate of spark recruitment as described by *Eq. 5*. To account for SCR we compute *R*_{scr}(*t*) using our lattice model to compute the spark rate attributable to SCR. Model parameters for our SCR model are given in Table 2 and have been adjusted so that the subcellular distribution of Ca, as described above, corresponds to that observed in the experiments in Ref. 41. To investigate the properties of SCR, we pace the cell at fixed cycle length (CL) for 100 beats followed by a 2-s pause. In Fig. 6, we plot the space-time distribution of active sarcomeres on our 1D lattice along with the current flux *J*_{rel} draining the SR and the Ca transient in the cell. Here we show the spatial distribution of Ca for five cells, where each cell is assigned a lattice of 150 coupled sarcomeres, and where the membrane voltage is modeled using the standard cable equation
*D* = 1.0 × 10^{-3} cm^{2}/ms is the effective diffusion coefficient of voltage in rabbit cardiac tissue, and where *I*_{ion} is formulated according to the Mahajan AP model. In Fig. 6 we see that, as the SR load is increased, the average latency to the first spontaneous Ca spark decreases. Furthermore, the number of SCRS increases with SR load as seen experimentally. As a result, the rate of spark recruitment increases as the SR load increases. This leads to a larger current *J*_{rel} draining the SR and a larger spontaneous Ca transient.

To validate the model, we have also studied the CL dependence of SCR. To study these properties we pace a cell at constant CL until steady state is reached (100 beats). We then cease pacing and measure the time between the last paced beat and the peak of the first spontaneous Ca transient (Fig. 7*A*), i.e., the Ca transient attributable to SCR activity. We then measured the maximum change of diastolic Ca as shown, which we denote by Δ*c*_{i}^{max}, which gives a measure of the amount of Ca released during the first episode of SCR. In Fig. 7*B* we plot the average time to peak following cessation of pacing as a function of CL. This average is computed by running 5,000 independent simulation runs. As shown, the average time decreases as the pacing rate is increased. In Fig. 7*C* we plot Δ*c*_{i}^{max} vs. CL showing that the amount of Ca released during SCR increases with decreasing CL. Both of these relationships can be explained by the fact the SR Ca content increases as the cell is paced faster (8). Thus SCR occurs earlier and is more synchronized as the pacing rate is increased.

We have also applied our coupled voltage-Ca model to study the properties of Ca-induced DADs. In Fig. 8*A* we plot a typical voltage deflection attributable to SCR. To characterize this excitation we measured the voltage difference δ*V* from the peak of the voltage to the resting potential. To determine the relationship between the Ca released and the corresponding voltage deflection, we paced the cell to steady state CL = 500 ms followed by a pause. We then measured δ*V* and Δ*c*_{i}^{max} for the first episode of SCR following the last paced AP. In Fig. 8*B* we plot δ*V* vs. Δ*c*_{i}^{max} measured for 5,000 independent simulation runs. As expected we find that δ*V* is roughly proportional to the corresponding change in the Ca transient Δ*c*_{i}^{max}.

#### Ectopic beat on a strip of cardiac tissue.

Our phenomenological model is roughly three orders of magnitude faster than our stochastic model of coupled CRUs and can therefore be used to simulate Ca-induced triggered activity in cardiac tissue. Here, we demonstrate the computational feasibility of tissue simulations by simulating the formation of an SCR-induced ectopic beat in a strip of cardiac tissue. To generate triggered activity we first modified the model parameters to create a “DAD-prone” cell model (see Fig. 9), where SCR is sufficient to induce a triggered AP, as opposed to a “normal” cell in which SCR typically does not induce enough inward current to generate a propagating beat. To accomplish this we decreased *I*_{K1} by 30% and then increased the initial SR load to 2,000 μM so that SCR during a pause following three paced beats typically induced an AP at the single cell level. We then formed a strip of cells of size 5 × 200 in the center of a 5 × 2,000 cable (Fig. 9). As shown, SCR activity in the DAD-prone region of tissue is sufficient to induce a propagating ectopic beat that originates from the DAD-prone region. Here, we show the voltage along the cable at the specified time intervals. In this simulation of 10^{4} cells, 1 s of real time was computed in a little more than 5 min of simulation time on a standard single processor workstation (3 GHz Intel Processor).

## DISCUSSION

In this study we have applied a detailed mathematical model of subcellular Ca to study SCR in cardiac myocytes. Our findings expand on recent experimental and theoretical work (41), showing that SCR occurs in discrete episodes where a few SCRS initiate propagating Ca waves, which annihilate as they collide with other waves or the cell boundary. In particular, we find that *1*) the timing of SCR is dictated by the first-latency distribution of the first few SCRS that occur in the cell and *2*) the rate of Ca released into the cell is determined by the number of SCRS per episode of SCR and the propagation speed of the triggered Ca waves. According to our detailed stochastic model the formation of SCRS is due to random fluctuation of RyRs within the ensemble of CRUs distributed in the cell. The timing of statistics of SCRS is governed by a Poisson distribution, which is characterized by the rate Γ of these spontaneous release events. Furthermore, this rate can be related to the spontaneous firing rate of a single CRU in the cell (γ). Assuming all CRUs are identical, it is straightforward to show that the average time for the first SCRS in a cell is
*t*_{CRU} is the average first-latency time for a single CRU, and where *n* is the number of CRUs in the cell. Now, in a typical ventricular myocyte we expect that there are *n* = 10^{4} dyadic junctions (1), so that the waiting time for an SCRS in the whole cell is roughly four orders of magnitude smaller than that of an isolated CRU. This result is relevant to arrhythmogenesis because triggered activity is likely to be dangerous when it occurs between beats following an AP, which is roughly on a time scale in the few 100-ms range. Thus for SCR to be dangerous then *t*_{CELL} ∼ 100 ms so that *t*_{CRU} ∼ 1 × 10^{6} ms, which implies that spontaneous activity at a single CRU might be extremely rare, yet a potentially arrhythmogenic episode of SCR is likely to occur between beats. Furthermore, because the number of CRUs in the myocyte is large, it is unlikely that the firing statistics of SCR are dominated by one or a few CRUs. This is because for a particular CRU (say the *k*^{th} CRU) to dictate the timing of SCR, then the average latency for that CRU should be roughly *t*_{k} ∼ *t*_{CELL}, so that the waiting time for that CRU should be ∼10^{4} times shorter than the average CRU in the cell. This is unlikely because we do not expect variations in timing statistics to vary to such a large degree in the cell. However, it may also be possible that the firing statistics of a single CRU is extremely sensitive to local properties such as RyR kinetics and the volume of the dyadic junction (5, 33). In this case, it may be possible that a single CRU can, on average, fire earlier than all other CRUs in the cell and therefore dictate the timing of SCR at the whole cell level. These features highlight the importance of treating the cell as an ensemble of stochastic components and, in particular, that Ca-mediated triggered activity is dictated by the statistics of the earliest events.

A key feature of our model is that SCR is extremely sensitive to SR load. In particular, to fit experimental data, we have taken the rate of spontaneous sparks in the cell to be a rapidly increasing function of the SR load. This fitting is consistent with our detailed stochastic model, which also predicts a strong sensitivity of spontaneous spark activation and SR load. In our model, this result is simply due to the higher gradient of Ca across the SR membrane, which ensures that, as the SR load is increased, a random RyR transition from a closed to open state is more likely to induce a Ca spark via CICR. In this study we did not include a Ca-binding site on the luminal side of the RyR channel, which has been suggested in several experimental studies (17, 39). However, in this case we expect that the spontaneous rate will be even more sensitive to SR load. This is because the added sensitivity attributable to a luminal sensor can only enhance the effect of the increased Ca concentration gradient across the SR membrane. However, it is important to realize that a luminal sensor is not essential to describe the strong dependence of SCR on the SR Ca load. A further important implication of a strong SR load dependence is that the timing variability of SCR, in a population of cells in tissue, will be sensitive to the time course of the SR Ca content. Thus, as argued previously (41), factors that influence the refilling rate of the SR load will influence the degree of synchrony of spontaneous Ca sparks. This result is important to arrhythmogenesis because triggered activity in tissue must arise from the coordination of SCR across a population of a few hundred to a thousand cells to generate enough outward current to depolarize cardiac tissue. Therefore, factors that increase the influence of the relative timing of SCR between cells will dictate whether a triggered beat will occur.

In this study we have presented a phenomenological model of SCR that is based on experimental measurements showing that SCR is attributable to multiple spontaneous Ca-release sources, which transition to Ca waves. In this model the timing distribution of SCR is accounted for by specifying the rate at which spontaneous Ca-release events occur at each sarcomere. In this way, we account for the timing statistics of the several hundred CRUs within a sarcomere without explicitly simulating the individual ion channels. Furthermore, our model formulation accounts for the spatial distribution of Ca within a cell, so that the time course of SCR can be represented as the summation of multiple spontaneous-release sources that initiate Ca waves. In this way, our phenomenological model can be used to simulate SCR within a cell in a computationally efficient manner while accounting for experimentally observed biological detail. Thus the essential parameters of the model are based directly on experimental measurements of subcellular Ca. Furthermore, we have coupled our phenomenological model of subcellular Ca with an experimentally based ionic model of the rabbit ventricular myocyte. To our knowledge, this is the first ionic model of a cardiac myocyte that accounts for the spatiotemporal dynamics of subcellular Ca during SCR. An important feature of this phenomenological model of SCR is that model parameters are determined directly from a detailed model in which the stochastic dynamics of individual ion channels are taken into account. Thus we can apply this model to study how mutations in Ca-cycling proteins influence SCR at the single cell level, multi-cell level, and in a large population of cells in cardiac tissue. This computational approach will allow us to address important questions on the nature of SCR-induced triggered activity. For instance, in a recent important study Cerrone et al. (4) demonstrated the emergence of focal activity in the purkinje fibers of CPVT mice. However, it is not completely understood why SCR-induced focal activity occurs preferentially in purkinje fibers. Also, the factors that govern the time course and morphology of SCR-induced focal activity in the complex geometry of the heart are not well understood. Our computational model can be implemented in whole heart simulations and therefore can be used to address some of these important questions relevant to the genesis of cardiac arrhythmias.

## GRANTS

This work was supported in part by the American Heart Association Grant 0830298N, and the National Heart, Lung, and Blood Institute grant R01HL101196.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

- Copyright © 2011 the American Physiological Society