Type-II ryanodine receptor channels (RYRs) play a fundamental role in intracellular Ca^{2+} 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 Ca^{2+} concentrations are made by considering 16 models, some of them previously reported in the literature. Results clearly show that single RYRs have Ca^{2+}-dependent gating modes. Moreover, our results demonstrate that single RYRs responding to a sudden change in Ca^{2+} 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 Ca^{2+} concentrations.

## Introduction

Type II calcium release RYR channels play a fundamental role in the intra cellular Ca^{2+} signaling dynamics of cardiac muscle cells that govern contractility of the heart. RYR channels are activated by a small Ca^{2+} influx through the surface membrane by a process known as Ca^{2+}-induced Ca^{2+} 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 Ca^{2+} concentration (Fabiato, 1985). Since then, the cytosolic RYR Ca^{2+} 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 Ca^{2+} 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 Ca^{2+}-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 Ca^{2+} release events known as Ca^{2+} 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 Ca^{2+}. The results obtained provide direct evidence for the existence of different Ca^{2+}-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 Ca^{2+} 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 Ca^{2+} 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 CsCH_{3}SO_{3}, 20 μM CaCl_{2}, 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 CsCH_{3}SO_{3}). 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 CsCH_{3}SO_{3} 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×10^{6} samples each) were made at three different steady-state free Ca^{2+} concentrations (1, 10, and 100 μM cis). Free Ca^{2+} concentration was verified using a Ca^{2+} electrode.

### Data Analysis

A representative single-channel recording at each Ca^{2+} 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* = {*O*_{1}, …, *O*_{v}} and *C* = {*C*_{1}, …, *C*_{w}}, with *v* + *w* = *n* and *O*_{i} and *C*_{i} 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 *p*_{ij}, grouped into a transition probability matrix *P* = [*p*_{ij}], for any *i*, *j* ϵ *O*, or *C*, the class conductance levels μ* _{C}* and μ

*, the class noise variances σ*

_{O}^{2}

*and σ*

_{C}^{2}

*, and finally the initial density for the hidden process Λ = [λ*

_{O}*],*

_{i}*i*ϵ{

*O*,

*C*}. For convenience hereafter we write

**θ**= (

*P*, μ

*, μ*

_{C}*, σ*

_{O}^{2}

*, σ*

_{C}^{2}

*, Λ).*

_{O}### 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

*y*

_{1},

*y*

_{2}, …

*y*

_{N}, and

**a particular realization**

*z**z*

_{1},

*z*

_{2}, …,

*z*

_{N}of the underlying process, i.e.,

**represents a sequence of states visited by the channel. Furthermore, let**

*z**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*(

*,*

**θ****) the likelihood function,**

*z*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

with

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

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. 2–4 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

*L*

**θ**

*,*

^{u}

*z*^{u}),

*u*= 1,2, …,

*m*, by running a Markov chain for a long time, such that in the limit as

*m*→ +∞, the pair (

**θ**

*,*

^{m}

*z*^{m}) 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},

*z*^{1}). 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 *M*_{3} 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 σ_{O2}^{2} and σ_{C2}^{2}. After a few steps the walk finally stays within a boundary that corresponds to a region of maximum posterior probability for (σ_{O2}^{2},σ_{C2}^{2}). This figure also presents the output generated for the hidden process shown as a randomly chosen set of idealized traces, *z*^{u}, 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.,

for some 1 < *b* < *m* and *G* = *m* − *b* + 1.

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

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* = [*k*_{ij}], 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 *k*_{ij} for *k*** _{ij}*, the elements of

*Q**. If the continuous process is sampled fast enough, i.e., if

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

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

and

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 Ca^{2+} conditions and each model (see Fig. 1). The initial hidden realization, *z*^{1}, was drawn from the initial density Λ^{1} and the initial probability matrix *P*^{1}. 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, (*z*^{u}), 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

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

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

This quantity is efficiently computed, once μ̅_{π}

**is partitioned into**

*y**s*independent segments of lengths

*N*

_{1},

*N*

_{2}, …,

*N*

_{s}and the likelihood for the

*i*th segment is denoted by

*L*

_{i}, then from Eq. 11 the

*BIC*takes the form

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

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*^{kδ} = exp(*Qk*δ) (see Eq. 19 of Colquhoun and Hawkes, 1995b). The open probability *P*_{o} as a function of *k*δ, is thus obtained via Eq. 14 by adding the entries in **η** that involve open states, i.e.

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 Ca^{2+} 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 Ca^{2+} (3.051 × 106 samples) with scheme *M*_{3} (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 Ca^{2+} 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 Ca^{2+} 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 *M*_{3}) and the worst one (scheme *M*_{15}). 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 *M*_{1} and *M*_{2}.

(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 *M*_{13}, share a common motif that includes states *O*_{2}, *O*_{3}, *C*_{4}, *C*_{5}, and *C*_{6}. Schemes, such as *M*_{5} and *M*_{6}, 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 *C*_{1}, *C*_{2}, and *C*_{3} in *M*_{4} and *M*_{5} 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 *M*_{3}.

### Q Matrix Estimates

Estimates for the *Q* matrices at each Ca^{2+} 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 *M*_{3} for each Ca^{2+} 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, σ_{O}_{2} and σ_{C}_{2}, respectively.

Simple inspection suggests that rate estimates for 1 and 100 μM Ca^{2+} have converged after ∼1,000 iterations. Convergence takes perhaps longer for the data at 10 μM Ca^{2+} (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 *C*_{6} at 1 μM Ca^{2+}. Rates leaving long-lived states with larger occupancy probabilities, such as *C*_{1} at 1 μM Ca^{2+}, 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 Ca^{2+}-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 Ca^{2+} induced changes. For this, the next section starts with a description of the RYR channel gating modes present at 1 μM Ca^{2+} and follows then by identifying the states involved in each mode.

Gating modes are clearly visible at 1 μM Ca^{2+}. Fig. 5

displays two representative single RYR channel records illustrating at least three gating modes. The first of them, labeled *G*_{l}, consists of long openings interrupted by brief closures. The second, named *G*_{z}, consists of brief openings followed by brief closures and finally the third mode, named *G*_{o}, is formed by brief openings followed by long shut periods. As shown in Fig. 5, open sojourns in *G*_{l} and *G*_{z} are grouped in clusters forming burst-like behavior. *G*_{l} has previously been identified in single RYR channel studies as a mode of high open probability named H gating mode, and *G*_{o} as a mode with low probability named L gating mode (Zaradniková and Zahradník, 1995, 1996; Armisén et al., 1996). *G*_{z} is not mentioned in these earlier references. The states participating in each of the *G*_{l}, *G*_{z}, and *G*_{o} modes and the closed periods between them may be identified directly by inspecting the last MCMC sampled realizations for the hidden process (*z*^{u}). 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 μ Ca^{2+} the channel spends most of its time at the *C*_{1} and *C*_{2} states, which account for the long closed sojourns between the *G*_{l} and *G*_{z} modes. From the *C*_{2} state, the channel may then enter the *O*_{3} state generating a brief opening. From the *O*_{3} state, it may enter the *C*_{4} state forming events corresponding to the *G*_{o} gating mode. Also from the *C*_{2} state, via *C*_{3}, the RYR channel may also (very rarely) enter the *O*_{1} sate generating a long lasting opening. These long openings may be followed by transitions toward the *O*_{2} state (brief openings) and subsequently to the *C*_{4} or *C*_{6} states (brief closings) forming bursts of channel activity characteristic of the *G*_{l} gating mode. From either the *O*_{1} or *O*_{2} states, *G*_{l} bursts are terminated either via *O*_{2}→*C*_{4}→*O*_{3}→*C*_{2} or *O*_{1}→*C*_{3}→*C*_{2}. The appearance of the *G*_{z} gating mode is associated with transitions of the form *O*_{3}→*C*_{4}↔*O*_{2}↔*C*_{6}, where the *C*_{6} state is responsible for the very fast closed characteristic of this mode. From the rates leaving *C*_{2} shown in Table III, there are 23 more transitions per unit of time toward the *C*_{3} state, the intermediate state for the initiation of a *G*_{l} gating mode bursts, than those for the *O*_{3} state that starts the *G*_{o} gating mode. Despite this, *G*_{l} gating mode bursts are quite rare when compared with the occurrence of the *G*_{o} gating mode. This is due to an extremely low rate for the transition from the *C*_{3} to *O*_{1} state. The appearance of the *G*_{z} or *G*_{o} modes is governed by the rates leaving the state *C*_{4} to the *O*_{2} and *O*_{3} states, respectively. Note that transitions to the *O*_{3} state are twice as frequent as those to the *O*_{2} state. This explains the apparent excess of the occurrence of *G*_{o} mode over *G*_{z} bursts. These results suggest the following gating mode frequency relations, *G*_{l }*G*_{z }*G*_{o}. At this Ca^{2+} 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 Ca^{2+} the channel tends to leave *C*_{1} and *C*_{2}, which are shorter than at 1 μM Ca^{2+}, and starts spending more time in states *C*_{3} and *C*_{4}, 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 Ca^{2+} concentration, the preferred open state is now *O*_{2}. This state is shorter than the most probable open state at 1 μM Ca^{2+} (*O*_{1}), which was responsible for long open sojourns during *G*_{l} bursts. This accounts for a considerable reduction in the total *G*_{l} open time. Bursts of type *G*_{z} are ∼9 times more frequent than sojourns in the *G*_{o} mode, whereas the proportion between the *G*_{l} and *G*_{o} gating modes is more equally distributed in contrast to the situation at 1 μM Ca^{2+}.

Finally, at 100 μM Ca^{2+} the channel now prefers to stay at *C*_{4}, a briefly lived closed state, and mostly frequents two intermediate brief open states, *O*_{2} and *O*_{3}. 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 Ca^{2+}-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 Ca^{2+} both open and closed densities show two clearly defined maxima. Fast events in the open class (Fig. 7 A) constituted by sojourns in *O*_{2} and *O*_{3} are about five times more frequent than slow sojourns in *O*_{1}. A high frequency of transitions to *O*_{3} accounts for most openings in the *G*_{o} mode. The closed class (Fig. 7 B) presents quite frequent fast and slow events corresponding to sojourns in *C*_{6} and *C*_{1}, *C*_{5}, respectively. The presence of multiple maxima in both classes at 10 μM Ca^{2+} is not as apparent as for 1 Ca^{2+}. This accounts for a partial disappearance of modal behavior. At this Ca^{2+} 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 Ca^{2+}. Furthermore, *O*_{3} is now the most frequent open state, being then closely followed by *O*_{2} and *O*_{1}. The most frequent states for the closed class (Fig. 7 D) are *C*_{6} and *C*_{4}. Provided the gating modes have still the same meaning as at 1 μM Ca^{2+}, both plots at 10 μM Ca^{2+} suggest the following mode frequency relations *G*_{l} < *G*_{o} < *G*_{z}. Finally at 100 μM Ca^{2+} (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 Ca^{2+}), the single-state frequency proportions are roughly the same as for 10 μM Ca^{2+}. These results suggest that a mayor change in the single state frequencies occurs only around 1–10 μM Ca^{2+}.

The univariate single-class densities suggest the existence of Ca^{2+}-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 Ca^{2+}. These modes appear to gradually merge into a single mode at progressively higher Ca^{2+} concentrations (10–100 μM Ca^{2+}). 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 Ca^{2+}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 *G*_{z} mode and is located at the bottom left corner of the plot. Another clearly identified maximum corresponds to long openings followed by brief closures (*G*_{l} 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 (*G*_{o} 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 Ca^{2+}. At 10 Ca^{2+} 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 Ca^{2+} (Fig. 8 F). However, correlations at 100 μM Ca^{2+} are much smaller (note contour level scales). This analysis presents direct evidence for RYR channel modal gating and reveals how it changes at different Ca^{2+} 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 Ca^{2+} 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 Ca^{2+}. 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 *C*_{6} (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 *G*_{z} gating mode. This is confirmed by comparison of the dependence-differences in Figs. 10 B and 8 B. The disappearance of *G*_{z} modifies the correlation pattern observed at 1 μM Ca^{2+}, 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 Ca^{2+} (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 *a*_{1}, *a*_{2}, …, *a*_{n−1}, *a*_{n} denote the states involved in the cycle of a particular gating mechanism. Gating is reversible if the transition rates in the cycle satisfy

For the cycle of scheme *M*_{3}, the clockwise product of the six rates involved at 1 μM Ca^{2+} is 6.9 × 10^{−6}. The product of the six rates in the opposite direction product is 6.6 × 10^{−6}. For 10 μM Ca^{2+} the products are 1.5 × 10^{−5} and 1.8 × 10^{−5} respectively. For 100 μM Ca^{2+} 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 Ca^{2+} Steps

The response of RYR channels to fast Ca^{2+} changes can be evaluated by computing the relaxation toward equilibrium of the probability of being in the open class, *P*_{o}. Here, we present the response of mechanism *M*_{3} to instantaneous changes in the Ca^{2+} concentration from 0 to 1, 10, or 100 μM Ca^{2+}. 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 Ca^{2+} steps and initial conditions considered, Fig. 12

shows that the *P*_{o} 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

**were then estimated and compared to their original values. To this end, two synthetic datasets were generated assuming two different Ca**

*z*^{2+}conditions. One synthetic dataset was generated using

*Q*matrix estimates for the

*M*

_{3}model at 1 μM Ca

^{2+}. The specific rate constants used for the

*M*

_{3}model are presented in Table III. The second synthetic dataset was generated using Q-matrix estimates for the

*M*

_{10}model at 10 μM Ca

^{2+}. The specific rate constants used for the

*M*

_{10}model are presented below (in ms

^{−1}).

Synthetic datasets with N = 2 × 10^{6} 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 σ^{2}* _{C}* = σ

^{2}

*= 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.*

_{O}For the synthetic dataset generated using the *M*_{3} model, the Gibbs sampler was run with four different gating schemes (*M*_{3}, *M*_{2}, *M*_{7}, and *M*_{14}). The actual *Q*-matrix samples for the *M*_{3} 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 *M*_{3}, *M*_{2}, *M*_{7}, and *M*_{14} 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 *M*_{3} > *M*_{2} > *M*_{7} > *M*_{14}. In the *M*_{3} model case, the estimates generated for μ_{C}, μ_{O}, and σ^{2}* _{C}* and σ

^{2}

*were 0.9979, 2.679, and 0.5961, respectively. The estimates for the rate constants (in ms*

_{O}^{−1}) are listed below.

For the synthetic dataset generated using the *M*_{10} model, the Gibbs sampler was also run for with four different gating schemes *M*_{1}, *M*_{3}, *M*_{10}, and *M*_{11}. The rate constant estimates (in ms^{−1}) are listed below.

The *BIC*/2 factors for the *M*_{1}, *M*_{3}, *M*_{10}, and *M*_{11} 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 *M*_{10} > *M*_{11} > *M*_{1} > *M*_{3}.

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 Ca^{2+} concentrations. Here, the reproducibility of the model selection is evaluated using real single RYR channel recordings at a single cytosolic Ca^{2+} 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 Ca^{2+} are presented in Table VI

. Estimates obtained from the analysis of another long (6.84 × 10^{6} samples) single RYR channel recording at 100 μM Ca^{2+} 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 (*M*_{2}, *M*_{3}, and *M*_{13}). 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 (*M*_{2}) 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 *M*_{4} > *M*_{12} > *M*_{9} > *M*_{10} > *M*_{16} >*M*_{14} > *M*_{15}. More importantly, the top five best-fit models were *M*_{2}, *M*_{3}, M_{5}, *M*_{6}, and *M*_{13} 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 Ca^{2+} 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. 6–9). 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 Ca^{2+} 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 Ca^{2+} 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 Ca^{2+} 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 Ca^{2+} gradient as the source of free energy that keeps the system out of time reversibility.

### Modal Gating and Adaptation

The steady-state Ca^{2+} 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 Ca^{2+} 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 Ca^{2+} elevation (Györke and Fill, 1993; Zaradnková et al., 1999). In some studies, channel activity peaks and then spontaneously decays following the fast Ca^{2+} elevation. The spontaneous decay was originally termed adaptation (Györke and Fill, 1993) and was particularly interesting because it occurred at Ca^{2+} concentrations that do not inactivate the RYR channel under steady-state conditions. The spontaneous decay (i.e., adaptation) has now been attributed to a Ca^{2+} 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 Ca^{2+} 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 Ca^{2+} elevation had a very fast (1 ms) overshoot at its leading edge. At first, this very fast Ca^{2+} overshoot was not considered to play a role in the much slower spontaneous decay (Velez et al., 1997). However, the possibility that the fast Ca^{2+} 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 Ca^{2+} overshoot. The analysis here supports the notion that this spontaneous decay is generated by a Ca^{2+}- and time-dependent modal gating shift (see below). The analysis here also suggests that a fast Ca^{2+} 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 Ca^{2+} 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 Ca^{2+} 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 Ca^{2+} activation rate depends on the initial conditions that exist before the fast Ca^{2+} change.

Several studies have detailed the Ca^{2+}-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 Ca^{2+}-dependent modal gating (see Fig. 8). The analyses here show that the Ca^{2+} dependence of modal gating is steepest at relatively low Ca^{2+} concentrations (≈1 μM). This corresponds to the Ca^{2+} 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 Ca^{2+}- and time-dependent modal gating shift. As initially proposed by Zaradniková and Zahradník (1995) the Ca^{2+}-dependent modal shift hypothesis rests essentially on two basic assumptions: (a) an increase of Ca^{2+} should favor transitions out of states with high *P*_{o}, but (b) at resting conditions Ca^{2+} 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 Ca^{2+}-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 Ca^{2+} 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 (*M*_{2}) 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 Ca^{2+} condition. The simultaneous estimation of a single common matrix should significantly improve its statistical properties; moreover, this should lead to the identification the Ca^{2+}-dependent transitions of a gating scheme. The latter would certainly clarify several aspects in the Ca^{2+}-dependent modal shift hypothesis. Finally, a single *Q* matrix with well-identified Ca^{2+}-dependent transitions would enclose all the necessary information for predicting the Ca^{2+}-dependent dynamic properties of single RYR channels. This includes the response to any imaginable Ca^{2+} 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 Ca^{2+}-regulatory environment that exists in living cells.

To make inferences by considering data at several Ca^{2+} concentrations simultaneously one has to be able introduce constraints for the rates, such that the involved states become well-defined occupancy sites for Ca^{2+}. For instance, one would generally like to set *k*_{ij} = α* _{ij}*[Ca

^{2+}], where α

*is the rate constant that results in the case [Ca*

_{ij}^{2+}] = 1 μM, if

*k*

_{ij}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,

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 Ca^{2+}-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

The online version of this article contains supplemental material.

*Abbreviations used in this paper: BIC*, Bayes information criterion; CICR, Ca^{2+}-induced Ca^{2+} release; HMM, hidden Markov model; MCMC, Markov chain Monte Carlo.