- Email: [email protected]
A comparison of three models of the inositol trisphosphate receptor
ARTICLE IN PRESS
Progress in Biophysics & Molecular Biology 85 (2004) 121–140
A comparison of three models of the inositol trisphosphate receptor J. Sneyda,*, M. Falckeb, J.-F. Dufourc, C. Foxa a
Department of Mathematics, University of Auckland, Auckland, New Zealand b Hahn Meitner Institute, SF5, Glienicker Str. 100, 14109 Berlin, Germany c University of Bern, Switzerland
Abstract The inositol (1,4,5)-trisphosphate receptor (IPR) plays a crucial role in calcium dynamics in a wide range of cell types, and is often a central feature in quantitative models of calcium oscillations and waves. We compare three mathematical models of the IPR, fitting each of them to the same data set to determine ranges for the parameter values. Each of the fits indicates that fast activation of the receptor, followed by slow inactivation, is an important feature of the model, and also that the speed of inositol trisphosphate ðIP3 Þ binding cannot necessarily be assumed to be faster than Ca2þ activation. In addition, the model which assumed saturating binding rates of Ca2þ to the IPR demonstrated the best fit. However, lack of convergence in the fitting procedure indicates that responses to step increases of ½Ca2þ and ½IP3 provide insufficient data to determine the parameters unambiguously in any of the models. r 2004 Elsevier Ltd. All rights reserved. Keywords: Calcium oscillations; Mathematical model; Markov chain Monte Carlo; Data fitting
1. Introduction In almost every cell type, the cytosolic concentration of free calcium ions ð½Ca2þ Þ exhibits complex dynamic behaviour, including oscillations and periodic waves (Sanderson et al., 1994; Clapham, 1995; Thomas et al., 1996; Berridge, 1997). Often, the physiological significance of such phenomena is not completely understood, but there is general acceptance that these oscillations and waves are one way in which cells can use Ca2þ as an intracellular signalling mechanism, while avoiding the toxic effects of prolonged raised ½Ca2þ : In electrically non-excitable cells, Ca2þ release is often mostly from internal stores, although influx from the outside can play an important modulatory role. Agonist stimulation results in the formation of inositol (1,4,5)-trisphosphate *Corresponding author. Fax: +64-9-3737-457. E-mail address: [email protected] (J. Sneyd). 0079-6107/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.pbiomolbio.2004.01.013
ARTICLE IN PRESS 122
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
ðIP3 Þ which diffuses through the cytoplasm and binds to the IP3 receptors (IPR) situated on the endoplasmic reticulum (ER). Subsequent release of Ca2þ from the ER can further modulate the open probability of the IPR, with cytosolic (and possibly lumenal) Ca2þ exerting both a positive and negative feedback effect on the IPR and on the formation of IP3 : It is these dynamic interactions, when coupled with Ca2þ removal processes and with other Ca2þ release pathways through, for instance, ryanodine receptors, that can result in oscillations and waves of ½Ca2þ : Thus, an understanding of IP3 receptors is central to a detailed understanding of Ca2þ dynamics in many cell types, and a number of IPR models have been constructed. The obvious questions to ask are, firstly, in what way have these models helped us understand the physiology or biophysics of the IPR, and, secondly, which is the best model? The first question is, we think, relatively easily answered. The models have contributed a great deal to our understanding of how the binding properties of the IPR underlie steady-state responses and adaptation, and have provided quantitative frameworks that have aided in the design of experiments. Not all of the useful models have been mathematically expressed. The qualitative, or diagrammatic, models proposed by Taylor (1998) (Adkins and Taylor, 1999; ! Taylor and Laude, 2002) and Hajnoczky and Thomas (1997) have led to a great deal of experimental work and been crucial in the development of our understanding of the IPR. Many of the ideas in those models have been incorporated into later mathematical models, particularly that of Sneyd and Dufour (2002). The second question is more difficult to answer. Indeed, it can only be satisfactorily answered by pointing out that (as has been said many times) all models are wrong, but some models are useful. Thus, the search for the ‘‘best’’ model is a futile one. Each IPR subtype behaves differently from the others, and its behaviour is modified by the environment of the IPR, making different models for different situations a necessity. Nevertheless, despite these complications, if we restrict ourselves to a less ambitious goal, models can be usefully compared. If we take a single data set, we can investigate how well a variety of models fit those data, and thus compare models with respect to the particular features of most importance in the data. Because of the effort involved, such detailed model comparisons are rarely done, but this does not make them any less important. We shall take a single data set (Dufour et al., 1997), and use it to determine the parameters in three models: the models of De Young and Keizer (1992), Sneyd and Dufour (2002) and Dawson et al. (2003). To determine the parameters, we performed parameter estimation within a Bayesian formalism computed via Markov chain Monte Carlo (MCMC) sampling from the posterior distribution for parameters. Examining the posterior not only gives ‘‘best fit’’ parameter values but also a quantitative idea of how well parameters are determined and how good a representation of the data the model is. The MCMC sampler was implemented using Metropolis–Hastings dynamics (Gilks et al., 1995).
2. The three models 2.1. The De Young–Keizer model In the early 1990s, two crucial features of the IPR became apparent. Firstly, it was shown (Iino, 1990; Parker and Ivorra, 1990; Bezprozvanny et al., 1991; Parys et al., 1992) that, at steady state,
ARTICLE IN PRESS J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
123
the open probability of the IPR is a bell-shaped function of ½Ca2þ ; with the maximum open probability occurring at about ½Ca2þ ¼ 300–500 nM: Secondly, it became clear that, in response to an increased ½Ca2þ ; the IPR responded in an adaptive manner, first activating, and subsequently inactivating (Parker and Yao, 1991; Finch et al., 1991; Iino and Endo, 1992; Yao and Parker, 1992; Iino and Tsukioka, 1994; Ilyin and Parker, 1994). These results were most easily interpreted by assuming that Ca2þ activates the IPR on a fast time scale, but then inactivates the IPR over a slower time scale. The first model of the IPR to take this into account was due to De Young and Keizer (1992), closely followed by the phenomenological model of Atri et al. (1993), and the sequential binding models of Othmer and Tang (1993) and of Bezprozvanny (1994). Here, we focus on the De Young–Keizer model, as, of these models, it has been the most influential. De Young and Keizer (1992) assumed that the IPR consists of three equivalent and independent subunits, all of which must be in a conducting state before the receptor allows Ca2þ flux. Each subunit has an IP3 binding site, an activating Ca2þ binding site, and an inactivating Ca2þ binding site, each of which can be either occupied or unoccupied, and thus each subunit can be in one of eight states. Each state of the subunit is labelled Sijk ; where i; j and k are equal to 0 or 1, with a 0 indicating the binding site is unoccupied and a 1 indicating it is occupied. The first index refers to the IP3 binding site, the second to the Ca2þ activation site, and the third to the Ca2þ inactivation site. This is illustrated in Fig. 1. Although a fully general model would include 24 rate constants, De Young and Keizer made a number of simplifying assumptions. The rate constants are assumed to be independent of whether or not activating Ca2þ is bound, and the rate of Ca2þ activation is assumed to be independent of whether or not IP3 is bound. Finally, Ca2þ inactivation is assumed to be independent of Ca2þ activation and of IP3 : This leaves only 10 rate constants, k71 ; y; k75 : The fraction of subunits in the state Sijk is denoted by xijk : The differential equations for these are based on mass action kinetics, and thus, for example, dx000 ¼ ðk1 x100 k1 px000 Þ þ ðk4 x001 k4 cx000 Þ þ ðk5 x010 k5 cx000 Þ; dt
ð1Þ
where p denotes ½IP3 and c denotes ½Ca2þ : The model assumes that the IP3 receptor passes Ca2þ current only when three subunits are in the state S110 (i.e., with one IP3 and one activating Ca2þ bound), and thus the open probability of the receptor is x3110 : This gives seven differential equations for the receptor states. k1p S110
S010
k2c
k-5 k-4
k-4
S011 k-3
k4c
k-5 k-2
k-2
k2c
S111
S101 k-3
k2c
k5c
k3p S001
k-5
S110
k-5
k5c
k3p S111
k4c
k5c S100
S000
k-1 k-2
k1p
k1c
k-5
Fig. 1. The binding diagram for the De Young–Keizer IPR model. c denotes ½Ca2þ ; and p denotes ½IP3 :
ARTICLE IN PRESS 124
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
Parameters of the model were chosen to obtain agreement with the steady-state data of Bezprozvanny et al. (1991). Nevertheless, the bell-shaped nature of the steady-state IPR open probability curve has very little, if anything, to do with the fact that the IPR can mediate oscillatory activity. As was originally pointed out by Keizer (see Sneyd et al., 1995, Fig. 2), during an oscillation the steady-state response is less important than the response to a changing ½Ca2þ : After all, the cell never ‘‘sees’’ a steady ½Ca2þ ; unless the time scale of the oscillation is very much slower than the time scale of the IPR kinetics. It is thus the fast activation, followed by the slow inactivation (k5 > k2 ) that is the crucial feature of the De Young–Keizer model and this separation of time scales is independent of the shape of the steady-state open probability curve. IPR models with monotonic steady-state curves support Ca2þ oscillations just as readily as IPR models with a bell-shaped curve (LeBeau et al., 1999). 2.2. The Sneyd–Dufour model Finch et al. (1991) were the first to develop a superfusion system in which the IPRs, isolated in microsomal vesicles, could be exposed to solutions in which the ½Ca2þ and ½IP3 were tightly controlled by a combination of fast flow and extensive buffering. This procedure gives a calcium clamp, similar to the voltage clamp used by Hodgkin and Huxley (1952). Since the work of Finch et al., other groups have used a superfusion calcium clamp (Dufour et al., 1997; Marchant et al., 1997; Marchant and Taylor, 1997, 1998; Adkins and Taylor, 1999) and now recordings with a time resolution of 9 ms are possible (Marchant and Taylor, 1997; Swatton and Taylor, 2002). Similar data were also obtained by Parker et al. (1996), although, because their results were obtained in vivo in Xenopus oocytes, they were not able to maintain a precise calcium clamp. Time-dependent data were used by Sneyd and Dufour (2002) in a model that incorporated several ideas suggested by experiments: sequential binding of IP3 and activating Ca2þ ; modulation of IP3 binding by Ca2þ ; a fit of the model to dynamical data instead of steady-state data only, time-dependent inactivation upon IP3 binding, and saturating binding rates of Ca2þ : The model scheme is shown in Fig. 2. The bare subunit, state R; can bind IP3 and move to state O: The binding rate depends on Ca2þ allowing for modulation of IP3 binding in a facilitating or inhibiting way depending on receptor type (Yoneshima et al., 1997; Cardy et al., 1997; Sienaert et al., 1997). Binding of Ca2þ to the activating site corresponds to a transition to state A: Since A can be reached from O only and not from R; IP3 and activating Ca2þ need to bind sequentially. Alternatively, the subunit could close to state S; corresponding to time-dependent inactivation upon IP3 binding. There are two inhibited states I1 and I2 : The state I1 has Ca2þ bound at the inhibiting site only whereas I2 also has Ca2þ bound at the activating site and IP3 : The receptor is open when four subunits are in state O or A or any combination of these states. However, the receptors in which some of the subunits are in the O state have a lower conductance. Hence, the IPR flux was chosen to be proportional to ðoO þ aAÞ4 with o in the range of 0.05–0.10 and a about 0.9. This accounts for the non-zero fluxes observed with only IP3 bound (Finch et al., 1991; Marchant and Taylor, 1997; Combettes et al., 1993; Dufour et al., 1997), and predicts that multiple conductance levels would be observed. Alternatively, to avoid having multiple conductance levels, it could be assumed that each state of the subunit flickers to and from the open configuration but that some states just have a higher probability of being in the open configuration.
ARTICLE IN PRESS J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
125
1(c) R
k-1+l-2
I1
-2(c)
2(c) p
4(c) O
-4(c)
5(c) A
I2 k-1+l-2
k-3
3(c) S
Fig. 2. Transition diagram of the receptor subunit in the Sneyd–Dufour model.
The functions in the model were derived by simplification of a more complex model, the full details of which are in Sneyd and Dufour (2002). This gives ðk1 L1 þ l2 Þc ; ð2Þ f1 ðcÞ ¼ L1 þ cð1 þ L1 =L3 Þ f2 ðcÞ ¼
k2 L3 þ l4 c ; L3 þ cð1 þ L3 =L1 Þ
f2 ðcÞ ¼
k2 þ l4 c ; 1 þ c=L5
ð3Þ ð4Þ
f3 ðcÞ ¼
k3 L5 ; c þ L5
ð5Þ
f4 ðcÞ ¼
ðk4 L5 þ l6 Þc ; L5 þ c
ð6Þ
f4 ðcÞ ¼ f1 ðcÞ ¼
L1 ðk4 þ l6 Þ ; L1 þ c
ðk1 L1 þ l2 Þc : L1 þ c
ð7Þ ð8Þ
One principal feature of the Sneyd–Dufour model is that the transition rates saturate as c increases. This allows for a limited range of activation speeds over a wide range of Ca2þ concentrations, a feature that is absent from the other two models. Saturating rates are obtained by a process essentially identical to the Michaelis–Menten model of enzyme kinetics (Sneyd and Dufour, 2002). The method is illustrated in Fig. 3. If we assume that A* and A% are in instantaneous
ARTICLE IN PRESS 126
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
k1c
~ A
|
k -1 l -1 l1c
l2
l -2
_ A
(c)
|
A k-1 + l -2
Fig. 3. Derivation of a saturating binding rate. If A% and A* are assumed to be in instantaneous equilibrium, the upper transition diagram reduces to the lower, as described in the text.
equilibrium, we get % cA* ¼ L1 A; from which it follows that L1 ¼ l1 =l1 : Thus, letting A ¼ A% þ A* we get dA ¼ ðk1 þ l2 ÞI fðcÞA; dt where cðk1 L1 þ l2 Þ : fðcÞ ¼ c þ L1
ð9Þ
ð10Þ
ð11Þ
Note that, in this model, the IPR can inactivate either in the absence of Ca2þ ; by going to state S; or in the absence of IP3 ; by going to state I1 ; both of which features are required by experimental data (Taylor, 1998). Furthermore, once IP3 is bound (and the receptor is in state O), there is an intrinsic competition between states A and S: If Ca2þ is abundant, the receptor will go preferentially to the activated state A; leading to a higher conductance (or higher open probability, depending on interpretation), but in conditions of low ½Ca2þ the receptor will move instead to the S state, which does not allow Ca2þ current. In a number of important respects, the Sneyd–Dufour model is a quantitative realisation of the qualitative model of Taylor (1998) (Adkins and Taylor, 1999; Marchant and Taylor, 1997, 1998; Taylor and Laude, 2002). However, one major difference is that the Sneyd–Dufour model does not assume that binding of IP3 protects the receptor from inactivation by Ca2þ (i.e., that there is no I2 state), as is proposed by Taylor for type 2 IPR (but not type 3). This proposal was not adopted by Sneyd and Dufour because with it they were unable to obtain acceptable fits to experimental data (unpublished results). One should not put undue weight on this negative result as it is possible that excellent fits do exist, just they have not been found yet. Another possibility is that ER depletion is playing a more important role than initially thought and thus one cannot expect IPR dynamics alone to fit the data well. To tell
ARTICLE IN PRESS J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
127
for sure if the I2 state is necessary will require detailed fits of the alternative models to data collected at the highest possible time resolution, and from types 2 and 3 IPR, as well as a detailed consideration of the possible effects of ER depletion and unstirred boundary layers close to the IPR. 2.3. The Dawson–Lea–Irvine model A different type of model, based upon the allosteric model of Monod et al. (1965), has been constructed by Dawson et al. (2003), with the goal of understanding adaptation and incremental responses by the IPR. They assumed that the IPR could exist in two conformations, R and T; R binding IP3 rapidly, but with low affinity, and T binding IP3 slowly, but with high affinity. Both R and T states can bind four IP3 molecules, and the only open state of the receptor is when the R conformation has four IP3 molecules bound. Finally, it is assumed that Ca2þ on the cytosolic side of the IPR can inactivate the receptor by binding to the R and T states, or activate the receptor by converting the open state, O1 ; to another open state, O2 ; that cannot be shifted across into the closed state. This leads to the state diagram shown in Fig. 4. The O1 to O2 transition gives positive feedback of Ca2þ on the IPR, while the R (or T) to I transition lets the receptor inactivate in the absence of IP3 but in the presence of Ca2þ ; as observed experimentally. The balance between the two conformational states provides fast inactivation and slower inactivation. Upon initial addition of IP3 ; the receptor is quickly pushed through to the O1 state (and, in the positive feedback process, to the O2 state). However, over a longer time scale, the receptor shifts across to the T conformation, pulling the receptor out of the O1 state, and leading to inactivation. To model Ca2þ fluxes, Dawson et al. connected their IPR model to a simplified model of Ca2þ handling, in which they kept track of ½Ca2þ in the ER. Ca2þ fluxes were determined by assuming that O1 has a different conductance than O2 ; and thus the flux was proportional to O1 þ 10O2 : Ca2þ flowing out of the ER was assumed to enter a local pool (with concentration denoted by c), flowing thence into the cytoplasm. It is the local concentration that affects the transition rates in the IPR, not the cytoplasmic concentration, and thus the model allows for the local build-up of Ca2þ around the mouth of the channel.
3. Comparing models 3.1. The data set and the posterior distribution Data were obtained by superfusion of 45 Ca2þ -loaded hepatic microsomes (Fig. 5). At t ¼ 0 s; IP3 and Ca2þ were added to the stimulus buffer at the indicated strengths. A combination of fast superfusion flow and high Ca2þ buffering meant that ½Ca2þ and ½IP3 were effectively held constant in the stimulus buffer, while Ca2þ flux through the IPR was measured as a function of time. Complete details can be found in Dufour et al. (1997). As can be clearly seen from the data in Fig. 5, if either ½Ca2þ or ½IP3 is stepped up and held fixed, the flux through the IPR first increases and then decreases, i.e., the IPR adapts to the maintained concentration. It has been argued that this adaptation is not due to depletion of microsomal Ca2þ ; but intrinsic to the receptor kinetics, a conclusion that is reinforced by similar results from permeabilised cells
ARTICLE IN PRESS 128
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140 I ck18
ck19 k-19
k-18 k1
T
R k-1
4 k5 p
4k9p
k -2
k 10
R1
k-6
T1
k -10
3 k5 p
3 k9 p
2 k-2
k11
2 k -6
T2
R2 k-11 2 k5 p
2k9p
3 k-2
3 k-6
k12 T3
R3
k5 p
k-12 k9 p
4 k-2
k 13 O1
4 k-6
C1
k-13 ck15
k-15
O2
Fig. 4. Transition diagram of the Dawson–Lea–Irvine model.
ARTICLE IN PRESS
Relative
45
Ca
2+
release
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140 0.30
129
[IP3] (µM)
[Ca ] = 0.4 µM 2+
0.3 1 3 5 10 30 50
0.25 0.20 0.15 0.10 0.05 0.00 0.0
0.5
1.0
1.5
2.0
0.25 0.20
Ca
0.10
45
0.15
Relative
2+
release
2+
[Ca ] (µM) 0.02 0.1 0.2 0.4 1 3 10 100
[IP3] = 10 µM
0.05 0.00 0.0
0.5
1.0 time (s)
1.5
2.0
Fig. 5. Data obtained by superfusion of 45 Ca2þ -loaded hepatic microsomes. Upper panel: at time t ¼ 0; 0:4 mM Ca2þ and the indicated IP3 concentrations were added to the stimulus buffer. Lower panel: at time t ¼ 0; 10 mM IP3 and the indicated Ca2þ concentrations were added to the stimulus buffer. Complete details of the experimental method are in Dufour et al. (1997).
(Swatton and Taylor, 2002). However, recent preliminary results (Falcke, unpublished work) raise the interesting possibility that depletion may be more important than previously thought. For now, we shall ignore this complication. From Fig. 5 it can be seen that the time courses were collected every 72 ms; a resolution too low to resolve the response peak for each time course. However, data collected at the higher resolution of 9 ms from type 2 IPR in permeabilised hepatocytes (Marchant et al., 1997; Adkins and Taylor, 1999; Swatton and Taylor, 2002) or from type 3 IPR in RINm5F cells (Swatton and Taylor, 2002) show that the times to peak are generally greater than 100 ms: Thus, the low time resolution of the data used here is unlikely to introduce severe systematic errors in the parameter estimation. Of course, fits to higher resolution data are preferable, but this is left for future work. For each model we constructed the posterior probability distribution of the parameters, f ; given the data, d; i.e., Prðf jdÞ; by using a Bayesian formalism. Thus Prðf jdÞpPrðdjf ÞPrðf Þ: Given any parameter set, we calculate the probability of the data given the parameters, Prðdjf Þ; by solving the model numerically using those parameters (with the IMSL routine DIVPAG) and
ARTICLE IN PRESS 130
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
assuming that the error at each data point was Gaussian distributed, with variance ð0:02Þ2 : We have no detailed knowledge about the statistical properties of the noise in the data, and thus, in the absence of any information about the error distribution, this was the simplest reasonable assumption. The value for the variance was chosen arbitrarily, but makes no essential difference to the fit, only to the absolute value of the log probability. Hence log Prðdjf Þ ¼
X jd*i di j2 ; 2ð0:02Þ2 i
ð12Þ
where d*i is the computed approximation at the ith data point, di is the ith data point, and where i runs over all the data points. The probability of the parameters, Prðf Þ; is the prior, which is chosen from prior knowledge of reasonable values for the parameters. We chose a non-negative, exponential prior. Thus, Prðf Þ ¼ 0; if any fi is negative, while otherwise Prðf Þ ¼ expðf =fmax Þ; where fmax is a vector of values that sets reasonable upper bounds on the parameter values. Different choices of fmax have little effect on the results; its major function is to ensure that the posterior distribution is normalisable. Once the posterior distribution is constructed, we wish to determine the mean of each parameter value. Alternatively, as is commonly done, we could maximise the posterior distribution, which corresponds to a regularised minimising the sum of squares in (12). However, the maximum of the distribution does not necessarily give the best parameter estimate, as the bulk of feasible parameter values could be nowhere near the maximising value for these skewed distributions. Thus, we choose to characterise the distribution by the parameter means. We determine parameter means and variations by MCMC sampling of the posterior distribution using the Metropolis–Hastings algorithm. For each run we take random initial conditions. It is important to note that, by sampling from the posterior distribution, we not only determine the mean of the distribution (and thus the mean closeness of the fit) but we also determine how much parameter variation is allowed by the data (and thus some idea of the sensitivity of the fit to the parameters). Parameters that converge to a distribution with a small variance are well defined by the fitting procedure, while parameters that converge to a distribution with a large variance have little effect on how well the model fits the data. This is illustrated in Fig. 6. Here, we do not perform an exhaustive analysis of the posterior distributions, merely point out some overall trends. Although, in general, the sampler did not display ergodic properties for the length of runs used, the resulting chains easily determined the feasible parameter ranges. When fitting these three models, we use the rate constants rather than their ratios, i.e., we fit using ki and ki rather than ki =ki : However, it sometimes happens that, in sampling from the posterior distribution, it becomes clear that only a ratio of parameters is determined, not the parameters themselves. This is illustrated in Fig. 7, which uses samples taken from a run using the Sneyd–Dufour model. It is clear that, although the run of samples for k4 and k4 do not converge well, the ratio k4 =k4 converges well and is thus determined by the fit. This indicates that this particular transition in the receptor model could be well approximated by a quasi-steady-state approximation. We have not tested each of the sample runs exhaustively for convergence of every possible parameter combination, but point out, in the appropriate captions, some of the convergent ratios we have found.
ARTICLE IN PRESS J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140 5000
6
4000 k -3
5 k2
131
4
3000 2000
3
1000
2
0 0
0
1000 2000 sample number/10
2000 sample number/10
Fig. 6. A typical run of samples from the posterior distribution generated by fitting the De Young–Keizer model to data, shown after the initial burn-in period. The left panel shows an example of a parameter that is unambiguously determined by the fit, while the right panel shows an example of a parameter that is not well determined by the fit.
2.5
k-4
2.0 1.5 1.0 0.5 0.0 0
2000
0
2000
4000
6000
80
k4
60 40 20 0 4000
6000
-3
k-4 /k4
30x10
20 10 0 0
500
1000
1500
2000
2500
sample number/10 Fig. 7. A typical run of samples from the posterior distribution generated by fitting the Sneyd–Dufour model to data, shown after the initial burn-in period. Although the run of samples for k4 and k4 do not converge well (top two panels), it is clear that variation in one is accompanied by variation in the other. Plotting their ratio (bottom panel) we see that the ratio k4 =k4 converges well and is thus determined by the fit.
ARTICLE IN PRESS 132
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
3.2. Fitting the De Young–Keizer model The first model we fit using this technique is the De Young–Keizer model, and results from a typical run of samples are shown in Fig. 6. Sample runs from two different parameters are included, each showing a different kind of behaviour. For many of the parameters, after an initial burn-in period caused by the randomly chosen starting point, the samples settle down to a distribution with a clear mean and variance, as shown in the sample run of k2 (left panel). We conclude that the mean value of k2 is 3.81, with standard deviation 0.18. However, the sample run for k3 (right panel) shows rather different behaviour, wandering from almost 0 to over 5000 without ever reaching a steady distribution. In addition, for this sample run, log Prðdjf Þ had mean 590, with standard deviation 2.4. We conclude from this that the data are not sufficient to determine k3 unambiguously; large variations in k3 make almost no difference to the posterior probability. Using the values given in Table 1 we plot the responses to Ca2þ or IP3 steps, corresponding to the data shown in Fig. 5. The model responses are shown in Fig. 8. The responses to steps of IP3 (Fig. 8A) agree well with the data, but the responses to steps of Ca2þ (Fig. 8B) become much too fast at higher ½Ca2þ : In fact, the responses become so fast that they are indistinguishable from the axes at this level of resolution. Interestingly, the steady-state open probability curves (Fig. 8C) move to the left with increasing ½IP3 ; the opposite of the behaviour when the original parameters are used. 3.3. Fitting the Sneyd–Dufour model Because of the greater number of parameters in this model, the data were less able to determine them unambiguously. Fewer than half the parameters converged during the MCMC sampling, and wide variations in many parameter values made almost no difference to the fit. A summary of the results is shown in Table 2, and typical solutions are shown in Fig. 9. In one typical sample run log Prðdjf Þ had mean 241, with standard deviation 2.3. From this perspective the Sneyd–Dufour model provides a better fit to the data than does the De Young–Keizer model. In addition, the Table 1 Parameters of the De Young–Keizer model, determined by fitting to kinetic data sd or range k1 ¼ 83:8 (400) k2 ¼ 3:81 (0.2) k3 ¼ 11 103 (400) k4 ¼ 0:15 (0.2) k5 ¼ 53:9 (20) k1 ¼ 114 (52) k2 ¼ 0:409 (0.21) k3 ¼ 2500 (377) k5 ¼ 4:52 (1.64)
sd ¼ 9:8 sd ¼ 0:18 7–16 103 0–0.35 sd ¼ 2:8 sd ¼ 20 sd ¼ 0:051 0–5000 sd ¼ 0:46
Original parameters of De Young and Keizer (1992) are shown in parentheses. If the parameter is determined unambiguously by the fit, the standard deviation is given. Otherwise, the typical range of values is given. k4 was determined by ensuring detailed balance and thus k1 k2 k3 k4 ¼ k1 k2 k3 k4 : All concentrations are in mM:
ARTICLE IN PRESS
relative release
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140 [IP3] = 1 µM 3 5 10 30 50
2+
[Ca ] = 0.4 µM
0.25
133
0.20 0.15 0.10 0.05 0.00 0.5
1.0 1.5 time (s)
(A)
relative release
0.30
2+
[Ca ] = 0.02 µM 0.1 0.2 0.4 1 3
[IP3] = 10 µM
0.25
2.0
0.20 0.15 0.10 0.05 0.00 0.0
0.5
1.0 1.5 time (s)
(B)
[IP3] = 0.02 0.2 1 2
0.010 open probability
2.0
0.008 0.006 0.004 0.002 0.000 0.01
(C)
0.1
1
10
100
2+
[Ca ] (µM)
Fig. 8. Responses of the De Young–Keizer model, using the best-fit parameters shown in Table 1. (A) Responses at the indicated ½Ca2þ ; using the same IP3 concentrations as in Fig. 5, upper panel. (B) Responses at the indicated ½IP3 ; using the same Ca2þ concentrations as in Fig. 5, lower panel. (C) The steady-state open probability curves at four different ½IP3 :
responses to step increases in ½Ca2þ (Fig. 9, middle panel) are not too fast at high ½Ca2þ ; as they are in the De Young–Keizer model. The steady-state open probability curves again shift to the left with increasing ½IP3 : 3.4. Fitting the Dawson–Lea–Irvine model When fitting the Dawson–Lea–Irvine model to the data we find, again, that the parameters are not well determined by the data, with large variations in parameters leading to little change in the
ARTICLE IN PRESS 134
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
Table 2 Parameters of the Sneyd–Dufour model, determined by fitting to kinetic data sd or range k1 ¼ 20:9 k2 ¼ 18:7 k3 ¼ 0:5 k4 ¼ 64 L1 ¼ 0:125 L3 ¼ 0:07 L5 ¼ 2346 k1 ¼ 1:2 k2 ¼ 0:8 k3 ¼ 241 k4 ¼ 1:5 l2 ¼ 0:36 l4 ¼ 4:2 l6 ¼ 50 898
sd ¼ 3:1 sd ¼ 3:1 0–2 24–97 sd ¼ 0:016 0.02–0.11 180–4400 sd ¼ 0:12 0–4 1–500 0.5–2.8 0–0.75 sd ¼ 0:6 10–140 000
If the parameter is determined unambiguously by the fit, the standard deviation is given. Otherwise, the range is given. The remaining parameters are determined by ensuring detailed balance, and thus l2 ¼ l2 k1 =ðk1 L1 Þ; l6 ¼ k4 l6 =ðk4 L5 Þ; l4 ¼ k2 l4 =ðk2 L3 Þ: All concentrations are in mM; time in seconds. Although neither k4 nor k4 were determined by the fit, the ratio k4 =k4 was, with a mean of 0.024 and standard deviation of 0.0024.
posterior (Table 3). However, the ratios k2 =k5 and k6 =k9 were determined unambiguously, leading to the prediction that these transitions could be eliminated by using a quasi-state-state approximation without sacrificing model accuracy. In a sample run log Prðdjf Þ had mean 760, with standard deviation 2.4. Typical solutions, and the steady-state curves, are shown in Fig. 10. For the given parameters, the steady-state curves have a shape quite unlike experimental data.
4. Discussion and conclusions Each of the three models we discuss was constructed with a particular aim in mind. The De Young–Keizer model parameters were chosen so as to reproduce the steady-state data of Bezprozvanny et al. (1991), and also so that receptor activation by Ca2þ would be faster than Ca2þ -induced inactivation. The principal motivation of the Dawson–Lea–Irvine model was to study how adaptation can occur in a receptor model. Only the Sneyd–Dufour model was originally constructed specifically to study the kinds of dynamic responses shown in Fig. 5. It is thus no great surprise that the Sneyd–Dufour model does the best job of fitting this particular data set. Nevertheless, by comparing these models it is clear that the use of simple mass-action binding of Ca2þ to the receptor results in kinetics that are much too fast, and response peaks that occur too early. We can thus conclude that saturation of the rate of Ca2þ binding is likely an important feature of the receptor kinetics. This conclusion holds even though the majority of the data points were in the tails of the responses, with relatively few around the peaks, and thus the
ARTICLE IN PRESS J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140 [IP3] = 1 µM 3 5 10 30 50
0.25 relative release
2+
[Ca ] = 0.4 µM
0.20
135
0.15 0.10 0.05 0.00 0.0
0.5
(A)
2.0
2+
[Ca ] = 0.02 µM 0.1 0.2 0.4 1 3
[IP3] = 10 µM
0.20 relative release
1.0 1.5 time (s)
0.15 0.10 0.05 0.00 0.0
0.5
(B)
1.0 1.5 time (s)
open probability
0.04
2.0
[IP3] = 0.02 0.2 1 2
0.03 0.02 0.01 0.00 0.01
(C)
0.1 1 2+ [Ca ] (µM)
10
100
Fig. 9. (A–C) Responses of the Sneyd–Dufour model, using the best-fit parameters shown in Table 2.
model fits will be biased towards giving agreement to the response tails. The best way to correct this would be to use a data set at a higher time resolution, but this is left for future work. In addition, the steady-state behaviour of the Dawson–Lea–Irvine model was not in good agreement with data, while the steady-state curves of the other two models exhibited the correct bell-shaped curves, but with peaks moving to the left as ½IP3 increased. Although the Sneyd– Dufour model gives the best quantitative fit to the selected data (i.e., with the highest mean likelihood), the Dawson–Lea–Irvine model the worst, it would be dangerous to draw conclusions as to which is the ‘‘best’’ model. For instance, the Sneyd-Dufour model paid little attention to adaptation of the IPR, while, as discussed below, there is uncertainty as to the correct steady-state
ARTICLE IN PRESS 136
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
Table 3 Parameters of the Dawson–Lea–Irvine model determined by the fit sd or range k1 ¼ 0:045 (1) k5 ¼ 1345 (1000) k9 ¼ 70 (100) k10 ¼ 0:033 (1) k12 ¼ 0:85 (1) k15 ¼ 19; 970 (100) k18 ¼ 45 (1) k19 ¼ 22 (10) k1 ¼ 173 (100) k2 ¼ 51; 522 (1000) k6 ¼ 24 (10) k11 ¼ 0:25 (1) k13 ¼ 2859 (10) k15 ¼ 12:25 (10) k18 ¼ 0:053 (0.1)
0.004–0.25 sd ¼ 367 17–115 104 –0.14 0.0001–2.2 sd ¼ 5; 531 sd ¼ 7 0.4–37 2–450 17,000–83,000 7–35 0.0003–0.7 sd ¼ 737 sd ¼ 1:14 106 0:12
Original parameters of Dawson et al. (2003) are shown in parentheses. Parameters that do not appear in this list were determined by detailed balance. For instance, k10 ¼ k10 k6 k1 k2 =ðk2 k1 k6 Þ: Although neither k5 nor k2 were determined by the fit, the ratio k2 =k5 was, with a mean of 33.3 and standard deviation of 2. The ratio k6 =k9 was also determined by the fit, with mean 0.5 and standard deviation 0.1.
data that should be used to constrain the fit. Nevertheless, it seems safe to conclude that the De Young–Keizer model exhibits responses that are much too fast at high ½Ca2þ ; and is thus not a good model in that regime, while the steady-state properties of Dawson–Lea–Irvine model need improvement. The best-fit parameters for the De Young–Keizer model have Ca2þ activation (k5 ) an order of magnitude faster than Ca2þ inactivation (k2 and k4 ), while binding of IP3 (k1 and k3 ) is faster still. Similarly, the best-fit parameters for the Sneyd–Dufour model have Ca2þ activation (governed principally by k4 and l6 ) much faster than Ca2þ inactivation (governed principally by k1 and l2 ), while IP3 binding (k2 ) is relatively slow. Thus, the major differences between the De Young– Keizer and the Sneyd–Dufour models are, firstly, the rate of response at high ½Ca2þ ; and, secondly, the relative rate of IP3 binding. In fact, contrary to a procedure that is often used, the Sneyd–Dufour model suggests that the binding of IP3 cannot be assumed to be fast enough to be eliminated from the equations. The best-fit parameters for the Dawson–Lea–Irvine model again have Ca2þ activation (k15 ) faster than Ca2þ inactivation (k18 ). Although binding of the first three IP3 molecules proceeds very quickly, binding of the fourth (k5 ) proceeds somewhat slower than Ca2þ activation. Thus, the Dawson–Lea–Irvine model also suggests that the kinetics of IP3 binding cannot be ignored in IPR models. When the best-fit parameters are used in the De Young–Keizer or Sneyd–Dufour model, the peak of the bell-shaped curve moves to the left with increasing ½IP3 : Since this is the opposite of the behaviour of the De Young–Keizer model when the original parameters are used, it emphasises the fact that this feature, being parameter-dependent, cannot be used to distinguish between models. In addition, it raises the important question of which steady-state data should be
ARTICLE IN PRESS J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140 [IP3] = 1 µM 3 5 10 30 50
relative release
2+
[Ca ] = 0.4 µM
0.30
137
0.20 0.10 0.00 0.0
0.5
(A)
relative release
0.25
1.0 1.5 time (s)
2.0
2+
[Ca ] = 0.02 µM 0.1 0.2 0.4 1 3
[IP3] = 10 µM
0.20 0.15 0.10 0.05 0.00 0.0
0.5
(B)
2.0
[IP3] = 0.2 2 10 20
0.06 open probability
1.0 1.5 time (s)
0.05 0.04 0.03 0.02 0.01 0.00 0.01
(C)
0.1 1 2+ [Ca ] (µM)
10
100
Fig. 10. (A–C) Responses of the Dawson–Lea–Irvine model, using the best-fit parameters shown in Table 3.
used to constrain the fitting procedure further. At present, there is no easy answer to this question. In each of the models discussed here, the steady-state open probabilities are low, less than 0.08. Although this agrees with the experimental data of Kaftan et al. (1997), Hagar et al. (1998) and Hagar and Ehrlich (2000), as well as being consistent with the time-dependent measurements of Finch et al. (1991), Swatton and Taylor (2002), and Marchant and Taylor (1997), it is in direct conflict with the data of Mak et al. (1998, 2000, 2001) and Ramos-Franco et al. (2000). Unfortunately, there is no agreement as to which set of data better represents in vivo behaviour. At this stage, the only thing modellers can do is wait until experimentalists manage to resolve these differences.
ARTICLE IN PRESS 138
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
Another obvious conclusion from this exercise in model comparison is that time-dependent responses to steps of ½Ca2þ and ½IP3 are insufficient to determine the model parameters unambiguously. In each of the models, large variations in model parameters make very little difference to the quality of the fit, resulting in a lack of convergence of the sampling procedure. The only way to improve this is to use additional data in the fit, or otherwise to constrain some of the parameters. Some obvious experimental data that have not been used in our fitting method are the time course of recovery from inactivation (for instance, see Fig. 3A of Swatton and Taylor, 2002), steady-state data (once agreement can be reached on what the correct data are), responses to steps of ½Ca2þ (½IP3 ) at a wider variety of ½IP3 (½Ca2þ ), and responses to double pulses of Ca2þ or IP3 : Ideally, the next generation of models will provide a good fit to this entire range of steady-state and kinetic data, exhibit the required degree of adaptation, incorporate possible vesicle or ER depletion, take into account the possible existence of microdomains around the channel mouth, and provide different parameter sets or model structures for each of the receptor sub-types. Construction of such a family of models, and fitting them to this much wider array of data will be a challenging problem.
Acknowledgements J. Sneyd was supported by the Marsden Fund of the Royal Society of New Zealand, and by the Mathematical Biosciences Institute of The Ohio State University. J.-F. Dufour was supported by Swiss National Found, Grant 3100A0-100513.
References Adkins, C.E., Taylor, C.W., 1999. Lateral inhibition of inositol 1,4,5-trisphosphate receptors by cytosolic Ca2þ : Curr. Biol. 9, 1115–1118. Atri, A., Amundson, J., Clapham, D., Sneyd, J., 1993. A single-pool model for intracellular calcium oscillations and waves in the Xenopus laevis oocyte. Biophys. J. 65, 1727–1739. Berridge, M.J., 1997. Elementary and global aspects of calcium signalling. J. Physiol. 499, 291–306. Bezprozvanny, I., 1994. Theoretical analysis of calcium wave propagation based on inositol (1,4,5)-trisphosphate ðInsP3 Þ receptor functional properties. Cell Calcium 16, 151–166. Bezprozvanny, I., Watras, J., Ehrlich, B.E., 1991. Bell-shaped calcium-response curves of Ins(1,4,5)P3 - and calciumgated channels from endoplasmic reticulum of cerebellum. Nature 351, 751–754. Cardy, T.J., Traynor, D., Taylor, C.W., 1997. Differential regulation of types-1 and -3 inositol trisphosphate receptors by cytosolic Ca2þ : Biochem. J. 328, 785–793. Clapham, D., 1995. Calcium signaling. Cell 80, 259–268. Combettes, L., Claret, M., Champeil, P., 1993. Calcium control on InsP3 -induced discharge of calcium from permeabilised hepatocyte pools. Cell Calcium 14, 279–292. Dawson, A.P., Lea, E.J., Irvine, R.F., 2003. Kinetic model of the inositol trisphosphate receptor that shows both steady-state and quantal patterns of Ca2þ release from intracellular stores. Biochem. J. 370, 621–629. De Young, G.W., Keizer, J., 1992. A single pool IP3 -receptor based model for agonist stimulated Ca2þ oscillations. Proc. Natl. Acad. Sci. USA 89, 9895–9899. Dufour, J.-F., Arias, I.M., Turner, T.J., 1997. Inositol 1,4,5-trisphosphate and calcium regulate the calcium channel function of the hepatic inositol 1,4,5-trisphosphate receptor. J. Biol. Chem. 272, 2675–2681.
ARTICLE IN PRESS J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
139
Finch, E.A., Turner, T.J., Goldin, S.M., 1991. Calcium as a coagonist of inositol 1,4,5-trisphosphate-induced calcium release. Science 252, 443–446. Gilks, W.R., Richardson, S., Spiegelhalter, D.J., 1995. Markov Chain Monte Carlo in Practice. CRC Press, Boca Raton, FL. Hagar, R.E., Ehrlich, B.E., 2000. Regulation of the type III InsP3 receptor by InsP3 and ATP. Biophys. J. 79, 271–278. Hagar, R.E., Burgstahler, A.D., Nathanson, M.H., Ehrlich, B.E., 1998. Type III InsP3 receptor channel stays open in the presence of increased calcium. Nature 396, 81–84. ! Hajnoczky, G., Thomas, A.P., 1997. Minimal requirements for calcium oscillations driven by the IP3 receptor. EMBO J. 16, 3533–3543. Hodgkin, A.L., Huxley, A.F., 1952. A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544. Iino, M., 1990. Biphasic Ca2þ dependence of inositol 1,4,5-trisphosphate-induced Ca2þ release in smooth muscle cells of the Guinea Pig Taenia Caeci. J. Gen. Physiol. 95, 1103–1122. Iino, M., Endo, M., 1992. Calcium-dependent immediate feedback control of inositol 1,4,5-trisphosphate-induced Ca2þ release. Nature 360, 76–78. Iino, M., Tsukioka, M., 1994. Feedback control of inositol trisphosphate signalling by calcium. Mol. Cell. Endocrinol. 98, 141–146. Ilyin, V., Parker, I., 1994. Role of cytosolic Ca2þ in inhibition of InsP3 -evoked Ca2þ release in Xenopus oocytes. J. Physiol. 477, 503–509. Kaftan, E.J., Ehrlich, B.E., Watras, J., 1997. Inositol 1,4,5-trisphosphate (InsP3 ) and calcium interact to increase the dynamic range of InsP3 receptor-dependent calcium signaling. J. Gen. Physiol. 110, 529–538. LeBeau, A.P., Yule, D.I., Groblewski, G.E., Sneyd, J., 1999. Agonist-dependent phosphorylation of the inositol 1,4,5trisphosphate receptor: a possible mechanism for agonist-specific calcium oscillations in pancreatic acinar cells. J. Gen. Physiol. 113, 851–871. Mak, D.-O.D., McBride, S., Foskett, J.K., 1998. Inositol 1,4,5-trisphosphate activation of inositol trisphosphate receptor Ca2þ channel by ligand tuning of Ca2þ inhibition. Proc. Natl. Acad. Sci. USA 95, 15821–15825. Mak, D.O., McBride, S., Raghuram, V., Yue, Y., Joseph, S.K., Foskett, J.K., 2000. Single-channel properties in endoplasmic reticulum membrane of recombinant type 3 inositol trisphosphate receptor. J. Gen. Physiol. 115, 241–256. Mak, D.-O.D., McBride, S., Foskett, J.K., 2001. Regulation by Ca2þ and inositol 1,4,5-trisphosphate (InsP3 ) of single recombinant type 3 InsP3 receptor channels: Ca2þ activation uniquely distinguishes types 1 and 3 InsP3 receptors. J. Gen. Physiol. 117, 435–446. Marchant, J.S., Taylor, C.W., 1997. Cooperative activation of IP3 receptors by sequential binding of IP3 and Ca2þ safeguards against spontaneous activity. Curr. Biol. 7, 510–518. Marchant, J.S., Taylor, C.W., 1998. Rapid activation and partial inactivation of inositol trisphosphate receptors by inositol trisphosphate. Biochemistry 37, 11524–11533. Marchant, J.S., Chang, Y.-T., Chung, S.-K., Irvine, R.F., Taylor, C.W., 1997. Rapid kinetic measurements of 45 Ca2þ mobilization reveal that Ins(2,4,5)P3 is a partial agonist at hepatic InsP3 receptors. Biochem. J. 321, 573–576. Monod, J., Wyman, J., Changeux, J.P., 1965. On the nature of allosteric transition: a plausible model. J. Mol. Biol. 12, 88–118. Othmer, H., Tang, Y., 1993. In: Othmer, H., Murry, J., Maini, P. (Eds.), Experimental and Theoretical Advances in Biological Pattern Formation. Plenum Press, London. Parker, I., Ivorra, I., 1990. Inhibition by Ca2þ of inositol trisphosphate-mediated Ca2þ liberation: a possible mechanism for oscillatory release of Ca2þ : Proc. Natl. Acad. Sci. USA 87, 260–264. Parker, I., Yao, Y., 1991. Regenerative release of calcium from functionally discrete subcellular stores by inositol trisphosphate. Proc. R. Soc. London B 246, 269–274. Parker, I., Yao, Y., Ilyin, V., 1996. Fast kinetics of calcium liberation induced in Xenopus oocytes by photoreleased inositol trisphosphate. Biophys. J. 70, 222–237. Parys, J.B., Sernett, S.W., DeLisle, S., Snyder, P.M., Welsh, M.J., Campbell, K.P., 1992. Isolation, characterization, and localization of the inositol 1,4,5-trisphosphate receptor protein in Xenopus laevis oocytes. J. Biol. Chem. 267, 18776–18782.
ARTICLE IN PRESS 140
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
Ramos-Franco, J., Bare, D., Caenepeel, S., Nani, A., Fill, M., Mignery, G., 2000. Single-channel function of recombinant type 2 inositol 1,4,5-trisphosphate receptor. Biophys. J. 79, 1388–1399. Sanderson, M.J., Charles, A.C., Boitano, S., Dirksen, E.R., 1994. Mechanisms and function of intercellular calcium signaling. Mol. Cell. Endocrinol. 98, 173–187. Sienaert, I., Missiaen, L., De Smedt, H., Parys, J.B., Sipma, H., Casteels, R., 1997. Molecular and functional evidence for multiple Ca2þ -binding domains in the type I inositol 1,4,5-trisphosphate receptor. J. Biol. Chem. 272, 25899–25906. Sneyd, J., Dufour, J.F., 2002. A dynamic model of the type-2 inositol trisphosphate receptor. Proc. Natl. Acad. Sci. USA 99, 2398–2403. Sneyd, J., Keizer, J., Sanderson, M.J., 1995. Mechanisms of calcium oscillations and waves: a quantitative analysis. FASEB J. 9, 1463–1472. Swatton, J.E., Taylor, C.W., 2002. Fast biphasic regulation of type 3 inositol trisphosphate receptors by cytosolic calcium. J. Biol. Chem. 277, 17571–17579. Taylor, C.W., 1998. Inositol trisphosphate receptors: Ca2þ -modulated intracellular Ca2þ channels. Biochim. Biophys. Acta 1436, 19–33. Taylor, C.W., Laude, A.J., 2002. IP3 receptors and their regulation by calmodulin and cytosolic Ca2þ : Cell Calcium 32, 321–334. ! Thomas, A.P., Bird, G.S.J., Hajnoczky, G., Robb-Gaspers, L.D., Putney, J.W.J., 1996. Spatial and temporal aspects of cellular calcium signaling. FASEB J. 10, 1505–1517. Yao, Y., Parker, I., 1992. Potentiation of inositol trisphosphate-induced Ca2þ mobilization in Xenopus oocytes by cytosolic Ca2þ : J. Physiol. 458, 319–338. Yoneshima, H., Miyawaki, A., Michikawa, T., Furuichi, T., Mikoshiba, K., 1997. Ca2þ differentially regulates the ligand-affinity states of type 1 and type 3 inositol 1,4,5-trisphosphate receptors. Biochem. J. 322, 591–596.
Progress in Biophysics & Molecular Biology 85 (2004) 121–140
A comparison of three models of the inositol trisphosphate receptor J. Sneyda,*, M. Falckeb, J.-F. Dufourc, C. Foxa a
Department of Mathematics, University of Auckland, Auckland, New Zealand b Hahn Meitner Institute, SF5, Glienicker Str. 100, 14109 Berlin, Germany c University of Bern, Switzerland
Abstract The inositol (1,4,5)-trisphosphate receptor (IPR) plays a crucial role in calcium dynamics in a wide range of cell types, and is often a central feature in quantitative models of calcium oscillations and waves. We compare three mathematical models of the IPR, fitting each of them to the same data set to determine ranges for the parameter values. Each of the fits indicates that fast activation of the receptor, followed by slow inactivation, is an important feature of the model, and also that the speed of inositol trisphosphate ðIP3 Þ binding cannot necessarily be assumed to be faster than Ca2þ activation. In addition, the model which assumed saturating binding rates of Ca2þ to the IPR demonstrated the best fit. However, lack of convergence in the fitting procedure indicates that responses to step increases of ½Ca2þ and ½IP3 provide insufficient data to determine the parameters unambiguously in any of the models. r 2004 Elsevier Ltd. All rights reserved. Keywords: Calcium oscillations; Mathematical model; Markov chain Monte Carlo; Data fitting
1. Introduction In almost every cell type, the cytosolic concentration of free calcium ions ð½Ca2þ Þ exhibits complex dynamic behaviour, including oscillations and periodic waves (Sanderson et al., 1994; Clapham, 1995; Thomas et al., 1996; Berridge, 1997). Often, the physiological significance of such phenomena is not completely understood, but there is general acceptance that these oscillations and waves are one way in which cells can use Ca2þ as an intracellular signalling mechanism, while avoiding the toxic effects of prolonged raised ½Ca2þ : In electrically non-excitable cells, Ca2þ release is often mostly from internal stores, although influx from the outside can play an important modulatory role. Agonist stimulation results in the formation of inositol (1,4,5)-trisphosphate *Corresponding author. Fax: +64-9-3737-457. E-mail address: [email protected] (J. Sneyd). 0079-6107/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.pbiomolbio.2004.01.013
ARTICLE IN PRESS 122
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
ðIP3 Þ which diffuses through the cytoplasm and binds to the IP3 receptors (IPR) situated on the endoplasmic reticulum (ER). Subsequent release of Ca2þ from the ER can further modulate the open probability of the IPR, with cytosolic (and possibly lumenal) Ca2þ exerting both a positive and negative feedback effect on the IPR and on the formation of IP3 : It is these dynamic interactions, when coupled with Ca2þ removal processes and with other Ca2þ release pathways through, for instance, ryanodine receptors, that can result in oscillations and waves of ½Ca2þ : Thus, an understanding of IP3 receptors is central to a detailed understanding of Ca2þ dynamics in many cell types, and a number of IPR models have been constructed. The obvious questions to ask are, firstly, in what way have these models helped us understand the physiology or biophysics of the IPR, and, secondly, which is the best model? The first question is, we think, relatively easily answered. The models have contributed a great deal to our understanding of how the binding properties of the IPR underlie steady-state responses and adaptation, and have provided quantitative frameworks that have aided in the design of experiments. Not all of the useful models have been mathematically expressed. The qualitative, or diagrammatic, models proposed by Taylor (1998) (Adkins and Taylor, 1999; ! Taylor and Laude, 2002) and Hajnoczky and Thomas (1997) have led to a great deal of experimental work and been crucial in the development of our understanding of the IPR. Many of the ideas in those models have been incorporated into later mathematical models, particularly that of Sneyd and Dufour (2002). The second question is more difficult to answer. Indeed, it can only be satisfactorily answered by pointing out that (as has been said many times) all models are wrong, but some models are useful. Thus, the search for the ‘‘best’’ model is a futile one. Each IPR subtype behaves differently from the others, and its behaviour is modified by the environment of the IPR, making different models for different situations a necessity. Nevertheless, despite these complications, if we restrict ourselves to a less ambitious goal, models can be usefully compared. If we take a single data set, we can investigate how well a variety of models fit those data, and thus compare models with respect to the particular features of most importance in the data. Because of the effort involved, such detailed model comparisons are rarely done, but this does not make them any less important. We shall take a single data set (Dufour et al., 1997), and use it to determine the parameters in three models: the models of De Young and Keizer (1992), Sneyd and Dufour (2002) and Dawson et al. (2003). To determine the parameters, we performed parameter estimation within a Bayesian formalism computed via Markov chain Monte Carlo (MCMC) sampling from the posterior distribution for parameters. Examining the posterior not only gives ‘‘best fit’’ parameter values but also a quantitative idea of how well parameters are determined and how good a representation of the data the model is. The MCMC sampler was implemented using Metropolis–Hastings dynamics (Gilks et al., 1995).
2. The three models 2.1. The De Young–Keizer model In the early 1990s, two crucial features of the IPR became apparent. Firstly, it was shown (Iino, 1990; Parker and Ivorra, 1990; Bezprozvanny et al., 1991; Parys et al., 1992) that, at steady state,
ARTICLE IN PRESS J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
123
the open probability of the IPR is a bell-shaped function of ½Ca2þ ; with the maximum open probability occurring at about ½Ca2þ ¼ 300–500 nM: Secondly, it became clear that, in response to an increased ½Ca2þ ; the IPR responded in an adaptive manner, first activating, and subsequently inactivating (Parker and Yao, 1991; Finch et al., 1991; Iino and Endo, 1992; Yao and Parker, 1992; Iino and Tsukioka, 1994; Ilyin and Parker, 1994). These results were most easily interpreted by assuming that Ca2þ activates the IPR on a fast time scale, but then inactivates the IPR over a slower time scale. The first model of the IPR to take this into account was due to De Young and Keizer (1992), closely followed by the phenomenological model of Atri et al. (1993), and the sequential binding models of Othmer and Tang (1993) and of Bezprozvanny (1994). Here, we focus on the De Young–Keizer model, as, of these models, it has been the most influential. De Young and Keizer (1992) assumed that the IPR consists of three equivalent and independent subunits, all of which must be in a conducting state before the receptor allows Ca2þ flux. Each subunit has an IP3 binding site, an activating Ca2þ binding site, and an inactivating Ca2þ binding site, each of which can be either occupied or unoccupied, and thus each subunit can be in one of eight states. Each state of the subunit is labelled Sijk ; where i; j and k are equal to 0 or 1, with a 0 indicating the binding site is unoccupied and a 1 indicating it is occupied. The first index refers to the IP3 binding site, the second to the Ca2þ activation site, and the third to the Ca2þ inactivation site. This is illustrated in Fig. 1. Although a fully general model would include 24 rate constants, De Young and Keizer made a number of simplifying assumptions. The rate constants are assumed to be independent of whether or not activating Ca2þ is bound, and the rate of Ca2þ activation is assumed to be independent of whether or not IP3 is bound. Finally, Ca2þ inactivation is assumed to be independent of Ca2þ activation and of IP3 : This leaves only 10 rate constants, k71 ; y; k75 : The fraction of subunits in the state Sijk is denoted by xijk : The differential equations for these are based on mass action kinetics, and thus, for example, dx000 ¼ ðk1 x100 k1 px000 Þ þ ðk4 x001 k4 cx000 Þ þ ðk5 x010 k5 cx000 Þ; dt
ð1Þ
where p denotes ½IP3 and c denotes ½Ca2þ : The model assumes that the IP3 receptor passes Ca2þ current only when three subunits are in the state S110 (i.e., with one IP3 and one activating Ca2þ bound), and thus the open probability of the receptor is x3110 : This gives seven differential equations for the receptor states. k1p S110
S010
k2c
k-5 k-4
k-4
S011 k-3
k4c
k-5 k-2
k-2
k2c
S111
S101 k-3
k2c
k5c
k3p S001
k-5
S110
k-5
k5c
k3p S111
k4c
k5c S100
S000
k-1 k-2
k1p
k1c
k-5
Fig. 1. The binding diagram for the De Young–Keizer IPR model. c denotes ½Ca2þ ; and p denotes ½IP3 :
ARTICLE IN PRESS 124
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
Parameters of the model were chosen to obtain agreement with the steady-state data of Bezprozvanny et al. (1991). Nevertheless, the bell-shaped nature of the steady-state IPR open probability curve has very little, if anything, to do with the fact that the IPR can mediate oscillatory activity. As was originally pointed out by Keizer (see Sneyd et al., 1995, Fig. 2), during an oscillation the steady-state response is less important than the response to a changing ½Ca2þ : After all, the cell never ‘‘sees’’ a steady ½Ca2þ ; unless the time scale of the oscillation is very much slower than the time scale of the IPR kinetics. It is thus the fast activation, followed by the slow inactivation (k5 > k2 ) that is the crucial feature of the De Young–Keizer model and this separation of time scales is independent of the shape of the steady-state open probability curve. IPR models with monotonic steady-state curves support Ca2þ oscillations just as readily as IPR models with a bell-shaped curve (LeBeau et al., 1999). 2.2. The Sneyd–Dufour model Finch et al. (1991) were the first to develop a superfusion system in which the IPRs, isolated in microsomal vesicles, could be exposed to solutions in which the ½Ca2þ and ½IP3 were tightly controlled by a combination of fast flow and extensive buffering. This procedure gives a calcium clamp, similar to the voltage clamp used by Hodgkin and Huxley (1952). Since the work of Finch et al., other groups have used a superfusion calcium clamp (Dufour et al., 1997; Marchant et al., 1997; Marchant and Taylor, 1997, 1998; Adkins and Taylor, 1999) and now recordings with a time resolution of 9 ms are possible (Marchant and Taylor, 1997; Swatton and Taylor, 2002). Similar data were also obtained by Parker et al. (1996), although, because their results were obtained in vivo in Xenopus oocytes, they were not able to maintain a precise calcium clamp. Time-dependent data were used by Sneyd and Dufour (2002) in a model that incorporated several ideas suggested by experiments: sequential binding of IP3 and activating Ca2þ ; modulation of IP3 binding by Ca2þ ; a fit of the model to dynamical data instead of steady-state data only, time-dependent inactivation upon IP3 binding, and saturating binding rates of Ca2þ : The model scheme is shown in Fig. 2. The bare subunit, state R; can bind IP3 and move to state O: The binding rate depends on Ca2þ allowing for modulation of IP3 binding in a facilitating or inhibiting way depending on receptor type (Yoneshima et al., 1997; Cardy et al., 1997; Sienaert et al., 1997). Binding of Ca2þ to the activating site corresponds to a transition to state A: Since A can be reached from O only and not from R; IP3 and activating Ca2þ need to bind sequentially. Alternatively, the subunit could close to state S; corresponding to time-dependent inactivation upon IP3 binding. There are two inhibited states I1 and I2 : The state I1 has Ca2þ bound at the inhibiting site only whereas I2 also has Ca2þ bound at the activating site and IP3 : The receptor is open when four subunits are in state O or A or any combination of these states. However, the receptors in which some of the subunits are in the O state have a lower conductance. Hence, the IPR flux was chosen to be proportional to ðoO þ aAÞ4 with o in the range of 0.05–0.10 and a about 0.9. This accounts for the non-zero fluxes observed with only IP3 bound (Finch et al., 1991; Marchant and Taylor, 1997; Combettes et al., 1993; Dufour et al., 1997), and predicts that multiple conductance levels would be observed. Alternatively, to avoid having multiple conductance levels, it could be assumed that each state of the subunit flickers to and from the open configuration but that some states just have a higher probability of being in the open configuration.
ARTICLE IN PRESS J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
125
1(c) R
k-1+l-2
I1
-2(c)
2(c) p
4(c) O
-4(c)
5(c) A
I2 k-1+l-2
k-3
3(c) S
Fig. 2. Transition diagram of the receptor subunit in the Sneyd–Dufour model.
The functions in the model were derived by simplification of a more complex model, the full details of which are in Sneyd and Dufour (2002). This gives ðk1 L1 þ l2 Þc ; ð2Þ f1 ðcÞ ¼ L1 þ cð1 þ L1 =L3 Þ f2 ðcÞ ¼
k2 L3 þ l4 c ; L3 þ cð1 þ L3 =L1 Þ
f2 ðcÞ ¼
k2 þ l4 c ; 1 þ c=L5
ð3Þ ð4Þ
f3 ðcÞ ¼
k3 L5 ; c þ L5
ð5Þ
f4 ðcÞ ¼
ðk4 L5 þ l6 Þc ; L5 þ c
ð6Þ
f4 ðcÞ ¼ f1 ðcÞ ¼
L1 ðk4 þ l6 Þ ; L1 þ c
ðk1 L1 þ l2 Þc : L1 þ c
ð7Þ ð8Þ
One principal feature of the Sneyd–Dufour model is that the transition rates saturate as c increases. This allows for a limited range of activation speeds over a wide range of Ca2þ concentrations, a feature that is absent from the other two models. Saturating rates are obtained by a process essentially identical to the Michaelis–Menten model of enzyme kinetics (Sneyd and Dufour, 2002). The method is illustrated in Fig. 3. If we assume that A* and A% are in instantaneous
ARTICLE IN PRESS 126
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
k1c
~ A
|
k -1 l -1 l1c
l2
l -2
_ A
(c)
|
A k-1 + l -2
Fig. 3. Derivation of a saturating binding rate. If A% and A* are assumed to be in instantaneous equilibrium, the upper transition diagram reduces to the lower, as described in the text.
equilibrium, we get % cA* ¼ L1 A; from which it follows that L1 ¼ l1 =l1 : Thus, letting A ¼ A% þ A* we get dA ¼ ðk1 þ l2 ÞI fðcÞA; dt where cðk1 L1 þ l2 Þ : fðcÞ ¼ c þ L1
ð9Þ
ð10Þ
ð11Þ
Note that, in this model, the IPR can inactivate either in the absence of Ca2þ ; by going to state S; or in the absence of IP3 ; by going to state I1 ; both of which features are required by experimental data (Taylor, 1998). Furthermore, once IP3 is bound (and the receptor is in state O), there is an intrinsic competition between states A and S: If Ca2þ is abundant, the receptor will go preferentially to the activated state A; leading to a higher conductance (or higher open probability, depending on interpretation), but in conditions of low ½Ca2þ the receptor will move instead to the S state, which does not allow Ca2þ current. In a number of important respects, the Sneyd–Dufour model is a quantitative realisation of the qualitative model of Taylor (1998) (Adkins and Taylor, 1999; Marchant and Taylor, 1997, 1998; Taylor and Laude, 2002). However, one major difference is that the Sneyd–Dufour model does not assume that binding of IP3 protects the receptor from inactivation by Ca2þ (i.e., that there is no I2 state), as is proposed by Taylor for type 2 IPR (but not type 3). This proposal was not adopted by Sneyd and Dufour because with it they were unable to obtain acceptable fits to experimental data (unpublished results). One should not put undue weight on this negative result as it is possible that excellent fits do exist, just they have not been found yet. Another possibility is that ER depletion is playing a more important role than initially thought and thus one cannot expect IPR dynamics alone to fit the data well. To tell
ARTICLE IN PRESS J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
127
for sure if the I2 state is necessary will require detailed fits of the alternative models to data collected at the highest possible time resolution, and from types 2 and 3 IPR, as well as a detailed consideration of the possible effects of ER depletion and unstirred boundary layers close to the IPR. 2.3. The Dawson–Lea–Irvine model A different type of model, based upon the allosteric model of Monod et al. (1965), has been constructed by Dawson et al. (2003), with the goal of understanding adaptation and incremental responses by the IPR. They assumed that the IPR could exist in two conformations, R and T; R binding IP3 rapidly, but with low affinity, and T binding IP3 slowly, but with high affinity. Both R and T states can bind four IP3 molecules, and the only open state of the receptor is when the R conformation has four IP3 molecules bound. Finally, it is assumed that Ca2þ on the cytosolic side of the IPR can inactivate the receptor by binding to the R and T states, or activate the receptor by converting the open state, O1 ; to another open state, O2 ; that cannot be shifted across into the closed state. This leads to the state diagram shown in Fig. 4. The O1 to O2 transition gives positive feedback of Ca2þ on the IPR, while the R (or T) to I transition lets the receptor inactivate in the absence of IP3 but in the presence of Ca2þ ; as observed experimentally. The balance between the two conformational states provides fast inactivation and slower inactivation. Upon initial addition of IP3 ; the receptor is quickly pushed through to the O1 state (and, in the positive feedback process, to the O2 state). However, over a longer time scale, the receptor shifts across to the T conformation, pulling the receptor out of the O1 state, and leading to inactivation. To model Ca2þ fluxes, Dawson et al. connected their IPR model to a simplified model of Ca2þ handling, in which they kept track of ½Ca2þ in the ER. Ca2þ fluxes were determined by assuming that O1 has a different conductance than O2 ; and thus the flux was proportional to O1 þ 10O2 : Ca2þ flowing out of the ER was assumed to enter a local pool (with concentration denoted by c), flowing thence into the cytoplasm. It is the local concentration that affects the transition rates in the IPR, not the cytoplasmic concentration, and thus the model allows for the local build-up of Ca2þ around the mouth of the channel.
3. Comparing models 3.1. The data set and the posterior distribution Data were obtained by superfusion of 45 Ca2þ -loaded hepatic microsomes (Fig. 5). At t ¼ 0 s; IP3 and Ca2þ were added to the stimulus buffer at the indicated strengths. A combination of fast superfusion flow and high Ca2þ buffering meant that ½Ca2þ and ½IP3 were effectively held constant in the stimulus buffer, while Ca2þ flux through the IPR was measured as a function of time. Complete details can be found in Dufour et al. (1997). As can be clearly seen from the data in Fig. 5, if either ½Ca2þ or ½IP3 is stepped up and held fixed, the flux through the IPR first increases and then decreases, i.e., the IPR adapts to the maintained concentration. It has been argued that this adaptation is not due to depletion of microsomal Ca2þ ; but intrinsic to the receptor kinetics, a conclusion that is reinforced by similar results from permeabilised cells
ARTICLE IN PRESS 128
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140 I ck18
ck19 k-19
k-18 k1
T
R k-1
4 k5 p
4k9p
k -2
k 10
R1
k-6
T1
k -10
3 k5 p
3 k9 p
2 k-2
k11
2 k -6
T2
R2 k-11 2 k5 p
2k9p
3 k-2
3 k-6
k12 T3
R3
k5 p
k-12 k9 p
4 k-2
k 13 O1
4 k-6
C1
k-13 ck15
k-15
O2
Fig. 4. Transition diagram of the Dawson–Lea–Irvine model.
ARTICLE IN PRESS
Relative
45
Ca
2+
release
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140 0.30
129
[IP3] (µM)
[Ca ] = 0.4 µM 2+
0.3 1 3 5 10 30 50
0.25 0.20 0.15 0.10 0.05 0.00 0.0
0.5
1.0
1.5
2.0
0.25 0.20
Ca
0.10
45
0.15
Relative
2+
release
2+
[Ca ] (µM) 0.02 0.1 0.2 0.4 1 3 10 100
[IP3] = 10 µM
0.05 0.00 0.0
0.5
1.0 time (s)
1.5
2.0
Fig. 5. Data obtained by superfusion of 45 Ca2þ -loaded hepatic microsomes. Upper panel: at time t ¼ 0; 0:4 mM Ca2þ and the indicated IP3 concentrations were added to the stimulus buffer. Lower panel: at time t ¼ 0; 10 mM IP3 and the indicated Ca2þ concentrations were added to the stimulus buffer. Complete details of the experimental method are in Dufour et al. (1997).
(Swatton and Taylor, 2002). However, recent preliminary results (Falcke, unpublished work) raise the interesting possibility that depletion may be more important than previously thought. For now, we shall ignore this complication. From Fig. 5 it can be seen that the time courses were collected every 72 ms; a resolution too low to resolve the response peak for each time course. However, data collected at the higher resolution of 9 ms from type 2 IPR in permeabilised hepatocytes (Marchant et al., 1997; Adkins and Taylor, 1999; Swatton and Taylor, 2002) or from type 3 IPR in RINm5F cells (Swatton and Taylor, 2002) show that the times to peak are generally greater than 100 ms: Thus, the low time resolution of the data used here is unlikely to introduce severe systematic errors in the parameter estimation. Of course, fits to higher resolution data are preferable, but this is left for future work. For each model we constructed the posterior probability distribution of the parameters, f ; given the data, d; i.e., Prðf jdÞ; by using a Bayesian formalism. Thus Prðf jdÞpPrðdjf ÞPrðf Þ: Given any parameter set, we calculate the probability of the data given the parameters, Prðdjf Þ; by solving the model numerically using those parameters (with the IMSL routine DIVPAG) and
ARTICLE IN PRESS 130
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
assuming that the error at each data point was Gaussian distributed, with variance ð0:02Þ2 : We have no detailed knowledge about the statistical properties of the noise in the data, and thus, in the absence of any information about the error distribution, this was the simplest reasonable assumption. The value for the variance was chosen arbitrarily, but makes no essential difference to the fit, only to the absolute value of the log probability. Hence log Prðdjf Þ ¼
X jd*i di j2 ; 2ð0:02Þ2 i
ð12Þ
where d*i is the computed approximation at the ith data point, di is the ith data point, and where i runs over all the data points. The probability of the parameters, Prðf Þ; is the prior, which is chosen from prior knowledge of reasonable values for the parameters. We chose a non-negative, exponential prior. Thus, Prðf Þ ¼ 0; if any fi is negative, while otherwise Prðf Þ ¼ expðf =fmax Þ; where fmax is a vector of values that sets reasonable upper bounds on the parameter values. Different choices of fmax have little effect on the results; its major function is to ensure that the posterior distribution is normalisable. Once the posterior distribution is constructed, we wish to determine the mean of each parameter value. Alternatively, as is commonly done, we could maximise the posterior distribution, which corresponds to a regularised minimising the sum of squares in (12). However, the maximum of the distribution does not necessarily give the best parameter estimate, as the bulk of feasible parameter values could be nowhere near the maximising value for these skewed distributions. Thus, we choose to characterise the distribution by the parameter means. We determine parameter means and variations by MCMC sampling of the posterior distribution using the Metropolis–Hastings algorithm. For each run we take random initial conditions. It is important to note that, by sampling from the posterior distribution, we not only determine the mean of the distribution (and thus the mean closeness of the fit) but we also determine how much parameter variation is allowed by the data (and thus some idea of the sensitivity of the fit to the parameters). Parameters that converge to a distribution with a small variance are well defined by the fitting procedure, while parameters that converge to a distribution with a large variance have little effect on how well the model fits the data. This is illustrated in Fig. 6. Here, we do not perform an exhaustive analysis of the posterior distributions, merely point out some overall trends. Although, in general, the sampler did not display ergodic properties for the length of runs used, the resulting chains easily determined the feasible parameter ranges. When fitting these three models, we use the rate constants rather than their ratios, i.e., we fit using ki and ki rather than ki =ki : However, it sometimes happens that, in sampling from the posterior distribution, it becomes clear that only a ratio of parameters is determined, not the parameters themselves. This is illustrated in Fig. 7, which uses samples taken from a run using the Sneyd–Dufour model. It is clear that, although the run of samples for k4 and k4 do not converge well, the ratio k4 =k4 converges well and is thus determined by the fit. This indicates that this particular transition in the receptor model could be well approximated by a quasi-steady-state approximation. We have not tested each of the sample runs exhaustively for convergence of every possible parameter combination, but point out, in the appropriate captions, some of the convergent ratios we have found.
ARTICLE IN PRESS J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140 5000
6
4000 k -3
5 k2
131
4
3000 2000
3
1000
2
0 0
0
1000 2000 sample number/10
2000 sample number/10
Fig. 6. A typical run of samples from the posterior distribution generated by fitting the De Young–Keizer model to data, shown after the initial burn-in period. The left panel shows an example of a parameter that is unambiguously determined by the fit, while the right panel shows an example of a parameter that is not well determined by the fit.
2.5
k-4
2.0 1.5 1.0 0.5 0.0 0
2000
0
2000
4000
6000
80
k4
60 40 20 0 4000
6000
-3
k-4 /k4
30x10
20 10 0 0
500
1000
1500
2000
2500
sample number/10 Fig. 7. A typical run of samples from the posterior distribution generated by fitting the Sneyd–Dufour model to data, shown after the initial burn-in period. Although the run of samples for k4 and k4 do not converge well (top two panels), it is clear that variation in one is accompanied by variation in the other. Plotting their ratio (bottom panel) we see that the ratio k4 =k4 converges well and is thus determined by the fit.
ARTICLE IN PRESS 132
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
3.2. Fitting the De Young–Keizer model The first model we fit using this technique is the De Young–Keizer model, and results from a typical run of samples are shown in Fig. 6. Sample runs from two different parameters are included, each showing a different kind of behaviour. For many of the parameters, after an initial burn-in period caused by the randomly chosen starting point, the samples settle down to a distribution with a clear mean and variance, as shown in the sample run of k2 (left panel). We conclude that the mean value of k2 is 3.81, with standard deviation 0.18. However, the sample run for k3 (right panel) shows rather different behaviour, wandering from almost 0 to over 5000 without ever reaching a steady distribution. In addition, for this sample run, log Prðdjf Þ had mean 590, with standard deviation 2.4. We conclude from this that the data are not sufficient to determine k3 unambiguously; large variations in k3 make almost no difference to the posterior probability. Using the values given in Table 1 we plot the responses to Ca2þ or IP3 steps, corresponding to the data shown in Fig. 5. The model responses are shown in Fig. 8. The responses to steps of IP3 (Fig. 8A) agree well with the data, but the responses to steps of Ca2þ (Fig. 8B) become much too fast at higher ½Ca2þ : In fact, the responses become so fast that they are indistinguishable from the axes at this level of resolution. Interestingly, the steady-state open probability curves (Fig. 8C) move to the left with increasing ½IP3 ; the opposite of the behaviour when the original parameters are used. 3.3. Fitting the Sneyd–Dufour model Because of the greater number of parameters in this model, the data were less able to determine them unambiguously. Fewer than half the parameters converged during the MCMC sampling, and wide variations in many parameter values made almost no difference to the fit. A summary of the results is shown in Table 2, and typical solutions are shown in Fig. 9. In one typical sample run log Prðdjf Þ had mean 241, with standard deviation 2.3. From this perspective the Sneyd–Dufour model provides a better fit to the data than does the De Young–Keizer model. In addition, the Table 1 Parameters of the De Young–Keizer model, determined by fitting to kinetic data sd or range k1 ¼ 83:8 (400) k2 ¼ 3:81 (0.2) k3 ¼ 11 103 (400) k4 ¼ 0:15 (0.2) k5 ¼ 53:9 (20) k1 ¼ 114 (52) k2 ¼ 0:409 (0.21) k3 ¼ 2500 (377) k5 ¼ 4:52 (1.64)
sd ¼ 9:8 sd ¼ 0:18 7–16 103 0–0.35 sd ¼ 2:8 sd ¼ 20 sd ¼ 0:051 0–5000 sd ¼ 0:46
Original parameters of De Young and Keizer (1992) are shown in parentheses. If the parameter is determined unambiguously by the fit, the standard deviation is given. Otherwise, the typical range of values is given. k4 was determined by ensuring detailed balance and thus k1 k2 k3 k4 ¼ k1 k2 k3 k4 : All concentrations are in mM:
ARTICLE IN PRESS
relative release
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140 [IP3] = 1 µM 3 5 10 30 50
2+
[Ca ] = 0.4 µM
0.25
133
0.20 0.15 0.10 0.05 0.00 0.5
1.0 1.5 time (s)
(A)
relative release
0.30
2+
[Ca ] = 0.02 µM 0.1 0.2 0.4 1 3
[IP3] = 10 µM
0.25
2.0
0.20 0.15 0.10 0.05 0.00 0.0
0.5
1.0 1.5 time (s)
(B)
[IP3] = 0.02 0.2 1 2
0.010 open probability
2.0
0.008 0.006 0.004 0.002 0.000 0.01
(C)
0.1
1
10
100
2+
[Ca ] (µM)
Fig. 8. Responses of the De Young–Keizer model, using the best-fit parameters shown in Table 1. (A) Responses at the indicated ½Ca2þ ; using the same IP3 concentrations as in Fig. 5, upper panel. (B) Responses at the indicated ½IP3 ; using the same Ca2þ concentrations as in Fig. 5, lower panel. (C) The steady-state open probability curves at four different ½IP3 :
responses to step increases in ½Ca2þ (Fig. 9, middle panel) are not too fast at high ½Ca2þ ; as they are in the De Young–Keizer model. The steady-state open probability curves again shift to the left with increasing ½IP3 : 3.4. Fitting the Dawson–Lea–Irvine model When fitting the Dawson–Lea–Irvine model to the data we find, again, that the parameters are not well determined by the data, with large variations in parameters leading to little change in the
ARTICLE IN PRESS 134
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
Table 2 Parameters of the Sneyd–Dufour model, determined by fitting to kinetic data sd or range k1 ¼ 20:9 k2 ¼ 18:7 k3 ¼ 0:5 k4 ¼ 64 L1 ¼ 0:125 L3 ¼ 0:07 L5 ¼ 2346 k1 ¼ 1:2 k2 ¼ 0:8 k3 ¼ 241 k4 ¼ 1:5 l2 ¼ 0:36 l4 ¼ 4:2 l6 ¼ 50 898
sd ¼ 3:1 sd ¼ 3:1 0–2 24–97 sd ¼ 0:016 0.02–0.11 180–4400 sd ¼ 0:12 0–4 1–500 0.5–2.8 0–0.75 sd ¼ 0:6 10–140 000
If the parameter is determined unambiguously by the fit, the standard deviation is given. Otherwise, the range is given. The remaining parameters are determined by ensuring detailed balance, and thus l2 ¼ l2 k1 =ðk1 L1 Þ; l6 ¼ k4 l6 =ðk4 L5 Þ; l4 ¼ k2 l4 =ðk2 L3 Þ: All concentrations are in mM; time in seconds. Although neither k4 nor k4 were determined by the fit, the ratio k4 =k4 was, with a mean of 0.024 and standard deviation of 0.0024.
posterior (Table 3). However, the ratios k2 =k5 and k6 =k9 were determined unambiguously, leading to the prediction that these transitions could be eliminated by using a quasi-state-state approximation without sacrificing model accuracy. In a sample run log Prðdjf Þ had mean 760, with standard deviation 2.4. Typical solutions, and the steady-state curves, are shown in Fig. 10. For the given parameters, the steady-state curves have a shape quite unlike experimental data.
4. Discussion and conclusions Each of the three models we discuss was constructed with a particular aim in mind. The De Young–Keizer model parameters were chosen so as to reproduce the steady-state data of Bezprozvanny et al. (1991), and also so that receptor activation by Ca2þ would be faster than Ca2þ -induced inactivation. The principal motivation of the Dawson–Lea–Irvine model was to study how adaptation can occur in a receptor model. Only the Sneyd–Dufour model was originally constructed specifically to study the kinds of dynamic responses shown in Fig. 5. It is thus no great surprise that the Sneyd–Dufour model does the best job of fitting this particular data set. Nevertheless, by comparing these models it is clear that the use of simple mass-action binding of Ca2þ to the receptor results in kinetics that are much too fast, and response peaks that occur too early. We can thus conclude that saturation of the rate of Ca2þ binding is likely an important feature of the receptor kinetics. This conclusion holds even though the majority of the data points were in the tails of the responses, with relatively few around the peaks, and thus the
ARTICLE IN PRESS J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140 [IP3] = 1 µM 3 5 10 30 50
0.25 relative release
2+
[Ca ] = 0.4 µM
0.20
135
0.15 0.10 0.05 0.00 0.0
0.5
(A)
2.0
2+
[Ca ] = 0.02 µM 0.1 0.2 0.4 1 3
[IP3] = 10 µM
0.20 relative release
1.0 1.5 time (s)
0.15 0.10 0.05 0.00 0.0
0.5
(B)
1.0 1.5 time (s)
open probability
0.04
2.0
[IP3] = 0.02 0.2 1 2
0.03 0.02 0.01 0.00 0.01
(C)
0.1 1 2+ [Ca ] (µM)
10
100
Fig. 9. (A–C) Responses of the Sneyd–Dufour model, using the best-fit parameters shown in Table 2.
model fits will be biased towards giving agreement to the response tails. The best way to correct this would be to use a data set at a higher time resolution, but this is left for future work. In addition, the steady-state behaviour of the Dawson–Lea–Irvine model was not in good agreement with data, while the steady-state curves of the other two models exhibited the correct bell-shaped curves, but with peaks moving to the left as ½IP3 increased. Although the Sneyd– Dufour model gives the best quantitative fit to the selected data (i.e., with the highest mean likelihood), the Dawson–Lea–Irvine model the worst, it would be dangerous to draw conclusions as to which is the ‘‘best’’ model. For instance, the Sneyd-Dufour model paid little attention to adaptation of the IPR, while, as discussed below, there is uncertainty as to the correct steady-state
ARTICLE IN PRESS 136
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
Table 3 Parameters of the Dawson–Lea–Irvine model determined by the fit sd or range k1 ¼ 0:045 (1) k5 ¼ 1345 (1000) k9 ¼ 70 (100) k10 ¼ 0:033 (1) k12 ¼ 0:85 (1) k15 ¼ 19; 970 (100) k18 ¼ 45 (1) k19 ¼ 22 (10) k1 ¼ 173 (100) k2 ¼ 51; 522 (1000) k6 ¼ 24 (10) k11 ¼ 0:25 (1) k13 ¼ 2859 (10) k15 ¼ 12:25 (10) k18 ¼ 0:053 (0.1)
0.004–0.25 sd ¼ 367 17–115 104 –0.14 0.0001–2.2 sd ¼ 5; 531 sd ¼ 7 0.4–37 2–450 17,000–83,000 7–35 0.0003–0.7 sd ¼ 737 sd ¼ 1:14 106 0:12
Original parameters of Dawson et al. (2003) are shown in parentheses. Parameters that do not appear in this list were determined by detailed balance. For instance, k10 ¼ k10 k6 k1 k2 =ðk2 k1 k6 Þ: Although neither k5 nor k2 were determined by the fit, the ratio k2 =k5 was, with a mean of 33.3 and standard deviation of 2. The ratio k6 =k9 was also determined by the fit, with mean 0.5 and standard deviation 0.1.
data that should be used to constrain the fit. Nevertheless, it seems safe to conclude that the De Young–Keizer model exhibits responses that are much too fast at high ½Ca2þ ; and is thus not a good model in that regime, while the steady-state properties of Dawson–Lea–Irvine model need improvement. The best-fit parameters for the De Young–Keizer model have Ca2þ activation (k5 ) an order of magnitude faster than Ca2þ inactivation (k2 and k4 ), while binding of IP3 (k1 and k3 ) is faster still. Similarly, the best-fit parameters for the Sneyd–Dufour model have Ca2þ activation (governed principally by k4 and l6 ) much faster than Ca2þ inactivation (governed principally by k1 and l2 ), while IP3 binding (k2 ) is relatively slow. Thus, the major differences between the De Young– Keizer and the Sneyd–Dufour models are, firstly, the rate of response at high ½Ca2þ ; and, secondly, the relative rate of IP3 binding. In fact, contrary to a procedure that is often used, the Sneyd–Dufour model suggests that the binding of IP3 cannot be assumed to be fast enough to be eliminated from the equations. The best-fit parameters for the Dawson–Lea–Irvine model again have Ca2þ activation (k15 ) faster than Ca2þ inactivation (k18 ). Although binding of the first three IP3 molecules proceeds very quickly, binding of the fourth (k5 ) proceeds somewhat slower than Ca2þ activation. Thus, the Dawson–Lea–Irvine model also suggests that the kinetics of IP3 binding cannot be ignored in IPR models. When the best-fit parameters are used in the De Young–Keizer or Sneyd–Dufour model, the peak of the bell-shaped curve moves to the left with increasing ½IP3 : Since this is the opposite of the behaviour of the De Young–Keizer model when the original parameters are used, it emphasises the fact that this feature, being parameter-dependent, cannot be used to distinguish between models. In addition, it raises the important question of which steady-state data should be
ARTICLE IN PRESS J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140 [IP3] = 1 µM 3 5 10 30 50
relative release
2+
[Ca ] = 0.4 µM
0.30
137
0.20 0.10 0.00 0.0
0.5
(A)
relative release
0.25
1.0 1.5 time (s)
2.0
2+
[Ca ] = 0.02 µM 0.1 0.2 0.4 1 3
[IP3] = 10 µM
0.20 0.15 0.10 0.05 0.00 0.0
0.5
(B)
2.0
[IP3] = 0.2 2 10 20
0.06 open probability
1.0 1.5 time (s)
0.05 0.04 0.03 0.02 0.01 0.00 0.01
(C)
0.1 1 2+ [Ca ] (µM)
10
100
Fig. 10. (A–C) Responses of the Dawson–Lea–Irvine model, using the best-fit parameters shown in Table 3.
used to constrain the fitting procedure further. At present, there is no easy answer to this question. In each of the models discussed here, the steady-state open probabilities are low, less than 0.08. Although this agrees with the experimental data of Kaftan et al. (1997), Hagar et al. (1998) and Hagar and Ehrlich (2000), as well as being consistent with the time-dependent measurements of Finch et al. (1991), Swatton and Taylor (2002), and Marchant and Taylor (1997), it is in direct conflict with the data of Mak et al. (1998, 2000, 2001) and Ramos-Franco et al. (2000). Unfortunately, there is no agreement as to which set of data better represents in vivo behaviour. At this stage, the only thing modellers can do is wait until experimentalists manage to resolve these differences.
ARTICLE IN PRESS 138
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
Another obvious conclusion from this exercise in model comparison is that time-dependent responses to steps of ½Ca2þ and ½IP3 are insufficient to determine the model parameters unambiguously. In each of the models, large variations in model parameters make very little difference to the quality of the fit, resulting in a lack of convergence of the sampling procedure. The only way to improve this is to use additional data in the fit, or otherwise to constrain some of the parameters. Some obvious experimental data that have not been used in our fitting method are the time course of recovery from inactivation (for instance, see Fig. 3A of Swatton and Taylor, 2002), steady-state data (once agreement can be reached on what the correct data are), responses to steps of ½Ca2þ (½IP3 ) at a wider variety of ½IP3 (½Ca2þ ), and responses to double pulses of Ca2þ or IP3 : Ideally, the next generation of models will provide a good fit to this entire range of steady-state and kinetic data, exhibit the required degree of adaptation, incorporate possible vesicle or ER depletion, take into account the possible existence of microdomains around the channel mouth, and provide different parameter sets or model structures for each of the receptor sub-types. Construction of such a family of models, and fitting them to this much wider array of data will be a challenging problem.
Acknowledgements J. Sneyd was supported by the Marsden Fund of the Royal Society of New Zealand, and by the Mathematical Biosciences Institute of The Ohio State University. J.-F. Dufour was supported by Swiss National Found, Grant 3100A0-100513.
References Adkins, C.E., Taylor, C.W., 1999. Lateral inhibition of inositol 1,4,5-trisphosphate receptors by cytosolic Ca2þ : Curr. Biol. 9, 1115–1118. Atri, A., Amundson, J., Clapham, D., Sneyd, J., 1993. A single-pool model for intracellular calcium oscillations and waves in the Xenopus laevis oocyte. Biophys. J. 65, 1727–1739. Berridge, M.J., 1997. Elementary and global aspects of calcium signalling. J. Physiol. 499, 291–306. Bezprozvanny, I., 1994. Theoretical analysis of calcium wave propagation based on inositol (1,4,5)-trisphosphate ðInsP3 Þ receptor functional properties. Cell Calcium 16, 151–166. Bezprozvanny, I., Watras, J., Ehrlich, B.E., 1991. Bell-shaped calcium-response curves of Ins(1,4,5)P3 - and calciumgated channels from endoplasmic reticulum of cerebellum. Nature 351, 751–754. Cardy, T.J., Traynor, D., Taylor, C.W., 1997. Differential regulation of types-1 and -3 inositol trisphosphate receptors by cytosolic Ca2þ : Biochem. J. 328, 785–793. Clapham, D., 1995. Calcium signaling. Cell 80, 259–268. Combettes, L., Claret, M., Champeil, P., 1993. Calcium control on InsP3 -induced discharge of calcium from permeabilised hepatocyte pools. Cell Calcium 14, 279–292. Dawson, A.P., Lea, E.J., Irvine, R.F., 2003. Kinetic model of the inositol trisphosphate receptor that shows both steady-state and quantal patterns of Ca2þ release from intracellular stores. Biochem. J. 370, 621–629. De Young, G.W., Keizer, J., 1992. A single pool IP3 -receptor based model for agonist stimulated Ca2þ oscillations. Proc. Natl. Acad. Sci. USA 89, 9895–9899. Dufour, J.-F., Arias, I.M., Turner, T.J., 1997. Inositol 1,4,5-trisphosphate and calcium regulate the calcium channel function of the hepatic inositol 1,4,5-trisphosphate receptor. J. Biol. Chem. 272, 2675–2681.
ARTICLE IN PRESS J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
139
Finch, E.A., Turner, T.J., Goldin, S.M., 1991. Calcium as a coagonist of inositol 1,4,5-trisphosphate-induced calcium release. Science 252, 443–446. Gilks, W.R., Richardson, S., Spiegelhalter, D.J., 1995. Markov Chain Monte Carlo in Practice. CRC Press, Boca Raton, FL. Hagar, R.E., Ehrlich, B.E., 2000. Regulation of the type III InsP3 receptor by InsP3 and ATP. Biophys. J. 79, 271–278. Hagar, R.E., Burgstahler, A.D., Nathanson, M.H., Ehrlich, B.E., 1998. Type III InsP3 receptor channel stays open in the presence of increased calcium. Nature 396, 81–84. ! Hajnoczky, G., Thomas, A.P., 1997. Minimal requirements for calcium oscillations driven by the IP3 receptor. EMBO J. 16, 3533–3543. Hodgkin, A.L., Huxley, A.F., 1952. A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544. Iino, M., 1990. Biphasic Ca2þ dependence of inositol 1,4,5-trisphosphate-induced Ca2þ release in smooth muscle cells of the Guinea Pig Taenia Caeci. J. Gen. Physiol. 95, 1103–1122. Iino, M., Endo, M., 1992. Calcium-dependent immediate feedback control of inositol 1,4,5-trisphosphate-induced Ca2þ release. Nature 360, 76–78. Iino, M., Tsukioka, M., 1994. Feedback control of inositol trisphosphate signalling by calcium. Mol. Cell. Endocrinol. 98, 141–146. Ilyin, V., Parker, I., 1994. Role of cytosolic Ca2þ in inhibition of InsP3 -evoked Ca2þ release in Xenopus oocytes. J. Physiol. 477, 503–509. Kaftan, E.J., Ehrlich, B.E., Watras, J., 1997. Inositol 1,4,5-trisphosphate (InsP3 ) and calcium interact to increase the dynamic range of InsP3 receptor-dependent calcium signaling. J. Gen. Physiol. 110, 529–538. LeBeau, A.P., Yule, D.I., Groblewski, G.E., Sneyd, J., 1999. Agonist-dependent phosphorylation of the inositol 1,4,5trisphosphate receptor: a possible mechanism for agonist-specific calcium oscillations in pancreatic acinar cells. J. Gen. Physiol. 113, 851–871. Mak, D.-O.D., McBride, S., Foskett, J.K., 1998. Inositol 1,4,5-trisphosphate activation of inositol trisphosphate receptor Ca2þ channel by ligand tuning of Ca2þ inhibition. Proc. Natl. Acad. Sci. USA 95, 15821–15825. Mak, D.O., McBride, S., Raghuram, V., Yue, Y., Joseph, S.K., Foskett, J.K., 2000. Single-channel properties in endoplasmic reticulum membrane of recombinant type 3 inositol trisphosphate receptor. J. Gen. Physiol. 115, 241–256. Mak, D.-O.D., McBride, S., Foskett, J.K., 2001. Regulation by Ca2þ and inositol 1,4,5-trisphosphate (InsP3 ) of single recombinant type 3 InsP3 receptor channels: Ca2þ activation uniquely distinguishes types 1 and 3 InsP3 receptors. J. Gen. Physiol. 117, 435–446. Marchant, J.S., Taylor, C.W., 1997. Cooperative activation of IP3 receptors by sequential binding of IP3 and Ca2þ safeguards against spontaneous activity. Curr. Biol. 7, 510–518. Marchant, J.S., Taylor, C.W., 1998. Rapid activation and partial inactivation of inositol trisphosphate receptors by inositol trisphosphate. Biochemistry 37, 11524–11533. Marchant, J.S., Chang, Y.-T., Chung, S.-K., Irvine, R.F., Taylor, C.W., 1997. Rapid kinetic measurements of 45 Ca2þ mobilization reveal that Ins(2,4,5)P3 is a partial agonist at hepatic InsP3 receptors. Biochem. J. 321, 573–576. Monod, J., Wyman, J., Changeux, J.P., 1965. On the nature of allosteric transition: a plausible model. J. Mol. Biol. 12, 88–118. Othmer, H., Tang, Y., 1993. In: Othmer, H., Murry, J., Maini, P. (Eds.), Experimental and Theoretical Advances in Biological Pattern Formation. Plenum Press, London. Parker, I., Ivorra, I., 1990. Inhibition by Ca2þ of inositol trisphosphate-mediated Ca2þ liberation: a possible mechanism for oscillatory release of Ca2þ : Proc. Natl. Acad. Sci. USA 87, 260–264. Parker, I., Yao, Y., 1991. Regenerative release of calcium from functionally discrete subcellular stores by inositol trisphosphate. Proc. R. Soc. London B 246, 269–274. Parker, I., Yao, Y., Ilyin, V., 1996. Fast kinetics of calcium liberation induced in Xenopus oocytes by photoreleased inositol trisphosphate. Biophys. J. 70, 222–237. Parys, J.B., Sernett, S.W., DeLisle, S., Snyder, P.M., Welsh, M.J., Campbell, K.P., 1992. Isolation, characterization, and localization of the inositol 1,4,5-trisphosphate receptor protein in Xenopus laevis oocytes. J. Biol. Chem. 267, 18776–18782.
ARTICLE IN PRESS 140
J. Sneyd et al. / Progress in Biophysics & Molecular Biology 85 (2004) 121–140
Ramos-Franco, J., Bare, D., Caenepeel, S., Nani, A., Fill, M., Mignery, G., 2000. Single-channel function of recombinant type 2 inositol 1,4,5-trisphosphate receptor. Biophys. J. 79, 1388–1399. Sanderson, M.J., Charles, A.C., Boitano, S., Dirksen, E.R., 1994. Mechanisms and function of intercellular calcium signaling. Mol. Cell. Endocrinol. 98, 173–187. Sienaert, I., Missiaen, L., De Smedt, H., Parys, J.B., Sipma, H., Casteels, R., 1997. Molecular and functional evidence for multiple Ca2þ -binding domains in the type I inositol 1,4,5-trisphosphate receptor. J. Biol. Chem. 272, 25899–25906. Sneyd, J., Dufour, J.F., 2002. A dynamic model of the type-2 inositol trisphosphate receptor. Proc. Natl. Acad. Sci. USA 99, 2398–2403. Sneyd, J., Keizer, J., Sanderson, M.J., 1995. Mechanisms of calcium oscillations and waves: a quantitative analysis. FASEB J. 9, 1463–1472. Swatton, J.E., Taylor, C.W., 2002. Fast biphasic regulation of type 3 inositol trisphosphate receptors by cytosolic calcium. J. Biol. Chem. 277, 17571–17579. Taylor, C.W., 1998. Inositol trisphosphate receptors: Ca2þ -modulated intracellular Ca2þ channels. Biochim. Biophys. Acta 1436, 19–33. Taylor, C.W., Laude, A.J., 2002. IP3 receptors and their regulation by calmodulin and cytosolic Ca2þ : Cell Calcium 32, 321–334. ! Thomas, A.P., Bird, G.S.J., Hajnoczky, G., Robb-Gaspers, L.D., Putney, J.W.J., 1996. Spatial and temporal aspects of cellular calcium signaling. FASEB J. 10, 1505–1517. Yao, Y., Parker, I., 1992. Potentiation of inositol trisphosphate-induced Ca2þ mobilization in Xenopus oocytes by cytosolic Ca2þ : J. Physiol. 458, 319–338. Yoneshima, H., Miyawaki, A., Michikawa, T., Furuichi, T., Mikoshiba, K., 1997. Ca2þ differentially regulates the ligand-affinity states of type 1 and type 3 inositol 1,4,5-trisphosphate receptors. Biochem. J. 322, 591–596.