Type-II ryanodine receptor channels (RYRs) play a fundamental role in intracellular Ca2+ dynamics in heart. The processes of activation, inactivation, and regulation of these channels have been the subject of intensive research and the focus of recent debates. Typically, approaches to understand these processes involve statistical analysis of single RYRs, involving signal restoration, model estimation, and selection. These tasks are usually performed by following rather phenomenological criteria that turn models into self-fulfilling prophecies. Here, a thorough statistical treatment is applied by modeling single RYRs using aggregated hidden Markov models. Inferences are made using Bayesian statistics and stochastic search methods known as Markov chain Monte Carlo. These methods allow extension of the temporal resolution of the analysis far beyond the limits of previous approaches and provide a direct measure of the uncertainties associated with every estimation step, together with a direct assessment of why and where a particular model fails. Analyses of single RYRs at several Ca2+ concentrations are made by considering 16 models, some of them previously reported in the literature. Results clearly show that single RYRs have Ca2+-dependent gating modes. Moreover, our results demonstrate that single RYRs responding to a sudden change in Ca2+ display adaptation kinetics. Interestingly, best ranked models predict microscopic reversibility when monovalent cations are used as the main permeating species. Finally, the extended bandwidth revealed the existence of novel fast buzz-mode at low Ca2+ concentrations.

Introduction

Type II calcium release RYR channels play a fundamental role in the intra cellular Ca2+ signaling dynamics of cardiac muscle cells that govern contractility of the heart. RYR channels are activated by a small Ca2+ influx through the surface membrane by a process known as Ca2+-induced Ca2+ release (CICR). The CICR process, an inherently self-regenerating process, is precisely controlled in cells. The regenerative feedback that counters the positive feedback of CICR is thought to depend on the cytosolic Ca2+ concentration (Fabiato, 1985). Since then, the cytosolic RYR Ca2+ regulation has been the subject of extensive research (Györke and Fill, 1993; Laver and Curtis, 1996; Schiefer et al., 1995; Sitsapesan et al., 1995; Zaradnková et al., 1999; Copello and Fill, 2002), and resent matter of debate (Fill et al., 2000; Lamb et al., 2000; Sitsapesan and Williams 2000).

Several models of single RYR channel gating have been developed to explain the RYRs role in the overall Ca2+ dynamics of cardiac muscle cells. These models, however, have been usually chosen and estimated by following rather phenomenological criteria in order to reproduce experimentally observed properties such as the Ca2+-dependent activation, inactivation (Zaradniková and Zahradník, 1995; Stern et al., 1999), adaptation (Cheng et al., 1995; Sachs et al., 1995; Keizer and Levine, 1996; Villaba-Galea et al., 1998; Fill et al., 2000), and modal gating—as well as the local RYR-mediated intracellular Ca2+ release events known as Ca2+ sparks (Jafri et al., 1998; Sobie et al., 2002; Villaba-Galea et al., 2002). Even in the presence of high quality data, model estimation and selection have only rarely been statistically addressed. For example, Saftenku et al. (2001) present a maximum likelihood approach that is subject to several drawbacks. The raw data is idealized by using standard threshold methods at a relatively limited bandwidth (i.e., 2 kHz). Second, models for particular visually identified gating modes are then separately estimated via maximum likelihood (Qin et al., 1996). This analysis heavily relies on the quality of the initial idealizations and the accuracy of the mode identification steps. This makes it subject to errors that are very difficult to statistically assess. Even more importantly, it does not reveal the connections between the identified modes and thus presents only a partial picture of the salient dynamics with limited predictive power.

This article presents a first step toward a statistical analysis of single type II RYRs at steady-state. The analysis provides some substantial benefits when compared with previously used approaches. First, the applicable bandwidth of the analysis is extended up to 10 kHz by the use of hidden Markov models (HMMs; see Michalek et al., 1999; Venkataramanan and Sigworth, 2002). Second, statistical inferences are made using a combination of Bayesian statistics and stochastic search methods (Markov chain Monte Carlo [MCMC]). These techniques have proven to be particularly well suited for complex estimation problems where the more commonly applied maximum likelihood–based approaches fail (Robert and Casella, 1999; Liu, 2001). The MCMC method not only provides a rich description of the modeling process but also a direct indication of where and why a particular model fails. More importantly, it allows us to directly assess the errors incurred at each estimation step, since all the estimates used throughout are endowed with a central limit theorem. The MCMC technique used here, namely a Gibbs sampler, is a generalization of the method described in Rosales et al. (2001) that allows us to consider the constraints imposed by the aggregation of states into conductance classes and the set of state interconnections. These techniques are similar in spirit to those of Ball et al. (1999) and Hodgson (1999); however, here the continuous gating of the channel is approximated by a discreetly sampled process, producing a less computationally intensive algorithm that is based on a solid probabilistic theory.

16 gating models were compared on this basis by analyzing three datasets containing the steady-state activity of single RYRs at 1, 10, and 100 μM Ca2+. The results obtained provide direct evidence for the existence of different Ca2+-dependent gating modes. The open/closed states associated with each mode were also identified. Additionally, gating schemes were used to predict the relaxation kinetics of single RYR channel response to Ca2+ step stimuli. The kinetic predictions for the top ranked schemes agree well with those that have been experimentally defined, showing that single RYRs respond transiently to sudden Ca2+ steps. Our results also demonstrate that gating follows detailed balance at the experimental conditions so far considered.

Materials And Methods

Microsome Preparation

Heavy microsomes enriched in type II RYRs obtained from canine ventricular cardiac muscle (Tate et al., 1985) were reconstituted in planar lipid bilayers (Miller and Racker, 1976). Briefly, the tissue was kept in a saline solution (154 mM NaCl, 10 mM Tris-malate, pH 6.8) at a temperature of 4°C before being chopped and then homogenized. Heavy microsomes were obtained by differential centrifugation and then kept at −70°C in saline solution with 300 mM sucrose until used.

Reconstitution and Recording

Planar lipid bilayers were formed across a 150 μm hole in a Delrin cup. Bilayers were obtained with a mixture (50 mg/ml in decane) of phosphatidylcholine (PC) and phosphatidylethanolamine (PE) (Avanti Polar Lipids, Inc.) in a 7:3 relation. Bilayers with an electrical capacitance of 200–400 pF were used in our experiments. The microsomes were added to one side of the bilayer defined as cis. The other side, defined as trans, was held at virtual ground. The standard solution was 20 mM CsCH3SO3, 20 μM CaCl2, 20 mM HEPES, pH 7.4. The fusion of the microsomes was promoted by mechanical stirring and an osmotic gradient (i.e., addition of 400 mM CsCH3SO3). The orientation of an incorporated RYR was such that their cytoplasmic side was in the cis compartment (i.e., defined by sidedness of ATP sensitivity). After channel incorporation, the trans CsCH3SO3 concentration was adjusted to 420 mM. The charge carrier was Cs+. All experiments were done at room temperature.

A patch clamp amplifier (Axopatch 200B; Axon Instruments, Inc.) was used to optimize single channel recording (Györke and Fill, 1993). The acquisition software (pClamp 6.0; Axon Instruments, Inc.) was run on a Pentium III computer controlling a 12-bit analogue/digital-digital/analogue converter (Axon Instruments, Inc.). Single-channel data were digitized at 100 kHz and filtered at 32 kHz. To reduce posterior computational effort, the data was digitally Gaussian filtered down to 10 kHz and sampled at 25 kHz. Long single RYR channel recordings (3×106 samples each) were made at three different steady-state free Ca2+ concentrations (1, 10, and 100 μM cis). Free Ca2+ concentration was verified using a Ca2+ electrode.

Data Analysis

A representative single-channel recording at each Ca2+ concentration was selected for analysis. The baseline was corrected to account for any drift and other irregularities. The files were analyzed separately with each of the 16 gating models shown in Fig. 1 

. The data was modeled by considering constrained HMMs and a combination of Bayesian statistics and MCMC methods. The most relevant aspects of these techniques together with a brief mention to further methods used are described bellow. A detailed account of the Bayes/MCMC framework is given in Rosales (2004).

Constrained HMMs

Observations were modeled by using constrained HMMs. Briefly, the gating dynamics of single RYR channels were assumed to follow a discreetly time sampled homogeneous Markov process with a finite state space E = {1, 2, …, n}. States were further aggregated into two conductance classes O = {O1, …, Ov} and C = {C1, …, Cw}, with v + w = n and Oi and Ci denoting single states, hereafter labeled open and closed. Transitions between the states of this process are not directly observed. The data arising from bilayer experiments is inherently corrupted by additive noise, assumed o be independent and Gaussian distributed as a first approximation. The constraints on the HMM parameter space considered here are those induced by the topology of the associated graph of a gating mechanism and those inherent to the association of states into C and O. The parameters involved are the transition probabilities pij, grouped into a transition probability matrix P = [pij], for any i, j ϵ O, or C, the class conductance levels μC and μO, the class noise variances σ2C and σ2O, and finally the initial density for the hidden process Λ = [λi], i ϵ{O, C}. For convenience hereafter we write θ = (P, μC, μO, σ2C, σ2O, Λ).

Bayesian Treatment

Within this framework, the identification of a gating model from the raw data leads to the following well posed statistical tasks: (A) signal idealization, (B) parameter estimation, and (C) model selection. Statistical inferences concerning C are considered by using an asymptotic form of penalized ratio test, namely the Bayesian information criterion, whereas tasks A and B are considered from a Bayesian perspective. Let y denote a set of N observations y1, y2, … yN, and z a particular realization z1, z2, …, zN of the underlying process, i.e., z represents a sequence of states visited by the channel. Furthermore, let Z stand for the space of all possible state sequences, and Θ a subset of the Euclidean space formed by the set of all possible values for θ. For each gating mechanism denote by L(θ, z) the likelihood function,

 
\[L\left(\mathbf{\mathit{{\theta}}},\mathbf{z}\right)=P\left(y,\mathbf{z}|\mathbf{\mathit{{\theta}}}\right)=P\left(y|\mathbf{z},\mathbf{\mathit{{\theta}}}\right)P\left(\mathbf{z}|\mathbf{\mathit{{\theta}}}\right)\mathrm{,}\]
(1)

defined as the joint probability for the occurrence of the data and a hidden realization, given a particular value for the parameters θ. Let ϕ(θ) denote the prior, a density that expresses belief or knowledge about the parameters before the observations have been examined. Then, following Bayes theorem, the joint posterior density for both unknowns θ and z, is

 
\[{\pi}\left(\mathbf{\mathit{{\theta}}},\mathbf{\mathit{z}}|\mathit{y}\right)=\mathit{L}\left(\mathbf{\mathit{{\theta}}},\mathbf{\mathit{z}}\right){\phi}\left(\mathbf{\mathit{{\theta}}}\right)\mathit{w}^{{-}\mathrm{1}},\]
(2)

with

 
\[w={{\sum}_{z{\in}Z}}{{\int}_{\mathrm{{\Theta}}}}L\left(\mathbf{\mathit{{\theta}}},\mathbf{\mathit{z}}\right){\phi}\left(\mathbf{\mathit{{\theta}}}\right)d\mathrm{\mathbf{\mathit{{\theta}}}.}\]

Estimates for model parameters, task B, are based on this density by integrating out z, for instance,

 
\[E_{\mathrm{{\pi}}}\left[\mathbf{\mathit{{\theta}}}\right]={{\sum}_{z{\in}Z}}{{\int}_{\mathrm{{\Theta}}}}\mathrm{\mathbf{\mathit{{\theta}}}{\pi}}\left(\mathbf{\mathit{{\theta}}},\mathbf{\mathit{z}}|\mathit{y}\right)d\mathrm{\mathbf{\mathit{{\theta}}},}\]
(3)
 
\[\mathrm{Var}_{\mathrm{{\pi}}}\left[\mathbf{\mathit{{\theta}}}\right]={{\sum}_{z{\in}Z}}{{\int}_{\mathrm{{\Theta}}}}\left(\mathbf{\mathit{{\theta}}}{-}E_{\mathrm{{\pi}}}\left[\mathbf{\mathit{{\theta}}}\right]\right)^{2}\mathrm{{\pi}}\left(\mathbf{\mathit{{\theta}}},\mathbf{\mathit{z}}|\mathit{y}\right)d\mathrm{\mathbf{\mathit{{\theta}}},}\]
(4)

constitute unbiased estimates for the model parameters and their associated variances. The same argument also applies for estimates of the hidden sequence z by integration of θ, which then represents a solution to task A.

MCMC

Although there are clear theoretical advantages for the use of posterior based inferences over maximum likelihood based approaches in this setting, none of the integrals in Eqs. 24 can be evaluated directly. This occurs even for the simplest gating mechanism and regardless of the choice of ϕ(θ), due the complex structure of the HMM marginal likelihood

\({\sum}_{\mathit{z}{\in}\mathit{Z}}\)
L
\(\left(\mathit{\mathbf{{\theta}},\mathbf{z}}\right)\)
. Thus, the required integrals have to be evaluated numerically. This task is performed here by using MCMC methods, more precisely by using the Gibbs sampler initially described by Rosales et al. (2001) and further modified by Rosales (2004) in order to accommodate class constraints. The basic principle behind any MCMC algorithm is to perform Monte Carlo integration by drawing samples of the required posterior in Eq. 2, and then evaluating sample averages to approximate the expectations in Eqs. 3 and 4. More precisely, in the current setting, the Gibbs sampler generates a sequence of random variables (θu, zu), u = 1,2, …, m, by running a Markov chain for a long time, such that in the limit as m→ +∞, the pair (θm, zm) is asymptotically distributed according to the target distribution π. In fact, for the constrained HMM outlined above, the Gibbs sampler converges toward π from any arbitrary starting point (θ1, z1). These and other regularity properties of the sampler have been formally established in Rosales (2004).

An example of the output generated by analyzing a small data segment with model M3 of Fig. 1 is shown in Fig. 2 

. This example shows how the sampler performs a random walk on two of the model parameters θ, namely on σO22 and σC22. After a few steps the walk finally stays within a boundary that corresponds to a region of maximum posterior probability for (σO22C22). This figure also presents the output generated for the hidden process shown as a randomly chosen set of idealized traces, zu, each one generated at a given iteration of the MCMC run. Well identifiable features are actually recovered in almost all iterations, whereas very brief events appear only in some of the sampled realizations.

Formally, due to convergence toward π, Monte Carlo estimates for the posterior mean and standard deviation of θ (Eqs. 3 and 4) are, respectively, the empirical averages of the obtained sequence (θu), i.e.,

 
\[\overline{\mathrm{{\mu}}}_{\mathrm{{\pi}}}\left(\mathbf{\mathit{{\theta}}}\right)=\frac{1}{G}{{\sum}_{u=b}^{m}}\mathbf{\mathit{{\theta}}}^{u}{\approx}E_{\mathrm{{\pi}}}\left[\mathbf{\mathit{{\theta}}}\right]\]
(5)
 
\[\overline{\mathrm{{\sigma}}}_{\mathrm{{\pi}}}\left(\mathbf{\mathit{{\theta}}}\right)=\left[\frac{1}{G}{{\sum}_{u=b}^{m}}\left(\mathbf{\mathit{{\theta}}}^{u}\right)^{2}{-}\overline{\mathrm{{\mu}}}_{\mathrm{{\pi}}}^{2}\left(\mathbf{\mathit{{\theta}}}\right)\right]^{{1}/{2}}{\approx}\mathrm{Var}_{\mathrm{{\pi}}}^{{1}/{2}}\left[\mathbf{\mathit{{\theta}}}\right]\mathrm{,}\]
(6)

for some 1 < b < m and G = mb + 1.

The use of these equations may be illustrated by the following example. Let Q* = [kij* for i,j ϵ{O,C} be a transition rate matrix for the process z estimated from the data and associated to P as

 
\[P=\mathrm{exp}\left(Q\mathrm{*{\delta}}\right)\mathrm{,}\]
(7)

with exp as the matrix exponential and δ as the sampling period of the acquisition system. It should be emphasized that Q* is not the infinitesimal generator, Q = [kij], of the hidden process, but rather a rate matrix associated to the limiting time resolution δ. Hereafter, however, for notational convenience we write Q for Q* and kij for k*ij, the elements of Q*. If the continuous process is sampled fast enough, i.e., if

 
\[\mathrm{{\delta}}^{{-}1}\ \mathrm{{<}{<}}\ sup_{i,j\mathrm{{\epsilon}}\left\{O,C\right\}}k_{ij}\mathrm{,}\]

then, following standard theory (see for instance Eq. 7 of Colquhoun and Hawkes, 1995a),

 
\[k_{ii}{\approx}{-}{{\sum}_{j{\neq}i}}{p_{ij}}/{\mathrm{{\delta}}}\mathrm{,}k_{ij}{\approx}\left({p_{ij}}/{\mathrm{{\delta}}}\right)\mathrm{.}\]
(8)

A posterior estimate for the observed transition rate matrix and its standard deviation follow immediately from these equations and Eqs. 5 and 6, namely

 
\[\begin{array}{l}\overline{\mathrm{{\mu}}}_{\mathrm{{\pi}}}\left(k_{ii}\right)={-}\frac{1}{\mathrm{{\delta}}G}{{\sum}_{u=b}^{m}}{{\sum}_{j{\neq}i}}p_{ij}^{u}\mathrm{,}\\\overline{\mathrm{{\mu}}}_{\mathrm{{\pi}}}\left(k_{ij}\right)=\frac{1}{\mathrm{{\delta}}G}{{\sum}_{u=b}^{m}}p_{ij}^{u}\mathrm{,}\end{array}\]
(9)

and

 
\[\begin{array}{l}\overline{\mathrm{{\sigma}}}_{\mathrm{{\pi}}}\left(k_{ii}\right)=\left[\frac{1}{\mathrm{{\delta}}^{2}G}{{\sum}_{u=b}^{m}}\left({{\sum}_{j{\neq}i}}p_{ij}^{u}\right)^{2}+\overline{\mathrm{{\mu}}}_{\mathrm{{\pi}}}\left(k_{ii}\right)^{2}\right]^{{1}/{2}}\mathrm{,}\\\overline{\mathrm{{\sigma}}}_{\mathrm{{\pi}}}\left(k_{ij}\right)=\left[\frac{1}{\mathrm{{\delta}}^{2}G}{{\sum}_{u=b}^{m}}\left(p_{ij}^{u}\right)^{2}{-}\overline{\mathrm{{\mu}}}_{\mathrm{{\pi}}}\left(k_{ij}\right)^{2}\right]^{{1}/{2}}\mathrm{.}\end{array}\]
(10)

Convergence of the sampler, and hence estimation of b in Eqs. 5 and 6, was assessed by the Raftery-Lewis and Heidelberger-Welch tests, both implemented in R (http://cran.r-project.org) or SPlus by the BOA program (Bayesian Output Analysis program; Version 1.0.0 for S-PLUS and R; http://www.public-health. uiowa.edu/boa/). MCMC runs were started at an arbitrary initial point θ1 ϵΘ for each model. The same values were used for any of the three Ca2+ conditions and each model (see Fig. 1). The initial hidden realization, z1, was drawn from the initial density Λ1 and the initial probability matrix P1. A weak informative prior for θ was considered (see Rosales et al., 2001).

Density Estimates

Posterior dwell-time density estimates based on the MCMC sampled realizations, (zu), for each model were computed by using optimal nonparametric kernel density estimation methods instead of standard log-binned histograms (Rosales et al., 2002). Where indicated, dwell-time densities were also calculated directly from the Q matrix estimates (Colquhoun and Hawkes, 1995b). Interclass state correlations were visually inspected by computing the dependency-difference plot introduced by Magleby and Song (1992). Estimates for the posterior parameter densities were also computed by using kernel methods and the MCMC output (θu). In this case estimates were obtained by following Bowman and Azzalini (1997).

Model Choice

Models are compared by using the Bayes information criterion (BIC), see for instance Gelfand and Dey (1994), which introduces a penalty according to the total number of free parameters and the length of the data record. For any mechanism this is given by

 
\[{-}2\mathrm{ln}\left[L\left(\mathbf{\mathit{{\theta}}}^{\mathrm{{\cdot}}}\right)\right]+2N_{\mathrm{{\theta}}}\mathrm{ln}\left(N\right)\mathrm{,}\]
(11)

with Nθ as the total number of associated free parameters and L as the marginal likelihood,

 
\[L\left(\mathbf{\mathit{{\theta}}}^{\mathrm{{\cdot}}}\right)={{\sum}_{z{\in}Z}}L\left(\mathbf{\mathit{{\theta}}}^{\mathrm{{\cdot}}},z\right)\mathrm{,}\]

evaluated at some θ· that maximizes its value. Here we have chosen the posterior mean estimate, i.e., θ· = μ̅π

\(\left(\mathit{\mathbf{{\theta}}}\right)\)
⁠. From Eq. 1,

 
\begin{eqnarray*}&&L\left(\mathbf{\mathit{{\theta}}}^{\mathrm{{\cdot}}}\right)=L\left(\overline{\mathrm{{\mu}}}_{\mathrm{{\pi}}}\left(\mathbf{\mathit{{\theta}}}\right)\right.\ \\&&\left.\ \right)={{\sum}_{z{\in}Z}}L\left(\overline{\mathrm{{\mu}}}_{\mathrm{{\pi}}}\left(\mathbf{\mathit{{\theta}}}\right),z\right)={{\sum}_{z{\in}Z}}P\left(y|z,\overline{\mathrm{{\mu}}}_{\mathrm{{\pi}}}\left(\mathbf{\mathit{{\theta}}}\right)\right)P\left(z|\overline{\mathrm{{\mu}}}_{\mathrm{{\pi}}}\left(\mathbf{\mathit{{\theta}}}\right)\right)\mathrm{.}\end{eqnarray*}
(12)

This quantity is efficiently computed, once μ̅π

\(\left(\mathit{\mathbf{{\theta}}}\right)\)
is obtained from the MCMC runs, by using Baum's forward procedure (Baum et al., 1970). Finally, if the data y is partitioned into s independent segments of lengths N1, N2, …, Ns and the likelihood for the ith segment is denoted by Li, then from Eq. 11 the BIC takes the form

 
\[BIC={-}2{{\sum}_{i=1}^{s}}\mathrm{ln}\left[L_{i}\left(\overline{\mathrm{{\mu}}}_{\mathrm{{\pi}}}\left(\mathbf{\mathit{{\theta}}}\right)\right)\right]+2N_{\mathrm{{\theta}}}{{\sum}_{i=1}^{s}}\mathrm{ln}\left(N_{i}\right)\mathrm{.}\]
(13)

A model is preferred over another if its BIC as defined by Eq. 13 is smaller.

Evolution Toward Equilibrium

Let η(kδ) = [ηi(kδ)], i ϵ{C,O}, be a n-dimensional row vector denoting the state occupation density at any time k=1,2, …, for δ > 0. Relaxation toward the equilibrium distribution, η(∞), is given by

 
\[\mathbf{{\eta}}\left(\mathrm{{\infty}}\right)=\begin{array}{l}\mathrm{lim}\\k{\rightarrow}\mathrm{{\infty}}\end{array}\mathbf{{\eta}}\left(k\mathrm{{\delta}}\right)=\begin{array}{l}\mathrm{lim}\\k{\rightarrow}\mathrm{{\infty}}\end{array}\mathbf{{\eta}}\left(0\right)P^{k\mathrm{{\delta}}}\mathrm{,}\]
(14)

with η(0) as an arbitrary initial distribution (see theorem 1.13 of Norris, 1997). This relation describes the evolution for the δ-discretisation of the continuous time process (see theorem 2.1.1 of Norris, 1997). In fact, the continuous time version for Eq. 14 is standard and may be found in the current context at several places by setting P = exp(Qkδ) (see Eq. 19 of Colquhoun and Hawkes, 1995b). The open probability Po as a function of kδ, is thus obtained via Eq. 14 by adding the entries in η that involve open states, i.e.

 
\[P_{o}\left(k\mathrm{{\delta}}\right)={{\sum}_{i{\in}O}}\mathrm{{\eta}}_{i}\left(k\mathrm{{\delta}}\right)\mathrm{.}\]
(15)

In this study, however, the scope of Eq. 15 is restricted to the evaluation of η(kδ) by using only one of the estimated Q at a time, representing one of the Ca2+ conditions considered experimentally.

Online Supplemental Material

The MCMC sampler used in this study was written in ANSI C and is freely available under the GNU license as online supplemental material . Code is distributed with documentation concerning compilation and running examples. Compilation under various platforms including SUN Sparks running Solaris 8, and several Pentiums with Linux and most Windows versions was made with gcc (http://gcc.gnu.org). Precompiled binaries may be available upon request by writing to the first author. Code and documentation for dwell-time kernel estimates is also available as online supplemental material. Kernel estimates for the parameter posterior marginals were computed from the MCMC output by using the sm library implemented in R and conforming to Bowman and Azzalini (1997). Simple plots were made with gnuplot (http://www.gnuplot.info) and others including surfaces with R. All Q matrix computations were made by using Maple 7 (Waterloo Maple, Inc.). When run on the dataset at 10 μM Ca2+ (3.051 × 106 samples) with scheme M3 (9 states and 29 free parameters), each iteration of the sampler took 17.8 s on a Dell Precision 330 with a 1.4 GHz Pentium 4 processor.

Results

Model Comparison

Results obtained by analyzing single RYR channel data collected at three different Ca2+ concentrations with the gating schemes shown in Fig. 1 are summarized in Table I 

. This table presents the values for the negative of the marginal log likelihood evaluated at the posterior mean estimates for the associated parameters at each Ca2+ concentration. Models are ordered from top to bottom according to their BIC ranking, as determined by Eq. 13. Simple inspection suggests that the information in the data strongly supports model complexity, as shown by the difference of 20,787.47 BIC units between the best ranked model (scheme M3) and the worst one (scheme M15). Complexity in this context is related to the total number of free parameters. The following general topological properties of the graph associated to a mechanism may be inferred from Table I:

(1) RYR channel gating schemes with three open states are preferred over those with two open states. Two open state schemes are preferred over those with one open state. Three open states would be consistent with the RYR channel experimental results published by Sitsapesan and Williams (1994) and Saftenku et al. (2001).

(2) RYR gating schemes with two communicating open states are preferred, see for example M1 and M2.

(3) RYR gating schemes with a hexagonal cycle of states have higher BIC values than those with smaller cycles.

(4) The best RYR gating schemes with three open states, except for scheme M13, share a common motif that includes states O2, O3, C4, C5, and C6. Schemes, such as M5 and M6, lacking this structure are ranked low.

(5) RYR gating schemes with more than two consecutive closed states not communicating with an open state at the left of the cycle such as C1, C2, and C3 in M4 and M5 are also heavily penalized.

(6) Table I also suggest the existence of more than four pathways between the open and the closed class. Interestingly, more than one pathway is necessary to explain modal RYR channel gating, an experimentally observed feature of single RYR channels (Zaradniková and Zahradník, 1995; Copello et al., 1997).

Results and predictions for the top ranked models with three open and six closed states are very similar and thus will not be detailed separately here. Instead, we will turn our attention below to the best ranked model, namely RYR gating scheme M3.

Q Matrix Estimates

Estimates for the Q matrices at each Ca2+ concentration constitute the basis for all further descriptions and predictions of the dynamics of single RYR channel behavior here. Fig. 3 

displays the sampled values of the transition rates between communicating states of mechanism M3 for each Ca2+ condition, plotted against the MCMC iteration. These rates were obtained from the sampled probabilities by using Eq. 8. This figure also presents the sampled values for the variance of the open and the closed levels, σO2 and σC2, respectively.

Simple inspection suggests that rate estimates for 1 and 100 μM Ca2+ have converged after ∼1,000 iterations. Convergence takes perhaps longer for the data at 10 μM Ca2+ (see Fig. 3, top right). Convergence at each concentration was quantitatively assessed by the Raftery-Lewis and Heidelberger-Welch tests (see Table II)

. The latter being more conservative than the first one indicates a convergence period lag of 1,500 iterations. Based on this information, the first 2,500 iterations were discarded and only the remaining 500 iterations were considered for further inferences.

The posterior mean and standard deviation for each rate, computed by using the last 500 samples via Eqs. 9 and 10, are shown in Table III 

. Few posterior marginal density estimates for some individual parameters are shown in Fig. 4 

. They were computed directly from the MCMC output according to Bowman and Azzalini (1997).

In general, the rates leaving short-lived states with small occupancy probability (see Tables V and IV below) have relatively high variances. An example is provided by C6 at 1 μM Ca2+. Rates leaving long-lived states with larger occupancy probabilities, such as C1 at 1 μM Ca2+, have relatively smaller variances.

The presence of multiple maxima in the posterior marginals (Fig. 4, A and C), and thus almost surely on the respective likelihood projections, should warn against the naive use of maximum likelihood methods for HMMs. This may include the use of standard deviations obtained from the likelihood curvature at “the maximum” (see Qin et al., 2000).

Single-state Ca2+-induced Changes

Straightforward computations with the Q matrix estimates lead to the associated stationary probabilities of being at each state and the expected time spent in each state. These quantities, shown in Tables IV and V 

, when considered together allow us to draw a detailed picture for the Ca2+ induced changes. For this, the next section starts with a description of the RYR channel gating modes present at 1 μM Ca2+ and follows then by identifying the states involved in each mode.

Gating modes are clearly visible at 1 μM Ca2+. Fig. 5 

displays two representative single RYR channel records illustrating at least three gating modes. The first of them, labeled Gl, consists of long openings interrupted by brief closures. The second, named Gz, consists of brief openings followed by brief closures and finally the third mode, named Go, is formed by brief openings followed by long shut periods. As shown in Fig. 5, open sojourns in Gl and Gz are grouped in clusters forming burst-like behavior. Gl has previously been identified in single RYR channel studies as a mode of high open probability named H gating mode, and Go as a mode with low probability named L gating mode (Zaradniková and Zahradník, 1995, 1996; Armisén et al., 1996). Gz is not mentioned in these earlier references. The states participating in each of the Gl, Gz, and Go modes and the closed periods between them may be identified directly by inspecting the last MCMC sampled realizations for the hidden process (zu). However, as mentioned above, this may also be done by considering the single state time constants and occupation probabilities (see Tables V and IV) obtained from the Q matrix estimates.

At 1 μ Ca2+ the channel spends most of its time at the C1 and C2 states, which account for the long closed sojourns between the Gl and Gz modes. From the C2 state, the channel may then enter the O3 state generating a brief opening. From the O3 state, it may enter the C4 state forming events corresponding to the Go gating mode. Also from the C2 state, via C3, the RYR channel may also (very rarely) enter the O1 sate generating a long lasting opening. These long openings may be followed by transitions toward the O2 state (brief openings) and subsequently to the C4 or C6 states (brief closings) forming bursts of channel activity characteristic of the Gl gating mode. From either the O1 or O2 states, Gl bursts are terminated either via O2C4O3C2 or O1C3C2. The appearance of the Gz gating mode is associated with transitions of the form O3C4O2C6, where the C6 state is responsible for the very fast closed characteristic of this mode. From the rates leaving C2 shown in Table III, there are 23 more transitions per unit of time toward the C3 state, the intermediate state for the initiation of a Gl gating mode bursts, than those for the O3 state that starts the Go gating mode. Despite this, Gl gating mode bursts are quite rare when compared with the occurrence of the Go gating mode. This is due to an extremely low rate for the transition from the C3 to O1 state. The appearance of the Gz or Go modes is governed by the rates leaving the state C4 to the O2 and O3 states, respectively. Note that transitions to the O3 state are twice as frequent as those to the O2 state. This explains the apparent excess of the occurrence of Go mode over Gz bursts. These results suggest the following gating mode frequency relations, Gl Gz Go. At this Ca2+ concentration the open probability for the channel is 0.06599. A schematic view of the modes described above and their participating states is summarized in Fig. 6 

.

At 10 μM Ca2+ the channel tends to leave C1 and C2, which are shorter than at 1 μM Ca2+, and starts spending more time in states C3 and C4, both intermediately lived. Consequently, the inter -burst time is reduced, and this makes burst-like behavior less apparent. The channel also spends more time in all three open states, elevating the open probability to 0.37983. At this Ca2+ concentration, the preferred open state is now O2. This state is shorter than the most probable open state at 1 μM Ca2+ (O1), which was responsible for long open sojourns during Gl bursts. This accounts for a considerable reduction in the total Gl open time. Bursts of type Gz are ∼9 times more frequent than sojourns in the Go mode, whereas the proportion between the Gl and Go gating modes is more equally distributed in contrast to the situation at 1 μM Ca2+.

Finally, at 100 μM Ca2+ the channel now prefers to stay at C4, a briefly lived closed state, and mostly frequents two intermediate brief open states, O2 and O3. This accounts for a decrease in the open probability, which is now 0.24098, and also to the disappearance of long openings and closings between bursts.

Dwell Time Density Estimates

Apart from the mean times and the probabilities at each state, further insights into the Ca2+-induced changes in single RYR channel activity may be gained by inspection of the relative frequencies of visits to particular states. These are shown in Fig. 7 

, which presents the log dwell-time densities for the sojourns in the open and closed classes together with their single-state components. Overall densities and components were computed by using the spectral decomposition of Q (Colquhoun and Hawkes, 1995b) and the estimates shown in Table III. At 1 μM Ca2+ both open and closed densities show two clearly defined maxima. Fast events in the open class (Fig. 7 A) constituted by sojourns in O2 and O3 are about five times more frequent than slow sojourns in O1. A high frequency of transitions to O3 accounts for most openings in the Go mode. The closed class (Fig. 7 B) presents quite frequent fast and slow events corresponding to sojourns in C6 and C1, C5, respectively. The presence of multiple maxima in both classes at 10 μM Ca2+ is not as apparent as for 1 Ca2+. This accounts for a partial disappearance of modal behavior. At this Ca2+ concentration (10 μM) the maximum for the open class is slightly shifted toward the right (Fig. 7 C) when compared with the mayor maximum at 1 Ca2+. Furthermore, O3 is now the most frequent open state, being then closely followed by O2 and O1. The most frequent states for the closed class (Fig. 7 D) are C6 and C4. Provided the gating modes have still the same meaning as at 1 μM Ca2+, both plots at 10 μM Ca2+ suggest the following mode frequency relations Gl < Go < Gz. Finally at 100 μM Ca2+ (Fig. 7, E and F) the densities for both classes have essentially a single well-defined maximum. This corresponds to a complete disappearance of modal behavior. At this concentration (100 μM Ca2+), the single-state frequency proportions are roughly the same as for 10 μM Ca2+. These results suggest that a mayor change in the single state frequencies occurs only around 1–10 μM Ca2+.

The univariate single-class densities suggest the existence of Ca2+-dependent RYR channel modal gating. However, direct evidence for this is only provided by the adjacent open-closed sojourn correlations. These are presented in Fig. 8 

, which displays a MCMC posterior estimate for the log open-closed dwell time densities and their derived dependency-difference plots. Dwell-time densities were computed by using only the last five MCMC sampled realizations of the hidden process. Use of more MCMC realizations did not appreciably improve these estimates. A gating mode corresponds (in principle) to the existence of a peak (or maximum) in the joint density estimate. Inspection of Fig. 8 (top) suggests the existence of two modes at 1 Ca2+. These modes appear to gradually merge into a single mode at progressively higher Ca2+ concentrations (10–100 μM Ca2+). Such a simple visual inspection of the joint density plots can be misleading. Specific correlations between communicating states may be shown instead via dependency-difference plots. At 1 μM Ca2+Fig. 8 B suggests the existence of at least three maxima (yellow to faint yellow contoured regions in the range of 0.0 to 0.04). The first of them corresponds to fast open-closed events associated to the Gz mode and is located at the bottom left corner of the plot. Another clearly identified maximum corresponds to long openings followed by brief closures (Gl mode). This is located in the region that starts at (2,0) at stretches horizontally to about (5,0). A third relatively minor maximum corresponds to intermediate openings followed by somewhat longer closings (Go mode). This one may be identified with a region that starts at (0,0) and stretches vertically up to about (0,4).

This analysis also suggests there is a major change in modal gating characteristics between 1 to 10 μM Ca2+. At 10 Ca2+ the dependency-difference plots show evidence of two maxima instead of three (Fig. 8 D). One of these corresponds to brief openings followed by intermediate length closures. The other corresponds to intermediate openings followed by short closures. This pattern is also evident at 100 μM Ca2+ (Fig. 8 F). However, correlations at 100 μM Ca2+ are much smaller (note contour level scales). This analysis presents direct evidence for RYR channel modal gating and reveals how it changes at different Ca2+ concentrations.

Filter and Threshold at 2 kHz

This section presents a standard filter and threshold analysis at 2 kHz. The purpose of this exercise is to illustrate the significance of the extended bandwidth analysis obtained by using the HMM approach. The data used for the MCMC analysis was first filtered down to 2 kHz with a digital Gaussian filter. The filtered data was then idealized by using a standard half amplitude threshold. The resulting dwell-time list was then finally used for dwell-time kernel density estimation. Fig. 9 

presents the single-class estimates, whereas Fig. 10 

shows the joint density for adjacent open-closed sojourns at each Ca2+ concentration.

Both figures show several differences when compared, respectively, to Figs. 7 and 8. Simple inspection shows that changes are quite dramatic at 1 μM Ca2+. Indeed, Fig. 9, A and B, shows that fast events in both the open and the closed classes have disappeared almost entirely at 2 kHz. As a consequence both densities are shifted to the right. This effect is more markedly appreciated in the closed class (Fig. 9 B), since at 10 kHz this density presents one of the fastest components identified with sojourns in state C6 (see Fig. 7 B and also Table V). Fig. 9 B also presents smaller probability values than Fig. 7 B for the slowest events. This effect is due to noise events that still cross the threshold at 2 kHz and break long shut periods. A comparison of the joint density in Fig. 10 A against the one in Fig. 8 A also suggests the disappearance of fast adjacent open-closed sojourns that form the Gz gating mode. This is confirmed by comparison of the dependence-differences in Figs. 10 B and 8 B. The disappearance of Gz modifies the correlation pattern observed at 1 μM Ca2+, which is actually inverted when compared against the one in Fig. 8 B. Fig. 10 B presents two regions of positive correlations that correspond to intermediary fast and slow open-closed sojourns, much like the pattern found for 10 and 100 Ca2+ (see Figs. 10, D and E, and 8, D and E).

Reversibility

Reversibility can be evaluated by the Kolmogorov's criterion for Markov chains with cycles (see Eq. 1.22 in Kelly, 1979) using the Q matrix estimates. Let a1, a2, …, an−1, an denote the states involved in the cycle of a particular gating mechanism. Gating is reversible if the transition rates in the cycle satisfy

 
\[k_{a_{1}a_{2}}\mathrm{{\ldots}}k_{a_{n{-}1}a_{n}}=k_{a_{2}a_{1}}\mathrm{{\ldots}}k_{a_{n}a_{n{-}1}}\mathrm{.}\]
(16)

For the cycle of scheme M3, the clockwise product of the six rates involved at 1 μM Ca2+ is 6.9 × 10−6. The product of the six rates in the opposite direction product is 6.6 × 10−6. For 10 μM Ca2+ the products are 1.5 × 10−5 and 1.8 × 10−5 respectively. For 100 μM Ca2+ these products are 2.1 × 10−3 and 2.3 × 10−3. These results present direct evidence for the reversibility of single RYR channel gating dynamics at the experimental conditions tested here. They are also consistent with Fig. 11 

. It should be emphasized that none of the cyclic models used so far considered reversibility constraints for the values of the transition probabilities pij. The results of these tests thus provide strong evidence in favor of time-reversible gating as a genuine feature of the data.

Response to Sudden Ca2+ Steps

The response of RYR channels to fast Ca2+ changes can be evaluated by computing the relaxation toward equilibrium of the probability of being in the open class, Po. Here, we present the response of mechanism M3 to instantaneous changes in the Ca2+ concentration from 0 to 1, 10, or 100 μM Ca2+. The predicted response of a single RYR channel was computed using Eq. 15, and considering several initial densities that represent possible candidates of the resting state. These were chosen by assigning different probability values to the closed states and 0 to any open state. For all the Ca2+ steps and initial conditions considered, Fig. 12 

shows that the Po rapidly peaks and the slowly relaxes toward a new equilibrium. This pattern of single RYR channel behavior has been experimentally observed (Györke and Fill, 1993; Valdivia et al., 1995). Closer inspection of the plots in Fig. 12 shows that the rising phase of the responses is strongly dependent on the resting condition.

MCMC Method Reliability

To test the reliability of the MCMC method, simulated single-channel data was obscured by adding random noise. θ and z were then estimated and compared to their original values. To this end, two synthetic datasets were generated assuming two different Ca2+ conditions. One synthetic dataset was generated using Q matrix estimates for the M3 model at 1 μM Ca2+. The specific rate constants used for the M3 model are presented in Table III. The second synthetic dataset was generated using Q-matrix estimates for the M10 model at 10 μM Ca2+. The specific rate constants used for the M10 model are presented below (in ms−1).

 
\[\begin{array}{l}\mathrm{k}_{\mathrm{C1C2}}=\mathrm{1.372,k}_{\mathrm{C2C1}}=\mathrm{0.759}\\\mathrm{k}_{\mathrm{C2O1}}=\mathrm{0.017,k}_{\mathrm{O1C2}}=\mathrm{0.025}\\\mathrm{k}_{\mathrm{C2O2}}=\mathrm{2.598,k}_{\mathrm{O2C2}}=\mathrm{0.329}\\\mathrm{k}_{\mathrm{C3O1}}=\mathrm{0.146,k}_{\mathrm{O1C3}}=\mathrm{1.362}\\\mathrm{k}_{\mathrm{C3O2}}=\mathrm{0.34,k}_{\mathrm{O2C3}}=\mathrm{0.025}\\\mathrm{k}_{\mathrm{C3C4}}=\mathrm{0.021,k}_{\mathrm{C4C3}}=\mathrm{0.18}\end{array}\]

Synthetic datasets with N = 2 × 106 samples were generated by sampling at δ = 0.04 ms. The following parameters were set to mimic those estimated from real data; μC = 0.9985, μO = 2.6784 and σ2C = σ2O = 0.5955. The synthetic dataset were then processed using the Gibbs sampler. The Gibbs sampler was run for 3,000 iterations with different gating schemes. Convergence in each case was diagnosed as described in Section 2.3.3. Only the last 1,000 iterations were used for estimation.

For the synthetic dataset generated using the M3 model, the Gibbs sampler was run with four different gating schemes (M3, M2, M7, and M14). The actual Q-matrix samples for the M3 model are shown in Fig. 13 

. This figure also contains results obtained from the real single channel data for comparison (also see Fig. 3).

The BIC/2 factors as given by Eq. 13 for M3, M2, M7, and M14 are, respectively, 3,227,097.713, 3,227,131.317, 3,227,155.629, and 3,798,851.394. This means that the rank order of best-fit models was M3 > M2 > M7 > M14. In the M3 model case, the estimates generated for μC, μO, and σ2C and σ2O were 0.9979, 2.679, and 0.5961, respectively. The estimates for the rate constants (in ms−1) are listed below.

 
\[\begin{array}{l}\mathrm{k}_{\mathrm{C1C2}}=\mathrm{1.801}+{\pm}\mathrm{0.193,k}_{\mathrm{C2C1}}=\mathrm{1.241}+{\pm}\mathrm{0.219}\\\mathrm{k}_{\mathrm{C2C3}}=\mathrm{0.570}+{\pm}\mathrm{0.127,k}_{\mathrm{C3C2}}=\mathrm{3.223}+{\pm}\mathrm{0.472}\\\mathrm{k}_{\mathrm{C2O3}}=\mathrm{0.041}+{\pm}\mathrm{0.003,k}_{\mathrm{O3C2}}=\mathrm{0.565}+{\pm}\mathrm{0.03}\\\mathrm{k}_{\mathrm{C3O1}}=\mathrm{0.001}+{\pm}\mathrm{0.001,k}_{\mathrm{O1C3}}=\mathrm{0.004}+{\pm}\mathrm{0.002}\\\mathrm{k}_{\mathrm{C4C5}}=\mathrm{2.41}+{\pm}\mathrm{0.957,k}_{\mathrm{C5C4}}=\mathrm{0.678}+{\pm}\mathrm{0.071}\\\mathrm{k}_{\mathrm{C4O2}}=\mathrm{0.248}+{\pm}\mathrm{0.077,k}_{\mathrm{O2C4}}=\mathrm{0.888}+{\pm}\mathrm{0.05}\\\mathrm{k}_{\mathrm{C4O3}}=\mathrm{0.524}+{\pm}\mathrm{0.209,k}_{\mathrm{O3C4}}=\mathrm{0.692}+{\pm}\mathrm{0.034}\\\mathrm{k}_{\mathrm{O1O2}}=\mathrm{0.039}+{\pm}\mathrm{0.004,k}_{\mathrm{O2O1}}=\mathrm{0.163}+{\pm}\mathrm{0.02}\\\mathrm{k}_{\mathrm{C6O2}}=\mathrm{20.927}+{\pm}\mathrm{0.47,k}_{\mathrm{O2C6}}=\mathrm{2.496}+{\pm}\mathrm{0.074}\end{array}\]

For the synthetic dataset generated using the M10 model, the Gibbs sampler was also run for with four different gating schemes M1, M3, M10, and M11. The rate constant estimates (in ms−1) are listed below.

 
\[\begin{array}{l}\mathrm{k}_{\mathrm{C1C2}}=\mathrm{1.31,k}_{\mathrm{C2C1}}=\mathrm{0.703}\\\mathrm{k}_{\mathrm{C2O1}}=\mathrm{0.021,k}_{\mathrm{O1C2}}=\mathrm{0.021}\\\mathrm{k}_{\mathrm{C2O2}}=\mathrm{2.59,k}_{\mathrm{O2C2}}=\mathrm{0.325}\\\mathrm{k}_{\mathrm{C3O1}}=\mathrm{0.121,k}_{\mathrm{O1C3}}=\mathrm{1.377}\\\mathrm{k}_{\mathrm{C3O2}}=\mathrm{0.34,k}_{\mathrm{O2C3}}=\mathrm{0.024}\\\mathrm{k}_{\mathrm{C3C4}}=\mathrm{0.022,k}_{\mathrm{C4C3}}=\mathrm{0.167}\end{array}\]

The BIC/2 factors for the M1, M3, M10, and M11 models are 2,571,612.52, 2,571,723.112, 223,01,276.171 and 345,101.05, respectively. This means that the rank order of best-fit model was M10 > M11 > M1 > M3.

In both cases, the estimated rate constants are similar to the original values used to generate the simulated datasets. Thus, the HMM-MCMC analysis performed reasonably well under noise conditions similar to those found experimentally. Further, the BIC values show that the method was also able to accurately choose the correct model.

Reproducibility of the Best-fit RYR Channel Model Prediction

Model selection was done using real single RYR channel data collected at three different cytosolic Ca2+ concentrations. Here, the reproducibility of the model selection is evaluated using real single RYR channel recordings at a single cytosolic Ca2+ concentration (100 μM). These recordings were collected from two different RYR channels. For each recording, the Gibbs sampler was used to generate estimates of rate constants and BIC values. Two complete rate prediction and BIC value sets from this analysis are provided here to illustrate the degree of channel-to-channel variability present. Rate estimates obtained from the analysis of the previously presented data collected at 100 μM Ca2+ are presented in Table VI 

. Estimates obtained from the analysis of another long (6.84 × 106 samples) single RYR channel recording at 100 μM Ca2+ are presented in Table VII 

. The experimental conditions in which both of these recording were done were exactly the same.

These tables contain the rate constant estimates and associated standard deviations for only three of the top ranked models (M2, M3, and M13). These rate estimates where computed by taking only the last 500 of 3,000 iterations of the Gibbs sampler. Convergence in each case was diagnosed as described in Section 2.3.3. Particular rate estimates obtained from the two channels are not identical. As expected, the natural biological variability between channels does impose variability in the absolute values of the predicted time constants. However, the relative magnitudes of the particular rate constants are relatively consistent between the channels. To illustrate the consistency in magnitude, the coefficient of variation between the rate constants estimations for all transitions in one of the top ranked models (M2) are plotted in Fig. 14 

.

The BIC values obtained from the different analysis sets were used to rank the models. The BIC values for each analysis set are presented in Tables I and VIII 

. In each case, the BIC values were then used to rank order the best-fit models for each analysis set. The best-fit rank orders for the two analysis sets are compared in Table VIII. The rank orders from both analysis sets contain the sequence M4 > M12 > M9 > M10 > M16 >M14 > M15. More importantly, the top five best-fit models were M2, M3, M5, M6, and M13 although not in exactly the same order. Note that the BIC value differences between these top five ranked models in most cases are quite small.

Discussion

Detailed knowledge about the dynamics of single RYR channel Ca2+ regulation is required to understand the role these molecules play in the local control of excitation-contraction coupling (Stern, 1992). Despite extensive research, many questions regarding the mechanisms that govern RYR channel gating remain unanswered. Several Markov models for the gating of single RYRs have been proposed. Typically, models have been selected and parameters estimated rather heuristically, usually by imposing a strong set of constraints that turn them into essentially a self-fulfilling prophecy (see Sachs et al., 1995; Zaradniková and Zahradník, 1996; Villaba-Galea et al., 1998; Sobie et al., 2002). These approaches end up simply showing that the selected model can reproduce a particular RYR channel attribute (e.g., modal gating, adaptation, and/or inactivation) with little (if any) statistical evaluation. The principal objective here is to infer single RYR properties directly from data by setting a minimum number of constraints on the model structure. By doing so, we have considered sixteen gating models, six of them previously reported in the literature. Statistical inferences were made from Bayesian perspective by using MCMC. The next sections summarize the implications of our findings and discuss their limitations.

Extended Bandwidth

The HMM approach applied here extends the analysis of single-channel data to relatively high bandwidths. This is particularly important for fast gating channels like the RYRs (Zaradniková et al., 1999). Here, the HMMs were applied to analyze single RYR channel data at 10 kHz. This analysis was directly compared with standard filter and threshold at 2 kHz. Although the same datasets were analyzed, there was a dramatic difference in outcome (see Figs. 69). Filtering at 2 kHz removed a considerable proportion of very fast events, shifted distributions, and obscured (in some cases eliminated) observed correlations between adjacent open-closed events. The filtering distortion was most marked at low Ca2+ levels (i.e., 1 μM) and generated a completely different pattern of single-channel gating. This is an important finding. It introduces new information concerning the physiological relevance of fast events of single RYR channels at low Ca2+ concentrations specifically.

Reversibility

Both the surface plots at Fig. 11 and the Kolmogorov's criterion for cycles for the Q matrix estimates show that RYR channels gate in detailed balance. It should be noted that this result is not in contradiction with recent evidence that supports irreversible gating when Ca2+ is the main diffusing species (Wang et al., 2002). Furthermore, this agrees and reinforces results of Rengifo et al. (2002) and Villaba-Galea et al. (2002), which propose the dissipation of the Ca2+ gradient as the source of free energy that keeps the system out of time reversibility.

Modal Gating and Adaptation

The steady-state Ca2+ sensitivity of single RYR channels is well established (Zaradniková and Zahradník, 1995, 1996; Armisén et al., 1996). The dynamics of single RYR channel Ca2+ regulation is not (Copello and Fill, 2002). It is clear from numerous studies that single RYR channels display rapid activation in response to a fast Ca2+ elevation (Györke and Fill, 1993; Zaradnková et al., 1999). In some studies, channel activity peaks and then spontaneously decays following the fast Ca2+ elevation. The spontaneous decay was originally termed adaptation (Györke and Fill, 1993) and was particularly interesting because it occurred at Ca2+ concentrations that do not inactivate the RYR channel under steady-state conditions. The spontaneous decay (i.e., adaptation) has now been attributed to a Ca2+ and time dependent modal gating shift (Zaradniková and Zahradník, 1995, 1996; Villaba-Galea et al., 1998; Györke, 1999; Zaradniková et al., 1999; Fill et al., 2000).

Interestingly, the highest ranked gating model selected from steady-state single-channel data predicts a spontaneous decay after a fast-step Ca2+ elevation (see Fig. 12). This is interesting because there has been vigorous debate about whether or not this is possible. In the original adaptation studies (Györke and Fill, 1993), the applied Ca2+ elevation had a very fast (1 ms) overshoot at its leading edge. At first, this very fast Ca2+ overshoot was not considered to play a role in the much slower spontaneous decay (Velez et al., 1997). However, the possibility that the fast Ca2+ overshoot may actually generate the slow spontaneous decay was quickly forwarded (Lamb and Stephenson, 1995; Lamb et al., 2000; Sitsapesan and Williams, 2000). Here, a spontaneous decay occurs in the absence of fast Ca2+ overshoot. The analysis here supports the notion that this spontaneous decay is generated by a Ca2+- and time-dependent modal gating shift (see below). The analysis here also suggests that a fast Ca2+ overshoot may indeed accelerate the initial transition into the open configuration, but will not have an impact on single RYR channel behavior after that.

As mentioned above, single RYR channels display rapid activation in response to a fast Ca2+ elevation (Zaradnková et al., 1999). There is no dispute concerning this RYR channel attribute. The model with highest BIC score, selected here from steady-state data, predicts a Ca2+ activation rate similar to those that have been experimentally defined (Györke and Fill, 1993; Valdivia et al., 1995; Velez et al., 1997). Not surprisingly, the Ca2+ activation rate depends on the initial conditions that exist before the fast Ca2+ change.

Several studies have detailed the Ca2+-dependent modal gating of single RYR channels (Armisén et al., 1996; Zaradniková and Zahradník, 1996; Zaradnková et al., 1999). Again, the highest ranked model predicts the existence of Ca2+-dependent modal gating (see Fig. 8). The analyses here show that the Ca2+ dependence of modal gating is steepest at relatively low Ca2+ concentrations (≈1 μM). This corresponds to the Ca2+ concentration range over which adaptation is experimentally observed (Fill et al., 2000). This also supports the notion that adaptation (i.e., the spontaneous decay) is due to a Ca2+- and time-dependent modal gating shift. As initially proposed by Zaradniková and Zahradník (1995) the Ca2+-dependent modal shift hypothesis rests essentially on two basic assumptions: (a) an increase of Ca2+ should favor transitions out of states with high Po, but (b) at resting conditions Ca2+ can only bind to transitions that lead to these states. These conditions still remain to be fully demonstrated. The analysis in Zaradniková et al. (1999) and the one presented here are consistent with the first point; however, none of them provides the means for the identification of the specific Ca2+-dependent transitions.

MCMC Method Reliability and Robustness of the Fit

The reliability of the MCMC method was tested by simulating two different synthetic datasets using the parameters fitted for two different models under two different cytosolic Ca2+ concentrations. Then noise was added to recordings simulated from the Q matrices and the HMM MCMC formalism was used for analyze the simulated data (Section 3.8). Excitingly, the analysis again correctly retrieves the model used to produce the data for the synthetic datasets.

To test the Robustness of our methods, two different channels under similar ionic conditions (100 μM free cytosolic Ca 2+ concentration) were analyzed using the HMM MCMC formalism. In Section 3.9 the model selection for two different channels and the rate constants are compared. Although the model selection rankings are note identical, the best first five models were the same for all the channels. It is not surprising that the ranking for different channels is not identical, mostly because the natural biological variability between channels. Nevertheless, BIC values for these top models are very similar. In addition, when the fitted Q matrices for the same model (M2) and two different datasets were compared (Fig. 14) the rate constants for each transition were not exactly the same but were in the same order. This difference could be explained again by the variability of the biological samples.

Future Directions

A key limitation of the approach applied here is the estimation of a separate Q matrix for each Ca2+ condition. The simultaneous estimation of a single common matrix should significantly improve its statistical properties; moreover, this should lead to the identification the Ca2+-dependent transitions of a gating scheme. The latter would certainly clarify several aspects in the Ca2+-dependent modal shift hypothesis. Finally, a single Q matrix with well-identified Ca2+-dependent transitions would enclose all the necessary information for predicting the Ca2+-dependent dynamic properties of single RYR channels. This includes the response to any imaginable Ca2+ waveform allowing comparison of predicted channel behavior with the myriad of experimental results that have been reported (Györke and Fill, 1993; Sitsapesan et al., 1995; Laver and Curtis, 1996; Zaradnková et al., 1999). More importantly, it would allow prediction of single RYR channels to the nonstationary Ca2+-regulatory environment that exists in living cells.

To make inferences by considering data at several Ca2+ concentrations simultaneously one has to be able introduce constraints for the rates, such that the involved states become well-defined occupancy sites for Ca2+. For instance, one would generally like to set kij = αij[Ca2+], where αij is the rate constant that results in the case [Ca2+] = 1 μM, if kij is to be expressed in μM × s−1. Unfortunately, these constraints do not translate in a simple manner to P and hence into the standard HMM formulation. The problem may be illustrated by considering an expansion of the exponential form which relates P and Q (see Eq. 7), namely,

 
\[P=\mathrm{exp}\left(\mathit{Q}{\delta}\right)=I+\mathit{Q}{\delta}+{\left(\mathit{Q}{\delta}\right)^{2}}/{2}+{\left(\mathit{Q}{\delta}\right)^{3}}/{3}+\mathrm{{\ldots},}\]
(17)

with I as the identity matrix. P may well be obtained from Q by simple matrix operations, but a simple change in one entry of Q takes a rather complicated algebraic form in P. Work in progress, however, shows that it is indeed possible to sample Q directly by using hybrid MCMC moves that include Metropolis-Hastings steps (Robert and Casella, 1999; Liu, 2001) inside the Gibbs sampler used here. This approach is highly promising, both in theory and in practice since it borrows many of the properties of the Gibbs sampler so far used.

In any case, unifying views of the Ca2+-regulatory dynamics of RYRs will certainly continue to stimulate the research of thorough quantitative analytical approaches and innovative experimental strategies.

Acknowledgments

We thank Dr. Carlos Villalba-Galea for helpful discussion. We also thank Ing. Mariano Reyes and Ing. Lorena Masine for their invaluable collaboration during the preparation of the manuscript.

Supported by TTUHSC SEED Grant to A. Escobar. and NIH-HL57832 to M. Fill.

Olaf S. Andersen served as editor.

References

References
Armisén, R., J. Sierralta, P. Vélez, D. Naranjo, and B.A. Súarez-Isla.
1996
. Modal gating in neuronal and skeletal muscle ryanodine-sensitive Ca2+ release channels.
Am. J. Physiol.
271
:
C144
–C153.
Ball, F.G., Y. Cai, J.B. Kadane, and A. O'Hagan.
1999
. Bayesian inference for ion channel gating mechanisms directly from single channel recordings, using Markov chain Monte Carlo.
Proc. Roy. Soc. Lond. A.
455
:
2879
–2932.
Baum, L.E., T. Petrie, G. Soules, and N. Weiss.
1970
. A maximization technique ocurring in the statistical analysis of probabilistic functions of Markov chains.
Ann. Math. Stat.
41
:
164
–171.
Bowman, A.W., and A. Azzalini. 1997. Applied Smoothing Techniques for Data Analysis. Oxford Statistical Science Series, 18. Clarendon Press, Oxford, UK.
Cheng, H., M. Fill, and H. Valdivia.
1995
. Models of Ca2+ release channels adaptation.
Science.
267
:
2009
–2010.
Colquhoun, D., and A.G. Hawkes. 1995a. Principles of the stochastic interpretation of ion-channel mechanisms. Single-channel Recording. B. Sackmann and E. Neher, editors. Plenum, New York. 397–482.
Colquhoun, D., and A.G. Hawkes. 1995b. A Q-Matrix cookbook. Single-channel Recording. B. Sackmann and E. Neher, editors. Plenum, New York. 589–633.
Copello, J.A., S. Barg, H. Onoue, and S. Fleischer.
1997
. Heterogeneity of Ca2+ gating of skeletal muscle and cardiac ryanodine receptors.
Biophys. J.
73
:
141
–156.
Copello, J.A., and M. Fill.
2002
. Ryanodine receptor calcium release channels.
Physiol. Rev.
82
:
893
–922.
Fabiato, A.
1985
. Time and calcium dependence of activation and inactivation of calcium-induced release of calcium from the sarcoplasmic reticulum of a skinned canine cardiac Purkinje cell.
J. Gen. Physiol.
85
:
247
–289.
Fill, M., A. Zahradnková, C.A. Villalba-Galea, I. Zahradnk, A.L. Escobar, and S. Györke.
2000
. Ryanodine receptor adaptation.
J. Gen. Physiol.
116
:
873
–802.
Gelfand, A.E., and D.K. Dey.
1994
. Bayesian model choice: asymptotics and exact calculations.
J. Roy. Stat. Soc. Ser. B.
56
:
501
–514.
Györke, S.
1999
. Ca2+ spark termination: inactivation and adaptation may be manifestations of the same mechanism.
J. Gen. Physiol.
114
:
163
–166.
Györke, S., and M. Fill.
1993
. Ryanodine receptor adaptation: control mechanism of Ca2+-induced Ca2+ release in heart.
Science.
260
:
807
–809.
Hodgson, M.E.A.
1999
. A Bayesian restoration of an ion channel signal.
J. Roy. Stat. Soc. B.
61
:
95
–114.
Jafri, M.S., J.J. Rice, and W.L. Winslow.
1998
. Cardiac Ca2+ dynamics: the roles of ryanodine receptor adaptation and sarcoplasmic reticulum load.
Biophys. J.
74
:
1149
–1168.
Keizer, J., and L. Levine.
1996
. Ryanodine receptor adaptation and CICR-dependent Ca2+ oscillations.
Biophys. J.
71
:
3477
–3487.
Kelly, F. 1979. Reversibility and Stochastic Networks. Wiley Series in Probability and Mathematical Statistics. John Wiley & Sons, New York.
Lamb, G.D., D.R. Laver, and D.G. Stephenson.
2000
. Questions about adaptation in ryanodine receptors.
J. Gen. Physiol.
116
:
883
–890.
Lamb, G.D., and D.G. Stephenson.
1995
. Activation of ryanodine receptors by flash photolisis of caged Ca2+.
Biophys. J.
68
:
946
–948.
Laver, D.R., and B.A. Curtis.
1996
. Response of ryanodine receptor channels to Ca2+ steps produced by rapid solution changes.
Biophys. J.
71
:
732
–741.
Liu, J.S. 2001. Monte Carlo Strategies in Scientific Computing. Springer Verlag, New York.
Magleby, K.L., and L. Song.
1992
. Dependency plots suggest the kinetic structure of ion channels.
Proc. R. Soc. Lond. B. Biol. Sci.
249
:
133
–142.
McManus, D.B., and K.L. Magleby.
1991
. Accounting for the Ca2+-dependent kinetics of single large-conductance Ca2+-activated K+ channels in rat skeletal muscle.
J. Physiol.
443
:
739
–777.
Michalek, S., H. Lerche, M. Wagner, N. Mitrovic, M. Schiebe, F. Lehmann-Horn, and J. Timmer.
1999
. On identification of Na+ channel gating schemes using moving-average filtered hidden Markov models.
Eur. Biophys. J.
28
:
605
–609.
Miller, C., and E. Racker.
1976
. Ca2+-induced fusion of fragmented sarcoplasmic reticulum with artificial planar bilayers.
J. Membr. Biol.
30
:
283
–300.
Norris, J.R. 1997. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge, UK.
Qin, F., A. Auerbach, and F. Sacks.
1996
. Estimating single-channel kinetic parameters from idealized patch-clamp data containing missed events.
Biophys. J.
70
:
264
–280.
Qin, F., A. Auerbach, and F. Sacks.
2000
. A direct optimization approach to hidden Markov modeling for single channel kinetics.
Biophys. J.
79
:
1915
–1927.
Rengifo, J., R. Rosales, A. González, H. Cheng, M.D. Stern, and E. Ros.
2002
. Intracellular Ca2+ release as irreversible Markov process.
Biophys. J.
83
:
2511
–2521.
Robert, C.P., and C. Casella. 1999. Monte Carlo Statistical Methods. Springer Verlag, New York.
Rosales, R., W.J. Fitzgerald, and S.B. Hladky.
2002
. Kernel estimates for one and two dimensional ion channel dwell-time densities.
Biophys. J.
82
:
29
–35.
Rosales, R., J.A. Stark, W.J. Fitzgerald, and S.B. Hladky.
2001
. Bayesian restoration of ion channel records using hidden Markov models.
Biophys. J.
80
:
1088
–1103.
Rosales, R.A.
2004
. MCMC for hidden markov models incorporating aggregation of states and filtering.
Bull. Math. Biol.
In press
.
Sachs, F., F. Qin, and P. Palade.
1995
. Models of Ca2+ release channel adaptation.
Science.
267
:
2010
–2011.
Saftenku, E., A.J. Williams, and R. Sitsapesan.
2001
. Markovian models of low and high activity levels of cardiac ryanodine receptors.
Biophys. J.
80
:
2727
–2741.
Schiefer, A., G. Meissner, and G. Isenberg.
1995
. Ca2+ activation and Ca2+-inactivation of canine reconstituted cardiac sarcoplasmic reticulum Ca2+-release channels.
J. Physiol.
489
:
337
–348.
Sitsapesan, R., and A.J. Williams.
1994
. Gating of the native and purified cardiac SR Ca2+-release channel with monovalent cations as permeant species.
Biophys. J.
67
:
1484
–1494.
Sitsapesan, R., and A.J. Williams.
2000
. Do inactivation mechanisms rather than adaptation hold the key to understanding ryanodine receptor channel gating?
J. Gen. Physiol.
116
:
867
–871.
Sitsapesan, R.A.P., R. Montgomery, and A.J. Williams.
1995
. Ca2+ activation of sheep cardiac SR Ca2+ release channel on physiologically relevant time course.
Circ. Res.
77
:
765
–772.
Sobie, E.A., K.W. Dilly, J. dos Santos Cruz, W.J. Lederer, and M.S. Jafri.
2002
. Termination of cardiac Ca2+ sparks: an investigative mathematical model of calcium-induced calcium release.
Biophys. J.
83
:
59
–78.
Song, L., and K.L. Magleby.
1994
. Testing for microscopic reversibility in the gating of maxi K+ channels using two-dimensional dwell-time distributions.
Biophys. J.
67
:
91
–104.
Stern, M.D.
1992
. Theory of excitation-contraction coupling in cardiac muscle.
Biophys. J.
63
:
497
–517.
Stern, M.D., L. Song, H. Cheng, J.S. Sham, H.T. Yang, K.R. Boheler, and E. Rios.
1999
. Local control models of cardiac excitation-contraction coupling. A possible role for allosteric interactions between ryanode receptors.
J. Gen. Physiol.
113
:
469
–489.
Tate, C.A., R.J. Bick, A. Chu, W.B. Van Winkle, and M.L. Entman.
1985
. Nucleotide specificity of cardiac sarcoplasmic reticulum. GTP-induced calcium accumulation and GTPase activity.
J. Biol. Chem.
260
:
9618
–9623.
Valdivia, H., J.H. Kaplan, G.C.R. Ellis-Davies, and W.J. Lederer.
1995
. Rapid adaptation of cardiac ryanodine receptors: modulation by Mg2+ and phosphorylation.
Science.
267
:
1997
–2000.
Velez, P., S. Györke, A. Escobar, J. Vergara, and M. Fill.
1997
. Adaptation of single cardiac ryanodine receptor channels.
Biophys. J.
72
:
691
–697.
Venkataramanan, L., and F.J. Sigworth.
2002
. Applying hidden Markov models to the analysis of single ion channel activity.
Biophys. J.
82
:
1930
–1942.
Villaba-Galea, C.A., R.A. Rosales, and A.L. Escobar.
2002
. RYR spatial organization is involved in Ca-sparks termination.
Biophys. J.
82
:
354
.
Villaba-Galea, C.A., B.A. Suarez-Isla, M. Fill, and A.L. Escobar.
1998
. Kinetic model for ryanodine receptor adaptation.
Biophys. J.
74
:
A58
.
Wang, S.Q., L.S. Song, L. Xu, G. Meissner, E.G. Lakatta, E. Ros, M.D. Stern, and H. Cheng.
2002
. Thermodynamically irreversible gating of ryanodine receptors in situ revealed by stereotyped duration of release in Ca2+ sparks.
Biophys. J.
83
:
242
–251.
Zaradniková, A., M. Dura, and S. Györke.
1999
. Modal gating transitions in cardiac ryanodine receptors during increases of Ca2+ concentration produce by photolysis of caged Ca2+.
Pflugers Arch
.
438
:
283
–288.
Zaradniková, A., and I. Zahradník.
1995
. Description of modal gating of the cardiac calcium release channel in planar lipid membranes.
Biophys. J.
69
:
1780
–1788.
Zaradniková, A., and I. Zahradník.
1996
. A minimal gating model for the cardiac calcium release channel.
Biophys. J.
71
:
2996
–3012.
Zaradnková, A., I. Zahradník, I. Györke, and S. Györke.
1999
. Rapid activation of the cardiac ryanodine receptor by submilisecond calcium stimuli.
J. Gen. Physiol.
114
:
787
–798.

The online version of this article contains supplemental material.

Abbreviations used in this paper: BIC, Bayes information criterion; CICR, Ca2+-induced Ca2+ release; HMM, hidden Markov model; MCMC, Markov chain Monte Carlo.