Pacific Symposium on Biocomputing 13:354-365(2008) September 11, 2007 13:8 WSPC - Pro ceedings Trim Size: 9in x 6in DeRemigio MARKOV CHAIN MODELS OF COUPLED INTRACELLULAR CALCIUM CHANNELS: KRONECKER STRUCTURED REPRESENTATIONS AND BENCHMARK STATIONARY DISTRIBUTION CALCULATIONS HILARY DEREMIGIO , PETER KEMPER , M. DREW LAMAR , and GREGORY D. SMITH Department of Applied Science, Department of Computer Science, The Col lege of Wil liam and Mary, Wil liamsburg, VA 23187 Mathematical models of calcium release sites derived from Markov chain models of intracellular calcium channels exhibit collective gating reminiscent of the experimentally observed phenomenon of stochastic calcium excitability (i.e., calcium puffs and sparks). We present a Kronecker structured representation for calcium release site models and perform benchmark stationary distribution calculations using numerical iterative solution techniques that leverage this structure. In this context we find multi-level methods and certain preconditioned pro jection methods superior to simple Gauss-Seidel type iterations. Response measures such as the number of channels in a particular state converge more quickly using these numerical iterative methods than occupation measures calculated via Monte Carlo simulation. 1. Intro duction The sto chastic gating of voltage- and ligand-gated ion channels in biological membranes that is observed by single channel recording techniques is often mo deled using discrete-state continuous-time Markov chains (CTMCs).1,2 While these single channel models can be relatively simple (e.g., two physico chemically distinct states) or complex (hundreds of states), most include only two conductance levels (closed and open). For example, a transition state diagram for a three-state calcium (Ca2+ )-regulated channel activated by sequential binding of two Ca2+ ions is given by + + ka c kb c C1 C2 O1 - - kb ka (1) Pacific Symposium on Biocomputing 13:354-365(2008) September 11, 2007 13:8 WSPC - Pro ceedings Trim Size: 9in x 6in DeRemigio + - where ki c and ki with i {a, b} are transition rates with units of reciprocal + time, ki is an asso ciation rate constant with units of conc-1 time-1 , and c is the [Ca2+ ] near the channel. If this local [Ca2+ ] is specified, the transitionstate diagram of the channel (1) defines a CTMC that takes on values in the state-space (C1 , C2 , O1 ). The experimentally observable conductance of this sto chastically gating channel is the aggregated process of transitions between the closed and open classes of states: C = {C1 , C2 } and O = {O1 }. The scientific literature developing stochastic models for the behavior of ion channels is largely focused on single channels or populations of independent channels. One notable exception is the work of Ball and colleagues analyzing interacting aggregated CTMC models of membrane patches containing several ion channels.3,4 A second example, the sub ject of this paper, are simulations of clusters of intracellular Ca2+ -regulated Ca2+ channels--inositol 1,4,5-trisphosphate receptors (IP3 Rs) and ryanodine receptors (RyRs)--lo cated on the surface of the endoplasmic reticulum or sarcoplasmic reticulum membrane--that give rise to localized intracellular [Ca2+ ] elevations known as Ca2+ puffs and sparks.5­9 1 10 5 0.5 0 0 100 200 300 400 500 0 1 2 3 4 5 6 7 8 9 10 Fig. 1. Left: Local [Ca2+ ] near 3 × 3 µm ER membrane with Ca2+ channels modeled as 0.05 pA point sources with positions randomly chosen from a uniform distribution on a disc of radius 2 µm. Buffered Ca2+ diffusion is modeled as in Ref. 10. Middle: Stochastic Ca2+ excitability reminiscent of Ca2+ puffs/sparks. Right: Probability distribution of the number of open channels leading to a puff/spark Score of 0.39. When Markov chain mo dels of Ca2+ -regulated Ca2+ channels such as (1) are coupled via a mathematical representation of buffered diffusion of intracellular Ca2+ , simulated Ca2+ release sites may exhibit the phenomenon of "sto chastic Ca2+ excitability" where the IP3 Rs or RyRs open and close in a concerted fashion10,11 (see Fig. 1 for representative simulation) . Such models are sto chastic automata networks (SANs) that involve a large number of functional transitions, that is, the transition probabilities of one automata (i.e., an individual channel) may depend on the local [Ca2+ ] and thus the state of the other channels. The experimentally observable quantity is ei- Pacific Symposium on Biocomputing 13:354-365(2008) September 11, 2007 13:8 WSPC - Pro ceedings Trim Size: 9in x 6in DeRemigio ther the lo cal [Ca2+ ] or the number of channels in the open class of states, NO (t) (see Fig. 1, middle). The relationship between single channel kinetics of Ca2+ -regulated channels and the emergent phenomenon of Ca2+ puffs and sparks is not well understo o d. However, if each release site configuration is known, several informative response measures can be determined from the steady-state probability distribution. For example, the so-called puff/spark Score 10 given by Var[fO ]/E[fO ] is the index of dispersion of the steady-state fraction of open channels, fO = NO /N (see Fig. 1, right). This response measure takes values between 0 and 1, and a puff/spark Score of greater than approximately 0.3 indicates the presence of Ca2+ excitability. However, Ca2+ release sites are composed of 5­250 channels and this leads to a state-space explosion that makes numerical calculation of the stationary distribution of mo del Ca2+ release sites difficult. 2. Formulation of Mo del In this paper we consider two single channel models: the three-state Ca2+ activated channel described above (1) and a six-state model that includes both fast Ca2+ activation and slow Ca2+ inactivation, processes that are important aspects of the dynamics of both IP3 Rs and RyRs. The six-state mo del assumes two identical channel subunits that both require Ca2+ binding to enter a permissive state and include a second Ca2+ -mediated transition into a long-lived non-permissive state (for transition state diagram and parameters see Ref. 12). In both the three- and six-state models, Ca2+ -mediated transitions out of open states can be accelerated due to the increase in local [Ca2+ ] when a Ca2+ -regulated Ca2+ channel is open.13,14 Assuming the formation and collapse of Ca2+ microdomains are fast compared to channel gating, we can denote the background and domain [Ca2+ ] experienced by the channel when closed and open as c and cd , respectively. With this assumption the three- and six-state single channel models take the form Q = K- + (c I + cd IO ) K+ where K- and K+ are M × M matrices - + that collect the unimolecular (ki ) and bimolecular (ki ) transition rates, IO = diag {eO }, and eO is a M × 1 vector indicating open states of the single channel mo del.10 In our model formulation, the interaction between channels is mediated through the buffered diffusion of intracellular Ca2+ (see Ref. 10 for a complete description). For the purposes of this paper we do not assume any particular cell type with known release site ultrastructure (e.g., cardiac myo cytes with channels arranged in a dyad) and instead Pacific Symposium on Biocomputing 13:354-365(2008) September 11, 2007 13:8 WSPC - Pro ceedings Trim Size: 9in x 6in DeRemigio consider that the N channels at the Ca2+ release site have positions chosen from a two-dimensional uniform distribution on a disc of radius 0.1­2.0 µm (i.e., constant surface density; see Fig. 1, left). When in the open state, each channel contributes to the landscape of [Ca2+ ] throughout the Ca2+ release site and influences the local [Ca2+ ] experienced by other channels. For simplicity we assume that the formation and collapse of individual peaks within the Ca2+ micro domain occur quickly compared to channel gating. We also assume the presence of a single high concentration Ca2+ buffer and the validity of superposing local [Ca2+ ] increases due to each of the N channels.15,16 Thus, channel interactions can be summarized by an N × N `coupling matrix' C = (cij ) that gives the increase over c experienced by channel j when channel i is open. 2.1. Instantaneous Coupling of Two Ca2+ -Regulated Ca2+ Channels In the case of two identical Ca2+ -regulated Ca2+ channels the interaction matrix takes the form a c d c12 C= c21 cd nd the expanded generator matrix is given by Q(2) = Q- + Q+ where Q- = K - I + I K - (2) (2) (2) (2) collects the unimolecular transition rates and denotes the Kronecker pro duct (see Ch. 9 in Ref. 17). The transition rates involving Ca2+ take the form Q+ = D1 (K+ I ) + D2 (I K+ ) , (2) (2) (2) (3) where the two terms represent Ca2+ -mediated transitions of each channel. (2) (2) The diagonal matrices D1 and D2 give the [Ca2+ ] experienced by channel 1 and 2, respectively, in every configuration of the release site, that is, D1 = diag {c (e e) + cd (eO e) + c21 (e eO )} = c (I I ) + cd (IO I ) + c21 (I IO ) and similarly for D2 . Using the Kronecker identities such as (I IO ) (I K+ ) = I IO K+ , Eq. 3 can be rearranged as Q+ = c K+ + cd (IO K+ I ) + c12 (IO K+ ) + c21 (K+ IO ) + cd (I IO K+ ) (4) (2) (2) (2) (2) Pacific Symposium on Biocomputing 13:354-365(2008) September 11, 2007 13:8 WSPC - Pro ceedings Trim Size: 9in x 6in DeRemigio where K+ = K+ I + I K+ . Combining Eqs. 2 and 4 and simplifying, Q(2) can be written compactly as Q(2) = Ad I + IO A12 + A21 IO + I Ad where Ad = K- + c K+ + cd IO K+ , and Aij = cij K+ . 2.2. Instantaneous Coupling of N Ca2+ -Regulated Ca2+ Channels In the case of N channels coupled at the Ca2+ release site, the expanded generator matrix--i.e., the SAN descriptor--is given by Q (N ) = Q - + Q + Q- Q+ (N ) (N ) (N ) (N ) (2) (5) (6) (7) ( 8) (9) = nN K- = N =1 nN cij I (n-1) K- I (N -n) Y1 1 ij Zij N · · · YiN Zij j =1 = c K + I O (N ) + i , j =1 Yin = j for i = n I otherwise n Zij = K + for j = n I otherwise (N ) where I (n) is an identity matrix of size M n and K+ = N=1 K+ . Comn bining Eqs. 7 and 8 and simplifying, Q(N ) can be written as N Q (N ) = A d for i = n I otherwise. i Xi1 · · · XiN j j IO for i = n Xin = Aij for j = n j I otherwise (10) , j =1 Xin = i and for i = j where Ad = K- + c K+ + cd IO K+ , and Aij = cij K+ . Note that all states of the expanded Markov chain Q(N ) are reachable, the matrices I , IO , Ad , Aij , and Xin are all M × M , and 2N 2 - N of the N 3 matrices denoted by j Xin are not identity matrices. j 3. Stationary Distribution Calculations The limiting probability distribution of a finite irreducible CTMC is the unique stationary distribution (N ) satisfying global balance,17 that is, (N ) Q (N ) = 0 sub ject to ( N ) e( N ) = 1 (11) Pacific Symposium on Biocomputing 13:354-365(2008) September 11, 2007 13:8 WSPC - Pro ceedings Trim Size: 9in x 6in DeRemigio where Q(N ) is the Ca2+ release site SAN descriptor (Eq. 10) and e(N ) is an M N × 1 column vector of ones. Although Monte Carlo simulation techniques such as Gillespie's Metho d18 can be implemented to estimate response measures such as the puff/spark Score, this is an inefficient approach when the convergence of the o ccupation measures to the limiting probability distribution is slow. This problem is compounded by the state-space explosion that o ccurs when the number of channels (N ) or number of states per channel (M ) is large (i.e., physiologically realistic). Both space requirements and quality of results can be addressed using the Kronecker representation (Eq. 10) and various iterative numerical methods to solve for (N ) . Many metho ds are available to solve Eq. 11 with different ranges of applicability (see Ref. 17 for review). For larger models, a variety of iterative metho ds are applicable including the method of Jacobi, and GaussSeidel, along with variants that use relaxation, e.g., Gauss-Seidel with relaxation (SOR). Such metho ds require space for iteration vectors and Q(N ) but usually converge quickly. More sophisticated pro jection methods--e.g., the generalized minimum residual method (GMRES) and the method of Arnoldi (ARNOLDI)--have better convergence properties but require more space. While the best metho d for a particular Markov chain is unclear in general, several options are available for exploration including the iterative metho ds described above that can be also enhanced with preconditioning, aggregation-disaggregation (AD), or Kronecker-specific multi-level (ML) metho ds that are inspired by multigrid and AD techniques. Unfortunately, we cannot acknowledge all relevant work on iterative methods due to limited space.19,20 A number of software tools are available that implement methods for Kronecker representations, and we selected the APNN toolbox21 and its numerical solution package Nsolve for its rich variety of numerical techniques for the steady state analysis of Markov chains. Nsolve provides more than 70 different metho ds and comes with an ASCII file format for a SAN descriptor easily interfaced with our MATLAB modeling environment. Nsolve mainly supports hierarchical Markovian models that include a trivial hierarchy with a single macrostate such as Eq. 10 as a special case (see Refs. 21­24). 4. Results In order to investigate which numerical techniques work best for the Kronecker representation of our Ca2+ release site models, we wrote a script for the matrix computing to ol MATLAB that takes a specific Ca2+ release site mo del--defined by K+ , K- , eO , c , and C --and produces the input Pacific Symposium on Biocomputing 13:354-365(2008) September 11, 2007 13:8 WSPC - Pro ceedings Trim Size: 9in x 6in DeRemigio files needed to interface with Nsolve. Using 10 three-state channels (1) we performed a preliminary study to determine which of the 70-plus numerical metho ds implemented in Nsolve were compatible with Eq. 10. 4.1. Benchmark Stationary Distribution Calculations Table 1 lists those solvers that converged in less than 20 minutes CPU time with a maximum residual less than 10-12 for one configuration of 10 three-state channels. For each method we report the maximum and sum of the residuals, the CPU and wall clock times (in seconds), and the total number of iterations performed. We find that traditional relaxation metho ds (e.g., JOR, RSOR) work well for this problem with 310 = 59, 049 states, but the addition of AD steps is not particularly helpful. AD steps do however greatly improve the performance of the GMRES solver and to a smaller extent the DQGMRES and ARNOLDI methods. The separable preconditioner (PRE) of Buchholz23 and the BSOR preconditioner are very effective and help to reduce solution times to less than 50 seconds for several pro jection methods. Among ML solvers, a JOR smoother gives the best results and dynamic (DYN) or cyclic (CYC) ordering is better than a fixed (FIX) order where V, W, or F indicate the type of cycle used.19,20 4.2. Problem Size and Method Performance In Sec. 4.1 we benchmarked the efficiency of several different algorithms that can be used to solve for the stationary distribution of Ca2+ release site models. To determine if this result depends strongly on problem size, we chose representatives of four classes of solvers (JOR, PRE ARNOLDI, BSOR BICGSTAB, and ML JOR F DYN) that worked well for release sites composed of 10 threestate channels (see Table 1). Using these four methods, Fig. 2 shows the wall clo ck time required for convergence of (N ) as a function of the number of channels (N ) for both the three- and six-state models (circles and squares, respectively). Because the N channels in each Ca2+ release site simulation have randomly chosen positions that may influence the time to convergence, Fig. 2 shows both the mean and standard deviation (error bars) of the wall clo ck time for five different release site configurations. Note that for each value of N in Fig. 2, the radius of each Ca2+ release site was chosen so that sto chastic Ca2+ excitability was observed. Due to irregular release site ultrastructure, these calculations can not be simplified using spatial symmetries. Figure 2 shows that the time until convergence is shorter when the Ca2+ Pacific Symposium on Biocomputing 13:354-365(2008) September 11, 2007 13:8 WSPC - Pro ceedings Trim Size: 9in x 6in DeRemigio Table 1. Benchmark calculations for 10 three-state channels computed using Linux PCs with dual core 3.8GHz EM64T Xeon processors and 8GB RAM solving Eq. 10. Solver JOR SOR RSOR JOR AD SOR AD DQGMRES ARNOLDI BICGSTAB GMRES AD DQGMRES AD ARNOLDI AD PRE POWER PRE GMRES PRE ARNOLDI PRE BICGSTAB BSOR BICGSTAB BSOR GMRES BSOR TFQMR PRE GMRES AD PRE ARNOLDI AD ML JOR V FIX ML JOR W FIX ML JOR F FIX ML JOR V CYC ML JOR W CYC ML JOR F CYC ML JOR V DYN ML JOR W DYN ML JOR F DYN Max Res 9.49E-13 9.49E-13 8.76E-13 9.44E-13 9.44E-13 9.87E-13 2.42E-13 8.66E-13 6.43E-13 1.03E-12 7.23E-13 9.37E-13 8.62E-15 8.62E-15 4.44E-16 8.22E-15 3.05E-13 1.83E-13 1.29E-13 4.32E-13 9.69E-13 9.12E-13 9.93E-13 8.35E-13 4.36E-13 6.76E-13 8.07E-13 2.81E-13 5.87E-13 Sum Res 5.16E-12 5.16E-12 2.40E-12 5.13E-12 5.13E-12 6.78E-10 4.04E-11 4.89E-11 3.61E-11 1.84E-10 7.60E-11 5.27E-12 3.73E-12 1.82E-12 2.49E-14 5.29E-13 7.73E-12 1.39E-12 1.52E-11 7.18E-12 3.54E-11 1.14E-10 1.01E-10 6.36E-12 5.41E-11 1.39E-11 6.09E-12 5.15E-11 1.68E-10 CPU 279 435 1190 415 413 490 214 146 88 184 109 246 45 26 28 19 20 17 36 27 105 156 146 42 26 18 58 14 15 Wall 279 436 1197 415 414 492 215 148 89 184 109 247 46 27 28 19 20 17 36 28 105 157 146 43 26 19 59 15 15 Iters 1840 1840 990 1550 1550 2940 1440 602 900 2008 1280 1670 180 160 188 52 49 48 140 140 372 326 330 168 38 56 152 38 46 release site is composed of three-state as opposed to six-state channels regardless of the numerical method used (compare circles to squares). Consistent with Table 1 we find that for large values of N the ML JOR F DYN (black) metho d requires the shortest amount of time, followed by BSOR BICGSTAB (dark gray), PRE ARNOLDI (light gray), and finally JOR (white). Though there are important differences in the speed of the four solvers, the wall clo ck time until convergence is approximately proportional to the number of states (M N ), that is, the slope of each line in Fig. 2 is nearly M = 3 or 6 depending on the single channel model used. We also experienced substantial differences in the amount of memory needed to run those solvers. While simple methods like JOR and SOR allocate space mainly for a few iteration vectors, Krylov subspace methods like Pacific Symposium on Biocomputing 13:354-365(2008) September 11, 2007 13:8 WSPC - Pro ceedings Trim Size: 9in x 6in DeRemigio 10 5 10 4 Wall Clock Time (sec) 10 3 10 2 10 1 10 0 10 -1 10 -2 4 5 6 7 8 9 10 11 12 N Fig. 2. Circles and error bars show the mean ± SD of wall clock time for five release site configurations of the three-state model (1) using: JOR (white), PRE ARNOLDI (light gray), + BSOR BICGSTAB (dark gray), and ML JOR F DYN (black). Three-state model parameters: ka -1 ms-1 , k - = 50 ms-1 , k + = 150 µM-1 ms-1 , k - = 1.5 ms-1 . Squares and = 1.5 µM a b b error bars give results for the six-state model (parameters as in Ref. 12). Calculations performed using 2.66 GHz Dual-Core Intel Xeon processors and 2 GB RAM. GMRES, DQGMRES and ARNOLDI use more vectors (20 in the default Nsolve configuration), and this can be prohibitive for large models. For pro jection metho ds that operate on a fixed and small set of vectors like TFQMR and BICGSTAB, we observe that the space for auxiliary data structures and vectors is on the order of 7­10 iteration vectors for these models. In general we find that the iterative numerical methods that incorporate pre-conditioning techniques are quite fast compared to more traditional relaxation techniques such as JOR. However, the power of pre-conditioning is only evident when problem size is less than some threshold that depends upon memory limitations. On the other hand, ML methods are constructed to take advantage of the Kronecker representation and to have very modest memory requirements. This is consistent with our experiments that indicate ML methods have the greatest potential to scale well with problem size, whether that be an increase in the number of channels (N ) or the number of states per channel (M ). 4.3. Comparison of Iterative Methods and Monte Carlo Simulation Although there may be problem size limitations, we expected that the stationary distribution of our Ca2+ release site models could be found more quickly using iterative methods than Monte Carlo simulation. This is confirmed in the convergence results of Fig. 3 using a release site composed of Pacific Symposium on Biocomputing 13:354-365(2008) September 11, 2007 13:8 WSPC - Pro ceedings Trim Size: 9in x 6in DeRemigio 10 10 10 10 2 0 -2 -4 Error 10 10 10 10 10 -6 -8 -10 -12 -14 0 1 2 3 10 10 10 10 Wall Time (sec) Fig. 3. Convergence of response measures for a release site composed of 10 three-state channels using ML JOR F DYN and Monte Carlo (fil led and open symbols, respectively). Circles and squares give 1- and -norms of the residual errors, upper pointing triangles give the relative error in the puff/spark S cor e for Monte Carlo (mean of 50 simulations shown) compared with the S cor e given by ML JOR F DYN upon convergence. Similarly, the lower pointing triangles give the relative error in the probability that all N channels are closed. Parameters as in Fig. 1. 10 three-state channels for both ML JOR F DYN (filled symbols) and Monte Carlo simulation (open symbols). We run a Monte Carlo simulation to estimate the stationary distribution and that estimate depends on the length of the simulation measured in seconds of wall clock time (our implementation averaged 1,260 transitions per second). The simulation starts with all N channels in state C1 --chosen because it is the most likely state at the background [Ca2+ ] (c ). Figure 3 shows the maximum and sum of 1and -norms of the residuals averaged over 50 simulations. As expected, the residuals asso ciated with the Monte Carlo simulations converge much slower than those obtained with ML JOR F DYN. Interestingly, Fig. 3 shows that even coarse response measures can be more quickly obtained using numerical iterative metho ds than Monte Carlo simulation. We find that the relative errors of the puff/spark Score (upwards pointing triangles) and the probability that all N channels were closed (downwards pointing triangles) obtained via Monte Carlo simulation did not converge significantly faster than the maximum residual error (open squares). 5. Conclusions We have presented a Kronecker structured representation for Ca2+ release sites composed of Ca2+ -regulated Ca2+ channels under the assumption that these channels interact instantaneously via the buffered diffusion of intra- Pacific Symposium on Biocomputing 13:354-365(2008) September 11, 2007 13:8 WSPC - Pro ceedings Trim Size: 9in x 6in DeRemigio cellular Ca2+ (Sec. 2). Because informative response measures such as the puff/spark Score can be determined if the steady-state probability of each release site configuration is known, we have identified numerical interative solution techniques that perform well in this biophysical context. The benchmark stationary distribution calculations presented here indicate significant performance differences among iterative solution methods. Multi-level metho ds provide excellent convergence with modest additional memory requirements for the Kronecker representation of our Ca2+ release site mo dels. When the available main memory permits, BSOR-preconditioned pro jection metho ds such as TFQMR and BICGSTAB are also effective, as is the metho d of Arnoldi combined with a simple preconditioner. In case of tight memory constraints, Jacobi and Gauss-Seidel iterations are also possible (but slower). When numerical iterative methods apply, they outperform our implementation of Monte Carlo simulation for estimates of response measures such as the puff/spark Score and the probability of a number of channels being in a particular state. Single channel mo dels of IP3 Rs and RyRs can be significantly more complicated than the three- and six-state models that are the focus of this manuscript. For example, the well-known DeYoung-Keizer IP3 R model includes four eight-state subunits per channel for a total of 330 distinguishable states.25 Because biophysically realistic Ca2+ release site simulations can involve tens or even hundreds of intracellular channels, we expect that the development of approximate methods for our SAN descriptor (Eq. 10) will be an important aspect of future work. Of course, some puff and spark statistics--such as puff/spark duration and inter-event interval distributions--cannot be determined from the Ca2+ release site stationary distribution. Consequently, it will be important to determine if transient analysis can also be accelerated by leveraging the Kronecker structure of Ca2+ release sites composed of instantanteously coupled Ca2+ -regulated Ca2+ channels. Furthermore, although the SAN conceptual framework and its asso ciated analysis techniques presented in this manuscript have focused solely on the emergent dynamics of Ca2+ release sites, it is also important to note that these techniques should be generally applicable to our understanding of signaling complexes of other kinds.26,27 Acknowledgments The authors thank Buchholz and Dayar for sharing their implementation of Nsolve. This material is based upon work supported by the National Science Foundation under Grants No. 0133132 and 0443843. Pacific Symposium on Biocomputing 13:354-365(2008) September 11, 2007 13:8 WSPC - Pro ceedings Trim Size: 9in x 6in DeRemigio References 1. D. Colquhoun and A. Hawkes, A Q-matrix cookb ook: how to write only one program to calculate the sigle-channel and macroscopic predictions for any kinetic mechanism, in Single-Channel Recording , eds. B. Sakmann and E. Neher (Plenum Press, New York, 1995) pp. 589­633. 2. G. Smith, Modeling the stochastic gating of ion channels, in Computational Cel l Biology , eds. C. Fall, E. Marland, J. Wagner and J. Tyson (SpringerVerlag, 2002) pp. 291­325. 3. F. Ball, R. Milne, I. Tame and G. Yeo, Advances in App Prob 29, 56 (1997). 4. F. Ball and G. Yeo, Methodology and Computing in App Prob 2, 93 (1999). 5. H. Cheng, W. Lederer and M. Cannell, Science 262, 740 (1993). 6. H. Cheng, M. Lederer, W. Lederer and M. Cannell, Am J Physiol 270, C148 (1996). 7. Y. Yao, J. Choi and I. Parker, J Physiol 482, 533 (1995). 8. I. Parker, J. Choi and Y. Yao, Cel l Calcium 20, 105 (1996). 9. M. Berridge, J Physiol (London) 499, 291 (1997). 10. V. Nguyen, R. Mathias and G. Smith, Bul l. Math. Biol. 67, 393 (2005). 11. S. Swillens, G. Dup ont, L. Comb ettes and P. Champ eil, Proc Natl Acad Sci USA 96, 13750(Nov 1999). 12. H. DeRemigio, P. Kemp er, M. LaMar and G. Smith, Technical Report WMCS-2007-06 (2007). 13. G. Smith, An extended DeYoung-Keizer-like IP3 receptor model that accounts for domain Ca2+ -mediated inactivation, in Recent Research Developments in Biophysical Chemistry, Vol. II , eds. C. Condat and A. Baruzzi (Research Signp ost, 2002). 14. I. Bezprozvanny, Cel l Calcium 16, 151 (1994). 15. M. Naraghi and E. Neher, J Neurosci 17, p. 6961(6973 1997). 16. G. Smith, L. Dai, R. Muira and A. Sherman, SIAM J Appl Math 61, 1816 (2001). 17. W. Stewart, Introduction to the Numerical Solution of Markov Chains (Princeton University Press, Princeton, 1994). 18. D. Gillespie, J Comp Phys 22, 403 (1976). 19. P. Buchholz and T. Dayar, Computing 73, 349 (2004). 20. P. Buchholz and T. Dayar, SIAM Matrix Analysis and App (to appear) (2007). 21. P. Buchholz and P. Kemp er, A toolb ox for the analysis of discrete event dynamic systems., in CAV. LNCS 1633 , 1999. 22. P. Buchholz and T. Dayar, SIAM J. Sci. Comput. 26, 1289 (2005). 23. P. Buchholz, Pro jection methods for the analysis of stochastic automata networks, in Numerical Solution of Markov Chains , eds. B. Plateau, W. Stewart and M. Silva (Prensas Univerversitarias de Zaragoza, 1999) pp. 149­168. 24. P. Buchholz, Prob in the Eng and Informational Sci 11, 229 (1997). 25. G. De Young and J. Keizer, Proc Natl Acad Sci USA 89, 9895 (1992). 26. J. Schlessinger, Cel l 103, 211 (2000). 27. H. Husi, M. A. Ward, J. S. Choudhary, W. P. Blackstock and S. G. Grant, Nat Neurosci 3, 661 (2000).