## Abstract

The inositol 1,4,5-trisphosphate (IP_{3}) receptor (IP_{3}R) channel is crucial for the generation and modulation of intracellular Ca^{2+} signals in animal cells. To gain insight into the complicated ligand regulation of this ubiquitous channel, we constructed a simple quantitative continuous-time Markov-chain model from the data. Our model accounts for most experimentally observed gating behaviors of single native IP_{3}R channels from insect Sf9 cells. Ligand (Ca^{2+} and IP_{3}) dependencies of channel activity established six main ligand-bound channel complexes, where a complex consists of one or more states with the same ligand stoichiometry and open or closed conformation. Channel gating in three distinct modes added one complex and indicated that three complexes gate in multiple modes. This also restricted the connectivity between channel complexes. Finally, latencies of channel responses to abrupt ligand concentration changes defined a model with specific network topology between 9 closed and 3 open states. The model with 28 parameters can closely reproduce the equilibrium gating statistics for all three gating modes over a broad range of ligand concentrations. It also captures the major features of channel response latency distributions. The model can generate falsifiable predictions of IP_{3}R channel gating behaviors and provide insights to both guide future experiment development and improve IP_{3}R channel gating analysis. Maximum likelihood estimates of the model parameters and of the parameters in the De Young–Keizer model yield strong statistical evidence in favor of our model. Our method is simple and easily applicable to the dynamics of other ion channels and molecules.

## INTRODUCTION

Ca^{2+} is one of the most important signaling ions in that it controls numerous cellular functions such as neurotransmitter release, muscle contraction, and nuclear transcription (Berridge et al., 1998). Ca^{2+} performs such a wide range of tasks through a hierarchy of signaling patterns ranging from single-channel events to waves sweeping over entire organs (Berridge et al., 1998). In many cell types, the inositol 1,4,5-trisphosphate (IP_{3}) receptor (IP_{3}R) Ca^{2+} release channel plays a central role in orchestrating the hierarchical Ca^{2+} patterns. Thus, it is crucial to understand the mechanisms regulating the IP_{3}R. Considerable efforts and expertise have gone into modeling the modulation of IP_{3}R channel gating by its main physiological ligands: IP_{3} and Ca^{2+} (De Young and Keizer, 1992; Sneyd and Dufour, 2002; Dawson et al., 2003; Mak et al., 2003; Shuai et al., 2007; Gin et al., 2009; Swaminathan et al., 2009). Most of these models can explain the Ca^{2+} dependence of the open probability (*P*_{o}) of IP_{3}R. Only the model in Mak et al. (2003) can also accurately reproduce the IP_{3} dependence of *P*_{o}, but it required 14 parameters to do so. One of us (Mak) has conducted several studies of the native IP_{3}R channel from insect Sf9 cells, which has a primary sequence most closely related to the type 1 IP_{3}R isoform in vertebrate cells (Foskett et al., 2007), to elucidate various behaviors of the IP_{3}R channel, including the ligand dependence of its *P*_{o} and mean open- and closed-channel dwell times (*τ*_{o} and *τ*_{c}, respectively), its modal gating behaviors in various ligand concentrations, and the transient response of the channel to rapid changes in ligand concentrations. To date, no model of the IP_{3}R has accounted for all these observations. Here, we propose such a model.

## MATERIALS AND METHODS

The many observations mentioned above have enabled us to develop a highly constrained mathematical model for the IP_{3}R channel gating kinetics, with the model structure and all of the parameters inferred directly from experimental data. The model consists of a Markov chain with 3 open states and 9 closed states, and requires 28 parameters. It reproduces the ligand-dependent equilibrium *P*_{o}, *τ*_{o}, and *τ*_{c} over wide ranges of cytoplasmic free Ca^{2+} () and IP_{3} () concentrations (Ionescu et al., 2006); the Ca^{2+} dependence of various modal gating features, including the relative prevalence (the probability of finding the channel in a particular gating mode) and lifetimes of the gating modes of the IP_{3}R channel, and channel *P*_{o} in each of the gating modes reported in Ionescu et al. (2007); as well as the distributions of the channel response latencies observed after various abrupt changes in the concentrations of one or both of its ligands (Mak et al., 2007).

We accomplished this working under the following major assumptions: (a) the model is a finite-state Markov chain, (b) the ligand-binding kinetics obey the law of mass action, (c) the chain is fully connected, (d) the channel must be able to unbind ligands, and (e) the model should be as simple as possible. We show in the supplemental text, section 1.1 (see also Fig. S1), that the gating and binding kinetics of IP_{3}R are consistent with the detailed balance hypothesis.

Structurally, the IP_{3}R channel is an oligomer. However, we do not assume that the IP_{3}R subunits in the channel are independent or that the ligand-binding sites in each subunit are independent. We make no a priori assumption about sequential or independent binding of IP_{3} and Ca^{2+} to the channel and therefore are consistent with experimental observations demonstrating that IP_{3} and Ca^{2+} bind to the channel cooperatively but with no specific sequential requirements (Mak et al., 2007), unlike many previous models that assumed sequential (Tang et al., 1996; Marchant and Taylor, 1997) or independent (Swillens et al., 1994; Tang et al., 1996; Hirose et al., 1998; Fraiman and Dawson, 2004) IP_{3} and Ca^{2+} binding.

Our one structural constraint is to search only over models with four IP_{3}-binding sites. Structure is at most an a priori weak constraint on dynamics. Beece et al. (1980) note that “the protein is not a rigid system in which a ligand moves in a fixed potential. Rather there is a strong mutual interaction between ligand and protein … A protein is not like a solid house into which the visitor (the ligand) enters by opening doors without changing the structure. Rather it is like a tent into which a cow strays.” Another difficulty with structurally motivated models is that they tend to have large numbers of states that cannot be deduced from patch-clamp recordings. Eigen’s “general allosteric model” has 55 states for a tetrameric channel in which each subunit can bind a single ligand and has two major conformations (Eigen, 1968). There is no a priori reason to believe that the IP_{3}R subunits have a particular number of conformations available to them. The IP_{3}R is known to have at least one IP_{3}-binding site and two Ca^{2+}-binding sites per monomer. The same reasoning that led to Eigen’s model would result in an IP_{3}R model with thousands of states. A large Markov chain will not have rates that can be learned from the data because of the large number of parameters that would need to be estimated from data. We believe that there are probably numerous small Markov chains that all yield similar likelihoods for the observed data. Based on the quality of the fit, one can use Markov chain selection criteria to at least select between specific categories of Markov chains. Here, we construct a model with cooperativity in the IP_{3}-binding dynamics. We refer to this model as the “cooperative model” (CM). We compare the likelihood of the data given the CM to the likelihood of the same data given the most frequently used De Young–Keizer model (DYKM) (De Young and Keizer, 1992). The DYKM yields a somewhat higher likelihood for the *P*_{o} data than the CM (although such was not the case for *Xenopus laevis* oocyte *P*_{o} data). We also performed maximum likelihood fits to dynamic data: latency measurements as well as stationary time series. Fits to the latency data only and the combined latency and current trace data indicated that the CM outperforms the DYKM. Thus, the DYKM does a very good job reproducing the *P*_{o} data from Sf9 cells. However, it does not compete well with the CM when the recent kinetic data are taken into account.

We use the theory of aggregated binary reversible Markov chains (Colquhoun and Hawkes, 1981; Fredkin et al. 1985. Proceedings of the Berkeley Conference in Honor of Jerzy Neyman and Jack Kiefer; Fredkin and Rice, 1986; Qin et al., 1997; Bruno et al., 2005) to derive a model via an iterative, data-driven approach. First, we consider the ligand dependencies of the equilibrium *P*_{o}, which provides evidence for the main ligand-binding stoichiometric complexes that must be included in a model Markov chain. A stoichiometric complex is specified by the number of Ca^{2+} ions (m) and IP_{3} molecules (n) bound to the channel, and by the gating state of the channel (open or closed). The C_{32} (O_{32}) complex refers to a closed (open) channel bound to three Ca^{2+} ions and two IP_{3} molecules. Such complexes can comprise more than one Markovian state (a state in a Markov chain from which the channel exits with time-independent rate constant(s) so that it has a mono-exponential lifetime distribution). The number of such states in a complex cannot be ascertained from ligand dependencies of the *P*_{o} data, but in the absence of data to the contrary, assumption (e) demands one state per complex. The parameters that can be inferred from *P*_{o} are the equilibrium occupancies of complexes and contain no information regarding network topology (supplemental text, section 1.2). Nor can the total number of complexes necessarily be inferred from *P*_{o}. At this stage, if the ligand dependencies of channel *P*_{o} do not require a particular complex, simplicity demands that the complex not be explicitly included in the chain. Then we consider the Ca^{2+} dependence of the modal gating behaviors: the relative prevalence of and the channel *P*_{o} in each of the modes. These data provide evidence for the inclusion of one more complex and suggest that some complexes can be in more than one gating state.

Furthermore, the Ca^{2+} dependence of the lifetimes of the gating modes provides some broad constraints on the connectivity between the various gating states. Once we have an estimate of the important states in the model and their occupancies, we deduce the full model with additional states and connections, and the rates related to the connections by fitting the response latencies of the channel to abrupt changes in ligand concentrations and current trace data. The full model that we constructed is shown in Fig. 1.

### Occupancy and flux

Instead of parameterizing our model with reaction rates, as done in most previous efforts to model the IP_{3}R channel gating, we use “state occupancy” and “flux” parameters as suggested in Yang et al. (2006). The two approaches are mathematically equivalent, but we find the occupancy–flux parameter approach simpler and more intuitive because it separates thermodynamic quantities (equilibrium occupancies) from kinetic quantities (reaction fluxes), or equivalently, it separates reaction rates from equilibrium constants. Occupancy and flux parameters are used also because occupancy determination is far more robust than state lifetime measurements. Failure to detect brief events (channel opening and closing or changing gating modes) can have a large impact on estimates of mean lifetimes of the states but does not affect state occupancy estimates significantly (supplemental text, section 1.3, and Fig. S2). Furthermore, the normalized occupancies are just rational functions of ligand. Thus, once the occupancies are learned by fitting the *P*_{o} function to the data (described below), they are fixed throughout the subsequent data-fitting process. This constraint is trivial to apply using occupancy parameters: we simply do not alter the occupancy parameters unless it can be done without significantly affecting channel *P*_{o}. Details of occupancy and flux parameters are given in the supplemental text, section 1.4.

### Low occupancy states

Another concept that we required is the notion of short-lived low occupancy states (Ullah et al., 2012). The *P*_{o} data suggest that there are complexes with no IP_{3} bound and complexes with four IP_{3} bound. Clearly there must be complexes with one, two, and three IP_{3} bound, but our analysis of channel *P*_{o} provides little evidence for them. We conclude that these states exist, but their occupancies are not high enough for them to be detected in the *P*_{o} data. Such low occupancy transition states are not included explicitly in the model and their occupancies are set to zero to satisfy the maximal simplicity assumption. Similarly, some complexes with no Ca^{2+} bound are linked directly to complexes with two Ca^{2+} bound in the model. However, the effects of these transition states are manifested in the particular functional form for the fluxes passing through them imposed by mass action and detailed balance. Proper treatment of low occupancy states is crucial, both for getting the rates right in the simplified model and for obtaining a model structure with learnable parameters. Applying such simplifications by excluding low occupancy states can result in a substantial reduction in the number of parameters required. More details regarding low occupancy states are given in the supplemental text, section 1.5. In the supplemental text, section 1.6, we derive the functional form of the transition rates between various model states (Fig. 1) and simplify the branches that involve low occupancy states by eliminating these states.

### Complex occupancies

We identify the relevant (high occupancy) complexes from the ligand dependencies of channel *P*_{o}. Relative to the reference unliganded closed complex C_{00}, the occupancies of the C_{mn} and O_{mn} complexes are proportional to , with occupancy parameters_{00} to each state (_{mn} and O_{mn} are respectively, where *Z*, the total occupancy over all states, is given bywith *m _{max}* and

*n*being the maximum number of Ca

_{max}^{2+}and IP

_{3}that the channel can bind. The equilibrium

*P*

_{o}function is:

Our initial *P*_{o}-based search resulted in a best-fit (according to the AIC score; Mcquarrie, 1998) model with five closed complexes (C_{00}, C_{04}, C_{24}, C_{32}, and C_{34}) and one open complex (O_{24}) (Fig. 1). We first used the Mathematica routine “NonlinearModelFit,” which does least squares and assumes a constant variance, meaning that it assumed the *P*_{o} function was Gaussian distributed with constant variance. We also assumed a variance proportional to *P*_{o} (1-*P*_{o}) and found little difference. Both distribution functions for the *P*_{o} data returned similar models. The fit to the *P*_{o} data are shown in Fig. 2. Although not essential to the *P*_{o} fit, the O_{14} complex is necessary to account for the Ca^{2+} dependence of channel *P*_{o} in the I gating mode (mode with intermediate *P*_{o}; see Modal occupancies). The C_{20} and C_{30} complexes are required to account for the latency distribution observed (see Model development using latencies), but the occupancies of these two complexes are so low that they do not contribute significantly to *P*_{o} or affect the quality of the *P*_{o} fits in Fig. 2. Theoretically, the *P*_{o} can be expressed as a function of the eight occupancy parameters of the nine complexes (occupancy of the reference C_{00} complex is 1). However, the terms 1, are all small relative to the other terms in *Z* for the values of used in the *P*_{o} measurements, and O_{14} complex is only required for the channel *P*_{o} in I mode. This implies that the *P*_{o} is effectively determined by the five occupancy parameters for the complexes with IP_{3} bound*P*_{o}.

Subsequent searches yielded other models with better AIC scores (for the fit to Sf9 *P*_{o} data) than the one we discuss here. However, we were unable to account for the modal gating statistics using these other models. One such model that gave a better AIC score on the Sf9 *P*_{o} data was DYKM. However, the CM performed better than the DYKM when fitted to the *P*_{o} data from *Xenopus* oocytes (see Comparison between the CM and the DYKM). The value of _{3} affinities along the path connecting C_{00} to C_{04} of ∼8 nM. The numerical values of the occupancy parameters are given in Table S1. In the remainder of this article, we construct a model based on these complexes.

### Modal occupancies

The Sf9 IP_{3}R channel exhibits three distinct gating patterns or modes: a quiescent L mode with low *P*_{o} in which brief openings are separated by long quiescent periods with long *τ*_{o}, an I mode with intermediate *P*_{o} in which the channel opens and closes rapidly with short *τ*_{c} and *τ*_{o}, and an H mode with high *P*_{o} in which the channel gates in bursts (Ionescu et al., 2007). All three modes are observed under all conditions for which the channel gates, and the channel spontaneously changes among all three modes, even under constant ligand conditions (Fig. 3 A). The gating characteristics of the channel (*P*_{o}, *τ*_{c}, and *τ*_{o}) remain remarkably consistent in each gating mode in all ligand conditions so the ligand dependencies of overall channel *P*_{o}, *τ*_{c}, and *τ*_{o} come from the ligand dependencies of the relative prevalence (or normalized occupancies), *π*^{M}, of the gating modes (M can be L for low mode, I for intermediate mode, and H for high mode) (Ionescu, et al., 2007). The key modal statistical data at our disposal are shown in Fig. 4. Each of the three plots shows nine data points corresponding to = 0.1, 1, and 89 µM ( = 10 µM for all data). The squares, bullets, and triangles in Fig. 4 represent the L, M, and H modes, respectively. To account for these observed modal gating behaviors, we have to extend our minimal model by splitting some complexes into two states when required, without adding new complexes other than O_{14}, which is required for the channel *P*_{o} in I mode.

Fig. 4 A shows the conditional open probabilities, ^{2+} dependence: rising from ∼0.2 at = 100 nM to 0.3 at = 1 µM, and then dropping back to 0.2 at = 89 µM.

Because of the ligand independence of _{24} complex into _{24} complex into

If the I mode contains only states with two Ca^{2+} and four IP_{3} bound, _{14} complex besides the C_{24} and O_{24} complexes are in the I mode. However, this will only give *π*^{I}) to increase with increasing (Fig. S4), which is not observed (Fig. 4 B). Thus, the peak in

Here, we have split the C_{04} complex into _{14}. Adding the O_{14} complex worsens the AIC score for the model fit to the *P*_{o} data (Fig. 2) by 1 but improves the quality of the fit by improving the log likelihood (see Eq. 5 below). Adding any other open states such as O_{04} or O_{34} worsens both the AIC and log-likelihood scores by a factor much higher than 1.

According to the AIC and BIC criteria (Mcquarrie, 1998), there is insufficient justification from experimental evidence to assume that the C_{34}, C_{32}, and C_{00} complexes consist of more than one state. For example, splitting the C_{34} complex into two states will add one more parameter to the model and will increase the AIC score by a factor of 2 (see the factor 2*k* in Eq. 5 below), without changing the *P*_{o} given by the model (Eq. 2). So we assign them as

Fig. 4 B shows the Ca^{2+} dependence of *π*^{M} = *Z*^{M}/*Z*. A fit of the *P*_{o} data in various and (curves in Fig. 2) followed by fits to the *π*^{M} data in different (curves in Fig. 4, A and B, respectively) using the Mathematica routine NonlinearModelFit (see Complex occupancies) determines the values of the occupancy parameters for the various states in the three gating modes.

By this point, the number of states in the model and the occupancy parameters have been determined. We will determine the connections between these states guided by the observed modal lifetimes and latency distributions. We will evaluate all possible connections and choose only those that lead to the observed behavior of the channel.

### Modal lifetimes

Some, but not all, connectivities between the various states can be inferred from the observed modal gating behaviors. Channel opening and closing events were observed as the channel remained in the same mode; therefore, there are direct connections between

The mean H mode lifetime, *τ*^{H}, has a maximum as a function of (Fig. 4 C, triangles). This suggested that the H mode is connected to *τ*^{H} at high . We also connect the H mode to *τ*^{H} decreases with decreasing . It should be noted that although moving from the H mode to the ^{2+} dissociation from the channel, the effective rate of the reaction is dependent (see Low occupancy states). The *τ*^{L} and *τ*^{I} in low . The long *τ*^{L} (Fig. 4 C, ∼ 3 s) indicates a low rate from

Finally, we connect the I and H mode open states via Ca^{2+}-independent rates to get the appropriate H mode lifetime under optimal ligand conditions. Further determination of the connectivity of the model requires consideration of channel response latencies (see Model development using latencies).

The remaining connections and the final values of the flux parameters will be determined below using the fit to the latency and time-series data.

### Model development using latencies

and on the cytosolic side of the IP_{3}R channel were abruptly changed in nuclear patch-clamp experiments in the cytoplasmic-side-out configuration using rapid perfusion solution exchange techniques to observe the latencies for the response of the channel to ligand concentration changes (Mak et al., 2007). was changed among low (<10 nM, which we will refer to as 0 for notational simplicity), optimal (2 µM) and inhibitory (300 µM) levels, and between 0 and high (10 µM) levels. (, ) were switched from nonactivating (2 µM, 0), (0, 10 µM), and (0, 0) to optimal (2 µM, 10 µM) to measure IP_{3} activation, Ca^{2+} activation, and Ca^{2+} and IP_{3} activation latencies, respectively. To measure corresponding deactivation latencies, (, ) were changed from optimal back to various nonactivating combinations. Vertebrate IP_{3}R channels exhibit spontaneous activity in (0, 0) (Mak et al., 2003). However, the insect Sf9 IP_{3}R does not exhibit such spontaneous activity (Mak et al., 2007), so the Ca^{2+} and Ca^{2+} and IP_{3} activation and deactivation latencies can be unambiguously determined. (, ) were changed from optimal to inhibitory (300 µM, 10 µM) and from (300 µM, 10 µM) to optimal to measure latencies of inhibition and recovery from inhibition, respectively. Latencies after activating or recovery solution switches are the durations between solution switches and the first observed openings of the channel in response to the changes, whereas latencies after deactivating or inhibitory solution switches are the durations from the switches to the last observed closings of the channel as its activity terminated. Histograms of these latencies were built up from repeated solution switches in multiple experiments as shown in Fig. 5.

Additional experiments were performed in which (, ) was switched from (0, 10 µM) to (300 µM, 10 µM) and back. In most of these switches, a burst of channel activity was observed before the channel was inhibited by high . But in 9 out of 103 switches, no channel open activity was observed. When was decreased from 300 µM to 0, the channel was transiently activated in only 6 out of 94 experiments (Mak et al., 2007). These observations indicate that _{24} complex.

Both the distributions of the inhibitory and inhibition recovery latencies after (, ) switching between (2 µM, 10 µM) and (300 µM, 10 µM) can be fitted with single-exponential curves with no deficit of short latencies (Mak et al., 2007). This suggests that the

The Ca^{2+} and IP_{3} activation latencies are short (Fig. 5 E), indicating that *τ*^{L}. Moreover, deficit in short latencies (<20 ms) for Ca^{2+} and IP_{3} activation and Ca^{2+} activation (Mak et al., 2007) suggests that

Furthermore, to satisfy assumption (d) and allow IP_{3} unbinding from the ^{2+} to unbind from _{3}, it must be connected to ^{2+} to unbind from

The three distinguishable peaks in the IP_{3} activation latency distribution (Fig. 5 A) require at least three pathways connecting the closed states with no IP_{3} bound to the _{3} activation.

Finally, to provide two pathways for the open states to reach ^{2+} deactivation (Fig. 5 D),

Next, we will determine the values of flux parameters involved in the transition rates between various states by fitting the model to the latency and time-series data. First, we set the initial values of flux parameters in all reactions, as described in the supplemental text, section 1.7. We will use these values as initial conditions for the global fit to the latency and time-series data described below. Although we estimate the initial values of flux parameters from observations, the initial guess proved unnecessary for the fit. The fit converged to the same final parameter values for random initial conditions.

### Maximum likelihood and AIC calculations

To determine various flux parameters, we performed maximum likelihood fits on current traces of channel gating under constant conditions of = 10 µM and = 0.1, 1, and 89 µM, and on the latency data. To accomplish this task, we used our own implementation and calculated the relevant “log-likelihood” functions for the CM in the open-source programming language “Octave.” We used the function “nelder_mead_min” in the optimization tool kit to minimize the likelihood scores (−log(likelihood(data))) of the latency and time-series data by varying the 18 free flux parameters (while holding the occupancy parameters fixed).

The (negative) of the logarithm of the likelihood function for the latency data are minus the sum of the logarithms of the likelihoods of all the latencies:

For the current traces, we calculated the log of the likelihood function:*π _{C}* is the initial probability of closed states being occupied at equilibrium,

*to*and

_{i}*tc*are the ith opening and closing in the time series, respectively, and

_{i}*u*

_{c}is an

*N*

_{c}(the number of close states) component vector of all 1’s.

*w*and

_{C}*w*are the equilibrium occupancies of the closed and open states, respectively. For the CM,

_{O}*N*= 9 and the number of open states,

_{C}*N*= 3. Thus

_{O}*Q*

_{CC},

*Q*

_{OO},

*Q*

_{CO}, and

*Q*

_{OC}are 9 × 9, 3 × 3, 9 × 3, and 3 × 9 matrices, respectively.

We did two different fits: one in which we fit to all the latency data and none of the time-series data, and a second in which we fit both the latency data and the time-series data. The total log likelihood of all data used in the fit was calculated as*N _{exp}* is the number of experiments, and data

_{i}is the dataset (latencies or time series) from experiment

*i*. There are eight latency experiments where each experiment consists of hundreds of trials and each trial represents a latency measurement. The time-series data consists of three experiments (each experiment results in one time series), one each for = 0.1, 1, and 89 µM. Finally, the AIC score is given by

*k*is the number of parameters in the model. The CM has 10 occupancy and 18 flux parameters (

*k*= 28).

To minimize the number of parameters, we assume that the effective rates for the reactions

### Online supplemental material

The supplemental text provides a more in-depth discussion and derivation of various statements made in the main text. In sections 1.1, 1.2, and 1.3, respectively, we show that the IP_{3}R obeys detailed balance (Fig. S1), the equilibrium occupancy data contain no information about the network connectivity, and missing events cause more error in *τ*_{o} and *τ*_{c} (Fig. S2). Sections 1.4, 1.5, 1.6, and 1.7 discuss the concepts of occupancy and probability flux, low occupancy states, simplification of the model, and parameter estimation, respectively, in detail. The derivation of mathematical expressions for *τ*_{o} and *τ*_{c}, and dwell-time distributions are given in sections 1.8 and 1.9. The experimental methods are described in section 1.10. Fig. S3 plots the ligand dependence of some rates in the model, and Fig. S4 proves that state _{3} activation latency has a minimum at intermediate value. The parameters and transition rates for the CM are given in Tables S1–S3, whereas those for the DYKM are given in Tables S4 and S5. The supplemental material is available at http://www.jgp.org/cgi/content/full/jgp.201110753/DC1.

## RESULTS

### Fitting modal lifetimes with the model

To reproduce the observed modal lifetimes (Fig. 4 C, *τ*^{M}) with our model, *τ*^{M} are expressed in terms of the occupancy and flux parameters in the model. In general, the lifetime *τ*^{X} of any aggregate X (complex, mode, or other combinations of Markovian states) in an aggregated Markov chain is given by*Z*^{X} is the unnormalized occupancy of aggregate X, and *J*_{X} is the unnormalized flux out of that aggregate. *J*_{X} is the sum of all the fluxes from all reactions from all Markovian states contained in the aggregate X to Markovian states not contained in X, so*r*_{SU} is the rate from a state S to a state U (supplemental text, section 1.6, and Table S2). The modal lifetimes given by the model are shown by the lines in Fig. 4 C.

### Fitting latency distributions with the model

To use the CM to account for the distributions of latencies observed in the rapid perfusion solution exchange experiments, we express the latencies in terms of occupancy and flux parameters. When the perfusion solution to which the channel is exposed is changed at time t = 0, the system is initially in an aggregate X containing *N*_{X} states. As time goes by, it will transition into another aggregate Y with *N*_{Y} states, which is disjoint from X. We assume that *N*_{X} + *N*_{Y} = *N*, the total number of states in the system. If *Q* is the generator matrix for the full system, we partition *Q* as*ij* element in *Q* is *J _{ij}*/

*Z*. At time t = 0, the occupancies of initial states (all states in X) are specified in the vector of initial normalized occupancies

_{i}*π*

_{(0)}. Let

*u*

_{X}and

*u*

_{Y}be the

*N*

_{X}and

*N*

_{Y}component vectors, respectively, in which all elements are 1, so

*π*

_{(0)}

*u*

_{X}= 1. The probability density function for the latency can be expressed as (Colquhoun and Hawkes, 1981; Fredkin et al. 1985. Proceedings of the Berkeley Conference in Honor of Jerzy Neyman and Jack Kiefer; Qin et al., 1997; Bruno et al., 2005):

*t*=

*t*

_{min}and

*t*

_{max}) in the model latency histograms

We apply this to our measurement of IP_{3} channel response latencies, taking care to choose states X and Y, as well as the initial occupancy vector *π*_{(0)} correctly, i.e., knowing the set of states in which the channel gates before and after the jump in the solution to which the channel is exposed. For activation and recovery from inhibition latencies, the choices of X and Y are more straightforward. The experiment starts with the channel being closed, and the latency is the time until the first channel opening as the channel enters an open state after the change of and/or . Thus, for activation and recovery from inhibition experiments, the aggregate X contains the nine closed states in the CM, tabulated in related matrices and vectors in this order: *π*_{(0)} is a vector having the initial normalized occupancies of the closed states. For the deactivation latencies, X is composed of

For IP_{3} activation in constant 2 µM Ca^{2+}, the logarithmically binned histogram for experimental latencies (Fig. 5 A, bars) has three peaks, the main one at ∼80 ms and two minor ones at ∼250 ms and 2 s (Mak et al., 2007), thus giving a mean IP_{3} activation latency of ∼500 ms, even though the majority of the activation latencies are shorter than 100 ms. This indicates that the Markovian model for ligand regulation of IP_{3}R channel has at least three distinct branches involved in IP_{3} activation (Mak et al., 2007). This is a natural feature of the CM. With = 2 µM,wherebecause only the _{3}. The _{3}, transiting through ^{2+}.

To fit the IP_{3} deactivation latency histogram (Fig. 5 B), we assume that the initial occupancies of the states are their equilibrium occupancies in = 10 µM and = 2 µM. The case of Ca^{2+} activation in constant (10 µM) is trickier and requires detailed consideration of the manner in which the experiments were performed. In an activation measurement, was switched from <10 nM to 2 µM. was not returned to <10 nM until ∼2 s after channel activity was detected to ensure that the channel was really activated. Similarly, in a deactivation measurement, ∼2 s was allowed after the last channel opening before was switched to 2 µM again, to confirm that the channel was really deactivated. We assume that for all Ca^{2+} activation measurements, the channel only occupies the ^{2+} activation experiments when the channel failed to go active. These failures to activate might have been caused by the IP_{3}R residing in state *Q*_{XX} and *Q*_{XY} were set in accord with = 2 µM and = 10 µM (Fig. 5 C).

The Ca^{2+} deactivation latency histogram has two distinct peaks (Fig. 5 D). The CM fails to capture this feature even though it has two distinct pathways of deactivation: one from

For Ca^{2+} and IP_{3} activation, the channel initially can only occupy the ^{2+} and IP_{3} activation and deactivation latency distributions that capture the general shape of the experimental distributions (Fig. 5, E and F).

For Ca^{2+} inhibition, *π*_{(0)} is set to the normalized occupancies of states ^{2+} inhibition latencies is unimpressive (Fig. 5 G). The observed inhibition latency distribution has at least two peaks, and the peaks are substantially sharper than the model distribution. The CM has only a single inhibition pathway that is characterized by a single flux parameter. Increasing the flux parameter to obtain a better fit for the inhibition latencies will shift the peak of the inhibition recovery latencies (Fig. 5 H) from the CM by the same amount in the same direction.

At this point, the inhibition pathways are opaque to us. Examination of experimental records reveals that after the majority (>80%) of switches from 2 to 300 µM, the channel quickly (∼30 ms) transits into a mode with low but nonzero *P*_{o}, followed by a transition into a different totally quiescent mode with *P*_{o} = 0. The inhibition latencies were evaluated as the durations the channel takes to enter the totally quiescent mode (Mak et al., 2007). Our present model does not distinguish the mode with low but nonzero *P*_{o} from the totally quiescent mode and therefore is not equipped to properly describe this gating behavior of the channel. Incorporating this feature in the model will require future experiments to gather enough data for modeling the kinetics of these two modes. We will then be able to devise a model distinguishing the two behaviors (with at least four gating modes) to more properly account for Ca^{2+} inhibition kinetics of the channel.

### Channel gating kinetics from the model

The mean open- and closed-channel durations (*τ*_{o} and *τ*_{c}) were calculated for various in = 10 µM (supplemental text, section 1.8). The values from the model agree reasonably well with experimental measurements over a broad range of (300 nM to 100 µM) (Fig. 6, A and B). The parameters were also used to generate simulated single-channel gating record exhibiting modal gating behavior that resembles experimental single IP_{3}R channel current record in 1 µM Ca^{2+} and 10 µM IP_{3} (Fig. 3 B). The open and closed dwell-time distributions from the model are computed as in the supplemental text, section 1.9, and are shown by red lines in Fig. 7. The distributions obtained from the observed time-series data are shown by the green bars in Fig. 7.

The supplemental text, section 1.9, also describes the derivation of the open and closed dwell-time distributions in H and I modes from the model. The dwell-time distributions from the model for = 1 µM in the I and H modes are given by the red lines in Fig. 8. The experimental data are presented by the green bars for comparison. The peaks near 1 ms in the observed closed-time distributions indicate that there are short-lived closed states available to both I and H modes (Fig. 8, B and D, respectively) that the model misses. There also seems to be a short-lived open state available to the I model (Fig. 8 A). Removing this discrepancy will require adding more states to the model. At this stage, we are not interested in complicating the model further. Nevertheless, our analysis about open and closed dwell-time distributions provides useful insight into the modal behavior of the channel. The discrepancy between the model and the observed dwell-time distributions in the I mode can also be removed by increasing the flux parameter involved in

### Comparison between the CM and the DYKM

The one model of note that gave a better AIC score for *P*_{o} data only is the DYKM. Here, we compare the maximum likelihood fit of the CM (Fig. 1) to the DYKM using time-series and latency data. We find that the CM has a higher likelihood of producing the observations than does the DYKM. The DYKM has only a single open state and cannot produce modal gating, so we do not try to extract modal statistics from it.

The DYKM is a 120-state model of the IP_{3}R. The DYKM treats the channel as consisting of three rather than four identical subunits. It assumes that the subunits are independent and that each subunit has three binding sites, one for IP_{3}, a Ca^{2+} activation, and Ca^{2+} inhibition. The subunit states can be labeled by three binary digits (000, 001, … 111), with the first, second, and third bits representing the occupancy of the IP_{3}, Ca^{2+} activation, and Ca^{2+} inhibition binding sites, respectively. The channel is open if all three subunits are in the (110) state and closed otherwise. The channel open probability can be written *K _{mno}* are the occupancy parameters for the subunit states. The full-channel model has 120 states of which 119 are closed. The DYKM requires 7 occupancy and 12 flux parameters, or, equivalently, 12 independent reaction rates.

The DYKM had an AIC score ∼10 better than the CM for the Sf9 *P*_{o} data. However, in the case of *Xenopus* *P*_{o} data, the DYKM had an AIC score of −92 as compared with −138 (better) for the CM. Furthermore, we show that *P*_{o} alone is inadequate to discriminate between models. To shed further light on the differences between the CM and the DYKM, we performed maximum likelihood fits on the latency data and current traces of channel gating using DYKM as described in Maximum likelihood and AIC calculations in Materials and methods. We varied the 12 free rate constants of the DYKM while holding the occupancy parameters fixed. The likelihood ratio for CM and DYKM, which is interpreted as the relative probabilities of the two models, is the exponential of the difference between their AIC scores:

For the DYKM, we took advantage of the fact that there is only a single open state, thus Eq. 3 becomes:*N _{C}* = 119 and

*N*= 1. Thus, w

_{O}_{O}is a scalar from which it follows that

*Q*is a 119 × 119 matrix,

_{CC}*Q*is scalar (equal to minus the inverse lifetime of the single open state), and

_{OO}*Q*and

_{CO}*Q*are 119 × 1 and 1 × 119 matrices, respectively.

_{OC}In both cases, fit to the latency data alone and fit to both the latency and time-series data (see Maximum likelihood and AIC calculations in Materials and methods), we calculated the likelihood score of the time-series data. In the first case, in which the time-series data were not used in the fit, the time-series data can be considered as “out-of-sample” data, with the differences in scores indicating which model does the best job of “predicting” the unseen data.

The results are summarized in Table 1, which demonstrate that there is strong statistical evidence for the CM. Adding more time-series data (up to 28 experiments for = 0.1, 1, and 89 µM each) further increases the probability for the CM as compared with DYKM.

The occupancy parameters for the DYKM are given in Table S4. The rate constants are given in Table S5. The fits to the latency distributions are shown in Fig. 9.

## DISCUSSION

### Deficiencies of the CM

In a famously titled section heading, George Box wrote “All models are wrong …” (Box, 1979). Ours is no exception. One of the major deficiencies of the CM is its lack of distinction between the gating mode in which the IP_{3}R channel gates with low (∼0.05) but nonzero *P*_{o} and the mode in which the channel is quiescent with *P*_{o} = 0 (Q mode). In the study of IP_{3}R modal gating behaviors under steady-state ligand conditions (Ionescu et al., 2007), it is difficult to tell if a channel remains closed for a long duration whether it stays in the L mode during that whole duration or has spontaneously entered into and exited from the Q mode during that time. However, reexamination of Ca^{2+} inhibition latency measurements reveals that the distinction between the L and Q mode is significant because the majority of such inhibition events proceed from channel in H mode, favored by the initial ligand conditions of (, ) = (2 µM, 10 µM), through the L mode to the Q mode, instead of going directly from the H mode into the Q mode. Thus, the latencies consist of two contributions: the time for the channel to transit from the H mode to the L mode and the time for the channel to enter the Q mode from the L mode.

Furthermore, the Ca^{2+} inhibition latency distribution clearly shows two peaks (at ∼70 and 900 ms). A natural explanation for this derived from the CM is that the short latencies are caused by channel transiting directly from *τ*^{I}, before they can be inhibited. However, examination of data records indicates that channel initially in the I mode when is changed enters the L mode very rapidly (∼3 ms). Thus, the long Ca^{2+} inhibition latencies are caused by the slow exit from the L mode into the Q mode. To properly model Ca^{2+} inhibition of IP_{3}R requires distinguishing the L mode from Q mode by including open O^{L} states for the channel in the L mode, which will increase the complexity of the model substantially. We believe that a more accurate description of the inhibition latencies is crucial for the understanding of how the channel functions and will be the subject of our future research.

The CM has a slightly higher number of short IP_{3} activation latencies (<10 ms) compared with the latency data. This may be a consequence of the exclusion of low occupancy states during IP_{3} binding (like the C_{i1}, C_{i2}, and C_{i3} states excluded from the _{3} at a rate of ∼0.4 ms^{−1}, substantially slower than the diffusion-limited rate of ∼20 ms^{−1} in = 10 µM (assuming a reaction radius of 1 nm and diffusion coefficient of 0.5 µm^{2} ms^{−1} for IP_{3}; Redner, 2001). This indicates that IP_{3} probably needs to enter its binding site in IP_{3}R with a particular orientation, as suggested by the structure of the site in the receptor with residues from different parts of the site interacting with different phosphate groups in IP_{3} through networks of hydrogen bonds (Bosanac et al., 2002). Alternately, this may be caused by slow conformation changes the channel has to undergo after each IP_{3} binding that prolongs the activation latencies.

With activation and deactivation latencies considered only for jumping between <10 nM and 2 µM and between 0 and 10 µM, there are insufficient data about the dynamic behaviors of the IP_{3}R channel with one to three IP_{3} or only one Ca^{2+} bound for the CM to include those states explicitly under the simplicity principle. The failure of the model to generate two resolved peaks for the Ca^{2+} deactivation latency distribution (Fig. 5 D) suggests that those processes probably involve multiple fast transitions through intermediate states, particularly for Ca^{2+} activation. Future latency measurements involving jumping between 0 and subsaturating levels (100 nM, say) and/or between <10 nM and suboptimal levels (300 nM, say) could provide more insights into these processes to improve the model. Finally, we cannot rule out the possibility of other simple models that might better describe our data.

### Strengths of the CM

The rest of Box’s section heading was “… some models are useful.” The CM is the result of the synthesis of a large body of data gathered in several studies of many different aspects of ligand regulations of IP_{3}R channel gating in a novel, data-driven minimalist approach. The model uses five to six complex occupancy parameters to capture accurately not only the Ca^{2+} but also the IP_{3} dependencies of the channel equilibrium *P*_{o}. The model does a fair job of accounting for the *P*_{o}, *τ*_{o}, *τ*_{c}, modal gating statistics, and the majority of latency distributions examined. Moreover, our model makes concrete predictions regarding the single-channel response to rapid changes in the cytosolic ligand concentrations of IP_{3} and Ca^{2+}, so its validity can be verified by future experiments. With jumps of (, ) from (2 µM, 0) to (2 µM, 10 µM), the channel mostly enters the H mode from ^{2+}-dependent transition to either ^{2+}-dependent _{3} activation latency should increase with increase in . Thus, our model predicts that the mean IP_{3} activation latency has a minimum of ≈0.26 s at moderate of ∼0.6 µM (Fig. S5).

The CM predicts that a higher number of longer Ca^{2+} activation latencies will be detected if the channel is exposed to < 10 nM longer than the 2 s used (Mak et al., 2007) before is changed to 2 µM. This is because the relatively Ca^{2+}-independent *τ*^{I} (∼0.4 s) observed suggests that the channel may not be fully equilibrated between the ^{2+} activation latencies.

As mentioned in Maximum likelihood and AIC calculations in Materials and methods, we minimize the number of parameters by assuming that the effective rates for the reactions from the states with no Ca^{2+} bound to the states with two Ca^{2+} bound are the same at = 0, 10 µM. Relaxing this constraint did not improve the fit to the data. This result indicates that the transition rate from the states with no Ca^{2+} bound to the states with two Ca^{2+} bound is independent of the number of IP_{3} bound. Whether this assertion is correct or not can be investigated further by performing more ligand jump experiments, where is jumped from 0 µM to suboptimal values, e.g., jumping (, ) from (0 µM , 0 µM) to (0.5 µM , 0 µM) and (0 µM , 10 µM) to (0.5 µM , 10 µM).

The CM, especially its deficiencies, provides novel insights that are not immediately obvious from available data. The failure of the CM to adequately account for the Ca^{2+} inhibition latency distribution reveals the importance of drawing the distinction between the L mode with low but nonzero *P*_{o} and the quiescent Q mode with *P*_{o} = 0 for improving our understanding of Ca^{2+} regulation of the channel, which is not obvious from the modal analysis in Ionescu et al. (2007).

Although the development of the CM is directed by experimental observations, it also provides insights to guide the designs of future experiments and data analysis protocols. Our model suggests that experiments in which is switched between 0 and subsaturating levels and between physiological resting to suboptimal levels can provide valuable information to improve the characterization of the kinetics of the channel with one to three IP_{3} bound, and channels with one Ca^{2+} bound. The model also reveals that the amount of time required for the channel to fully equilibrate in one combination of (, ) before the ligand conditions can be changed in latency measurements, and that varying the time the channel stays in one ligand combination before perfusion solution switches can provide valuable information on kinetic rates. Furthermore, analysis of Ca^{2+} inhibition latencies should include measuring the duration between the time when is changed and the time when the channel enters the L mode, and the duration when the channel stays in the L mode until its activity terminates as it enters the Q mode. We suspect that there are significant surprises in store regarding the mechanistic interpretation of inhibition.

Although the aim here was to model the IP_{3}R Ca^{2+} channel, our methods can be easily applied to other channels and single-molecule dynamics.

## Acknowledgments

We would like to acknowledge the resources and support provided by J. Kevin Foskett and many useful conversations with Bill Bruno, Andy Fraser, Nick Hengartner, Ben McMahon, Paul Fenimore, Hans Frauenfelder, Joel Berenzden, Peter Hraber, Ian Parker, Jian-Wei Shuai, Silvina Ponce Dawson, Tony Auerbach, David Colquhoun, Fred Sachs, and Horia Vais.

This work was supported by National Institutes of Health (grant 5RO1GM065830-08).

Sharona E. Gordon served as editor.

## Footnotes

- Abbreviations used in this paper:
- CM
- cooperative model
- DYKM
- De Young–Keizer model
- IP
_{3} - inositol 1,4,5-trisphosphate
- IP
_{3}R - IP
_{3}receptor

- Submitted: 13 December 2011
- Accepted: 12 July 2012

This article is distributed under the terms of an Attribution–Noncommercial–Share Alike–No Mirror Sites license for the first six months after the publication date (see http://www.rupress.org/terms). After six months it is available under a Creative Commons License (Attribution–Noncommercial–Share Alike 3.0 Unported license, as described at http://creativecommons.org/licenses/by-nc-sa/3.0/).