Novel Ca²⁺ wave mechanisms in cardiac myocytes revealed by multiscale Ca²⁺ release model

The internal structure of cardiac muscle cells (cardiomyocytes) is well adapted for transducing surface action potentials to the release of Ca²⁺ from intracellular stores and consequent myofilament contraction (excitation–contraction coupling). Ventricular myocytes are roughly cylindrical with a diameter of ∼20 μm and a length of 100–150 μm. They are delineated by a sarcolemma that forms deep invaginations into the cell interior called transverse tubules (T-tubules). Electron microscope and super-resolution optical images reveal that these T-tubules weave around the myofilaments to form transverse planar lattice structures every 1.7–1.9 μm along the cell (Bers, 2001). They also show specialized intracellular synapses (dyad junctions) between the T-tubules and the sarcoplasmic reticulum (SR, the intracellular Ca²⁺ store). The dyad junction is comprised of a 14-nm cytoplasmic cleft between the membranes of the T-tubule and a specialized subcompartment of the SR called the terminal SR (TSR), which is joined via a nexus to a longitudinal SR (LSR) network running along and around the myofilaments (Brochet et al., 2005). The dyad junctions occur along the T-tubule with a mean edge-to-edge spacing of ∼150 nm (Hou et al., 2015).

During excitation–contraction coupling, the action potential opens voltage-dependent, L-type Ca²⁺ channels in the sarcolemma and T-tubules, which causes Ca²⁺ to enter the dyad cleft. This releases Ca²⁺ from the TSR via clusters of Ca²⁺-activated ryanodine receptor Ca²⁺ channels (type 2 ryanodine receptors, RyR2s) in the TSR membrane (Bers, 2002). The subsequent release of calcium from the TSR further increases [Ca²⁺] in the cleft, which strongly reinforces local RyR2 activation, a process called Ca²⁺-induced Ca²⁺ release (CICR) (Fabiato, 1983). The combined action of Ca²⁺ release from many such dyad junctions produces a global increase in cytoplasmic [Ca²⁺], which is the trigger for myocyte contraction (systole). The resulting depletion of the SR terminates Ca²⁺ release, and Ca²⁺ is either sequestered back into the SR by the ATP-powered Ca²⁺ pumps (SERCA2a) in the LSR membrane or extruded from the cell by the Na/Ca exchanger (NCX) in the sarcolemma, allowing muscle relaxation (diastole) (Bers, 2002). This results in a Ca²⁺ transient in which cytoplasmic [Ca²⁺] increases from ∼0.1 to ∼1 μM and where LSR [Ca²⁺] decreases from ∼1 to 0.3 mM. In ventricular myocytes, the action potential normally drives the cyclic release and uptake of Ca²⁺ by the SR (Bers, 2002). However, pathological overactivation of Ca²⁺ release by altered RyR2 function or Ca²⁺ store overload can lead to propagating waves of SR Ca²⁺ release (Berlin et al., 1989) that initiate independently of the action potential leading to heart arrhythmias and sudden cardiac death (Belevych et al., 2012; Blayney and Lai, 2009).

Integrating these cellular Ca²⁺ release phenomena with the well-known Ca²⁺ activation of the RyR2 that underlies CICR has proved to be an ongoing challenge. CICR provides a strong positive feedback mechanism that, once initiated, should drive near-total Ca²⁺ depletion of the SR, independently of the sarcolemma membrane potential. To explain the observed graded response of Ca²⁺ release to sarcolemma membrane potential and the incomplete depletion of the SR, Stern (1992) proposed the “local control” model whereby CICR is moderated by the spatial separation of discrete Ca²⁺ release sites (see review by Cannell and Kong [2012]). Indeed, such a moderating process was demonstrated by Cheng et al. (1993) who observed discrete Ca²⁺ release events as Ca²⁺ sparks by virtue of their appearance as flashes of light (∼2 μm diameter) in the presence of fluorescent Ca²⁺ indicators. Ca²⁺ sparks are localized Ca²⁺ release events that occur at single-dyad junctions and are thought to be the basic quanta of global Ca²⁺ release phenomena such as Ca²⁺ waves and transients (Cheng et al., 1996). Their fluorescence time course has an exponential rising phase lasting 10–20 ms, which roughly corresponds to the time that the RyR2s are open and releasing Ca²⁺ (Kong et al., 2013). When the RyR2s close and Ca²⁺ release is terminated, the Ca²⁺ spark fluorescence declines to baseline within 50–100 ms. Although it was quickly understood that RyR2 opening during a Ca²⁺ spark is strongly reinforced by CICR between individual RyR2s in the dyad, it was not until much later that Ca²⁺ spark termination could be explained by a computational model incorporating experimentally measured properties of RyR2 gating and dyad architecture (Cannell et al., 2013; Laver et al., 2013). This computational model and a concurrent experimental investigation (Guo et al., 2012) both concluded that Ca²⁺ release is strongly reinforced by CICR between RyR2s in the dyad and cessation of release is due to decay of the CICR process by Ca²⁺ depletion of the TSR; a process dubbed as either induction decay (Cannell et al., 2013) or pernicious attrition (Guo et al., 2012).

Ca²⁺ sparks and waves also occur in chemically skinned cardiomyocytes that lack an effective sarcolemma (Lukyanenko and Gyorke, 1999), thus demonstrating that they emerge from properties of the SR and not the sarcolemma. Although it is assumed that Ca²⁺ waves are the summation of Ca²⁺ sparks that are coupled by CICR between neighboring dyads (Cheng et al., 1996), the details of this process remain unclear and the development of a convincing Ca²⁺ wave model has proved challenging (Izu et al., 2013). For example, RyR2 sensitivity to cytoplasmic Ca²⁺ (half-activation by 20–50 μM Ca²⁺ [Gyorke and Gyorke, 1998; Laver and Honen, 2008; Walweel et al., 2014]) is too low to explain Ca²⁺ wave propagation between dyads when wave amplitudes are typically only a few μM. Models that compensate by increasing RyR2 Ca²⁺ sensitivity generate unrealistically high spark frequencies and initiation of myriad wavelets rather than a single propagating wave (Izu et al., 2013). Other models solve this problem by incorporating store overload–induced Ca²⁺ release (SOICR) into their RyR2 gating models; a phenomenon in which SR Ca²⁺ release drastically increases when the SR [Ca²⁺] exceeds 0.5 mM (Jiang et al., 2004; Shannon et al., 2000a). These Ca²⁺ wave models frequently adopt a fourth-power (or more [Restrepo et al., 2008]) dependence for RyR2 sensitivity to SR [Ca²⁺] (Walker et al., 2017; Walker et al., 2014; Xie et al., 2019), which virtually amounts to a step change in RyR2 activity at a threshold SR [Ca²⁺]. Thus, RyR2s are made quiescent between waves when the SR load is low and yet very permissive to wave propagation when the SR load is high. The calcium release properties of RyR2s have been intensely studied and are very well characterized and demonstrate that the RyR2 can be activated by Ca²⁺ in both the SR and the cytoplasm (Laver, 2007; Sitsapesan and Williams, 1994). However, these studies reveal that the RyR2 has a smaller, more graded response to SR [Ca²⁺] (Gyorke and Gyorke, 1998; Laver and Honen, 2008). Over the SR [Ca²⁺] range 0.1–1 mM (in the presence of physiological luminal [Mg²⁺] of 1 mM), the RyR2 open probability increased only threefold (Walweel et al., 2014) and then only when cytoplasmic [Ca²⁺] was below 1 μM (Laver and Honen, 2008). At higher cytoplasmic concentrations, SR [Ca²⁺] had no effect on RyR2 activity. Thus, it is still not clear how Ca²⁺ signaling phenomena such as SOICR and Ca²⁺ waves arise in cells.

Here, we develop a model that successfully reproduces store load–dependent Ca²⁺ release and Ca²⁺ waves that use experimentally verified Ca²⁺-dependent rates of RyR2 opening and closing. The model is based on our previous induction-decay model of the dyad junction that explained the time course of Ca²⁺ sparks, Ca²⁺ blinks, and Ca²⁺ spark restitution (Cannell et al., 2013; Laver et al., 2013). We use parallel computing to calculate Ca²⁺ release in regions of the cell containing 10³–10⁵ dyad junctions that are distributed in array formations informed by super-resolution micrographs (Hou et al., 2015).

Details of the numerical methods are given by Vysma et al. (2023). Briefly, the model cell volume was divided into cubic voxels with sides of 0.25 μm containing either LSR or LSR plus a Ca²⁺-releasing dyad. Ca²⁺ diffusion and binding to Ca²⁺ buffers within and between voxels were calculated by solving Fick’s equation, which required the solution of a system of 62 ordinary differential equations (described in the Supplementary text at the bottom of the PDF.) using Diagonally Iterated fully Implicit Runge–Kutta (DIIRK). DIIRK offered a 10-fold speed improvement and removed out-of-memory exceptions for large models when used as a replacement for MATLAB’s ode15s solver (Vysma et al., 2023). Each voxel was given a parallel thread on either the CPU (two Intel Xeon Gold 5218 processors with a total of 64 hyper-threads) or GPU devices (four Nvidia Quadro RTX 6000 graphics cards). DIIRK uses a variable number of iterations depending on the difficulty of finding a solution. To efficiently use heterogeneous computation hardware (GPUs and CPUs), we developed a predictor function that estimates the number of iterations required for a voxel given the cytoplasmic calcium concentration. This heuristic predictor function was then used to dynamically assign voxels with higher expected computation costs to more powerful computation hardware. Model parameters, constants, and SR morphology were set within MATLAB, which then “just in time” compiles the model into C and OpenCL applications for each simulation. MATLAB collects the completed simulation data into a format suitable for presentation.

Since Ca²⁺ sparks and waves are similar in intact cells and chemically skinned cardiomyocytes (Lukyanenko and Gyorke, 1999), it is likely that the sarcolemma does not play a significant role in the initiation and propagation of Ca²⁺ waves. Therefore, we explored mechanisms of spark and wave generation in the simpler skinned cardiomyocyte model that excludes the sarcolemma. The dyad model is derived from a 3D spatial model that accurately simulated the activation and termination of Ca²⁺ sparks (Cannell et al., 2013; Laver et al., 2013) where Ca²⁺ activation arose from CICR between RyR2s within a dyad and termination arose from induction decay; a decay in CICR due to TSR depletion. The simplified dyad compartment model used here aims to capture both CICR and induction-decay mechanisms while being computationally feasible in a model cell containing many thousands of such dyads.

The dyad cleft is modeled by a 14-nm planar gap between the T-tubule and TSR (Fig. 1 A). The TSR is divided into four concentric interconnected pools. The RyR2 sits in a 5 × 5 square cluster with a nearest neighbor separation of 31 nm, which occupies the three inner pools. We expect that diffusion of Ca²⁺ within the RyR2 cluster is slowed by the tortuosity of diffusion paths around the RyR2 proteins. We calculated the fractional reduction in the diffusion rate (tortuosity factor = 0.5) by simulating diffusion in a 2D space containing a representative RyR2 cluster taken from Asghari et al. (2014) (Fig. S1).

The total Ca²⁺ flux through a RyR2 cluster is proportional to their unitary Ca²⁺ current, open probability, and their number in the cluster. The unitary current for a range of LSR [Ca²⁺] under physiological ionic conditions was measured by single-channel recording and is shown in Fig. S2. The RyR2 open probability in the dyad is calculated using a hybrid stochastic–deterministic method using experimentally determined Ca²⁺ dependencies of RyR2 opening and closing rates from sheep, human, and rat (Fig. S3 and Table S3). Ca²⁺ release is initiated by a RyR2 opening, modeled as a discrete-time stochastic process generated by the Monte Carlo technique at 20-μs intervals. Although any RyR2 in the cluster may open first, the initial current was assigned to the central, smallest pool, containing only 1 RyR, for CICR to rapidly (∼1 ms) activate the remaining RyR2s that are modeled using a deterministic process requiring the solution of differential equations.

The TSR portion of the dyad is connected via a nexus to the LSR (Fig. 1 B). Depletion of TSR and spark termination occurs when TSR refill via the nexus is slower than Ca²⁺ release through the RyR2 cluster. Also, the nexus controls the rate of TSR refill following a Ca²⁺ spark and hence the rate of spark restitution. In the absence of any details of nexus geometry, we estimated the nexus diffusion using an approach we adopted previously (Laver et al., 2013). The nexus Ca²⁺ diffusion constant was set as a fraction of the diffusion between voxels (LSR–TSR nexus factor, N_f), which was adjusted so that the rate of restitution of the spark amplitude following a spark at the same release site (see Fig. S4 B) was comparable to published values (Zima et al., 2008b).

The LSR is distributed uniformly throughout the cell, occupying 3.5% of the cytoplasmic volume (Bers, 2001) such that Ca²⁺ can diffuse between the LSR in neighboring voxels. The diffusion coefficient for Ca²⁺ in the LSR was chosen to match the fluorescence recovery after localized Ca²⁺ depletion of the LSR measured by Wu and Bers (2006) (Fig. S5). The model includes the Ca²⁺ transport properties of SERCA2a taken from Shannon et al. (2000b).

The cellular distribution of the Ca²⁺ release sites in the model is based on super-resolution micrographs of cardiac cells (see Fig. 1 of Hou et al. [2015]). Ca²⁺ release sites are distributed in a 1-μm planar grid at each Z-line, and Z-lines are spaced at 1.75-μm intervals along the entire cell (Fig. 1, C and D). RyR2 clusters contained a mean of 25 RyR2s, and the resulting RyR2 cluster density in the model is 4 μm⁻³, equating to 100 RyR2 μm⁻³ (estimates in the literature range from 18 to 130 μm⁻³ [Hayashi et al., 2009; Hou et al., 2015]). We considered two alternatives for dyad sizes in the model: (1) all dyads were of equal size and contained 25 RyR2s, and (2) dyads of equal size were randomly allocated a variable RyR2 cluster size ranging between 10 and 120, exponentially distributed with a mean of 25.

ChatGPT assisted in condensing the abstract.

Supplemental figures describe the experimental data underlying key model parameters (Figs. S1, S2, S3, and S5), model simulations of [Ca²⁺] in the dyad cleft, TSR and LSR during Ca²⁺ sparks and embers (Figs. S4 and S6), spatial and temporal distributions of Ca²⁺ sparks and embers during a Ca²⁺ wave (Fig. S7), z-t scans showing the effect of Ca²⁺ embers on wave initiation and propagation (Fig. S8), the effect of cytoplasmic Ca²⁺ buffering on Ca²⁺ waves (Fig. S9), parameter sensitivity analysis of the model (Fig. S10), and the effect of RyR2 gating kinetics on Ca²⁺ waves (Figs. S11 and S12). Table S1 shows diffusion and buffering coefficients. Table S2 shows structural parameters. Table S3 shows RyR2 Ca²⁺ activation parameters. Table S4 shows miscellaneous model parameters. Supplemental text at the end of the PDF provides further information on the experimental methods and their analysis.

The hybrid stochastic–deterministic method for calculating RyR2 activity during a single Ca²⁺ spark described here captures the main features of dyad [Ca²⁺] transients seen with the more computationally intensive, fully stochastic method used previously (Cannell et al., 2013; Laver et al., 2013). This, together with the DIIRK ODE solver and parallel computing approach, could carry out a single-spark simulation comparable to our previous studies with a 100-fold speed improvement. The hybrid method produces a slightly faster rise and decay in dyad [Ca²⁺] and a deeper faster depletion of the TSR (Fig. S4 A). With both methods, the initial rise in dyad [Ca²⁺] results from a regenerative phase of CICR between neighboring volume elements in the dyad that persists for 5 ms. TSR [Ca²⁺] declines to ∼10% of resting levels within 10 ms, which reduces the Ca²⁺ flux through RyR2s and hence dyad [Ca²⁺] to a point where RyR2s start to close, further reducing the Ca²⁺ flux in an induction-decay cycle that continues until Ca²⁺ release ceases all together after 30 ms. The TSR then slowly refills with Ca²⁺ that diffuses from the LSR via the nexus leading to Ca²⁺ spark restitution over the next few hundred ms (Fig. S4 B).

Ca²⁺ sparks generated by various RyR2 cluster sizes had similar morphology, but, as expected, occurred with frequencies in proportion to cluster number (Fig. 2 A). A representative simulated line-scan image of a Ca²⁺ spark cell is shown in Fig. 2 B for [Ca²⁺]. Morphology distributions and frequencies were gathered from those with dF/Fo amplitudes above 0.5. The amplitude of fluorescent transient signals varied from near zero up to dF/Fo ∼1.2. Smaller amplitude signals originated from Ca²⁺ sparks away from the scan line. Ca²⁺ spark frequency increased from 10 to 40 s⁻¹(100 μm)⁻¹ as LSR [Ca²⁺] increased from 0.4 to 0.75 mM (Fig. 2 C) primarily because of the luminal Ca²⁺ dependence of the RyR2 opening rate. Although reported Ca²⁺ spark frequencies vary considerably between studies, the relative dependence on LSR [Ca²⁺] is qualitatively like that seen by Zima et al. (2010). Histograms of Ca²⁺ spark amplitude, time-to-peak, width, and duration (LSR [Ca²⁺] = 0.6 mM, Fig. 2, D–G) when compared to experimental distributions (Picht et al., 2007) reveal similar spark amplitudes and durations but ∼2-fold faster time-to-peak and ∼1.5-fold lower width.

We investigated the effect of increasing LSR [Ca²⁺] on Ca²⁺ release during a single Ca²⁺ spark and found that there was a threshold LSR [Ca²⁺] of ∼0.95 mM at which induction decay became substantially delayed, leading to an extended Ca²⁺ release (Fig. 3 A and Fig. S6); a phenomenon dubbed Ca²⁺ embers (Zhou et al., 2003). Above this threshold, Ca²⁺ diffusion from the LSR to TSR was able to keep up with Ca²⁺ release through RyR2 and the TSR failed to deplete sufficiently to initiate the induction-decay cycle. Ca²⁺ release continued until LSR [Ca²⁺] decreased to ∼0.7 mM after 1 s. The threshold LSR [Ca²⁺] for Ca²⁺ embers decreased as the external [Ca²⁺] increased (Fig. 3 B), indicating that higher cytoplasmic [Ca²⁺] promoted Ca²⁺ embers.

The Ca²⁺ flux during individual sparks and embers exhibited a similar initial transient lasting ∼20 ms, but only embers had an extended plateau with ∼5% of the amplitude of the initial transient (Fig. S6 E, solid line). The cytoplasmic [Ca²⁺] increase that occurred at the nearest neighboring dyad reached 2 μM during the transient and 100 nM during the plateau (Fig. 3 C, peak 2). During the triggering of two adjacent dyads, the extended plateau fluxes had a summative effect on the cytoplasmic [Ca²⁺] increase (Fig. 3 C, peak 1) that was unlikely to occur with brief transients. We calculated the cumulative probability that an ember would trigger Ca²⁺ release from a neighboring 25 RyR2 cluster (Fig. 3 D). During the flux transient (common to both embers and sparks), this probability rose to 0.13 in the first 20 ms, and during the plateau flux, it increased substantially more over the next few hundred milliseconds.

To investigate how Ca²⁺ release at individual dyads might contribute to global Ca²⁺ release, we considered a permeabilized cell volume, 16 × 16 × 55 μm, like that shown in Fig. 1 D. The computation speed was substantially increased by relying on fourfold symmetry and doing calculations in only one cell quadrant with reflecting intracellular boundaries (i.e., no diffusion; see Fig. 4 A) plus a reflecting boundary at one end. Unless otherwise stated, we doubled the K_max of SERCA2a to 0.4 mMs⁻¹ from the commonly used value (see Shannon et al. [2000b]) to reduce wave periods and hence simulation durations. The model parameters are given in Tables S1, S2, S3, and S4. The computed fluorescence (⁠$∝$ [CaFluo-4], in the x–z plane) from spontaneous Ca²⁺ waves (propagating from right to left) from sheep/human RyR2 is shown in Fig. 4 B for the case of exponentially distributed RyR2 cluster sizes. Waves are shown at various time points, starting 200 ms before wave initiation. Ca²⁺ sparks are seen as cyan circular patches that appear mainly before the wave front.

Fig. 5, A and B, shows the [Ca²⁺] in the cytoplasm and LSR during waves generated from sheep/human in dyads with uniform RyR2 cluster sizes (black and red lines) or exponentially distributed cluster sizes (green lines) and that the waves were about the same. In the case of uniform clusters, [Ca²⁺] profiles are shown near the cell axis (black lines) and cell boundary (red lines). During the Ca²⁺ wave peak, Ca²⁺ diffuses out of the cell along its concentration gradient between the cell interior (8 μM) and the external bath (300 nM, Fig. 5 A, dashed line) so that at 1 μm from the cell boundary the wave peak is 3 μM. During the wave trough, Ca²⁺ diffuses into the cell along its concentration gradient between the external bath and the cell interior (50 nM) so that near the cell boundary the trough is 150 nM. Although the LSR is not directly connected with the external bath, it supports a transverse [Ca²⁺] gradient due to differences in SERCA2a-driven Ca²⁺ uptake (Fig. 5 B). During the trough, Ca²⁺ uptake into the LSR is larger, and Ca²⁺ loading is faster, near the cell surface because the SERCA2a substrate concentration there is higher. Ca²⁺ loading continued until LSR [Ca²⁺] reached ∼1 mM at which point the Ca²⁺ wave initiated and unloaded the LSR. Ca²⁺ waves appeared to initiate near the cell surface (see Fig. 4 B at 3.1 s) because Ca²⁺ loading of the LSR there is faster.

Ca²⁺ transients from waves generated by rat RyR2s are shown in Fig. 5, C and D. Rat RyR2s showed a threefold lower sensitivity to cytoplasmic Ca²⁺ compared with those from sheep/human (Fig. S3). These waves could only be sustained with reduced Ca²⁺ buffering (i.e., 30 μM Fluo-4 compared with 100 μM for sheep/human). Ca²⁺ waves in rat had higher LSR [Ca²⁺] threshold, larger cytoplasmic [Ca²⁺] amplitude, and longer periods than those in sheep/human. Waves had propagation velocities of 60 μm/s (sheep/human) and 100 μm/s (rat).

Fig. 6, A and B, compares the z-t line scans of Ca²⁺ waves generated from uniform RyR2 clusters (A) and exponentially distributed RyR2 cluster sizes (B). Embers were very dim in these z-t line scans due to the large dynamic range of brightness required to present the Ca²⁺ waves. The display of embers was enhanced in Fig. 6 C by plotting the logarithm of relative fluorescence, revealing that Ca²⁺ waves were preceded by Ca²⁺ embers that occurred both near and far from the wave origin (e.g., see arrows in Fig. 6 C). Waves generated by exponentially distributed cluster sizes were like uniform clusters except the former showed a twofold higher occurrence of embers ahead of the wave (Fig. S6 F). Variations in RyR2 cluster size led to a distribution in the threshold LSR [Ca²⁺] for triggering embers whereby outliers with low thresholds produced more embers further ahead of the wave.

We investigated the roles of Ca²⁺ sparks and embers in the initiation and propagation of a Ca²⁺ wave (the one simulated in Fig. 6 A) by plotting their locations at 100-ms intervals during the Ca²⁺ wave (Fig. S7, representative times shown in Fig. 7). Embers (red dots) were defined as Ca²⁺ release events that failed to terminate within 50 ms as opposed to Ca²⁺ sparks (blue dots). During the Ca²⁺ wave trough (Fig. 7, 2.4 s), Ca²⁺ sparks were uniformly distributed throughout the cell, whereas Ca²⁺ embers tended to locate the region of highest LSR [Ca²⁺] near the cell surface (x = y = 8 μm), especially along the corners (see Fig. S7 for LSR [Ca²⁺] distributions). At 2.7 s, embers become more numerous than sparks at z > 50 μm mainly due to increased LSR [Ca²⁺], while cytoplasmic [Ca²⁺] remains near baseline levels. At 3.0 s, cytoplasmic [Ca²⁺] at z = 53 μm increases above 1 μM, which should also decrease the LSR [Ca²⁺] threshold for embers by ∼10% (Fig. 3 B), providing an additional positive reinforcement for ember generation despite the falling LSR [Ca²⁺]. From 3.0 to 3.4 s, the cytoplasmic [Ca²⁺] distribution morphs into a 5-μM peak at the cell interior at z = 47 μm and depletion of the LSR returns the predominance of sparks over embers at higher z. From 3.4 s onward, CICR produces a propagating wave in cytoplasmic [Ca²⁺] where the Ca²⁺ embers are strongly stimulated at the wave peak and behind the wave front, and absent afterward (3.8 s).

We investigated the role of Ca²⁺ sparks and embers in wave initiation and propagation by calculating the Ca²⁺ fluxes from all the dyads in two x-y planes, one near the point of wave initiation at z = 46 μm (Fig. 8, A and C) and the other that experiences the propagating wave at z = 21 μm (Fig. 8 B). The total Ca²⁺ flux time courses in each x-y plane due to transient (from both sparks and embers) and plateau phases (embers only) are plotted separately in Fig. 8. Sparks occurred throughout the displayed interval (2–4 s) with a surge in spark frequency during the wave, whereas embers only started to occur at 2.7 s when LSR [Ca²⁺] reached 0.9 mM and ceased once LSR [Ca²⁺] had declined to 0.3 mM. Just before the wave begins (3.0–3.4 s, Fig. 8 A), the total flux from transients was about twice as high as the steady plateau flux. When the wave occurs (3.4–3.6 s), the transient flux increased to roughly 10 times higher than the plateau flux. During the wave, the plateau flux itself was about 50% higher compared with the period before the wave. Further from the wave’s point of origin (Fig. 8 B), the ratio of transient to plateau flux was like that near the origin. However, in this region, the plateau flux during the wave was four times greater than it was before the wave.

In Fig. 8 C, we did the same simulation up to the 3-s time point and subsequently suppressed plateau Ca²⁺ fluxes without significantly altering the frequency of initial transients by restricting diffusion between the LSR and TSR. This maneuver not only prevented wave initiation but also quenched an established Ca²⁺ wave within 200 ms despite the resulting increase in the SR load (Fig. S8), demonstrating that the Ca²⁺ flux plateau (i.e., embers) is a prerequisite for wave initiation and propagation.

The ability of Ca²⁺ release at one dyad to activate its neighbors (the essence of CICR) is dependent on the Ca²⁺ efflux from the dyad and the ability of Ca²⁺ to diffuse between dyads, which, in turn, is strongly influenced by Ca²⁺ buffering power. The buffering power of the cytoplasm is known to decrease with increasing cytoplasmic [Ca²⁺] (e.g., see a review by Smith and Eisner [2019]), and the buffering power of most Ca²⁺ buffers in our model decreases 10-fold between 100 nM and 1 μM and 100-fold at peak [Ca²⁺] (Fig. S9 A). We examined the role of this decreased buffering power in Ca²⁺ wave propagation by simulating Ca²⁺ release in the presence of a fictitious Ca²⁺ buffer with normal buffering power at 100 nM cytoplasmic [Ca²⁺] that remained constant up to 10 μM cytoplasmic [Ca²⁺] (Fig. S9 A, dashed line). Under these conditions, we could not sustain a propagating Ca²⁺ wave; instead, Ca²⁺ release appeared as a collection of stochastic, long-lasting release events (embers, Fig. S9, B and C). Therefore, the role of reduced cytoplasmic buffering power in CICR is important for the generation of a propagating wave.

Therefore, both Ca²⁺ embers and CICR are necessary to produce a Ca²⁺ wave in this model. Even though embers contributed only 30% to the total flux, they had a disproportionate ability to create a cumulative increase in cytoplasmic [Ca²⁺] (e.g., Fig. 3 C). Once cytoplasmic [Ca²⁺] nears 0.5 μM, CICR begins to recruit additional embers and sparks, providing the positive feedback in Ca²⁺ release needed to sustain the Ca²⁺ wave. Thus, our model predicts that embers are required to prime the cytoplasmic milieu for CICR as illustrated in Fig. 9 (the priming period is indicated by the black bar in Fig. 8). In the lead-up to a Ca²⁺ wave, Ca²⁺ uptake by the SERCA2a outpaces Ca²⁺ leak via Ca²⁺ sparks, resulting in Ca²⁺ loading of the LSR. At an LSR [Ca²⁺] threshold (in this case, 0.9 mM), embers start occurring that produce the increase in cytoplasmic [Ca²⁺] and decrease in cytoplasmic buffering power needed to initiate CICR. Moreover, Fig. S8 B shows that priming of the cytoplasm by embers throughout the cell is essential for sustaining CICR during the propagation of an established Ca²⁺ wave. At the wave peak, cytoplasmic [Ca²⁺] is high enough to drive the LSR [Ca²⁺] ember threshold down to 0.3 mM (Fig. 3 B), thus allowing increased ember frequency and plateau Ca²⁺ flux, even while LSR [Ca²⁺] decreases below the threshold for wave initiation. Once the Ca²⁺ wave is established, there is a more pronounced peak in the total plateau Ca²⁺ flux at peak cytoplasmic [Ca²⁺] (Fig. 8 B, blue, 3.7 s), indicating that cytoplasmic [Ca²⁺] is generating more embers. Thus, while embers initiate the wave, the established wave initiates more embers.

The process of wave termination is a negative feedback process that begins when flux due to transients drops 80% and plateaus drop 100% due to a depleting SR. Cytoplasmic [Ca²⁺] peaks and begins to fall as the flux balance shifts between Ca²⁺ release, diffusion into the surrounding cell and bath, and uptake by the LSR. As cytoplasmic [Ca²⁺] decreases, increased cytoplasmic buffering power decreases Ca²⁺ diffusion between neighboring dyads, which together with their decreased Ca²⁺ efflux attenuates CICR in a negative feedback loop that continues until cytoplasmic [Ca²⁺] returns to basal levels (Fig. 9). During this process, the LSR [Ca²⁺] threshold for embers returns to higher levels, which eliminates embers until the LSR refills.

Generally, we found that 2 μM was the smallest cytoplasmic [Ca²⁺] amplitude that could sustain Ca²⁺ waves with sheep and human RyR2s. For example, decreasing Ca²⁺ store capacity by decreasing its main Ca²⁺ buffer concentration (TSR [CSQ]) from 30 to 20 mM (Fig. 10), rather than merely decreasing wave amplitude as extrapolated from our sensitivity analysis (see below), caused the cell to exhibit wavelets in the case of uniform RyR2 clusters (Fig. 10 C) or to settle to a steady state for nonuniform clusters after two brief oscillations in the first 3 s in which cytoplasmic [Ca²⁺] = 0.3 μM and LSR [Ca²⁺] = 0.8 mM (Fig. 10, A and B, red line; and Fig. 10 D). Fig. 10 E shows the contributions of Ca²⁺ flux transients and plateaus to Ca²⁺ release for the case of nonuniform RyR2 clusters. As LSR [Ca²⁺] reaches a threshold, embers appear as seen by a rise in the plateau flux (blue), but in this case, it does not produce the surge in transient fluxes (red) normally seen with CICR. Hence, the Ca²⁺ flux is insufficient to create enough CICR or to reduce LSR [Ca²⁺] to below the threshold for embers so that the plateau current does not terminate. Thus, our model predicts that in cases of reduced Ca²⁺ release, sustained Ca²⁺ waves fail to occur because both positive and negative feedback mechanisms are insufficient to maintain system oscillations.

Due to the enormous computational size of the model, the sensitivity analysis used a one-factor-at-a-time method to determine a ratio of the relative change in the model output with respect to relative changes in each model input. Fig. 11 shows the first of two analyses where the model outputs were the Ca²⁺ wave period, peak cytoplasmic [Ca²⁺] during a wave (cyt_max), the minimum cytoplasmic [Ca²⁺] between waves (cyt_min), the peak LSR [Ca²⁺] just before wave initiation (LSR_max), and the minimum LSR [Ca²⁺] immediately following wave termination (LSR_min). Model input parameters were grouped in the following categories: (1) Ca²⁺ diffusion, including Ca²⁺ diffusion constant, its relative permeability across the perforated external membrane, and [Ca²⁺] in the external bath; (2) the dimensions of the model compartments, RyR2 separation, and mean cluster size; (3) RyR2 unitary current and Ca²⁺ activation kinetics; (4) SERCA2a activity; and (5) Ca²⁺ buffering. The resulting heat map matrix of model sensitivities in Fig. 11 shows that wave period, LSR_max, and cyt_max were most sensitive to the size of cell compartments and RyR2 distribution and moderately sensitive to SERCA2a activity, RyR2 gating parameters that determine CICR within and between dyads, and Ca²⁺ buffering in the cytoplasm by ATP and the TSR by CSQ.

In the second sensitivity analysis, the model output was the threshold LSR [Ca²⁺] for ember formation at a single release site using the same model inputs as our first analysis. We found a positive correlation between the effect of LSR [Ca²⁺] on ember formation and wave initiation in 22 out of the 31 model parameters inspected (Fig. S10). This is consistent with the transition of Ca²⁺ sparks to Ca²⁺ embers being a requirement for Ca²⁺ wave initiation seen in Fig. 8. Thus, our analyses support two governing principles of Ca²⁺ wave initiation: (1) Ca²⁺ embers occur when [Ca²⁺] (hence CICR) in the dyad cleft remains too high to initiate the induction-decay cycle that causes Ca²⁺ spark termination, and (2) Ca²⁺ waves initiate at the LSR [Ca²⁺] threshold where Ca²⁺ sparks transform into Ca²⁺ embers. In light of these principles, we explain (below) how the model parameters influence Ca²⁺ wave properties.

LSR_max is effectively the LSR [Ca²⁺] that can sustain sufficient CICR in the dyad cleft to prevent the induction-decay cycle. Thus, LSR_max is lowered by factors that increase CICR in the dyad cleft such as increasing RyR2 sensitivity to cytoplasmic Ca²⁺ or increasing RyR2 unitary current. We note that the RyR2 closing rate has a threefold stronger effect on LSR_max and ember threshold than the opening rate (Fig. S10, c.f. #13 and #14). LSR_max is also lowered by decreasing Ca²⁺ efflux from the dyad cleft via more constrained dyad geometry or decreasing Ca²⁺ buffering by ATP, T-tubules, and SR lipids, which also tends to maintain a higher dyad [Ca²⁺]. (The negative surface potential on the T-tubule strongly attracts Ca²⁺, increasing the Ca²⁺ buffering capacity of T-tubule lipids.) Increasing Ca²⁺ diffusion between the LSR and TSR or increasing TSR Ca²⁺ buffering by CSQ hinders TSR depletion, thus increasing [Ca²⁺] release and dyad [Ca²⁺]. Finally, the increase in LSR_max with increasing SERCA2a activity is because of its ability to lower cytoplasmic [Ca²⁺] between waves, hence reducing the cytoplasmic stimulus of RyR2s.

LSR_min was sensitive to EC₅₀ for RyR2 activation, TSR [CSQ], diffusion between the LSR and TSR, and diffusion in the dyad cleft (notably, dyad gap and tortuosity factor). Otherwise, LSR_min followed a similar pattern of sensitivities as LSR_max albeit to a much lesser degree.

The wave period is influenced by most of the factors that influence LSR_max because it reflects the time taken for a given level of SERCA2a activity to load the LSR to a threshold level. By the same argument, the wave period is reduced by increasing SERCA2a activity. In addition, the wave period is sensitive to bath [Ca²⁺], the permeability of the external membrane to Ca²⁺, and the Ca²⁺ diffusion constant. This is because these factors limit SERCA2a activity by limiting the rate of Ca²⁺ diffusion from the external bath to SERCA2a.

Cyt_max is positively correlated with factors determining the Ca²⁺ store load at wave initiation (SR Ca²⁺ buffering capacity and LSR_max) and negatively correlated with the cytoplasmic buffering by the fixed buffers, troponin (low-affinity site) and CaM, and relatively insensitive to troponin (high-affinity site), myosin, and SERCA2a. Interestingly, the diffusible buffers, ATP and Fluo-4, have the opposite effect on cyt_max to the fixed buffers because diffusible buffers also act as Ca²⁺ carriers to increase Ca²⁺ diffusion in the cytoplasm and external bath, hence increasing SERCA2a activity and LSR_max.

Cyt_min was mainly sensitive to factors influencing the transport of Ca²⁺ into the cytoplasm such as bath [Ca²⁺], surface membrane permeability, and Ca²⁺ diffusion constant, as well as SERCA2a activity, which determines transport of Ca²⁺ out of the cytoplasm. In addition, [ATP] had a positive correlation with cyt_min for the same reasons as it did for cyt_max.

Variables that had no significant impact on any aspect of the wave simulations were: (1) the EC₅₀ for the RyR2 closing rate, which determines the Ca²⁺ dependence of the RyR2 closing rate; (2) [SR lipid], which is the concentration of Ca²⁺ binding lipids in the SR membrane; and (3) cytoplasmic Ca²⁺ buffering by myosin, SERCA2a, CaM, and high-affinity troponin sites and CSQ in the LSR.

Our sensitivity analysis revealed differential effects of the RyR2 opening and closing rate on ember and wave initiation. We investigated the relative importance of RyR2 gating rates and open probability to Ca²⁺ wave properties by measuring the effect on Ca²⁺ waves of a uniform threefold reduction in RyR2 opening and closing rates over the full range of cytoplasmic [Ca²⁺] (Fig. S11. E and F, c.f. blue and green lines), thus reducing RyR2 gating rates without altering channel open probability. The resulting simulations produced Ca²⁺ waves with ∼10% longer period and 10% higher LSR [Ca²⁺] threshold for wave initiation (Fig. S11, A and B, c.f. blue and green lines).

The reason for these changes lay in the differential effects of the opening and closing rate on CICR within dyads that underlies ember formation and CICR between dyads that together initiate waves (see Fig. 9). Slowing the closing rate resulted in a greater reduction in the LSR [Ca²⁺] threshold for ember formation compared with the increase caused by slowing the opening rate (see Fig. S10). The net effect was to decrease this threshold so that embers initiated earlier in the wave cycle, allowing them to generate a larger, longer plateau flux and hence a stronger priming phase (Fig. S11 D, green line). However, stronger priming of the cytoplasm was offset by the slower RyR2 opening rate, which delayed CICR between dyads and wave initiation, consistent with the reduced transient flux (Fig. S11 C, green line).

We also investigated how RyR2 activation at diastolic [Ca²⁺] may affect Ca²⁺ waves by reducing the RyR2 opening rate by approximately threefold at diastolic [Ca²⁺] via a threefold reduction in the parameter G_f, in Eq. 14 in the Supplemental text at the end of the PDF (Fig. S11, red lines). This caused a 5% increase in the wave period and amplitude, and a 10% higher LSR [Ca²⁺] threshold for wave initiation. The reduction in the opening rate caused a relatively small (0.05 mM) decrease in the LSR [Ca²⁺] threshold for ember formation but delayed the CICR phase so that wave initiation occurred at a higher LSR [Ca²⁺]. Once cytoplasmic [Ca²⁺] increased above diastolic levels, reestablishment of the normal RyR2 opening rate along with the higher LSR [Ca²⁺] produced a larger transient flux and wave amplitude (Fig. S11 C, c.f. red and blue lines).

Finally, we investigated how inhibition of RyR2 activity via decreasing channel open durations (open channel block) or increasing channel closed durations (closed channel block) might alter Ca²⁺ waves. We simulated the open channel block by a threefold increase in the RyR2 closing rate and the closed channel block by a threefold decrease in the opening rate (Fig. S12). The open channel block caused a 20% increase in the wave period, amplitude, and LSR [Ca²⁺] threshold, whereas the closed channel block increased all these by 40%. With the open channel block (red lines), the faster RyR2 closing rate increased the LSR [Ca²⁺] threshold for embers, causing them to initiate later in the wave cycle, but once embers did initiate, the higher LSR [Ca²⁺] increased the plateau flux, increasing CICR (note the increased transient flux in Fig. S12 C), thus rapidly stimulating wave initiation. With the closed channel block (green lines), the slower opening rate increased the LSR [Ca²⁺] threshold for embers, but only by a third of that for the open channel block. However, the slower opening rate also decreased CICR, which restrained wave initiation until a higher LSR [Ca²⁺] than seen with the open channel block, thus producing longer wave periods and amplitudes than for the open channel block.

We show that it is possible to generate realistic simulations of Ca²⁺ waves using experimentally determined RyR2 gating kinetics without the need to introduce exaggerated RyR2 channel activation by SR luminal Ca²⁺. We find that sustained, cumulative increases in cytoplasmic [Ca²⁺] occur when brief localized Ca²⁺ release events (Ca²⁺ sparks) become prolonged by orders of magnitude (becoming Ca²⁺ embers) above a threshold SR [Ca²⁺] when CICR between RyR2s within the dyad does not decay to allow normal Ca²⁺ spark termination (Cannell et al., 2013). The resulting small, sustained increase in cytoplasmic [Ca²⁺], along with an associated fall in cytoplasmic buffering power, primes the cytoplasm throughout the cell for CICR between dyads that produces Ca²⁺ waves following a scheme illustrated in Fig. 9. In our model, the progression from sustained CICR within dyads (Ca²⁺ ember formation) to CICR between dyads is both essential for initiation and propagation of Ca²⁺ waves. We predict that the SR [Ca²⁺] threshold for Ca²⁺ embers underpins the [Ca²⁺] threshold for Ca²⁺ waves and SOICR.

Although our seven-compartment dyad model is a considerable simplification of our previous 3D dyad model for Ca²⁺ sparks (Cannell et al., 2013; Laver et al., 2013), the simplified model could still approximate the time course of Ca²⁺ spark activation, termination, and restoration (Fig. S4). Interestingly, the Ca²⁺ spark full width at half maximum values we obtain are 30% smaller than normally seen (Picht et al., 2007), which is a problem common to Ca²⁺ release models that may be due to an error in the common assumption that diffusion obeys Fick’s law (Tan et al., 2007). Like our previous model, the Ca²⁺ spark activation phase resulted from regenerative CICR between RyR2s in the dyad cleft and spark termination resulted from an induction-decay cycle in which Ca²⁺ depletion of the TSR reduces the RyR2 Ca²⁺ flux and hence dyad [Ca²⁺] to a point where RyR2s start to deactivate, further reducing the Ca²⁺ flux. Guo et al. (2012) presented experimental evidence for such a mechanism for spark termination which they called pernicious attrition. Our model also predicted that when the [Ca²⁺] in the LSR exceeds a certain threshold, the TSR fails to deplete sufficiently to initiate the induction-decay cycle that terminates a spark. In this case, replenishment of the TSR by Ca²⁺ diffusion from the LSR is sufficient to sustain a quasi- steady-state release of Ca²⁺ from the dyad, which lasts until the LSR [Ca²⁺] is sufficiently decreased.

Although sustained (>50 ms) Ca²⁺ release occurred quite frequently in our simulations in the lead-up to a Ca²⁺ wave (Fig. S6 F and Fig. S7), the Ca²⁺ embers that result from these sustained releases were detected less frequently in line scans such as in Fig. 6. This is because of the low amplitude of embers compared with sparks and waves, go undetected when the release sites are displaced a mere 0.5 μm from the scan line. (To clearly distinguish embers in our simulations, we utilized a logarithmic fluorescence scale, as illustrated in Fig. 6 C.)

The properties of Ca²⁺ embers produced by our model are in accord with experiments. Long-lasting, localized Ca²⁺ release events (Ca²⁺ embers) were first seen in skeletal muscle where RyR1s were modified by Imperatoxin A (Shtifman et al., 2000) and in cardiomyocytes where RyR2s were either partially blocked by tetracaine or modified by ryanodine or FK506 (Sobie et al., 2002; Zima et al., 2008a). More recent studies observed embers in the absence of RyR2 modifiers. Zima et al. (2008b) reported that embers were promoted at higher LSR [Ca²⁺] and in regions where dyads had a greater connection (i.e., faster diffusion) with the LSR. Brochet et al. (2011) found that Ca²⁺ release underlying Ca²⁺ sparks in rabbit ventricular myocytes had a brief, large initial transient followed by a sustained, low-flux component; both components contributing similar to the total Ca²⁺ release. A cardiomyocyte model by Sato et al. (2016) produced embers lasting >1 s in response to the same factors that we also found important for ember generation, such as reducing the number of RyR2s in the dyad and/or their open probability, increasing Ca²⁺ diffusion between the TSR and LSR, and reducing Ca²⁺ diffusion in the cytosol. Importantly, and in agreement with our simulations, they showed that increased LSR [Ca²⁺] promoted embers.

The model produced waves with velocities that ranged from 60 to 100 μM/s and with amplitudes as small as 2 μM, in line with experimental observations (Bers, 2002; Lipp and Niggli, 1993). Our waves tended to initiate at cell/bath boundaries (at x = 8 μm, y = 8 μm, and z = 52 μm) where the LSR would load more quickly and reach the threshold first. Therefore, waves in our simulations initiated at the end of the cell near z = 52 μm. Experimental observations often show that waves originate in central cell regions. Our model suggests that this could arise from cellular inhomogeneities, which cause local LSR [Ca²⁺] to reach the ember threshold earlier in these regions. This phenomenon may result from a more permeable surface membrane, increased LSR Ca²⁺ uptake, lower ember threshold, or other factors depicted in Fig. 11 that locally reduce the Ca²⁺ wave period.

Despite their modest Ca²⁺ flux, embers are essential precursors to Ca²⁺ waves in our model (Fig. 7 and Fig. 8) and that without embers, Ca²⁺ waves fail to initiate or propagate (Fig. S8). Our analysis of Ca²⁺ release from individual dyads (Fig. 3) shows that the sustained plateau flux of embers is much more likely than the initial transients of embers and sparks to cause CICR at neighboring release sites and produce cumulative increases in local cytoplasmic [Ca²⁺] needed to prime the cytoplasm for waves. Embers are clearly apparent in our line-scan simulations of Ca²⁺ waves at various points along the cell (Fig. 6), but in the wave initiation region, they are difficult to distinguish due to cluttering of many ember events, which may explain why embers have not been reported in association with wave initiation.

Our sensitivity analysis (Fig. 11 and Fig. S10) showed a strong positive correlation with 22 out of 31 model parameters between their effects on LSR [Ca²⁺] thresholds for promoting Ca²⁺ embers and Ca²⁺ waves. Thus, our model predicts that changes in cell properties that lower the LSR [Ca²⁺] threshold for embers are also likely to promote spontaneous Ca²⁺ waves and hence arrhythmias. Therefore, the threshold LSR [Ca²⁺] for SOICR depends on RyR2 sensitivity not only to SR luminal Ca²⁺ as previously thought (Jiang et al., 2004), but also to other factors such as dyad structure (see below), Ca²⁺ buffering in the cytoplasm by ATP and in the TSR by CSQ, and RyR2 sensitivity to cytoplasmic Ca²⁺. Indeed, recent experiments (Kurebayashi et al., 2022) show that SOICR is much more dependent on RyR2 sensitivity to cytoplasmic [Ca²⁺] than to SR [Ca²⁺] in accord with our model predictions.

Although our simulations show that Ca²⁺ embers are an essential precursor to waves, they do not always produce Ca²⁺ waves. It has been long recognized that increased Ca²⁺ buffering in the cytoplasm retards Ca²⁺ wave propagation (Gyorke and Gyorke, 1998) because it reduces CICR between adjacent release sites. A more nuanced buffering action on Ca²⁺ signaling occurs when cytoplasmic Ca²⁺ buffering power is reduced at higher [Ca²⁺] (e.g., see a review by Smith and Eisner [2019]) as shown in Fig. S9 A. Under these conditions, our simulations generated Ca²⁺ embers and wavelets but not waves (Fig. S9 C), demonstrating that declining buffering power at high [Ca²⁺] is a necessary component of the positive feedback provided by CICR that allows embers to initiate waves. More generally, where the initial peak cytoplasmic [Ca²⁺] is small, such as in the case of reduced CSQ in the TSR, cells may adopt a steady state in which the SR is highly permeable to Ca²⁺, but waves do not occur (Fig. 10). This is consistent with observations that skinned cardiomyocytes sometimes do not exhibit Ca²⁺ waves, especially so in cell lines where CSQ is lacking (B.C. Knollmann, personal communication).

Experimental studies indicate that heart rhythm depends more on RyR2 opening and closing rates than it does on RyR2 open probability per se. During the discovery of the first effective treatment for ventricular arrhythmias arising from the overactive RyR2 in catecholaminergic polymorphic ventricular tachycardia (CPVT), it was found that reducing RyR2 open probability alone was insufficient to prevent CPVT. Flecainide, which reduced RyR2 open durations (considered an open channel blocker), was antiarrhythmic, whereas tetracaine, another RyR2 inhibitor that increased closed durations (a closed channel blocker), increased the Ca²⁺ wave amplitude and was proarrhythmic (Watanabe et al., 2009). In permeabilized cardiomyocytes, tetracaine and flecainide increased the period of Ca²⁺ waves. However, tetracaine alone increased the LSR content and wave amplitude (Hilliard et al., 2010).

These different effects of RyR2 opening and closing rates on Ca²⁺ waves can be understood within our model. Simulated open and closed channel RyR2 block in Fig. S12 show that the closed channel block produces increases in peak store load and wave amplitude that are larger than seen with the open channel block. These differences arise from the property that the RyR2 closing rate has a threefold bigger effect on the LSR [Ca²⁺] ember threshold than opening rates, whereas the opening rate is the main driver of CICR between dyads.

Carvajal et al. (2024) reported that the species (rabbit, human) with shorter systolic RyR2 open and closed times have a lower tendency to form delayed after depolarizations (DADs) and a lower SR Ca²⁺ content than species with the longer systolic open and closed times and higher diastolic open probability (mouse, rat) even though all species had similar RyR2 open probabilities at systolic (10 μM) cytoplasmic [Ca²⁺]. We indeed find that shorter RyR2 open and closed times lead to lower peak SR Ca²⁺ content (Fig. S11 B, c.f. blue and green lines). Although our simulations of Carvajal et al. (2024) findings explain how RyR2 open and closing rates in different species can produce different SR Ca²⁺ handling, just how these changes lead to different propensities for DADs requires a cell model that includes the surface membrane.

Our sensitivity analysis predicts that Ca²⁺ wave properties depend on SR structural elements such as RyR2 separation, dyad surface area, TSR thickness, TSR/T-tubule separation, LSR–TSR connectivity, and Ca²⁺ diffusion tortuosity. (Tortuosity in the dyad cleft would depend on the distribution of proteins that limit diffusion paths.) Nanoscale dyad morphology is quite heterogeneous. Electron tomography of ultrarapid high-pressure frozen cardiomyocyte samples from healthy hearts revealed variations between dyads in dyad volume (twofold range) and TSR/T-tubule separation (3–20 nm) (Rog-Zielinska et al., 2021). The size of RyR2 clusters follows an exponential distribution (Hou et al., 2015), and RyR2 spatial distribution within clusters also varies (Asghari et al., 2014). Some of these dyad-to-dyad variations in the RyR2 cluster size have been included in our model. It has been proposed that such variations in cluster size should aid wave initiation on the assumption that waves are initiated by sparks (Xie et al., 2019). However, in our model, waves are initiated by embers, so cluster size variations only produce a minor effect on wave simulations as seen in line-scan images of Ca²⁺ waves (Fig. 6).

Confocal and super-resolution optical microscopy of cardiomyocytes from failing hearts has found that dysfunctional Ca²⁺ release is associated with remodeling of T-tubule structures, colocalization of structural proteins, and distribution and sizes of RyR2 clusters (see a review by Jayasinghe et al. [2018]). Previous models have addressed the influence of ion channel clustering on Ca²⁺ sparks (Iaparov et al., 2021; Walker et al., 2014) and excitation–contraction coupling (Cosi et al., 2019; Tanskanen et al., 2007). Our model predicts that proarrhythmic Ca²⁺ waves will occur at a lower LSR [Ca²⁺] threshold, by structural changes that constrain Ca²⁺ diffusion in the dyad such as increasing tortuosity of diffusion paths or dyad area or decreasing the dyad gap between the TSR and T-tubule. A lower LSR [Ca²⁺] threshold can also be achieved by closer average spacing of RyR2s in clusters and increased connectivity between the TSR and LSR.

Given the high sensitivity of SR Ca²⁺ release to the dyad structure, care should be taken when applying insights from this model to Ca²⁺ release in other cell types (e.g., HEK293) where SR organization is quite different. Also, the accuracy of the model is limited by our current knowledge about important dyad structures and their connectivity with the LSR, which awaits experimental data. Our simplified model of the Ca²⁺ release site structure, RyR2 distribution within clusters, and the hybrid stochastic–deterministic approximation for RyR2 gating restricts the scope of the model to mechanisms by which multiple Ca²⁺ release sites interact to produce Ca²⁺ waves. The deterministic calculation of RyR2 activity during Ca²⁺ release does not capture the stochastic aspects of Ca²⁺ spark transients such as the possibility of Ca²⁺ quarks in cases of insufficient CICR to generate a full Ca²⁺ spark nor does it capture the effects of fragmented RyR2 arrays on Ca²⁺ release investigated in previous Ca²⁺ spark models (Cannell et al., 2013; Cosi et al., 2019; Iaparov et al., 2021). Also, we adopted a simplified two-state RyR2 gating scheme based on single-channel recordings of RyR2s under steady conditions. However, RyR2s are known to have more complex Mg²⁺-dependent gating kinetics (Laver and Honen, 2008; Zahradníková and Zahradník, 2012), which, in the rapidly changing ion concentrations during sparks and waves, could generate a different Ca²⁺ release time course to the steady-state model used here (Iaparov et al., 2022). Despite these simplifications, our model still captured the overall time course of Ca²⁺ spark activation by CICR, termination by induction decay, and restitution by TSR refill.

Since our model does not include an active surface membrane, it mainly offers insights into spontaneous Ca²⁺ release phenomena generated by the SR alone such as Ca²⁺ sparks, embers, and waves. This study demonstrates that Ca²⁺ wave formation, driven by increasing LSR [Ca²⁺], involves a progression from sustained CICR within dyads (Ca²⁺ ember formation) to CICR between dyads. This provides a plausible alternative to the traditional mechanism of CICR between dyads alone in generating SOICR and Ca²⁺ waves. CICR both within and between dyads arises from the complex interplay of RyR2 Ca²⁺-dependent kinetics, organelle structure, and Ca²⁺ diffusion and buffering within the cell. We anticipate these and other novel factors that regulate SOICR in our simulations will guide the development of more complete cell models and future experimental investigations into cardiac excitation–contraction coupling and cardiac rhythm.

Data are available in the published article and its online supplemental materials.

David A. Eisner served as an editor.

This work was funded by the NSW Health Infrastructure grant through the Hunter Medical Research Institute to D.R. Laver. This research was supported by Australian Government Research Training Program Scholarship to M. Vysma.

Author contributions: M. Vysma: investigation, methodology, software, and writing—review and editing. J.S. Welsh: conceptualization, data curation, formal analysis, methodology, software, and supervision. D.R. Laver: conceptualization, formal analysis, investigation, methodology, project administration, software, supervision, visualization, and Writing—original draft, review, and editing.