Ryanodine Receptor Allosteric Coupling and the Dynamics of Calcium Sparks

A two-state Ca²⁺-activated RyR model

Stochastic models of single channel gating often take the form of continuous-time discrete-state Markov chains (for review, see (42,43)). For example, the state-transition diagram for a two-state Ca²⁺-activated RyR model is defined as

(1)

where k⁺c^η and k⁻ are transition rates with units of time⁻¹, k⁺ is an association rate constant with units conc^−η time⁻¹, η is the cooperativity of Ca²⁺ binding (usually chosen to be η = 2), and c is the local [Ca²⁺]. If c(t) is specified, then Eq. 1 defines a discrete-state continuous-time stochastic process, S(t), with the state space S ∈ {𝒞, 𝒪}. When the local [Ca²⁺] is not time-varying—for example, a fixed background [Ca²⁺] that we denote as c_∞—then Eq. 1 corresponds to the well-known telegraph process with infinitesimal generator or Q-matrix (42,44) given by

(2)

Each off-diagonal element of Eq. 2 is the probability per unit time of a transition from state i to state j,

and the diagonal elements are selected to ensure that the row sums of Q are zero This condition ensures conservation of probability where

is the element in the i^th row and j^th column of the matrix exponential. For example, during a small time step Δt, the probability that a channel initially in state i makes a transition into state j is approximated by p_ij ≈ [I + QΔt]_ij, and in this case it is clear that is required for Note that all of the statistical properties of the two-state channel model diagrammed in Eq. 1 can be calculated from the Q-matrix (Eq. 2), and that this matrix can be decomposed as

(3)

where the matrices

collect the dissociation and association rate constants, respectively.

Collective gating of RyR clusters

In a natural extension of the single channel modeling approach, a model Ca²⁺ release site composed of N channels is the vector-valued Markov chain, where S_n(t) is the state of channel n at time t (45). We will denote release site configurations as a vector where i_n is the state of channel n. The transition rate from release site configuration to denoted by q_ij,

(4)

is nonzero if the origin (i) and destination (j) release site configurations are identical with the exception of one channel—that is, i_n = j_n for all n ≠ n′ where 1 ≤ n′(i, j) ≤ N is the index of the channel changing state—and the i_n′ → j_n′ transition is included in the single-channel model.

More formally, the transition rates q_ij for a release site composed of N identical Ca²⁺-regulated channels (Eq. 2) are given by

(5a)

(5b)

where either K ⁻[i_n′, j_n′] or K⁺[i_n′, j_n′] is the rate constant for the transition being made (only one of which is nonzero) and c(i, j) is the relevant [Ca²⁺], that is, the concentration experienced by channel n′(i, j) in the origin configuration i. In the following section, we show how c(i, j) depends on the mathematical representation of the release site ultrastructure and buffered Ca²⁺ diffusion.

Although it may not be practical to do so for large release sites, the infinitesimal generator matrix, Q = (q_ij), for a model Ca²⁺ release site can be constructed by enumerating transition rates according to Eq. 5 and selecting the diagonal elements q_ii to ensure the rows sum to zero.

Release site ultrastructure and the Ca²⁺ microdomain

Because Ca²⁺-activated RyRs experience coupling mediated by the buffered diffusion of intracellular Ca²⁺, the model includes a mathematical representation for the landscape of local [Ca²⁺] near the Ca²⁺ release site (the so-called Ca²⁺ microdomain) required to specify c(i, j) in Eq. 5b. For simplicity, we assume channels are instantaneously coupled via the Ca²⁺ microdomain (30,31)—that is, the formation and collapse of the local peaks in the Ca²⁺ profile are fast compared to the closed and open dwell times of the channels—and we assume the validity of linearly superposing local [Ca²⁺] increases due to individual channels at the release site (25,27). We also assume that all channels are localized on a planar section of SR membrane (z = 0) so that the position of the pore of channel n can be written as

Assuming a single high concentration Ca²⁺ buffer and using the steady-state solution of the linearized equations for the buffered diffusion of intracellular Ca²⁺ (25,26), the increase in [Ca²⁺] above background at position is given by

(6a)

where

(6b)

(6c)

and

(6d)

In these equations, the sum is over all channels at the release site, σ_n is the source amplitude of channel n (number of Ca²⁺ ions per unit time); D_c and D_b are the diffusion coefficients for free Ca²⁺ and the Ca²⁺ buffer, respectively; is the buffer association rate constant; is the buffer dissociation rate constant, and [B]_T is the total concentration of the Ca²⁺ buffer. Assuming all RyRs have identical source amplitudes,

(7a)

and

(7b)

where i_Ca is the unitary current of each channel, 2 is the valence of Ca²⁺, and F is Faraday's constant.

While Eqs. 6 and 7 define the [Ca²⁺] at any position r for a given release site ultrastructure, {r_n}, it is helpful to summarize channel-to-channel Ca²⁺ interactions using an N × N coupling matrix C = (c_nm) that provides the increase in [Ca²⁺] over the background (c_∞) experienced by channel m when channel n is open. If specifies the position of the Ca²⁺ regulatory site for channel m located a small distance r_d above the channel pore, then

(8)

Using this expression we can determine the Ca²⁺ concentrations needed to specify the rates of Ca²⁺-mediated transitions in Eqs. 5a and 5b, that is,

(9a)

where

(9b)

n′(i, j) is the index of the channel changing state, and i_n is the state of channel n.

Fig. 2 A uses Eqs. 6 and 7 and calmodulin-like buffer parameters (see Table 1) to calculate the Ca²⁺ microdomain near a cluster of N = 25 open RyRs organized on a Cartesian lattice (Fig. 1 B). The strength of Ca²⁺ interactions at the release site can be modified by changing any of the parameters in Eqs. 6 and 7, including the channel source amplitude, buffer parameters, or the diffusion constant for free Ca²⁺. For example, Fig. 2 B shows that increasing the total buffer concentration ([B]_T) decreases the local [Ca²⁺] experienced by the RyRs. Similarly, Fig. 2 C shows that the Ca²⁺ coupling strength defined as the average of the off-diagonal elements of the coupling matrix,

(10)

is a decreasing function of the total buffer concentration [B]_T for any fixed unitary current i_Ca and an increasing function of i_Ca for any fixed [B]_T. Note that the unitary current of RyRs in vivo has been estimated to be <0.6 pA and as low as 0.07 pA in the presence of physiological concentrations of Mg²⁺ (46–48). Because the model does not explicitly include localized depletion of luminal Ca²⁺, a phenomenon that is expected to reduce the effective unitary current of RyRs in vivo, our standard parameter set includes a unitary current of 0.04 pA (see Table 1).

TABLE 1.

Default parameters used in Ca²⁺ release site simulations for both the full model and the mean-field reduction (when applicable)

Parameter	Value	Unit	Description
Single channel parameters
k+	0.04	μM−ηms−1	Association rate constant
k−	1	ms−1	Dissociation rate constant
c∞	0.1	μM	Background [Ca2+]
η	2		Cooperativity of Ca2+ binding
iCa	0.04	pA	Effective unitary current
rd	30	nm	Pore to regulatory site distance
Buffer parameters
	100	μM−1 s−1	Association rate constant
	38	s−1	Dissociation rate constant
Dc	250	μm2 s−1	Ca2+ diffusion coefficient
Db	32	μm2 s−1	Buffer diffusion coefficient

Open in a new tab

Single-channel kinetic parameters are selected for a dissociation constant or K_d = 5 μM (70). Buffer parameters correspond to calmodulin (27,82). Although the exact location of the Ca²⁺-regulatory site is unknown, the pore-to-regulatory site distance is consistent with cryo-electron microscopy data that suggests the RyR oligomer has a large 29 × 29 × 12 nm cytoplasmic assembly and a transmembrane assembly that protrudes 7 nm from the center of the cytoplasmic assembly (16,83).

Allosteric interactions between physically coupled channels

Following the methodology presented in Stern et al. (36), the RyR cluster model with Ca²⁺-mediated coupling is extended to include allosteric interactions between neighboring channels. We begin by defining dimensionless free energies of interaction ɛ_ij (units of k_BT) that specify the change in free energy experienced by a channel in state j when allosterically coupled to a channel in state i. For convenience we collect these interaction energies in an M × M matrix ℰ where M is the number of states in the single-channel model and ɛ_ij = ɛ_ji (i ≠ j) to satisfy the requirement of thermodynamic reversibility. For the two-state single-channel model considered in this article,

(11)

where Because allosteric interactions require physical contact between neighboring RyRs, the model formulation includes a symmetric N × N adjacency matrix defined as

(12)

where a_nn = 0 because channels do not experience allosteric interactions with themselves. The nonzero elements of A are chosen consistent with release site ultrastructure (e.g., dotted lines in Fig. 1 B).

To include the effect of allosteric coupling in the Ca²⁺ release site model, the total allosteric energy experienced by channel n′(i, j) in the origin and destination configurations of an i → j transition are calculated as

(13)

where the sum is over all N channels, a_nn′ are elements of A, and and are entries of ℰ. The difference between these total allosteric energies (γ_j − γ_i) is used to modify the equilibrium constant of the i → j transition, that is,

(14a)

(14b)

and

(14c)

where and denote unmodified rates calculated using Eq. 5 and the parameters 0 ≤ ν_ij ≤ 1 and ν_ji = 1 −ν_ij (36) partition contributions due to allosteric coupling between the forward (q_ij) and reverse (q_ji) rates. While ν_ij and ν_ji can potentially have different values for every transition i → j, we assume transition rates involving the association of Ca²⁺ are diffusion-limited. Thus, transition rates for release site configuration changes where channels make 𝒞 → 𝒪 transitions are assigned ν = 0. Conversely, ν = 1 for all other configuration changes where channels make 𝒪 → 𝒞 transitions.

Ca²⁺ and allosteric coupling at a three RyR cluster

To clarify the model formulation, transition rate expressions corresponding to the example configuration changes shown in Fig. 3 A are written below. These configuration changes involve a triangular cluster of three two-state RyRs experiencing Ca²⁺ coupling and nearest-neighbor allosteric interactions. The corresponding Ca²⁺ coupling matrix and allosteric adjacency matrix are

(15)

respectively, where the c_nm are determined using Eq. 8. In each panel of Fig. 3 A, the total allosteric energy experienced by the RyR changing state (labeled with asterisks) is calculated for both the origin (i) and destination (j) configurations using Eq. 13.

The i → j configuration changes shown in Fig. 3 A each involve an RyR making a Ca²⁺-mediated 𝒞 → 𝒪 transition (see Eq. 1) at rate q_ij that is a function of c(i,j), that is, the [Ca²⁺] experienced by the channels changing state (Eq. 9). Let us number the RyRs in a counterclockwise fashion beginning with the channel changing state. For the 𝒞𝒞𝒞 → 𝒪𝒞𝒞 configuration change shown in Fig. 3 Aa, c(i,j) = c_∞ because all channels are closed in the origin configuration 𝒞𝒞𝒞. For the 𝒞𝒞𝒪 → 𝒪𝒞𝒪 configuration change, c(i,j) = c_∞ + c₃₁ because channel 3 is open in configuration 𝒞𝒞𝒪 (Fig. 3 Ab). Similarly, for the 𝒞𝒪𝒪 → 𝒪𝒪𝒪 configuration change, (Fig. 3 Ac). Having determined the appropriate [Ca²⁺] concentrations, Eq. 5b is needed to calculate the transition rates:

(16a)

(16b)

and

(16c)

Because it is assumed that configuration changes involving the binding of Ca²⁺ are diffusion limited, these rates are not modified due to allosteric interactions (i.e., ν_ij = 0).

Conversely, j → i configuration changes shown in Fig. 3 A involve channels making unimolecular 𝒪 → 𝒞 transitions at the base rate q_ji = k⁻ that is modified by the change in allosteric interaction energy experienced by the channel changing state. Using ν_ji = 1, the rates for j → i configuration changes are given by (Eq. 14c)

(17a)

(17b)

and

(17c)

Note that in these transition rate expressions, the elements of the allosteric interaction energy matrix occur as the differences ɛ_𝒞𝒞 − ɛ_𝒞𝒪 and ɛ_𝒪𝒞 − ɛ_𝒪𝒪. This is true regardless of the number of channels, and we may without loss of generality fix ɛ_𝒞𝒪 = ɛ_𝒪𝒞 = 0. That is, we will probe the effects of allosteric interactions on Ca²⁺ release site dynamics by varying only the change in free energy due to allosterically interacting closed-closed (ɛ_𝒞𝒞) and open-open (ɛ_𝒪𝒪) channel pairs. Because we are primarily concerned with the effects of allosteric interactions that promote synchronous gating, we assume allosteric interactions stabilize closed-closed and/or open-open channel pairs (i.e., ɛ_𝒞𝒞 ≤ 0 and ɛ_𝒪𝒪 ≤ 0). For simplicity, we focus on three allosteric coupling paradigms in which allosteric interactions stabilize

1. Closed-closed channel pairs (ɛ_𝒞𝒞 = < 0, ɛ_𝒪𝒪 = 0).

2. Open-open channel pairs (ɛ_𝒪𝒪 = 0, ɛ_𝒪𝒪 < 0).

3. Both closed-closed and open-open channel pairs in a balanced fashion (ɛ_𝒪𝒪 = ɛ_𝒪𝒪 < 0).

The simulations shown in Fig. 3, B–E, demonstrate how synchronizing allosteric interactions included in these three ways affect the dynamics of the synchronous gating of the three RyRs. Simulations are carried out using the exact numerical method presented in Appendix A and, for simplicity, the configuration of the RyRs is summarized by plotting only the number of open channels (N_𝒪) as a function of time. Interestingly, Fig. 3 B demonstrates that synchronizing allosteric interactions are not required (ɛ_𝒞𝒞 = ɛ_𝒪𝒪 = 0) for channels to exhibit synchronous gating. Rather, channels may exhibit coupled gating that is mediated entirely via the buffered diffusion of local Ca²⁺ as long as the average Ca²⁺ coupling strength is sufficient (c_* = 0.75 μM) (31). Shaded bars in the left panel of Fig. 3 B show the steady-state probability distribution for the number of open RyRs (N_𝒪) directly calculated from the relevant Q-matrix as described in Appendix B. The disagreement between these results and the open bars, showing a binomial distribution with the same mean, is a signature of the cooperative gating of these RyRs.

While Fig. 3 B demonstrates that the synchronous gating of channels can be mediated entirely via Ca²⁺, Fig. 3, C–E, show how synchronizing allosteric interactions affect the dynamics of coupled gating. For example, Fig. 3 C demonstrates that when closed channel pairs are stabilized (ɛ_𝒞𝒞 = −0.8, ɛ_𝒪𝒪 = 0), the steady-state probability of having zero open channels (N_𝒪 = 0) increases while the probability of N_𝒪 = 3 decreases relative to Fig. 3 B. Conversely, Fig. 3 D illustrates that when allosteric interactions stabilize open channel pairs (ɛ_𝒞𝒞 = 0, ɛ_𝒪𝒪 = −0.8), the probability of having a maximally activated release site (N_𝒪 = 3) increases. In Fig. 3 E allosteric interactions stabilize closed-closed and open-open channel pairs in a balanced fashion (ɛ_𝒞𝒞 = ɛ_𝒪𝒪 = −0.8) and the probability of both N_𝒪 = 0 and N_𝒪 = 3 increases while the probability of N_𝒪 = 1 and N_𝒪 = 2 decreases compared to Fig. 3 B.

Effects of Ca²⁺ and allosteric coupling strength on spontaneous sparks

The previous section demonstrated how the dynamics of coupled RyR gating may depend on synchronizing allosteric interactions that stabilize closed channel pairs, open channel pairs, or both in a balanced fashion. In this section, release sites composed of 25 nearest-neighbor coupled RyRs organized on a Cartesian lattice (see Fig. 1 B) are used to investigate how Ca²⁺ spark generation and termination depend on both the strength of coupling mediated by the Ca²⁺ microdomain and the strength of synchronizing allosteric interactions introduced in one of these three ways. Note that nearest-neighbor coupling implies that each channel experiences allosteric interactions with 2–4 other channels, while increases in the Ca²⁺ microdomain due to open RyRs are experienced by all channels.

Fig. 4 A shows a simulation in which the strength of allosteric interactions (ɛ_𝒞𝒞 = −0.2, ɛ_𝒪𝒪 = 0) and Ca²⁺ coupling (c_* = 0.55 μM) are selected to illustrate the phenomenon of stochastic Ca²⁺ excitability reminiscent of spontaneous Ca²⁺ sparks. While the channels at the release site are closed most of the time (N_𝒪 < 5), on occasion the RyRs simultaneously open (N_𝒪 ≈ 25). Fig. 5 shows that the sparks observed in Fig. 4 A are sensitive to changes in the strength of allosteric interactions that stabilize closed channel pairs. For example, the release site is tonically active when allosteric interactions are not included in simulations (ɛ_𝒞𝒞 = ɛ_𝒪𝒪 = 0, Fig. 5 A). On the other hand, sparks fail to initiate when the strength of allosteric interactions that stabilize closed channel pairs is greater than in Fig. 4 A (ɛ_𝒞𝒞 = −0.4, ɛ_𝒪𝒪 = 0, Fig. 5 B).

A response measure that is strongly correlated with the presence of sparks in Monte Carlo simulations is the so-called Ca²⁺ spark Score introduced in (31). The Score is defined as the index of dispersion of the fraction of open channels (f_O = N_𝒪/N) and is given by

(18)

Score values >0.3 are indicative of spark-like excitability in stochastic Ca²⁺ release site simulations (30,31). For example, using the observed probability distribution for the number of open channels at the release site estimated from a long Monte Carlo simulation as described in Appendix B (Fig. 4 B), the Score corresponding to the simulation shown in Fig. 4 is a high value of 0.71. Conversely, the tonically active release site shown in Fig. 5 A has a low Score of 0.013 because E[N_𝒪] is large. The quiescent release site shown in Fig. 5 B also has a low Score of 0.082.

While Figs. 4 and 5 demonstrated that Ca²⁺ sparks are sensitive to changes in the strength of allosteric interactions that stabilize closed channel pairs, Fig. 6 A shows that sparks observed in simulations of a 25 RyR release site are also sensitive to the Ca²⁺ coupling strength. For example, triangles show the Score (reported as the mean ± SD of 10 Monte Carlo simulations) as a function of c_* when the strength of allosteric interactions that stabilize closed channel pairs is ɛ_𝒞𝒞 = −0.2 as in Fig. 4. The Ca²⁺ coupling strength (c_*) is systematically varied by increasing or decreasing the total buffer concentration ([B]_T). Note that sparks are observed in simulations (Score > 0.3) over a range of Ca²⁺ coupling strengths but are not observed (Score < 0.3) in simulations that use c_* < 0.4 μM because the Ca²⁺ coupling strength is insufficient to initiate sparks. Similarly, Score < 0.3 when c_* > 0.7 μM because the Ca²⁺ coupling strength is too large to allow spark termination. Fig. 6 A also shows that the optimal Ca²⁺ coupling strength—that is, the c_* resulting in the highest Score—is sensitive to the strength of allosteric interactions that stabilize closed channels. Indeed, comparing circles (ɛ_𝒞𝒞 = 0) and squares (ɛ_𝒞𝒞 = −0.4) to triangles (ɛ_𝒞𝒞 = −0.2), we notice that the optimal c_* is an increasing function of the magnitude of ɛ_𝒞𝒞. In comparison, Fig. 6 B demonstrates that as the strength of allosteric interactions that stabilize open channel pairs increases, the optimal c_* decreases. On the other hand, Fig. 6 C shows that increasing the strength of allosteric interactions that stabilize both closed-closed and open-open channel pairs in a balanced fashion has little effect on the optimal value of c_*.

Fig. 6 demonstrates that sparks depend on c_*, ɛ_𝒪𝒪, and ɛ_𝒞𝒞 in a complicated manner. For example, sparks that are eliminated as c_* increases may be recovered by increasing the strength of allosteric interactions that stabilize closed channel pairs (ɛ_𝒞𝒞) or by decreasing the strength of allosteric interactions that stabilize open channel pairs (ɛ_𝒪𝒪). On the other hand, sparks that are eliminated as c_* decreases may be recovered by decreasing the magnitude of ɛ_𝒞𝒞 or increasing the magnitude of ɛ_𝒪𝒪. Note that for all three types of allosteric interactions there are Ca²⁺ coupling strengths (c_*) for which stronger interactions lead to more robust sparks. Indeed, summary plots in Fig. 7 A show that the Score at these optimal c_* values is a monotonically increasing function of the strength of allosteric interactions. Interestingly, the Score is enhanced the most when both closed-closed and open-open channel pairs are increasingly stabilized in a balanced fashion (circles).

In Fig. 7 B the sensitivity of sparks to the Ca²⁺ coupling strength is quantified using the full width at half-maximum (FWHM) of cubic spline fits to the results of Fig. 6 (dashed lines); a larger FWHM implies less sensitivity to changes in c_*. The triangles of Fig. 7 B show that sparks are less sensitive to variations in c_* as the strength of allosteric interactions that stabilize closed channel pairs increases. Conversely, the squares show that sparks are more sensitive to c_* as the strength of allosteric interactions that stabilize open channel pairs increases. The circles show that increasing the strength of allosteric interactions that stabilize both closed-closed and open-open channel pairs in a balanced fashion has little effect on the FWHM.

The effect of washing out allosteric interactions on spark statistics

In the previous section we showed how the presence or absence of Ca²⁺ sparks depends on both the strength of Ca²⁺ coupling (c_*) and the strength of stabilizing allosteric interactions (ɛ_𝒞𝒞 and ɛ_𝒪𝒪). Next, we investigate how spark statistics (duration, interspark interval, and frequency) are affected by washing-out stabilizing allosteric interactions, that is, we study how these spark statistics change as an increasing fraction of nearest-neighbor allosteric couplings are removed. Many experimental studies show that genetic deficiencies in, and the pharmacological washout of, the FK-binding proteins that mediate allosteric interactions lead to cardiac arrhythmias and changes in spark dynamics (49–51). The following simulations aim to clarify how these experimental results may be interpreted as evidence for allosteric interactions that stabilize closed channel pairs, open channel pairs, or both (see Discussion).

The shaded bars in Fig. 8 A are probability distributions of spark duration and interspark interval estimated from multiple spark simulations (the mean is indicated by shaded triangles). As in Fig. 4, twenty-five RyRs experience nearest-neighbor allosteric interactions that stabilize closed channel pairs (ɛ_𝒞𝒞 = −0.2, ɛ_𝒪𝒪 = 0), and the Ca²⁺ coupling strength is selected to ensure a high Score (c_* = 0.58 μM). Spark duration is defined as the period beginning when one-fifth of the channels at the release site open (N_𝒪 = 4 → 5) and ending when all channels close (N_𝒪 = 0), thus excluding small sparks from the calculation. Interspark interval is the time between the end of a spark and the beginning of the subsequent spark.

For comparison, open bars in Fig. 8 A are the spark duration and interspark interval distributions after one-fifth of the nearest-neighbor allosteric couplings are selected at random and eliminated from the simulations. Notice that this washout of allosteric interactions that stabilize closed channel pairs has the effect of increasing the expected spark duration and decreasing the expected interspark interval (compare open and shaded triangles). On the other hand, Fig. 8 B shows that when allosteric interactions stabilize open channel pairs (ɛ_𝒞𝒞 = 0, ɛ_𝒞𝒞 = −0.2, and c_* = 0.40 μM), removing one-fifth of these couplings decreases the expected spark duration with little change to the interspark interval. When both closed-closed and open-open channel pairs are stabilized in a balanced fashion (ɛ_𝒞𝒞 = ɛ_𝒪𝒪 = −0.2, c_* = 0.48 μM), washout of allosteric couplings decreases interspark interval but has little effect on spark duration (Fig. 8 C).

To further probe the effects of washing out allosteric interactions, Fig. 9, A and B, show the mean and standard deviation of spark duration and interspark interval plotted against the fraction of allosteric couplings removed from simulations (denoted as φ). Similar to Fig. 8, allosteric interactions are included to stabilize closed channel pairs (triangles), open channel pairs (squares), or both in a balanced fashion (circles). In each case, the Ca²⁺ coupling strength of c_* = 0.58, 0.40, and 0.48 μM, respectively, is selected to maximize the Ca²⁺ spark Score before the washout of synchronizing allosteric interactions (φ = 0); thus, the Score is always a decreasing function of φ (not shown). When the squares and circles of Fig. 9, A and B, are recalculated using c_* = 0.58 μM, qualitatively similar results are obtained.

The results shown in Fig. 9, A and B, are consistent with those shown in Fig. 8. When allosteric interactions that stabilize closed channel pairs are washed out (increasing φ), spark duration increases and interspark interval decreases (triangles). When allosteric interactions that stabilize open channel pairs are washed out, spark duration decreases but interspark interval is largely unchanged (squares). When both closed-closed and open-open channel pairs are stabilized in a balanced fashion, washout causes interspark interval to decrease but spark duration is unchanged (circles). Notice that the standard deviations (bars) of spark statistics are approximately proportional to the means regardless of the degree of washout.

Because the allosteric couplings washed out in the simulations of Fig. 9, A and B, are randomly selected, there are many realizations of the allosteric adjacency matrix A consistent with any nonzero φ. To show the effects of variations in allosteric connectivity on spark dynamics, multiple symbols plotted at each value of φ show the mean spark duration and interspark interval using different realizations of A. The proximity of these symbols to each other at any given value of φ indicates that the dynamics of Ca²⁺ sparks—as measured by duration and interspark interval—are largely insensitive to these variations in allosteric connectivity.

Fig. 9 C shows the spark frequency—defined as the reciprocal of the sum of the mean spark duration and interspark interval—plotted against φ for the three allosteric coupling paradigms. When allosteric interactions that stabilize closed channel pairs are washed out, spark frequency increases but ultimately decreases (triangles). When allosteric interactions stabilize open channel pairs (squares), spark frequency is a nearly constant function of φ. When both closed-closed and open-open channel pairs are stabilized in a balanced fashion, spark frequency increases during washout (circles).

A mean-field RyR cluster model

In previous sections, we used Monte Carlo simulations to study how both the strength of Ca²⁺ coupling and stabilizing allosteric interactions contribute to the dynamics of sparks. Much of the complexity of these simulations is due to the spatially explicit account of channel-to-channel coupling represented by the Ca²⁺ coupling matrix C and the allosteric adjacency matrix A. To facilitate parameter studies of the effects of allosteric coupling on spark statistics, this section presents a mean-field approximation applicable to a cluster of two-state RyRs coupled via the buffered diffusion of Ca²⁺ and nearest-neighbor allosteric interactions.

The mean-field approximation is perhaps best introduced by considering a simplified Ca²⁺ coupling matrix that takes the form (31)

(19)

where the identical off-diagonal elements (c_*) are the average of the N(N – 1) off-diagonal elements of the original Ca²⁺ coupling matrix C (Eq. 10). (The diagonal elements c_d that represent domain Ca²⁺ are inconsequential to simulations involving clusters of RyRs with no Ca²⁺-mediated transition out of an open state.) Consider also an allosteric adjacency matrix that takes a similar simplified form,

(20)

where 0 ≤ a_* ≤ 1 is the average allosteric connectivity calculated from the off-diagonal elements of the original allosteric adjacency matrix A = (a_nm),

(21)

Note that it is not possible to choose a release site ultrastructure so that is equal to C with N > 3 channels on a planar membrane. Likewise, will not be equal to A unless allosteric coupling is all-to-all, a situation not consistent with RyR clusters in which the extent of interchannel physical coupling is limited to nearest neighbors. Nevertheless, in simulations performed using and the RyRs are indistinguishable and the [Ca²⁺] and allosteric interaction energy experienced by channels depends only on the number of open and closed RyRs at the release site. Importantly, simulations using and satisfy a lumpability condition that allows all release site configurations with the same number of channels in each state to be agglomerated without further approximation (31,52). This yields a contracted Markov chain with state-transition diagram

(22)

where the state of the system S(t) ∈ {0, 1, …, N} is the number of open channels N_𝒪 at the release site and q_ij is the rate of the N_𝒪 = i → j transition (see below). The number of closed channels is given by N_𝒞 = N − N_𝒪.

Equation 22 describes a birth-death process with boundaries with skip-free transitions that increase (N_𝒪 = n → n + 1) or decrease (N_𝒪 = n → n − 1) the number of open channels at the release site. The N + 1 by N + 1 generator matrix corresponding to Eq. 22 is tridiagonal,

(23)

with diagonal elements (⋄) selected to ensure row sums of zero. The birth rate (q_{n, n+1}) for the n → n + 1 transition is given by

(24)

where (N – n) is the number of closed channels at the release site that may potentially open and c_∞ + N_𝒪c_* is the [Ca²⁺] experienced by all RyRs. While our assumption that the binding of Ca²⁺ is diffusion-limited leads to birth rates that are not dependent on allosteric energies, the death rates are modified due to allosteric interactions and are given by

(25)

where n is the number of open channels at the release site that may potentially close, the coefficients (n – 1) and (N – n) are the number of open and closed neighbors, and ɛ_𝒪𝒞 − ɛ_𝒪𝒪 and ɛ_𝒞𝒞 − ɛ_𝒞𝒪 are the differences in free energies experienced by a transitioning channel due to allosteric couplings with neighboring open and closed channels. Because we have without loss of generality set ɛ_𝒞𝒪 = ɛ_𝒪𝒞 = 0, Eq. 25 simplifies to

(26)

Note that the mean-field RyR cluster model has only nine parameters: N, k⁺, k⁻, η, c_∞, c_*, ɛ_𝒞𝒞, ɛ_𝒪𝒪, and a_*.

Representative mean-field simulations

Fig. 10 A shows representative simulations of 25 mean-field coupled RyRs arranged according to the strength of Ca²⁺ coupling (c_*) and allosteric interactions (ɛ_𝒞𝒞) used. These allosteric interactions stabilize closed channel pairs (ɛ_𝒪𝒪 = 0) and the average allosteric connectivity is a_* = 0.13, as calculated using the adjacency matrix for 25 nearest-neighbor coupled RyRs organized on a Cartesian lattice (see Fig. 1). Notice that sparks are only observed on the diagonal panels of Fig. 10 A, indicating that increased c_* can be compensated for by more negative ɛ_𝒞𝒞. Release sites are tonically active when c_* is large and represents weak allosteric interactions (upper right panels), while release sites are quiescent when c_* is small and represents strong allosteric interactions (lower left panels). These mean-field results are consistent with simulations that use the full model when allosteric interactions stabilize closed channel pairs (Figs. 4 and 5, and Fig. 6 A). Mean-field simulations that include allosteric interactions that stabilize open channel pairs or both closed-closed and open-open channel pairs in a balanced fashion (not shown) are also consistent with the full model (Fig. 6, B and C).

The panels of Fig. 10 B show the birth rates (q_{n, n+1}) used in each column of Fig. 10 A plotted as a function of the number of open channels at the release site (n = N_𝒪). Note that while the q_{n, n+1} are small when n is either small or large, the birth rates are accelerated for intermediate n, and this acceleration is enhanced as c_* increases. The panels of Fig. 10 C show the death rates (q_{n, n–1}) used in the simulations of each row of Fig. 10 A plotted as a function of n. Notice that when allosteric interactions are not included in simulations (top panel, ɛ_𝒞𝒞 = ɛ_𝒪𝒪 = 0), the death rates q_{n, n–1} are a linear increasing function of n. However, as the magnitude of ɛ_𝒞𝒞 increases, q_{n, n–1} is accelerated for all values of n < N with the most significant acceleration at intermediate n. While Fig. 10 C shows how the death rates change with the strength of allosteric interactions that stabilize closed channel pairs, qualitatively different changes in the death rates are observed when allosteric interactions stabilize open channel pairs, or both closed-closed and open-open channel pairs in a balanced fashion. For example, when ɛ_𝒞𝒞 = 0, the death rates q_{n, n–1} decrease for all n > 1 as ɛ_𝒪𝒪 becomes more negative. On the other hand, the birth rates q_{n, n–1} increase for small n but decreases for large n when both ɛ_𝒞𝒞 and ɛ_𝒪𝒪 become more negative (results not shown).

Comparison of mean-field approximation and full model

In the previous section, we demonstrated mean-field simulations may exhibit stochastic Ca²⁺ excitability reminiscent of Ca²⁺ sparks. Similar to full model simulations, these sparks are sensitive to variations of the Ca²⁺ coupling strength (c_*) and the allosteric coupling strengths (ɛ_𝒞𝒞,ɛ_𝒪𝒪) used. In this section we validate the mean-field approximation by comparing the Ca²⁺ spark Score estimated from Monte Carlo simulations of the full model to the Score calculated directly from the Q-matrix of the corresponding mean-field model. In this comparison, the c_* and a_* of the mean-field model are calculated from the C and A of the full model, and the parameters of the single-channel models used are identical.

The symbols in Fig. 11 A plot the Score (mean ± SD of 10 trials) of Monte Carlo simulations using the full model as a function of the Ca²⁺ coupling strength (c_*) for release sites of different sizes (N) when allosteric interactions stabilize closed channel pairs (ɛ_𝒞𝒞 = −0.2, ɛ_𝒪𝒪 = 0). The dashed lines show the Score calculated using Q of the corresponding mean-field approximations. Both full and reduced models demonstrate that the optimal Ca²⁺ coupling strength, that is, the c_* that yields the highest Score, decreases as a function of N. Moreover, the range of c_* values that result in sparks (Score > 0.3) decreases as N increases. This inverse relationship between the optimal c_* and the release site size N, and the increase in the sensitivity of sparks to variations in c_* as N increases, are also observed when allosteric interactions stabilize open channel pairs (ɛ_𝒞𝒞 = 0, ɛ_𝒪𝒪 = −0.2) or both closed and open channel pairs in a balanced fashion (ɛ_𝒞𝒞 = ɛ_𝒪𝒪 = −0.2) (not shown).

Although the Score obtained using the full model and the mean-field approximation agree qualitatively (Fig. 11 A), the optimal c_* and the maximum Score for any given value of N show quantitative differences that becomes more evident with large N. Fig. 11 B shows that the Score (open circles) of simulations that use mean-field Ca²⁺ coupling () and nearest-neighbor allosteric coupling (A) are similar to mean-field model results (dashed line). Similarly, the Score (solid circles) of simulations that use the full Ca²⁺ coupling matrix (C) and mean-field allosteric interactions () show improved agreement with full model results (open triangles). These results suggest that the differences between the full model and the mean-field approximation are a consequence of the assumption of mean-field Ca²⁺ coupling and not the assumption of mean-field allosteric coupling.

Effect of allosteric coupling on Ca²⁺ spark statistics

The reduced state space of the mean-field approximation (N + 1) as opposed to the full model (2^N) greatly facilitates the calculation of spark statistics. For release site size of N = 25, the 2^N × 2^N Q-matrices of the full model exceed the memory limitations of modern workstations; consequently, the probability distribution for N_𝒪 and the Score must be estimated from Monte Carlo simulation. Because the N + 1 × N + 1 Q-matrices of the mean-field approximation are comparatively small, direct matrix analytic methods can be used to calculate these response measures (Appendix B) as well as spark statistics such as duration, interspark interval, and frequency (Appendix C).

In this matrix analytic approach it is convenient to reduce the number of parameters of the mean-field model via nondimensionalization. Accordingly, we express Ca²⁺ concentrations in units of the dissociation constant of Ca²⁺ binding (K where K^η = k⁻/k⁺) and denote the nondimensional Ca²⁺ coupling strength and background [Ca²⁺] as and respectively. Substituting and into Eq. 23 and expressing time in units of the reciprocal of the dissociation rate constant (1/k⁻), we arrive at the dimensionless generator matrix After nondimensionalizing, the nine parameters of the mean-field model (N, η, ɛ_𝒞𝒞, ɛ_𝒪𝒪, k⁺, k⁻, c_∞, c_*, and a_*) are reduced to seven dimensionless parameters (N, η, ɛ_𝒞𝒞, ɛ_𝒪𝒪, and a_*).

Using the for 25 mean-field coupled RyRs, Fig. 12 shows spark duration, interspark interval, and spark frequency (grayscale) as a function of the strength of dimensionless Ca²⁺ coupling () and allosteric interactions that stabilize closed (ɛ_𝒞𝒞) and open (ɛ_𝒪𝒪) channel pairs. Each panel explores a slice of this three-dimensional parameter space indicated by the shaded region of the cubes shown at the left. These correspond to allosteric interactions that stabilize closed channel pairs (Fig. 12 A, ɛ_𝒞𝒞 < 0, ɛ_𝒪𝒪 = 0), open channel pairs (Fig. 12 B, ɛ_𝒞𝒞 = 0, ɛ_𝒪𝒪 < 0), and both in a balance fashion (Fig. 12 C, ɛ_𝒞𝒞 = ɛ_𝒪𝒪 < 0). Spark statistics are only shown when sparks are present (Score > 0.3).

Note that similar to simulations using the full model (Figs. 6 and 7), the magnitude and range of values that result in sparks increase as the strength of allosteric interactions that stabilize closed channel pairs increases (Fig. 12 A) and decreases as the strength of allosteric interactions that stabilize open channel pairs increases (Fig. 12 B). The magnitude and range of values that result in sparks does not vary significantly as the magnitude of ɛ_𝒞𝒞 and ɛ_𝒪𝒪 increases in a balanced fashion (Fig. 12 C). Regardless of how stabilizing allosteric interactions are introduced, spark duration and interspark interval are increasing and decreasing functions of respectively. In Fig. 12, A–C, spark duration increases and interspark interval decreases in such a manner that spark frequency at first increases but ultimately decreases as a function of

While similar changes of spark statistics are seen as increases regardless of how stabilizing allosteric interactions are included, qualitatively different changes are observed in Fig. 12, A–C, as the strength of allosteric interactions increases. Fig. 12 A shows that increasing the strength of allosteric interactions that stabilize closed channel pairs decreases spark duration and increases interspark interval. Fig. 12 B shows that increasing the strength of allosteric interactions that stabilize open channel pairs increases spark duration but has little effect on interspark interval. Fig. 12 C shows that increasing the strength of allosteric interactions that stabilize both closed-closed and open-open channel pairs in a balanced fashion decreases interspark interval, while spark duration is largely unaffected. While in Fig. 12 C spark frequency is a decreasing function of the strength of allosteric interactions, in Fig. 12, A and B, spark frequency may increase, decrease, or both, depending on the coupling strength Fig. 12, A–C, is qualitatively unchanged when the dimensionless background [Ca²⁺] () is doubled or halved (not shown).

Allosteric coupling and Ca²⁺ spark generation and termination

A significant result of this study is the observation that synchronizing allosteric interactions always promote Ca²⁺ sparks (i.e., result in a higher Score) for some value of the strength of Ca²⁺ coupling (c_*), regardless of whether synchronizing allosteric interactions stabilize closed channel pairs, open channel pairs, or both (see Figs. 6 and 7). When the strength of Ca²⁺ coupling is sufficiently large to preclude termination of simulated sparks, allosteric interactions that stabilize closed channel pairs can promote spark termination. Similarly, allosteric interactions that stabilize open channel pairs facilitate spark initiation when Ca²⁺ coupling is too weak to mediate stochastic Ca²⁺ excitability. Sparks are less sensitive to variations in c_* when the strength of allosteric interactions that stabilize closed channel pairs is increased, and more sensitive to c_* when the strength of allosteric interactions that stabilize open channel pairs is increased.

Allosteric coupling washout, cardiac dysfunction, and Ca²⁺ spark statistics

A substantial body of experimental evidence demonstrates that normal cardiac function requires the association of the 12.6 kDa FK506 binding protein FKBP12.6 to the RyR channel complex (55–58). For example, pharmacological or exercise-induced PKA hyperphosphorylation of RyRs has been shown to substantially dissociate FKBP12.6 from RyRs and has been linked to increased frequency of ventricular arrhythmias and sudden cardiac death (58,59). In addition, the absence of FKBP12.6 in knockout mice has been associated with increased systolic [Ca²⁺] and cardiac hypertrophy (60).

While the connection between FKBP12.6 depletion and cardiac dysfunction is not clearly established, evidence that FK-binding proteins are responsible for coupled gating of RyRs suggests that organ-level failure may be inherited from defects in the collective gating of channels leading to irregularities in the dynamics of Ca²⁺ sparks. In striated (skeletal and cardiac) and smooth muscle, both the frequency and duration of spontaneous sparks increase upon knockout of genes encoding relevant FK-binding proteins or treatment with FK506 or rapamycin, two drugs that physically and/or functionally dissociate FK-binding proteins from RyRs (17,49,50,60–64). Conversely, overexpression of FKBP12.6 has been shown to decrease spark frequency (51). Interestingly, these experimentally observed changes in spark duration and frequency are consistent with simulated washout of allosteric interactions that stabilize closed-closed channel pairs or both closed-closed and open-open channel pairs, but inconsistent with simulations involving the washout of allosteric interactions that stabilize only open-open channel pairs (Figs. 8 and 9). While in principle these different types of allosteric coupling could leave a signature in the distribution of spark durations, this does not appear to be the case for the minimal two-state RyR model used here (Fig. 8). While these simulations aim to clarify how changes in spark statistics due to pharmacological washout of the accessory proteins mediating allosteric coupling may indicate the type of synchronizing allosteric interactions exhibited by physically coupled RyRs, it is unclear the degree to which the results will generalize to more complicated and realistic RyR models (see below).

The mean-field approximation for allosteric interactions

The mean-field approximation formulated in this study is applicable to a cluster of RyRs coupled via both Ca²⁺ and allosteric interactions. Although this reduced model has a drastically contracted state space compared to full model simulations, the mean-field coupled RyRs exhibit Ca²⁺ sparks that are qualitatively similar to sparks of the full model (Fig. 11). However, for mean-field simulations involving a fixed number of channels and fixed allosteric coupling parameters, the Ca²⁺ coupling strength (c_*) that results in the highest Score is slightly elevated compared to the optimal c_* of corresponding full model simulations. This difference becomes more evident as the number of channels at release sites increases (Fig. 11), and may be a consequence of the spatial spread of activation or the clustering of open channels in full model simulations.

The mean-field reduction formulated here is analogous to the sticky cluster model of Sobie et al. (40) where the coupled gating of RyRs is represented by multiplying the 𝒞 → 𝒪 and 𝒪 → 𝒞 transition rates by cooperativity factors (χ_𝒪 and χ_𝒞) that depend on the number of open and closed channels in the cluster. For example, in Sobie et al. (40) the death rates are given by where

(27)

and the scaling factor k_coop sets the strength of RyR coupling. By inspecting the death rates presented in this article (Eq. 26), one finds that the cooperativity factor in the mean-field model is

(28)

which when expressed in terms of is

(29)

Note that Eq. 27 is an increasing function of N_𝒞, consistent with Eq. 29, when ɛ_𝒪𝒪 + ɛ_𝒞𝒞 < 0, as in most of the simulations presented here. However, Eq. 29 is a nonlinear function of N_𝒞 (Eq. 27 is linear), and the scaling factor for the strength of allosteric coupling (a_*) enters Eq. 29 differently than k_coop in Eq. 27. Furthermore, χ′_c = 1 when N_𝒞 = 0 and ɛ_𝒪𝒪 = 0 (and when and ɛ_𝒞𝒞 = 0) regardless of the strength of allosteric coupling (not so for χ_c in Eq. 27). While Eq. 27 has only one free parameter (k_coop), we would recommend using Eq. 29 because the three parameters (a_*, ɛ_𝒪𝒪, ɛ_𝒞𝒞) are not post hoc additions to an N + 1 state model, but rather derived from the microscopic parameters of the 2^N state Ca²⁺ release site that is reduced to N + 1 states using the mean-field approximation. Equation 29 has the additional advantage of being able to incorporate synchronizing (or desynchronizing) allosteric interactions that stabilize (or destabilize) closed channel pairs, open channel pairs, or both. Perhaps most importantly, the comparatively small state space of mean-field coupled RyR clusters could be used to mitigate against the difficulties inherent in realistic multiscale modeling of cardiac myocyte excitation-contraction coupling (65–68).

Generalizing the mean-field approximation

Although the single-channel model used in this article includes only two states (closed and open), the mean-field approximation can be applied to clusters of channels with more complicated single-channel dynamics that include mechanisms suspected to contribute to Ca²⁺ spark dynamics in situ such as luminal regulation, Ca²⁺-dependent inactivation, or adaptation (15,69,70). For N M-state channels there are n-choose-k{N + M – 1, N} states in the mean-field approximation, each of which can be expressed as a vector of the form (N₁, N₂, ···, N_M) where N_m is the number of channels in state m, 1 ≤ m ≤ M, and If the current state of the release site is (N₁, N₂, ···, N_M) and a channel makes an i → j transition, the transition rate is N_ik_ijχ_ij and the appropriate cooperativity factor is

(30)

where ν_ij is the previously encountered coefficient that partitions allosteric coupling between the forward and reverse transitions (0 ≤ ν_j ≤ 1 and ν_ji = 1 – ν_ij), and δ_ki is the Krönecker delta function defined by

(31)

Limitations of the model

A potential limitation of this study is the assumption of instantaneous coupling via the local [Ca²⁺]. Theoretical studies of two-state Ca²⁺-activated channels coupled via a time-dependent Ca²⁺ microdomain demonstrate that the time constant of Ca²⁺ domain formation and collapse can affect the dynamics of puffs and sparks (29,32). For example, slow domain formation can make the triggering of sparks less likely while slow domain collapse can prohibit the termination of Ca²⁺ release events. On the other hand, Ca²⁺ release via clusters of RyRs in ventricular myocytes occurs within dyadic clefts, spatially restricted regions of the cytosol located between the sarcolemma of T-tubules and the sarcoplasmic reticulum membrane (4,71,72). Theoretical studies indicate that the time constant of Ca²⁺ domain formation decreases as the volume of a dyad decreases and may be <1 ms (73,74), while the decay of elevated [Ca²⁺] to background levels after termination of release may require 10s of milliseconds due to low affinity binding sites on the cytosolic face of the sarcolemma (74). Thus, in the context of Ca²⁺ release via RyR clusters in ventricular myocytes, the assumption of instantaneous coupling is more justified during the rising phase of Ca²⁺ release events than during the falling phase. Our prior work (29,30,32) suggests that this feature of the modeling formalism will increase the likelihood that all the open RyRs will close simultaneously, a mechanism referred to as stochastic attrition (15,71).

While the analysis of this article is simplified by assuming instantaneous Ca²⁺ coupling and a minimal two-state RyR, the lack of an explicit mechanism for spark termination—e.g., depletion of luminal Ca²⁺, Ca²⁺-dependent inactivation, or adaptation—results in sparks that terminate exclusively via stochastic attrition. Consequently, sparks of physiologically realistic durations are only observed over a finite range of Ca²⁺ coupling strengths, even when allosteric interactions are included. While allosteric interactions that stabilize closed channel pairs may potentiate spark termination via stochastic attrition when the Ca²⁺ coupling strength (c_*) is elevated (resulting in sparks that are less sensitive to c_*; see Figs. 6 A and 7 A), stabilizing allosteric interactions between closed channels do not result in robust termination of sparks at all Ca²⁺ coupling strengths. Taken as a whole, our simulations demonstrate that allosteric interactions may facilitate spark generation, and are often sufficient for spark termination in the absence of another mechanism such as depletion of luminal Ca²⁺ or Ca²⁺-dependent inactivation. When the strength of Ca²⁺ coupling is not optimal, the strength of allosteric coupling can usually be adjusted to yield robust Ca²⁺ sparks (Fig. 12). On the other hand, for fixed allosteric coupling parameters, the range of Ca²⁺ coupling strengths leading to robust sparks was never observed to be greater than 25% of the optimal Ca²⁺ coupling strength.

While many buffers with various binding kinetics, affinities, and diffusion constants contribute to the landscape of [Ca²⁺] in vivo, the mathematical representation of the Ca²⁺ microdomain used in this article assumes a single Ca²⁺ buffer. Because the single-channel model does not include mechanisms that would promote spark termination, high buffer concentrations are required to achieve Ca²⁺ coupling strengths that allow sparks to spontaneously terminate via stochastic attrition. For example, when the RyR is modeled with dissociation constant K_d = 5 μM and unitary current of i_Ca = 0.04 pA, simulations that do not include allosteric interactions require [B]_T ≈ 1.2 mM to achieve the optimal Ca²⁺ coupling strength of c_* ≈ 0.48 μM. As shown in Fig. 2 C, this coupling strength can be obtained using a variety of different values for [B]_T or i_Ca; as expected, simulations using lower buffer concentrations with lower unitary current yield results that are similar to Fig. 6. When allosteric interactions stabilizing closed channel pair are included (ɛ_𝒞𝒞 = −0.4, squares of Fig. 6 A), the optimal coupling strength of c_* ≈ 0.71 μM corresponds to a total buffer concentration of [B]_T ≈ 570 μM. Utilization of complex RyR gating schemes and explicit modeling of the depletion of luminal Ca²⁺ would likely decrease the total buffer concentration required for spark termination.

Perhaps the most significant limitation of this study is that the degree to which the results will generalize to more complicated and realistic RyR models is unknown. This concern is ever present when minimal single-channel models that reproduce select features of Ca²⁺-regulation are used to study the collective gating that gives rise to Ca²⁺ sparks (28–38). Although beyond the scope of this article, it might be possible to extend inference methods commonly used in conjunction with single-channel recording (75–77) to the collective gating of mean-field coupled intracellular channels. In this way, experimentally observed statistics of sparks (e.g., the shape of the distribution of spark durations and interspark intervals) might be used to distinguish between channels that are coupled via local [Ca²⁺], allosteric interactions, or both.

For now, the generalization of our results to other single-channel models can only be addressed on a case-by-case basis. For example, in the absence of allosteric interactions, instantaneously coupled two-state RyRs do not exhibit Ca²⁺ sparks unless the cooperativity of Ca²⁺ binding is two or more (η ≥ 2) (30,31). Similarly, a preliminary survey of all possible three-state single-channel models that include unimolecular Ca²⁺ binding suggests that multiple Ca²⁺-binding transitions are required for sparks (not shown). However, when stabilizing allosteric interactions are included, cooperative Ca²⁺ binding is no longer required, that is, η = 1 can yield robust sparks (not shown).

Our validation of the mean-field approach to modeling allosteric interactions suggests that studies utilizing more realistic RyR models could be performed using the coupling factors (Eq. 30) that are derived here for the first time. Our attempts at this further analysis include simulations of mean-field coupled three-state RyRs that include a long-lived closed state (R),

(32)

These simulations demonstrate that both Ca²⁺-dependent and Ca²⁺-independent inactivation often reduce the sensitivity of sparks to variations in the coupling strength (78). In preliminary studies we have found that stabilizing allosteric interactions can further extend the range of c_* values that result in robust sparks (not shown). However, it remains to be determined whether the statistics of Ca²⁺ sparks can ever be used to rule out allosteric coupling as a synchronization mechanism.