Estimating the probabilities of rare arrhythmic events in multiscale computational models of cardiac cells and tissue

Ectopic heartbeats can trigger reentrant arrhythmias, leading to ventricular fibrillation and sudden cardiac death. Such events have been attributed to perturbed Ca2+ handling in cardiac myocytes leading to spontaneous Ca2+ release and delayed afterdepolarizations (DADs). However, the ways in which perturbation of specific molecular mechanisms alters the probability of ectopic beats is not understood. We present a multiscale model of cardiac tissue incorporating a biophysically detailed three-dimensional model of the ventricular myocyte. This model reproduces realistic Ca2+ waves and DADs driven by stochastic Ca2+ release channel (RyR) gating and is used to study mechanisms of DAD variability. In agreement with previous experimental and modeling studies, key factors influencing the distribution of DAD amplitude and timing include cytosolic and sarcoplasmic reticulum Ca2+ concentrations, inwardly rectifying potassium current (IK1) density, and gap junction conductance. The cardiac tissue model is used to investigate how random RyR gating gives rise to probabilistic triggered activity in a one-dimensional myocyte tissue model. A novel spatial-average filtering method for estimating the probability of extreme (i.e. rare, high-amplitude) stochastic events from a limited set of spontaneous Ca2+ release profiles is presented. These events occur when randomly organized clusters of cells exhibit synchronized, high amplitude Ca2+ release flux. It is shown how reduced IK1 density and gap junction coupling, as observed in heart failure, increase the probability of extreme DADs by multiple orders of magnitude. This method enables prediction of arrhythmia likelihood and its modulation by alterations of other cellular mechanisms.


Introduction
In cardiac myocytes, dyads are sites where the junctional sarcoplasmic reticulum (JSR) membrane closely approaches (~15 nm) invaginations of the cell membrane known as transverse tubules (TTs). Voltage-sensitive L-type calcium (Ca 2+ ) channels (LCCs) are preferentially localized to the TT membrane of the dyad, where they closely appose Ca 2+ -binding Ca 2+ -release channels known as ryanodine receptors (RyRs) in the dyad JSR membrane. Depolarization of the cell membrane during an action potential (AP) increases LCC open probability, generating a flux of Ca 2+ ions into the dyad. The resulting local increases of dyad Ca 2+ concentration ([Ca 2+ ] d ) increase RyR open probability, which when open allow flux of Ca 2+ ions from the JSR into the dyad. This process, known as Ca 2+ -induced Ca 2+ release (CICR), causes brief, spatially localized release events known as a Ca 2+ sparks [1]. Due to their synchronization, these Ca 2+ sparks cause a cell-wide rise in cytosolic [Ca 2+ ] ([Ca 2+ ] i ), leading to myofilament activation and force generation. This process, known as excitation-contraction coupling (ECC), is central to the function of the myocyte [2].
Ca 2+ sparks can also occur randomly at a single release site when the spontaneous opening of a single RyR triggers the CICR process [3]. Under conditions promoting cellular Ca 2+ overload, Ca 2+ sparks are more likely to trigger RyRs at nearby release sites, thereby generating propagating Ca 2+ waves [4]. Spontaneous Ca 2+ release events generate an inward current via the Na + /Ca 2+ exchanger (NCX), which transports 3 Na + ions into the cell for every Ca 2+ ion extruded, and additionally in canine myocytes via the Ca 2+ -activated chloride channel [5]. In diastole, this produces a net inward current, resulting in an elevation of cell membrane potential (V) known as a delayed afterdepolarization (DAD) [6]. DADs of sufficient amplitude can lead to the activation of the fast inward Na + current (I Na ) and trigger a premature AP. Gap junctions joining adjacent cells then conduct the aberrant AP across the myocardial syncytium. Such ectopic events in the heart can induce reentrant ventricular arrhythmias that lead to sudden cardiac death [7]. Furthermore, the propensity for spontaneous Ca 2+ release is increased in heart diseases such as heart failure, hypertrophy, and some forms of long-QT syndrome, which are associated with increased risk for sudden cardiac death. Therefore, understanding Ca 2+ dynamics in ventricular myocytes and the Ca 2+ handling instability that arises under pathological conditions is fundamental to our understanding of cardiac arrhythmogenesis.
Experimental studies have observed triggered activity under conditions evoking spontaneous Ca 2+ release in myocardial wedge preparations [8] and whole heart [9,10]. These studies show that the likelihood of observing ectopic foci is correlated with the degree of Ca 2+ loading.
In isolated myocytes, Ca 2+ waves are observed when the sarcoplasmic reticulum (SR) Ca 2+ load achieves a critical level [11]. However, electrotonic coupling in tissue attenuates DAD amplitude by diverting inward current to adjacent myocytes through the gap junctions. Wasserstrom et al. reported that with increasing SR Ca 2+ load, spontaneous Ca 2+ waves exhibited greater synchrony following cessation of rapid pacing in intact heart [10]. Synchronous DADs result in smaller spatial gradients in membrane potential, less loss of depolarizing current into neighboring cells, and therefore larger DAD amplitude.
An elegant theoretical study by Chen et al. analytically investigated the probability of triggered events in a 1D fiber as a first passage time problem using a minimal model of Ca 2+ release and membrane currents [12]. They showed that the expected time to a triggered event decreased according to a power law of the number of cells in the fiber, and that this effect depended on the balance of gap junction, NCX, and inwardly-rectifying potassium current (I K1 ) conductance. However, to be tractable, the model employed simplified models of the membrane currents. Here we present an approach that retains mechanistically realistic characterization of membrane current properties while at the same time is computationally tractable, permitting the evaluation of a sufficiently large number of stochastic simulations to estimate rare event probabilities.
In this study, we present a multicellular model of cardiac tissue which incorporates a stochastic biophysically detailed spatial model of the ventricular myocyte as its fundamental building block. The single cell model is used to study the roles of stochastic RyR gating and Ca 2+ wave dynamics on the statistical distribution of DADs under pathophysiological conditions. The cellular Ca 2+ load and density of I K1 are shown to be two important factors influencing the mean and variance of DAD amplitude and timing. We then develop a one-dimensional (1D) myocyte tissue model comprised of these cells, and present a method for estimating the probability of rare events that are "surrogates" for ectopic beats. These surrogate events are defined as the occurrence of (rare) high-amplitude DADs with membrane potential exceeding a threshold value V T . We develop a computationally efficient approach for estimating the probability of these threshold crossings. This method enables us to estimate how changes in model parameters influence the probability of these surrogate events. As an example, the effects of reductions of I K1 density and gap junction coupling on the distribution of DAD amplitude and therefore the probability of threshold crossings are demonstrated, providing quantitative insight into how electrophysiological remodeling affects the probability of potentially arrhythmogenic ectopic beats.

Ventricular myocyte model
We have developed a three-dimensional (3D) spatial model of a single myocyte based on the (non-spatial) Greenstein-Winslow canine ventricular myocyte model [13]. To enable reproduction of realistic Ca 2+ waves and DADs, the original model was adapted to include spatial Ca 2+ diffusion on a rectangular lattice of 25,000 Ca 2+ release sites (Fig 1) distributed within the cell. The cell was divided into 25 and 20 lattice points in the two transverse directions and 50 lattice points in the longitudinal direction. Release sites were spaced 1 and 2 μm in the transverse and longitudinal directions, respectively [14]. The time constant for longitudinal Ca 2+ diffusion was twice (i.e. 2x slower) that for the transverse direction such that the model exhibited symmetric Ca 2+ wave propagation [15]. This difference in diffusion rates arises from differences in diffusion along versus across TTs.
A sub-membrane (SM) release site compartment (Fig 1A) was added to describe the volume under the TT membrane where the Ca 2+ concentration ([Ca 2+ ]) is elevated during Ca 2+ sparks and cell-wide Ca 2+ release [16]. Detailed imaging studies of Ca 2+ release sites have revealed that RyR clusters exhibit edge-to-edge spacing of less than 100 nm [17,18]. This suggests that neighboring sites may be functionally coupled through local Ca 2+ diffusion. Therefore Ca 2+ diffusion between SM compartments was implemented to reflect Ca 2+ transport across steep [Ca 2+ ] gradients on the periphery of the release site during release [19]. The SM compartment was modeled as a cylinder encircling the TT membrane with inner radius 100 nm, outer radius 180 nm, and 1 μm axis. It was assumed that 50% of NCX are located in the TT membrane of the SM compartment, and the remaining 50% are in the sarcolemmal membrane of the cytosolic compartment [20,21]. Ca 2+ in the SM is buffered by calmodulin and sarcolemmal binding sites. Ca 2+ transport rates from the SM to the cytosol and between SM compartments were constrained to yield a realistic Ca 2+ wave threshold (~100-150 μmol/[L cytosol]) [22] and velocity (50-100 μm/s) [4] in the baseline model.
Each spatial site, with location defined by coordinate (i,j,k), is represented by a set of ordinary differential equations describing local Ca 2+ transport ( Fig 1A). [Ca 2+ ] is assumed to be uniform within the local JSR ([Ca 2+ ] JSR,i,j,k ), dyadic subspace ([Ca 2+ ] d,i,j,k ), and SM ([Ca 2+ ] SM,i,j,k ) compartments. Spatial Ca 2+ diffusion is modeled as transport between SM compartments of adjacent release sites in the 3D lattice ( Fig 1B). Model equations and parameters are given in S1 Equations and S1 Each release site contains a set of 48 RyRs and 8 LCCs that gate stochastically according to Markov chain models. The LCC model is as described in the Greenstein-Winslow model [13], with adjustments made to the rate of Ca 2+ -dependent inactivation (see S1 Text). Briefly, the LCC inactivation rate is a saturating function of [Ca 2+ ] SS , which was necessary to reproduce inactivation kinetics consistent with the original model. Note, however, that as a result AP duration does not substantially decrease at reduced SR Ca 2+ loads as exhibited previously [13]. RyR gating is described by a minimal two-state Markov model based on the work of Williams et al. [23], described in detail in Walker et al. [19]. Briefly, mean open time of each channel is 2 ms, and the opening rate is given by where k + = 1.107 × 10 −4 ms -1 μM -η is the opening rate constant, η = 2.1 is the Ca 2+ Hill coefficient, and ϕ is a [Ca] JSR,i,j,k -dependent regulation term given by

Fiber model
The 3D cell model was incorporated into a tissue-scale model of a 1D fiber of myocytes. Fig 1C  depicts the multiple biological scales represented in the model, from single ion channels to the multicellular fiber. The 3D cell model was augmented with a current carried by the gap junctions at either end of each cell (I gap ). The current from cell i into an adjacent cell i+1 is given by: where g gap is the gap junction conductance, which was adjusted to yield a conduction velocity of 55 cm/s. The membrane potential in the fiber was solved using the Crank-Nicolson method [24] with 50 μs time steps. Operator splitting was used to explicitly solve each cell model using an embedded adaptive Runge-Kutta method described previously [13].

Beta-adrenergic stimulation
Sympathetic stimulation of the heart occurs through β-adrenergic receptor activation, which activates intracellular signaling pathways, most notably Protein Kinase A (PKA) and Ca 2+ /calmodulin-activated protein kinase II (CaMKII), that increase contractility [25]. β-adrenergic stimulation is also known to be pro-arrhythmic, and can contribute to spontaneous Ca 2+ release [26] in pathological conditions. Cell model parameters were modified to reflect the effects of acute β-adrenergic stimulation. LCC open probability was increased [27] by changing the fraction of gating LCCs from 25% to 60% and setting 3-5% of the channels to gate in a high-activity mode in which the mean open time was increased from 0.5 to 5.8 ms [28,29]. Enhanced activation of inward currents was implemented for I Kr using modifications described previously [29] and for I Ks by shifting the voltage-dependence of activation by -35 mV and increasing conductance by 40% [30]. SR Ca 2+ loading was facilitated by reducing SERCA pump K d for [Ca 2+ ] i by 50% [31]. RyR opening rate was increased by either 50% to reflect increased SR Ca 2+ leak observed in experimental studies [32] or 400% to reproduce pathological behavior after ouabain overdose [33]. Unless otherwise noted, these conditions were applied to all simulations in this study.

Model properties
In order to reproduce protocols designed to measure ECC properties, membrane potential was stepped to varying test potentials for 200 ms from a holding potential of -80 mV. The dependence of normalized peak RyR and LCC Ca 2+ fluxes on the test potential and corresponding ECC gain values are similar to those observed experimentally [34] (Fig 2A and 2B). The peak normalized RyR flux curve is right-shifted with respect to that of the LCC flux curve. The ECC gain is 13.8 at 0 mV and decreases monotonically with increasingly depolarized test potentials. The lack of high gain at potentials below -10 mV is a result of consolidating the four distinct dyadic subspace compartments utilized in the previous Greenstein-Winslow canine ventricular myocyte model [13] [35], with an AP duration (APD), defined as the time to reach 90% repolarization, of approximately 320 ms (Fig 2C and 2D). The model reaches a stable steady state after~10 seconds when paced at 1 Hz (S2 Fig). The NCX current during an action potential (S1 Fig) is in agreement with experimental and theoretical studies [36][37][38]. Reducing I K1 density by 50% prolongs the action potential to 469 ms due to the reduction in repolarizing current. Note that this is shown for the first paced beat with identical initial conditions as the baseline model, as the model fails to repolarize under these conditions after continued pacing (see Discussion). Simulating the effect of β-adrenergic stimulation increases the amplitude of the AP plateau as well as [Ca 2+ ] i transient amplitude and decay rate in addition to decreasing the APD to~255 ms [29]. Additionally reducing I K1 density by 50% resulted in a marked increase in APD to~300 ms at steady state. This prolonged AP increased Ca 2+ loading, resulting in marginally higher systolic [Ca 2+ ] i .
[Ca 2+ ] NSR was clamped at increasing values to test the relationship between SR Ca 2+ load and leak. The model exhibits an exponential leak-load relationship that is similar to experimental estimates [39,40] (Fig 2E). Spontaneous Ca 2+ waves form at a threshold SR Ca 2+ load, at which wave fronts of propagating Ca 2+ sparks emanate from random regions of high spark activity. Fig 2F shows a plot of SR Ca 2+ leak along a cross-section through the center of the cell during a representative Ca 2+ wave under β-adrenergic stimulation, this wave occurred at a lower SR Ca 2+ load (82 μmol/L cytosol) compared to under baseline conditions due to greater RyR Ca 2+ sensitivity. The wave shape and velocity of 68 μm/s are similar to those observed in experimental studies [4].

Delayed afterdepolarizations during pacing
Liu et al. demonstrated that ouabain overdose causes accumulation of [Na + ] i , leading to Ca 2+ overload and DADs [33]. In addition, the authors showed that the production of reactive oxygen species, which are known to oxidize RyRs [41] and CaMKII [42], both of which enhance RyR activity, contributed to DAD generation. To induce Ca 2+ overload, model parameters were modified to simulate β-adrenergic stimulation and to reflect the conditions in Liu et al. [33] by inhibiting the Na + /K + ATPase by 90%, which leads to accumulation of intracellular Na + (>20 mM).
To demonstrate the emergence of DADs in the model, we simulated a protocol in which the RyR opening rate and [Na + ] i were controlled over time. β-adrenergic stimulation was applied and the cell was paced at 1 Hz. First, [Na + ] i was fixed to 15 mM and the RyR opening rate was ramped from 1.5x to 5x that of baseline over t = 0 to 5 s. Intermittent Ca 2+ waves, triggered by overload of JSR Ca 2+ , cause spontaneous [Ca 2+ ] i transients, which began to occur (Fig 3A). RyR sensitization initially causes an increase in [Ca 2+ ] i transient peak from 1.73 to 1.93 μM, followed by a decrease (Fig 3B) as the Ca 2+ wave threshold decreases ( Fig 3C) and cellular Ca 2+ is extruded by NCX. Many of these events begin just prior to (<100 ms before) the next stimulus (Fig 3A inset). Prominent DADs occur at t = 4 and 6 s. This causes lower Ca 2+ transient amplitude in the following beats due to the lost Ca 2+ stores. The threshold SR load for Ca 2+ waves was reduced to 80 μmol/[L cytosol] compared to the baseline model threshold of 140 μmol/[L cytosol] (see Fig 2E). This is due to high [Na + ] i and the greater RyR opening rate, which increases Ca 2+ spark frequency and results in Ca 2+ wave nucleation and propagation.
To induce triggered APs, [Na + ] i was ramped from 15 to 23 mM over t = 10 to 15 s. This causes Ca 2+ waves to form more readily by reducing extrusion of Ca 2+ via NCX in the SM compartment, thus elevating [Ca 2+ ] d at sites adjacent to Ca 2+ sparks and increasing the probability of Ca 2+ spark propagation. Note that elevated diastolic [Ca 2+ ] i has been implicated in DAD formation in experimental studies [8,43]. Diastolic [Ca 2+ ] i was at first~140 nM when the DADs are sub-threshold. Higher [Na + ] i decreases the delay until the DADs, resulting in higher diastolic [Ca 2+ ] i . A DAD of sufficient amplitude to activate I Na occurs at t = 18.7 s and triggers a spontaneous AP. Triggered APs result in greater spontaneous [Ca 2+ ] i transients due to activation of LCCs, and the cells exhibit elevated diastolic [Ca 2+ ] i in the range of~470-520 nM. This causes a reduction in the SR Ca 2+ load threshold for spontaneous release to 63 μmol/ [L cytosol] due to the resulting increase in RyR opening rate and Ca 2+ spark frequency. Note that under these conditions, triggered APs exhibit pacemaker-like automaticity, occurring every~570 ms after cessation of pacing (S3 Fig) but stop when [Na + ] i is lowered below~20 mM. This behavior is consistent with spontaneous contractions observed in mouse heart in the presence of ouabain [44].
The protocol also produced changes in AP duration (APD) (S3 Fig). The first beat has an APD of 235 ms, 27% shorter than under normal conditions due reduced inward NCX current in the presence of elevated [Na + ] i , [45]. Upon RyR sensitization, APD decreases to~215 ms due to decreased inward NCX current accompanying the smaller [Ca 2+ ] i transients. Further increasing [Na + ] i to 23 mM reduced APD to 150 ms. A previous study showed that APD increases following spontaneous Ca 2+ release due to slowed Ca 2+ -dependent inactivation of the LCCs [46]. The model does not reproduce this behavior because of the LCC Markov chain has a saturating dependence on [Ca 2+ ] d and is therefore not sensitive to Ca 2+ load (see Discussion). Rather, the lower [Ca 2+ ] i transient results in less inward NCX current, thus reducing time to repolarization.
These results demonstrate how the model reproduces SR Ca 2+ overload as a driver of DADs under pathophysiological conditions. Elevated RyR sensitivity and [Na + ] i accumulation led to DADs of sufficient amplitude to trigger action potentials during the diastolic intervals. In addition, these results illustrate the interplay between [Ca 2+ ] i and SR load dynamics during spontaneous Ca 2+ release.

Effect of SR Ca 2+ load on DADs
The relationship between SR Ca 2+ load and spontaneous Ca 2+ release was investigated. Simulations were run using initial conditions reflecting the cell state just prior to the moment when SR Ca 2+ load reaches the Ca 2+ wave threshold following an AP. Initial [Ca 2+ ] i was set to 150 nM, similar to the level during the late decay phase of a cytosolic Ca 2+ transient. Fig 4A shows DADs occurring at the five different values of initial SR Ca 2+ load shown in Fig 4B. At the highest SR Ca 2+ load, DAD amplitude is large enough to trigger an AP, as shown in simulation (v) in Fig 4A. Elevating SR Ca 2+ load reduces the delay until the spontaneous release event, consistent with the observations of Wasserstrom et al. [10]. The increase in DAD amplitude is consistent with a study by Schlotthauer and Bers, who demonstrated increased amplitude of caffeine-induced DADs at higher SR Ca 2+ loads [47]. Fig 4C shows volume renderings of [Ca 2+ ] SM at three time points in each simulation. The number of Ca 2+ wave nucleation sites (N nuc ) was defined as the number of independently formed wave front epicenters and was estimated by inspection of the volume renderings (S1 Movie). N nuc generally increases with SR Ca 2+ load, in agreement with experimental studies in intact heart [9,10]. Therefore, the increase in SR Ca 2+ load also increases RyR Ca 2+ release flux (J RyR ) by enhancing the synchrony of RyR opening and number of nucleation sites.

Ensemble properties of DADs
We hypothesized that stochastic gating of the RyRs drives variability in Ca 2+ wave dynamics and thus DAD amplitude and timing. To test this, five independent realizations were generated, each of which had identical initial conditions similar to simulation (i) from Fig 4. The pseudorandom number generator seed was varied among the realizations in order to produce independent patterns of RyR gating. Fig 5A shows the resulting DADs, which exhibit marked variability in timing and amplitude. Time of occurrence of the DAD peak varies from 520 to 1209 ms and peak amplitudes range from 2.3 to 6.2 mV. These DADs appear qualitatively similar to experimental observations in rat myocytes, with delays ranging by~1 s and amplitudẽ 2-fold, although the SR Ca 2+ load was not reported [6]. The area under the curve of each DAD, measured relative to the resting potential, varied from 0.79 in (iii) to 1.18 in (ii) (50% greater). This roughly correlated with DAD amplitude, with the exception of the prolonged DAD (iv), which had the second greatest area under the curve despite having the lowest amplitude. Thus substantial DAD variability can be attributed to the stochastic nature of RyR gating. Spontaneous Ca 2+ release generates DADs by driving an inward current through NCX [6]. Imaging, electrophysiological, and modeling [38] studies have all suggested that NCX senses a [Ca 2+ ] SM that is greater than [Ca 2+ ] i because it can be localized near the release sites [16,20,21]. Therefore, the driving force for inward NCX current is likely to be determined, in part, by the aggregate number of concurrent Ca 2+ sparks occurring across the cell. This is consistent with a study showing that peak Ca 2+ release flux is strongly correlated with the likelihood of ectopic activity [48]. An ensemble of 98 independent simulations were performed to determine how peak [Ca 2+ ] i and J RyR are related to DAD amplitude. Fig 5B and 5C show there is a strong linear correlation between the peak membrane potential during the DAD (V max ) and the maximum of [Ca 2+ ] i during the spontaneous [Ca 2+ ] i transient (r 2 = 0.950). There is an even stronger relationship between V max and the peak J RyR value (r 2 = 0.998). This confirms that the inward NCX current is primarily driven by J RyR via the resulting rise in [Ca 2+ ] SM .
Note that J RyR reflects the combined Ca 2+ release flux associated with all Ca 2+ sparks occurring at any time. It was therefore expected that the variability in DADs was the result of spatiotemporal variations in Ca 2+ wave dynamics. Fig 5D and S2 Movie show volume renderings of the Ca 2+ waves in each simulation. Waves emanate from nucleation sites of high Ca 2+ spark activity. For example, in the simulation (iv) of Fig 5D at 700 ms, a large cluster of Ca 2+ sparks is visible at the left end of the cell, and by 900 ms a wave front of Ca 2+ is seen propagating radially away from this cluster site. The random nature of nucleation site locations results in DAD variability across simulations. Simulations i and ii are both associated with the highest-amplitude DADs as well as the greatest number of nucleation sites, as shown in Fig 5D. These nucleation sites are spaced widely across the cell, resulting in independent propagating wave fronts of high [Ca 2+ ] SM that drive inward I NCX . The remaining three simulations have lower amplitudes due to fewer nucleation sites. Note that in simulation iv, two separate Ca 2+ waves form over 100 ms apart, resulting in a prolonged low-amplitude DAD with two peaks. These results are consistent with the strong correlation between V max and maximum J RyR , which reflects the timing and pattern of Ca 2+ wave formation.

Dependence of DAD distribution on Ca 2+ load and I K1 density
In this section, the statistical relationships between Ca 2+ loading and DAD amplitude and timing in ensemble simulations are investigated. Fig 6A shows variability in sub-threshold DAD respectively. This suggests that there is an increase in DAD amplitude variability at higher SR Ca 2+ loads. Note also that the width of the shaded regions decreases as SR load increased, reflecting increased DAD synchrony.
Recall that diastolic [Ca 2+ ] i plays a critical role in determining the SR Ca 2+ wave threshold during pacing (Fig 3). The effect of increasing [Ca 2+ ] i on the DAD distribution with SR Ca 2+ load held constant was next tested (Fig 6B). Elevating [Ca 2+ ] i yields identical effects as seen when increasing SR Ca 2+ load, with both reducing the delay and increasing the amplitude of DADs.
The inward rectifier K + current, I K1 , is the primary membrane current that stabilizes V at the cell's resting potential and plays a critical role in protecting the cell from triggered APs. Down-regulation of I K1 is associated with ventricular arrhythmias in diseases such as heart failure [49], Andersen's syndrome [50], and long QT syndrome [51]. Fig 6C shows SR and is consistent with the notion that total cell Ca 2+ plays a major role in setting the DAD distribution, as observed in other models [52].
Second, reducing I K1 by 50% causes a marked increase in the average and SD of V max . In this case, V max was more sensitive to inward I NCX during the DAD because I K1 rectification reduces its outward current upon membrane depolarization which increases the contribution of I NCX to the net membrane current, and this imbalance between the currents becomes greater with reduction in I K1 density. This phenomenon is illustrated in a phase plane plot of example simulations at baseline and 50% I K1 densities exhibiting both sub-and supra-threshold DADs (Fig 6E). For a given NCX current, V rises substantially higher at 50% I K1 density compared to the baseline model. The rightward shift of the trajectories also illustrates thẽ 50% reduction in peak NCX current required to trigger an AP (~3.4 to~1.7 pA/pF). Therefore, DAD amplitudes are greater for a given I NCX amplitude, causing an increase in the mean. Similarly, amplitude SD is also increased. At the lowest Ca 2+ load tested, 50% IK1 reduction caused~3x greater SD, which is mainly accounted for by the fact that the amplitude is on aver-age~2.6x greater.
DAD delay and synchrony is also strongly correlated with total cell Ca 2+ , as measured by the distribution of the time until the DAD peak occurred (Fig 6F). Increasing initial total cell Ca 2+ reduces DAD delay and increases synchrony, similar to the observations of Wasserstrom et al. [10]. In the baseline model, the standard deviation of the delay was 178 ms at the lowest SR load tested. This variability is similar to that measured by Wasserstrom et al. (162 ms) in whole rat heart paced at 1 Hz and elevated [Ca] o = 7.0 mM. DAD delay does not considerably change with 50% I K1 reduction because it is determined primarily by the timing of Ca 2+ wave formation, which is not affected by I K1 density.
With increasing initial total cell Ca 2+ , in the form of increased initial [Ca 2+ ] i or [Ca 2+ ] SR , the higher DAD amplitudes makes it more likely that the threshold membrane potential for triggering an AP (~-55 mV) will be reached. The probability that a DAD-triggered beat occurred is shown in Fig 6G. Increasing cell Ca 2+ results in a steep increase in trigger probability. This switch-like behavior is due to the DAD amplitude reaching threshold AP-triggering potential with high certainty. Reducing I K1 by 50% reduces the critical Ca 2+ load at which APs are triggered, consistent with the observed increase in DAD amplitude. The effect is significant: at normal I K1 density zero of 96 cells exhibit triggered APs when initial [Ca 2+ ] is 5.8 fmol, Rare arrhythmic events in computational models of cardiac tissue but reducing I K1 density by 50% increases the probability to~0.5. Thus I K1 density plays a critical role in modulating the Ca 2+ load at which arrhythmia-triggering events occur.

Probabilistic triggered activity in a paced fiber of myocytes
In the previous section, it was shown that triggered APs occur with a probability that depends on Ca 2+ load and that varying I K1 density shifts the Ca 2+ load dependence. Simulations were performed to test whether the model could produce probabilistic triggered activity in a 1D fiber of myocytes during pacing, where electrotonic loading of any cell by adjacent cells becomes important. Fig 7A shows the membrane potential of a fiber paced at 0.5 Hz under conditions similar to those in Fig 3 with elevated [Na + ] i set to 19 mM, as may be observed in heart failure [53]. Note also that [Na + ] sensed by membrane channels during contraction may be higher than [Na + ] i in the bulk cytosol due to local buildup of Na + imported by NCX [54].
To reflect a state of pathological remodeling, I K1 density and gap junction conductance were each reduced by 50% as observed in HF [49,55]. To limit boundary effects, [Na + ] i was set to 10 mM in the outer twenty-four cells near each end of the cable in order to prevent DAD generation. DADs, appearing as faint bands of depolarization between the paced beats, and reach V max values between approximately -70 and -60 mV. This is consistent with experimental observations of synchronized spontaneous Ca 2+ release causing DADs in intact heart following rapid pacing [9].
The model exhibited considerable variability in V max along the length of the fiber. To investigate the source of this variability, the state of the model at time t restart immediately after the third paced beat was recorded. A set of three independent realizations were run, each starting from this same model initial condition but with simulations performed using different initial pseudorandom number generator seeds. These realizations therefore differ in the particular pattern of LCC and RyR channel gating. Fig 7B shows: (i) the DAD from the original simulation in Fig 7A at 4.9 s; and (ii)-(iv) the three additional simulations initialized at t restart . Each simulation exhibits substantial differences in DAD amplitude along the fiber. In panel (iv), a spontaneous traveling action potential wave is observed. The wave originates from a region of locally high DAD amplitude near cell #50 and emanates in either direction along the fiber.
These results illustrate that arrhythmic events can occur probabilistically in tissue. By restarting the simulations immediately prior to the DADs, it became evident that the variability in DAD amplitude was due to variations in the stochastic events that occurred in this brief time window. Taken together with the results from Figs 5 and 6, the variability in DAD amplitude in the fiber can therefore be attributed to the underlying randomness of Ca 2+ wave dynamics and thus stochastic RyR gating.

Roles of Ca 2+ loading, I K1 density, and g gap in fiber DADs
Conditions reflecting pathological remodeling influenced the distribution of V max in the fiber.  Fig 6D). This reduction in variability is due to electrotonic coupling through inter-cellular gap junctional conductance, g gap , which attenuates spatial gradients in V. Fig 8B shows similar simulations where I K1 density was reduced by 50%. This results in greater fluctuations of V max over a range of 8.0 mV (SD 1.5 mV) compared to baseline. This is consistent with the increase in V max SD from 1.7 to 6.1 mV when I K1 density was decreased to 50% in isolated cells (see Fig 6D) and also reflects a substantial reduction in V max variability due to electrotonic coupling in the fiber. Reducing g gap by 50% in addition to the 50% I K1 reduction further amplifies V max spatial fluctuations, resulting in a range of 10.2 mV (SD 2.2 mV, Fig 8C). The increase in variability arises from the reduced electrotonic load experienced by each cell. These results demonstrate how perturbations of I K1 and g gap increase the likelihood of observing large DADs.

Overview of method for estimating rare event probabilities
The DAD results presented thus far indicate that spatial fluctuations in V max can result in a triggered beat emanating from a cluster of cells in which V max exceeds the threshold voltage (-55 mV). For a fiber of a given length, the probability of this event depends upon both the mean V max and the likelihood of a deviation from the mean with sufficient amplitude to reach threshold. Under conditions where the mean V max is far below threshold, however, it is unclear whether it would be plausible to observe a sufficiently large fluctuation that causes a triggered beat. We sought to characterize the likelihood of such events by estimating the upper tail of the V max distribution.
In Fig 5, we showed that V max was strongly correlated with J max in DAD simulations of isolated cells. In tissue, however, electrotonic coupling attenuates DAD amplitude by drawing current to neighboring cells via gap junctions. However, if spontaneous Ca 2+ release occurs synchronously in a cluster of cells, the potential gradient between cells is smaller, thus reducing the gap junction current and increasing DAD amplitude. Therefore the V max of a given myocyte in tissue depends on the relative timing and amplitude of Ca 2+ release in nearby cells.
We hypothesize that extreme deviations in V max are caused by rare spatially clustered synchronized Ca 2+ release events. However, estimating the probability of extremely rare events (e.g. 1 in 10 6 ) by direct simulation is computationally prohibitive using the full biophysical model since it would require potentially millions of simulations. Therefore, a method was developed for estimating the probability of such events using the output from a limited set of simulations.
As shown in the previous sections, stochastic RyR activity gives rise to random Ca 2+ wave patterns and thus variable profiles of J RyR . Therefore J RyR is a stochastic process that is dependent on complex microscopic events, which are computationally intensive to sample. The method leverages the fact that J RyR has a relatively simple distribution by resampling J RyR profiles from a set of simulations that is sufficiently large to approximate the distribution of J RyR . Multiple realizations of release events (i.e., the J RyR values from each separate simulation of the full cable model) are generated using a resampling method in which cell positions along the cable are shuffled to produce independent, distinct realizations. The key to this approach is the fact that spontaneous Ca 2+ release events are decoupled across the cells in the fiber (see Discussion). In principle, in using this approach one could develop a modified tissue model in which J RyR in each cell is fixed to one of the resampled J RyR profiles rather than being simulated, thus greatly reducing the computational burden. However, this would still require immense computational power to simulate the remaining differential equations when estimating the probability of extreme events. The second part of the method further increases computational efficiency by fitting a linear spatial-averaging model to predict V max. from J RyR alone. This permits the rapid estimation of the greatest V max in a 1D tissue model in millions of simulations and thus can estimate the probability that the cable will reach a proarrhythmic threshold potential.
The following sections describe this method and validate its accuracy for estimating DAD probabilities. We then apply this new technique to show how myocyte properties within a fiber, using I K1 density and g gap as examples, affect the likelihood of rare large amplitude DADs. ] i was initialized to 300 nM, and the right column illustrates the steps used in the filtering method. The first step in the method is to apply a uniform averaging filter to J RyR at each point in time to obtain a spatially smoothed profile

Filtering method for estimating V max from J RyR
where x refers to cell index, t refers to time, and W is the width of the filter (an odd integer). The rational for doing this is that membrane potential fluctuations of the cell at position x is influenced not only by the way its own complement of NCX responding to the local J RyR , but is also influenced by membrane potential fluctuations produced in response to J RyR and NCX activity in neighboring cells. The filter width W will therefore depend on the strength of gap junction coupling. For each cell, the maximum J f RyR value over all time was computed as The value of W maximizing the correlation coefficient between V max and J f max was calculated at each different value of gap junction coupling conductance used. J f max was then normalized to obtain estimates of V max using the formula where μ V and σ V are estimates of the mean and SD respectively of V max , and μ J and σ J are estimates of the mean and SD respectively of J f max . These estimates are calculated from a single full 1D tissue simulation containing~450 cells. Note how in the example of Fig 9A the estimated voltage profile V f closely resembles that of V. This approach therefore leverages the strong linear correlation between J RyR and V max (Fig 5B and 5C) and the statistical independence of Ca 2+ release events to dramatically reduce the number of simulations needed to estimate the probability distribution of membrane potential fluctuations.
The filtering method was applied to simulations of DADs with baseline conditions, with 50% I K1 , and with both 50% I K1 and 50% g gap (Fig 9B). In the non-baseline conditions, initial [Ca 2+ ] i was adjusted from 300 to 220 nM so that μ V would be approximately equal to that of the baseline. The width W of the smoothing filter is dependent on the fiber model parameters and therefore was optimized separately for each case. Table 1 shows the parameter fits for each condition. The resulting fits of V f max to V max have high correlation coefficient values (r 2 ! 0.90). All values of μ V fall within a narrow range from -77.0 mV to -79.1 mV, while the values of μ J in the two non-baseline conditions are less than half of the baseline value due to their lower Ca 2+ loads.
The increase in V max variability in the two pathological conditions is reflected in the parameters of the filtering model. It can be shown that the quantity S V = (σ V /σ J )/W scales with the Illustration of the filtering method for an example fiber simulation. A uniformly-weighted spatial smoothing filter of width W was applied to the spatiotemporal J RyR profile. The maximum values of the filtered profile were then normalized to obtain an estimate of the voltage profile (see Eq 6). (b) The filtering method was applied to three fiber simulations using the conditions indicated above each plot. The filtering method estimate (red) accurately reproduced V max from the original simulation (black). Parameters are listed in Table 1 SD of V f max (see S1 Text). Recall that 50% I K1 reduction increased variability of V in the cell model (Fig 6). For this condition in the fiber, S V is 28% larger compared to baseline, primarily reflecting the higher value of σ V . Imposing 50% g gap results in a filter width of only 27 cells compared to 43 and 49 in the other cases due to local decoupling of cells when gap junction conductance is decreased. This resulted in an S V that is 87% larger than baseline and 46% larger than with 50% I K1 alone.
These results show that the filtering method accurately estimates fluctuations in V max based on the J RyR profile. Moreover, the new empirical relationships derived here between J RyR and V max provide an approach to yield new insight and intuition into how changes in cellular properties and environment, such as reductions in outward currents and electrotonic coupling analyzed here, alter and in this case enhance spatial fluctuations in V max .

Resampling method for estimating rare event probabilities
The filtering method enables studies of properties of V max without the need to simulate the full biophysical model. Studying the statistical properties of V max requires the generation of large numbers of realizations. To do this, independent samples of V max were generated by shuffling cell positions in the J RyR profile, as illustrated in Fig 10A. Note that boundary effects the outer 24 cells on either end of the fiber are initialized to normal SR Ca 2+ loads and are not included in the shuffle. The filter method was then applied to the shuffled J RyR profile to estimate V max . An example V max profile obtained using this method is shown on the right. Note that the fluctuations about the mean μ V are qualitatively similar to those of the original simulation. In Fig  7, it was shown how a triggered propagating ectopic beat originated from a region of the fiber that experienced extreme DADs. Unfortunately it is impossible to perform the millions of cable simulations using the full biophysical model that are needed to study the probability of generating ectopic beats as model parameters are varied. However, our shuffling and filtering method does enable us to study the probability of a surrogate event of interest for each fiber simulation, with the surrogate event defined as the greatest value of V max realized along the fiber, referred to here as V peak .
The shuffling of cell positions makes two assumptions. The first is that the J RyR profiles of the cells are decoupled stochastic processes. The second assumption is that the fiber contains a sufficient number of cells such that the true distribution of J RyR is well represented by the collection obtained from a single simulation.
Previous work has shown that sub-threshold membrane depolarization induces Ca 2+ sparks and waves [56]. It is unclear whether local depolarizations in the membrane potential caused by spontaneous release affects release in neighboring cells, which would invalidate the independence assumption. This was tested by computing the peak amplitude and peak time of [Ca 2+ ] i in the fiber simulation represented by the red trace in Fig 8C. The absolute difference in the peak amplitude and peak time were then computed for all adjacent cell pairs and all cell pairs separated by a 50-cell distance along the fiber. Assuming that Ca 2+ release in 50 th neighbors are decoupled, if release in adjacent cells were coupled then one would expect the distributions of peak amplitude and timing to differ from those of the 50 th neighbor pairs. In particular, if local Ca 2+ release events were more synchronized, one would expect the difference in timing to be smaller. However, no substantial differences in Ca 2+ release peak or timing were observed (S4 Fig), as these distributions also were not significantly different according to non-parametric Kruskal-Wallis tests (p = 0.77, 0.66). Therefore the coupling effect of adjacent cells does not substantially affect spontaneous Ca 2+ release, which can thus be treated as a decoupled process in each cell in the fiber.
Validation of these assumptions also requires that the distribution of V peak generated using this method match that obtained when using the detailed biophysical model. This distribution was estimated from 1,000 independent realizations. Simulating the~500-cell fiber with 25,000 release sites in each cell is computationally prohibitive. For the purposes of the validation, the number of release sites was reduced to 2,500 and the fiber length to 96 cells. The values of μ V and σ V were computed using the V max values from all cells in all fibers. Fig 10B compares the true distribution of V peak computed using the full biophysical model to that obtained from applying only the filtering method to the simulated J RyR profiles (without shuffling). While this distribution exhibits a very small bias of~-0.05 mV, it is not significantly different from the true distribution on the basis of a Kruskal-Wallis test after correcting for the bias by subtracting their means (p = 0.89). The shuffling method was then validated by resampling from the J RyR s in all 1,000 fibers to produce 1,000 96-cell fiber realizations and computing V peak with the filtering method. While the resulting distribution also exhibits a very small bias of~-0.08 mV, it is not significantly different from the true distribution after adjusting the means (p = 0.69) (Fig 10B). The upper tails of the distributions containing the extreme DADs of interest are also similar (Fig 10C). To validate the second assumption that the sampling population of J RyR is sufficiently large, these tests were repeated using a subset of 5 of the 1,000 fibers to compute μ V , σ V , and the population of resampled J RyR 's (bias +0.05 mV, p = 0.68). Therefore, the method is accurate using a total of~500 simulated cells. In summary, the method outlined here is a computationally efficient approach that permits the rapid estimation of the V peak distribution using output from a single 500-cell fiber simulation.

Prediction of rare events
Fig 10D plots the probability that V peak -μ V exceeds a given potential estimated from an ensemble of 10 6 realizations. Reducing I K1 by 50% shifts the distribution tail to greater amplitude DADs compared to the baseline model. The probability of observing a V peak 3 mV greater than μ V (-77.1 mV) is 9.9 × 10 −5 in the baseline model compared to 0.096 with 50% reduction in I K1 , a nearly 1000-fold increase. The most extreme event observed with 50% I K1 was 6.1 mV greater than μ V compared to 3.4 mV in baseline, corresponding to DAD amplitudes 41% and 21% above average, respectively. Reducing g gap by 50% further increases the likelihood of large-amplitude DADs, with the probability of a 3 mV DAD reaching 0.50, a 5000-fold increase over baseline, and the most extreme event observed at 7.5 mV, corresponding to a DAD amplitude 53% higher than average. These results demonstrate how changes ion channel expression or function can alter the probability of supra-threshold membrane potential fluctuations. In this case, reductions in I K1 and g gap considerably widen the distribution of V peak and increase the probability of occurrence of larger V peak values by multiple orders of magnitude.
To illustrate the nature of these rare events, the realizations exhibiting the greatest V peak values were examined in each condition. Fig 11A plots V max -μ V . Because V max at any position is determined by an average of the Ca 2+ release profiles over a window of cells, each rare depolarization occurs at a cluster of cells whose width corresponds to the window size. Fig 11B shows the underlying J RyR profiles and the filter window centered on the V peak . In each case, the extreme event occurs at a cluster of cells within the window where J RyR tends to be greater and more synchronized than in the rest of the fiber. This result demonstrate that rare (occurring 1 in 10 6 beats in a 500-cell fiber) proarrhythmic events result from spontaneous Ca 2+ release exhibiting high (1) flux magnitude, (2) synchronization, and (3) spatial clustering. None of the 999,999 other simulations result in such an extreme event due to a lack of one or more of these properties. Note that because the window size is smaller in the case of 50% g gap , fewer cells need to fulfill these three criteria and thus the probability of extreme events is much higher.

Discussion
In this study we have presented a biophysically detailed stochastic computational model of the ventricular cardiac myocyte describing spatial Ca 2+ diffusion between release sites and incorporated it into a tissue-scale model to study the mechanisms and statistical properties of DADs. Loading of SR Ca 2+ is known to cause spontaneous Ca 2+ release [6,47,57]. In this model, Ca 2+ waves were generated when passive diffusion of Ca 2+ between release sites caused Ca 2+ sparks to propagate across the cell. In agreement with experimental studies, RyR sensitivity modulated the threshold SR Ca 2+ load at which this instability arises [58]. Rather than explicitly modeling RyR regulation mechanisms such as interaction with calmodulin [59], phosphorylation by activated CaMKII [60] and PKA [61], allosteric channel decoupling due to PKA-dependent dissociation of FKBP12.6 [62], and oxidation by reactive oxygen species [33,41], RyR opening rate was scaled such that the model reproduced experimentally observed triggered APs at 1 Hz pacing [33]. [Ca 2+ ] JSR -dependent regulation of the RyRs increases their sensitivity to cytosolic [Ca 2+ ] [63,64], but its role in spontaneous Ca 2+ release is controversial. In our model of RyR gating, it was assumed that this mechanism has little effect on RyR open probability when [Ca 2+ ] JSR < 1 mM [65,66] and therefore played a negligible role in dynamically regulating the RyRs in this study.
In our investigation, we found that pacing the baseline model at 1 Hz with 50% I K1 inhibition (in the absence of β-adrenergic stimulation) caused increased APD and cellular Ca 2+ loading over time. This resulted in greater inward NCX current during the repolarization phase, thus further increasing APD and Ca 2+ loading. This positive-feedback loop led to EADs and ultimately failure to repolarize. It is possible that the saturating behavior of LCC Ca 2+ -dependent activation in this model prevents sufficient LCC inhibition required for repolarization in the presence of elevated SR Ca 2+ 3). Elevated diastolic [Ca 2+ ] i increased RyR opening rate and thus perpetuated Ca 2+ wave formation at lower SR Ca 2+ loads. It also caused more rapid loading of the SR to induce overload. Cellular Ca 2+ loading also increased the amplitude of DADs due to the greater number and concurrence of Ca 2+ wave nucleation sites, which is in agreement with experimental studies [9,10]. Consistent with our results, Wasserstrom et al. reported that cellular Ca 2+ loading reduced DAD delay and increased synchrony associated this with greater likelihood of triggered activity [10].
The formation of Ca 2+ sparks and waves may be influenced by additional factors not considered in this study. These include intrinsic properties such as release site heterogeneity, including dyad geometry, RyR cluster morphology [19,68], local SR connectivity [69], and release site spacing [70,71]. Additionally, cell-to-cell heterogeneity would also affect DADs in tissue. Further work is needed to quantify the contributions of these factors to DAD variability.
Triggered events were investigated using a fiber model. Experimental study of the nature of DAD-induced arrhythmias is difficult due the limited temporal and spatial resolutions of live multiplex fluorescence imaging of tissue preparations. Xie et al. used a computational tissue model to study how the critical size of an isolated cluster of cells exhibiting identical spontaneous Ca 2+ release events required to trigger an AP depended on a variety of tissue geometries and was reduced under pathophysiological conditions [72]. They further showed that DAD amplitude and the critical cell mass depended sensitively on the amplitude of the spontaneous [Ca 2+ ] i transient. Here we have investigated stochastic variability in DAD amplitude in a long fiber of Ca 2+ -overloaded cells. We showed that large DADs occurred probabilistically due to random patterns of RyR gating and the formation of Ca 2+ waves that gave rise to high-amplitude, synchronized Ca 2+ release flux in a cluster of cells. Such events trigger arrhythmias by causing ectopic beats, inducing regional conduction block via Na + channel inactivation [73], and increasing variability of repolarization [46]. Heterogeneity of cell types and intercellular variability of Ca 2+ handling may also play an important role in determining the likelihood and location of triggered foci [10] but exploration of these factors remains beyond the scope of this study. Furthermore, the potential effects of pro-arrhythmic beat-to-beat AP variability [74] or EADs [75], which could affect Ca 2+ loading and DAD timing, were not examined.
The role of I K1 , which acts to stabilize the resting membrane potential, affected the DAD distribution in isolated cells. Loss of I K1 function has been associated with arrhythmogenesis in diseases such as heart failure [49], Andersen's syndrome [50], and long QT syndrome [51]. Maruyama et al. showed that the Ca 2+ -membrane voltage gain, defined as the ratio of DAD amplitude to spontaneous [Ca 2+ ] i transient amplitude, increased following I K1 suppression [76]. Consistent with this finding, the model exhibited a considerable increase in DAD amplitude in both isolated cells and in the fiber when I K1 density was reduced by 50%. Most notably, I K1 suppression increased DAD amplitude variability, dramatically increasing (by 1000-fold) the probability of occurrence of potentially arrhythmogenic large-amplitude DADs. Note that the model does not include Ca 2+ -dependent regulation of I K1 , as this mechanism has been controversial [77][78][79]. While a recent study by Nagy et al. [80] suggests elevating [Ca 2+ ] i at a fixed concentration upregulates I K1 via an increase in CaMKII activity, there remains insufficient data to constrain a model of this mechanism. Future models could be informed by further studying dynamic changes in CaMKII activity and I K1 in response to elevated [Ca 2+ ] i . Reduction of gap junction conductance, another pathological feature of diseases such as HF [55], further increased the variability of DAD amplitude in the fiber by reducing the spatial scale of electrotonic coupling. These findings are consistent with a previous modeling study that showed that fewer contiguous cells are required to exhibit DADs to produce a triggered beat under such conditions [72].
The relationship between I K1 and/or g gap to the probability of occurrence of extreme DADs demonstrated here is but one of many possible examples of how alterations of cellular mechanisms or environment can modulate propensity for arrhythmia. For example, HF also causes extensive remodeling of the TTs and dyads, leading to reduced CICR efficiency [81,82]. The model can be used to directly probe the effects of corresponding perturbations such as the number of release sites, RyR cluster size, and dyad volume on intracellular Ca 2+ cycling, the formation of spontaneous Ca 2+ sparks and waves, and rare arrhythmic events. Thus the novel method presented here provides a general framework for prediction of arrhythmia likelihoods in response to specific aspects of disease-related physiological remodeling.
A significant contribution of this work is that the emergence of sudden arrhythmias can be causally linked to stochastic molecular events. A computationally efficient method was developed to estimate the probability of extreme DADs. An important assumption of this method is that the spontaneous Ca 2+ release events in neighboring cells are decoupled. It is known that membrane depolarization increases the frequency of Ca 2+ waves by reducing NCX-mediated Ca 2+ efflux and thus promoting Ca 2+ waves due to increased intracellular [Ca 2+ ] [56]. Therefore, a release event in one cell that causes a local increase in V may hasten spontaneous Ca 2+ release in its neighbors. Validation of our method suggests that this coupling phenomenon has negligible impact on the V max distribution because the effect is weak compared to the intrinsic variability of spontaneous Ca 2+ release. Note that the model does not account for gap junction mediated intercellular Ca 2+ diffusion, resulting in cell-to-cell transfer of Ca 2+ waves, though the prevalence of such events remains unclear [9,[83][84][85][86][87].
There are two important conclusions that emerge from using this method. First, variability of the inward current due to stochastic RyR gating causes random patterns of Ca 2+ wave dynamics and results in substantial DAD variability at the tissue scale, particularly in the pathological states tested where I K1 and g gap were reduced. Second, while one can imagine a case where a contiguous cluster of cells exhibit large synchronized spontaneous Ca 2+ release, the probability distribution of such events has not been well characterized. In a 496-cell fiber with reduced I K1 and gap junction coupling, the largest DAD amplitude observed in an ensemble of 10 6 realizations had amplitude~50% greater than the mean DAD amplitude. For such a fiber paced at 1 Hz and exhibiting a DAD after every beat, one could expect to observe such an event approximately once every 11 days. Thus extreme DADs, while quite rare, are still possible over relevant time frames. Further work is needed to estimate the probability of such events in whole heart, as 3D tissue is~1-2 orders of magnitude less likely to exhibit triggered beats due to the increased electrotonic coupling [72] but also contain a much greater number of cells than tested here. Nevertheless, the results presented here suggest that variability due to stochastic molecular events play a large role in the initiation of cardiac arrhythmias and sudden cardiac death.  Fig 8C, a 496-cell fiber with 50% I K1 and 50% g gap . The absolute difference in peak [Ca 2+ ] i was computed for each pair of adjacent (dC 1 ) and 50 th neighbor (dC 50 ) cells. The absolute difference in the times of the peaks was also computed (dT 1 , dT 50 ). Histograms and QQ plots of the differences in (a) time of peak and (b) peak [Ca 2+ ] i . Both QQ plots exhibit linear trends, implying that the distributions are indeed alike. Therefore, the timing and amplitude of spontaneous Ca 2+ release of adjacent cells do not differ substantially from those of distant cells. (EPS) S1 Equations. Release site Ca 2+ transport equations. (DOCX) S1 Table. Release site Ca 2+ transport parameters.