Comparison of models for IP3 receptor kinetics using stochastic simulations.

Inositol 1,4,5-trisphosphate receptor (IP3R) is a ubiquitous intracellular calcium (Ca(2+)) channel which has a major role in controlling Ca(2+) levels in neurons. A variety of computational models have been developed to describe the kinetic function of IP3R under different conditions. In the field of computational neuroscience, it is of great interest to apply the existing models of IP3R when modeling local Ca(2+) transients in dendrites or overall Ca(2+) dynamics in large neuronal models. The goal of this study was to evaluate existing IP3R models, based on electrophysiological data. This was done in order to be able to suggest suitable models for neuronal modeling. Altogether four models (Othmer and Tang, 1993; Dawson et al., 2003; Fraiman and Dawson, 2004; Doi et al., 2005) were selected for a more detailed comparison. The selection was based on the computational efficiency of the models and the type of experimental data that was used in developing the model. The kinetics of all four models were simulated by stochastic means, using the simulation software STEPS, which implements the Gillespie stochastic simulation algorithm. The results show major differences in the statistical properties of model functionality. Of the four compared models, the one by Fraiman and Dawson (2004) proved most satisfactory in producing the specific features of experimental findings reported in literature. To our knowledge, the present study is the first detailed evaluation of IP3R models using stochastic simulation methods, thus providing an important setting for constructing a new, realistic model of IP3R channel kinetics for compartmental modeling of neuronal functions. We conclude that the kinetics of IP3R with different concentrations of Ca(2+) and IP3 should be more carefully addressed when new models for IP3R are developed.

Introduction

Inositol 1,4,5-trisphosphate receptor (IP3R) is a ligand-gated calcium (Ca2+) release channel typically expressed on the endoplasmic reticulum (ER) in neurons and many other cell types. It has a major role in intracellular Ca2+ dynamics which, in turn, is involved in many cellular processes such as muscle contraction, neurotransmitter release, vesicle secretion, fertilization, gene transcription, immunity, and apoptosis. In neurons, dynamical changes in Ca2+ concentration ([Ca2+]) are involved, among others, in neuroplasticity and development (see recent reviews [1], [2]), and in neurodegeneration (see [3], [4]). Transient, repetitive changes in cytosolic Ca2+ concentration are crucial for synapse modification and plasticity, including long-term potentiation (LTP) and long-term depression (LTD) [5]–[8]. These phenomena constitute the biological basis for learning and memory formation in the brain [8], [9]. Particularly in the cerebellum, IP3Rs are relatively highly expressed in Purkinje cells [10]. Ca2+ release from ER has been shown to be a key mediator of cerebellar LTD [11].

The inositol 1,4,5-trisphosphate receptor is a tetrameric receptor-channel, consisting of four sub-units. In total, three different genes (ITPR1, ITPR2, and ITPR3) encode three different types (1, 2, and 3) of IP3R and their splice variants from which homo- or heterotetramers can form [12]. IP3R is activated and opened by both IP3 and Ca2+. Ca2+ can also act as the inhibitor of IP3R in higher concentrations. IP3 is produced from phosphatidylinositol 4,5-bisphosphate (PIP) by phospholipase C (PLC). After a cell is stimulated (for example by glutamate in neurons) certain G protein- or tyrosine kinase-linked receptors are activated. These, in turn, can activate PLC. ER acts as a Ca2+ store, and while open, IP3R can release Ca2+ from ER lumen to the cytosol. Transient rises or oscillations in Ca2+ concentration can then activate various enzymes and even induce changes in the transcriptional level. IP3Rs are known to be responsible for the phenomenon called Ca2+-induced Ca2+ release (CICR), in addition to ryanodine receptors (RyRs) [13], [14].

In order to develop models for ion channels and receptors detailed data on the structure and function of the modeled entity is required. The function of IP3R has been studied with electrophysiological techniques. However, since IP3Rs are prevalently located on the endoplasmic reticulum of a cell, performing the recordings is not straightforward. The first recordings performed on IP3Rs involved isolated microsomes from smooth muscle cells incorporated into artificial lipid bilayer [15]. Later, the same technique has been used, for example, for IP3R in canine cerebellum [16]–[20], in mouse cerebellum [21], and in HEK cells [22] (IP3R recombinantly expressed). IP3Rs have also been recorded from the plasma membranes of DT40 cells [23] (IP3R endogenously expressed (native)) and DT40-3KO cells [24], [25] (stably expressed IP3R construct, native IP3R ablated). Since the nuclear membrane is a continuation of the ER, IP3Rs have also been recorded from isolated nuclei of Xenopus oocytes (for example [26] (recombinantly expressed and native IP3Rs), Purkinje neurons and granule cells [27], [28] (IP3R endogenously expressed), and DT40 cells [23], [29]. These kind of data are of great value when developing a model for ion channel kinetics. However, the electrophysiological raw data on IP3R is not available in any of the publicly available databases, but its statistics is described in publications. For example, the dependence of open probability on cytosolic Ca2+ or IP3 concentrations is given ([16], [19], [20], [29], [30]). In some cases, the open and closed time distributions [18], [20], [22] or mean open time [17], [18], [22], [29] are also reported. In an ideal case, the raw data would be publicly available in a database and a modeler could extract all needed statistical measures out of the data or use the raw data for automated estimation of model parameter values.

In addition to electrophysiological measurements, Ca2+ imaging and radioactive assays have also been used to study the behavior of IP3R in vitro. For example, Fujiwara et al. [31] analyzed the kinetics of Ca2+ release via IP3R in controlled cytoplasmic environment in permeabilized cerebellar Purkinje cells. In addition, superfusion and 45Ca2+ release assay (radioactive assay) have been used for studying the Ca2+ release and inhibition of IP3R by Ca2+ in hepatic microsomes [32]–[34]. These kind of studies give more detailed information on the IP3R regulation by IP3and Ca2+ and their affinities than electrophysiological studies. In some cases, the data obtained from Ca2+ imaging studies or from radioactive assays has been used in modeling studies, for example Fujiwara et al. [31] by Doi et al. [32] and Dufour et al. [35] by Sneyd et al. [36].

In order to reach a better understanding of the dynamical behavior of IP3R, as well as its involvement in various cellular processes, it is of interest to build models of IP3R. Computational models are important for understanding the time evolution, dynamics, and regulation of ion channels and intracellular proteins and enzymes [37], [38]. Several models have previously been proposed to describe the behavior of IP3R (for a comprehensive review, see, for example [39]). There are models presented for different types of IP3R (type 1, 2, and 3) [12] in different animals, tissues and cells (for example Xenopus oocyte [40], cerebellar cells [41], pancreatic acinar cells [42], and hepatic cells [32]). The first and most well-known model is the one by De Young and Keizer [41]. Some models for IP3R have been compared either analytically or by means of simulation [36], [43]–[45], and later reviewed [39], [46].

The majority of the existing models is deterministic. Deterministic approaches, however, do not give biologically valid results and are not always capable of modeling the random behavior observed with small numbers of molecules [47]–[50]. Stochastic modeling is therefore more and more used for describing the dynamics of a biochemical system. The stochastic approach is always valid whenever the deterministic approach is valid, but when the deterministic is not, the stochastic might sometimes be valid [51]. Most commonly, deterministic methods and, in some cases, analytical methods are used to investigate the properties of IP3R models (see, for example [43] or [52]). More rarely, stochastic methods are applied [53], [54], even though it is known that the behavior of ion channels is stochastic.

Despite the wealth of IP3R models the selection of a specific model for describing IP3R related calcium dynamics or signaling is not straightforward. The models are seldom generic in nature and capable of describing all possible data obtained for a specific IP3R or cell type. The reason for this is that the models are developed for some specific purpose, describe the behavior only in certain experimental conditions, or the dynamics are not fully analyzed to validate the model. This can be due to the limited access to experimental data. We therefore wanted to study the dynamics of existing models in detail and to specifically address their suitability in the context of complex neuronal models. In this work, the interest is set on the type 1 IP3R because it is most commonly expressed in neurons [10]. After a preliminary study, we chose four models [35], [55]–[57] for a more detailed analysis and comparison. Other models did not meet our criteria. The chosen models were originally developed by using data either from IP3R in canine cerebellum or type 1 IP3R. As the selected models are biophysically realistic and based on the law of mass action, they can be implemented to the stochastic simulation tool STEPS [58], [59] used in this study. Additionally, we decided to concentrate on computationally inexpensive IP3R models so that it would be possible to integrate them as part of larger model for calcium dynamics or synaptic plasticity. We validated the functionality of the models by comparing the statistical behavior of IP3R channel kinetics (open probability curves, mean open times, and open and closed time distributions) to the equivalent obtained by electrophysiological recordings from IP3Rs expressed in neurons.

Our results show firstly, that the behavior of the studied models varies in similar simulation conditions and, secondly, some models show quite unrealistic kinetic behavior. We therefore conclude that the kinetics of IP3R (open and closed times and the open probability) with different concentrations of both Ca2+ and IP3 should be more carefully addressed when new models for IP3R are developed.

Materials and Methods

In our present work, after a preliminary review on existing IP3R models, we selected four models [35], [55]–[57] for comparison. The selection was based on the following criteria: (1) relative simpicity (i.e. the model should have less than 20 states), (2) development based on data obtained from neuronal or type 1 IP3R, and (3) basis in the law of mass action (the reactions include binding and unbinding reaction and state transitions). As our ultimate goal is to find a model that can be an integral part of a larger model for Ca2+ dynamics or synaptic plasticity in neurons, it is an advantage to have a structurally simple model. The selected models are based on the law of mass action and can thus be implemented into the stochastic simulators such as STEPS [58].

Models

The model of Othmer and Tang

The model of Othmer and Tang [55] is one of the earliest and small- scaled models regarding the number of states. There is only four states, since the binding order of Ca2+ or IP3 is not free, but sequential, opposite to the models of De Young and Keizer [41] or Bezprozvanny and Ehrlich [18]. Othmer and Tang [55] assume that IP3 has to bind to its binding site before Ca2+ can bind and the channel can open, as well as the activating Ca2+ has to bind to its site before the inhibition by Ca2+ can occur. The schematic representation of the model of Othmer and Tang [55] in Figure 1A and the parameter values in Table 1 were used in this study. The model of Othmer and Tang [55] has been used before as a part of a larger model for calcium dynamics, for example, by Mishra and Bhalla [60].

Figure 1

Schematic representation of the states and transitions of the IP3R models.

(A) Othmer and Tang [55] (forward direction of a reaction is to the right) (B) Doi et al. [35] (forward direction of a reaction is to the right or up), (C) Fraiman and Dawson [57] (forward direction of a reaction is to the right or down) (D) Dawson et al. [56] (forward direction of a reaction is the to the direction of binding a ligand or in the plain state transitions from left to the right).

Table 1

Rate constants for IP3R model of Othmer and Tang [55].

Reaction	kf	kb
r1
r2
r3

r1 to r3 refer to reactions represented in Figure 1A.

The model of Dawson et al

Dawson et al. [56] built a model for IP3R, using a RyR model by Sachs et al. [61] as their starting point, to understand the adaptive and incremental behavior of IP3R. The model of Dawson et al. [56] is applicable to type 1 and 2 IP3Rs and with some modification to type 3. Dawson et al. [56] assume that IP3R has two conformations, R and P. The conformation R can bind four IP3molecules rapidly, but with low affinity, to reach an open state. The conformation P, on the other hand, slowly binds four IP3 molecules, but with high affinity, to reach a closed state where it is thereafter possible to reach the open state. In this work, Scheme 2 from the original paper was used with two exceptions: the flux through an open channel (reactions 14 and 16 in the original paper) and the diffusion of released Ca2+ (reaction 17) were not taken into account in order to make the model comparable with other models. This does not have an effect on the actual channel kinetics of the receptor as the removed reactions deal with Ca2+ flux and diffusion. Moreover, we used constant Ca2+ concentration and the simulated reactions happened in well-mixed system and in the present work only the kinetics of the IP3R, not Ca2+ dynamics was studied. We used the the model presented in Figure 1D and the parameter values given in Table 2.

Table 2

Rate constants for IP3R model of Dawson et al. [56].

Reaction	kf	kb	Reaction	kf	kb
r1			r9
r2			r10
r3			r11
r4			r12
r5			r13
r6			r15
r7			r18
r8			r19

r1 to r19 refer to reactions presented in Figure 1D.

The model of Fraiman and Dawson

The IP3R model of Fraiman and Dawson [57] was originally built to study the effects of different Ca2+ concentrations inside the ER to the kinetics of IP3R. It is the only model included in the present study that has a Ca2+ binding site inside the ER in addition to the cytosolic binding sites. The state scheme of the model of Fraiman and Dawson [57] is presented in Figure 1C and the parameter values used in this work are in Table 3.

Table 3

Rate constants for IP3R model of Fraiman and Dawson [57], taken from [67].

Reaction	kf	kb
r1
r2
r3
r4
r5
r6
r7
r8
r9
r10
r11
r12
r13
r14

r1 to r14 refer to reactions represented in Figure 1C.

Originally, six states, O, O, O, P, P, and P, were considered open. However, it has been experimentally shown that IP3R needs IP3 to reach a stable open conformation [33], [62]. For this reason, we neglected three of the original open states (i.e., they were considered closed) in the present work and only states O, O, and O were considered open. In addition, in the original publication [57], the rate constant of the transition from A10 to A00 is defined as ‘detailed balance’, with no given numerical value. In our study, it was mandatory to have a numerical value for the parameter and thus we fixed the parameter by testing three values with open probability simulations (data not shown). The parameter values of 0 s−1 and 200 s−1 produced identical results which were in accordance with the results in the original publication [57], while the value of 2000 s−1 slightly upraised the left side of the open probability curve. Based on these simulations we chose the value of 200 s−1 for the transition from A10 to A00 (reaction 7, k) and concluded that it was in the range of what was originally used.

The model of Doi et al

The IP3R model of Doi et al. [35] was originally published as part of a larger model for Ca2+ dynamics in the cerebellar Purkinje cell spine to investigate the role of IP3Rs as a coincidence detector of two input signals. Doi et al. [35] constructed their model based on a conceptual model of Adkins and Taylor [34]. Doi et al. [35] used experimental data by Khodakhah and Ogden [63], Marchant and Taylor [33], and Fujiwara et al. [31] to define the structure and kinetics of the model and experimental data by Bezprozvanny et al. [16] to test how well the model can reproduce the bell-shaped curve. A schematic representation of the model is presented in Figure 1B and the rate constants for each reaction in Table 4. In the model of Doi et al. [35], IP3R has seven states and the receptor needs to bind both IP3 and Ca2+ to open and thus provide Ca2+ flux from ER lumen to cytosol. In this model, IP3R has one open state, RIC.

Table 4

Rate constants for IP3R model of Doi et al. [35].

Reaction	kf	kb
r1
r2
r3
r4
r5
r6

r1 to r6 refer to reactions represented in Figure 1B.

Simulations and data analysis

In the present study, the simulations were designed to reproduce the data produced in experimental electrophysiological measurements from neuronal IP3Rs. We used stochastic simulation approaches since deterministic approaches were not applicable due to the stochastic nature of ion channel gating. The simulated data was compared with experimental data available in literature. The four selected models were implemented according to the information presented in the original publications with some exceptions presented in the section ‘Models’. Our work does not include parameter estimation (as, for example, [36]) since raw data on channel kinetics of IP3Rs in neurons is not publicly available.

In this work, STEPS (STochastic Engine for Pathway Simulation) ([58], [59]; http://steps.sourceforge.net/) version 1.1.2 was used for simulation. With STEPS, it is possible to perform full stochastic simulation of reactions and diffusion of molecules in three dimensions and also deterministic simulations. For stochastic simulations, STEPS uses the stochastic simulation algorithm (SSA) described by Gillespie [64]. The model scripts are available at ModelDB (http://senselab.med.yale.edu/ModelDB/).

In our simulations, we assumed a well-mixed system. Our models had two compartments, cytosol and ER lumen, each having volume of 0.1 fl and a surface, ER, between them. The IP3R was placed on the surface and the cytosolic concentrations of Ca2+ and IP3 were kept constant in the simulations to mimic the buffered conditions in patch-clamp recording.

The simulations were run on a stand-alone Linux computer. For open probability curves, simulations were repeated, depending on the model, 750–12 000 times and averaged over the repetitions for each data point. To produce one such curve, the simulations lasted from an hour to several hours. Simulations for open and closed time distributions were run once for 10–5000 s to obtain sufficient number of events to get statistically significant results. These computations took from less than a second to a couple of seconds each. Analysis of the simulated data was performed and the figures were drawn with MATLAB [65].

Results

We compared four kinetic models previously developed for IP3 receptor function by simulating them with the Gillespie stochastic simulation algorithm of STEPS simulator. The comparison was done by analyzing the steady-state behavior, such as the open probability, open and closed time distributions, and the mean open and closed time. Here we show that the behavior of the models varies and some models behave somewhat unrealistically.

Open probability

It has been experimentally shown that the open probability () of IP3R is dependent on the cytosolic Ca2+ concentration and that the dependence is bell-shaped [16]. We repeated similar experiments by computational means and tested whether the selected four models are capable of expressing the bell-shaped curve. All the models except the model of Dawson et al. [56] produced the bell-shaped curve (see Figure 2A). Instead, the model of Dawson et al. [56] (blue in Figure 2A) produced an s-shaped curve similarly as in a previous comparison study by Sneyd et al. [36]. The model of Othmer and Tang [55] (green in Figure 2A) reaches the highest ( = 0.33) at cytosolic Ca2+ concentration around 80 nM. The model of Doi et al. [35] (magenta in Figure 2A) and the model of Fraiman and Dawson[57] (red in Figure 2A) reach the highest (  = 0.15 and  = 0.38, respectively) around [Ca2+]  =  300 nM, which is closest to the experimentally obtained values ([Ca2+]  = 250 nM by Bezprozvanny et al. [16] and [Ca2+]  =  200 nM by Kaznacheyeva et al. [22]). The absolute value of obtained in simulations cannot be directly compared to the experimental data, because Bezprozvanny et al. [16] and Kaznacheyeva et al. [22] report only normalized values, not absolute values, for .

Figure 2

Open probability of IP3R as a function of (A) cytosolic Ca2+ concentration (IP3  = 10 M) and (B) cytosolic IP3 concentration (Ca2+  = 0.25 M).

Green: Othmer and Tang [55], Blue: Dawson et al. [56], Red: Fraiman and Dawson [57], Magenta: Doi et al. [35].

The open probability of IP3R is also dependent on cytosolic IP3 concentration (see for example [17], [27], [29]). The open probability curves of the models obtained in simulations are shown in Figure 2B. All the models except the model of Dawson et al. [56] (blue in Figure 2B) follow the s-shape that is reported in experimental studies [17], [27], [29]. In their study on IP3Rs on Purkinje cell nuclear membrane, Marchenko et al. [27] have shown that the stays close to 0 until IP3 concentration reaches 0.3 M and keeps rising until IP3 concentration is 3 M ([Ca2+]  =  .25 M). Watras et al. [17] have shown that the rise starts when IP3 concentration is 0.03 M and settles after 1 M. The in models of Dawson et al. [56] (blue in Figure 2B) and Doi et al. [35] (magenta in Figure 2B) starts rising approximately at the same IP3 concentration as in [27], but the elevation does not stop at the right concentrations. In the models of Othmer and Tang [55] (green in Figure 2B) and Fraiman and Dawson [57] (red in Figure 2B), starts rising one or two orders of magnitude too low when compared to the experimental results.

Kaftan et al. [19] have shown in their experiments on cerebellar IP3R that the bell-shaped Ca2+-dependence curve moves upward and to the right when IP3 concentration is increased. They used IP3 concentration values of 0.02, 0.2, 2, and 180 M. We used the same concentrations, in addition to their fivefold values, except 180 M in our simulation for all the models (results in Figure 3). The model of Othmer and Tang [55] (Figure 3A) shows a shift upward and to the left, the model of Dawson et al. [56] (Figure 3B) upward, and the models of Fraiman and Dawson [57] (Figure 3C) and Doi et al. [35] (Figure 3D) upward and slightly to the left when IP3 concentration increases. Similar trend has also been shown for the model of Othmer and Tang [55] by Diambra and Guisoni [66] and Tang et al. [43]. None of the models reproduced the results presented by Kaftan et al. [19].

Figure 3

Open probability of IP3R as a function of cytosolic Ca2+ concentration in different IP3 concentrations.

(A) Othmer and Tang [55] (B) Dawson et al. [56] (C) Fraiman and Dawson [57] (D) Doi et al. [35].

Mean open and closed times and distributions of open and closed times

Bezprozvanny and Ehrlich [18] reported that the mean open time of canine cerebellar IP3R is 2.9 0.2 ms and Kaznacheyeva et al. [22] that the mean open time of wild-type rat cerebellar IP3R is 4.20.5 ms and that the open and closed times have exponential distributions (black dashed line in Figure 4E and 4K) in certain experimental conditions (lipid bilayer experiments, [IP3]  =  2 M, [Ca2+]  =  0.2 M). We simulated the selected models in these same conditions (Sim 1, results in Table 5 and Figure 4A–F) and, in order to take into account the affinity difference [31], with five times greater IP3 concentration (Sim 2, results in Table 5 and Figure 4G–L). The mean open times of the model of Fraiman and Dawson [57] are 2.5 ms (Sim 1) and 2.6 ms (Sim 2). These values are close to the experimentally obtained values. The mean open times obtained with the other models are an order of magnitude smaller (0.5 ms for Dawson et al. [56] and Doi et al. [35]) or significantly greater (460 ms, Othmer and Tang [55]). None of open time distributions of the selected models (Figures 4A-C and 4G-I) follow the experimental distribution by Kaznacheyeva et al. [22] fully, but all give, however, the exponential shape (see Figures 4B and 4K). The open time distribution of the model of Fraiman and Dawson [57] is the closest to experimentally [22] obtained distribution (see Figures 4B, 4H). The same applies also to the closed time distributions (see Figures 4E, 4K).

Figure 4

Distribution of IP3R open and closed times for all the selected models obtained in simulation conditions Sim 1 (A–F) and Sim 2 (G–L).

(A) Open time distributions of all the models in conditions Sim 1, (B) Enlarged from A, (C) Enlarged from B, (D) Closed time distributions of all the models conditions Sim 1, (E) Enlarged from D, (F) Enlarged from E, (G) Open time distributions of all the models conditions Sim 2, (H) Enlarged from G, (I) Enlarged from H, (J) Closed time distributions of all the models conditions Sim 2, (K) Enlarged from J, (L) Enlarged from K. Experimental data is from [22]. In simulation conditions Sim 1 [Ca2+]  = 0.2 M, [IP3]  =  2 M and Sim 2 [Ca2+]  = 0.2 M, [IP3]  = 10 M (as shown in Table 6).

Table 5

Mean open and closed times of IP3R of the selected models.

	Model	mean open time (ms)	mean closed time (ms)	n	simulation time (s)
Sim 1	Othmer and Tang	451.19±423.06	12892563	1068	1 800
	Dawson et al.	0.59±5.46	10.33120.24	1797	20
	Fraiman and Dawson	2.45±2.52	4.0111.22	1535	10
	Doi et al.	0.470.46	11.2138.08	1711	20
Sim 2	Othmer and Tang	463.55463.96	12902793	1045	1 800
	Dawson et al.	0.524.70	9.38167.55	1897	20
	Fraiman and Dawson	2.572.76	4.9219.84	1391	10
	Doi et al.	0.470.46	3.7623.11	1501	6
Sim 3	Othmer and Tang	510.08526.46	1047±2074	1927	3 000
	Dawson et al.	0.465.90	8.83±97.56	2004	20
	Fraiman and Dawson	2.482.64	5.38±2.64	1293	10
	Doi et al.	0.530.53	19.00±29.60	1024	20
Sim 4	Othmer and Tang	509.68525.59	958.50±2073	2044	3 000
	Dawson et al.	0.6510.64	5.21±100.02	2063	10
	Fraiman and Dawson	2.512.67	5.50±15.55	1249	10
	Doi et al.	0.510.50	4.72±0.50	1866	10
Sim 5	Othmer and Tang	598.32598.68	3356±3384	1263	5 000
	Dawson et al.	0.250.25	12.60±129.80	1161	10
	Fraiman and Dawson	2.472.60	25.18±88.87	1446	40
	Doi et al.	0.470.47	107.37±123.06	1854	200
Sim 6	Othmer and Tang	596.98592.01	2712±2709	1509	5 000
	Dawson et al.	0.250.26	9.27±163.25	2098	20
	Fraiman and Dawson	2.492.61	27.54±95.62	1331	40
	Doi et al.	0.460.46	27.92±47.50	1407	40

The different simulation conditions (Sim 1 – Sim 6) are presented in Table 6.

Table 6

Ca2+ and IP3 concentration used in different simulations for open and closed time distributions.

	[Ca2+] (μM)	[IP3] (μM)
Sim 1	0.2	2
Sim 2	0.2	10
Sim 3	0.1	2
Sim 4	0.1	10
Sim 5	0.01	2
Sim 6	0.01	10

The simulations were done in the same conditions as wet-lab experiments [20], [22] and with five times greater IP3 concentration in order to take into account the affinity difference between in vivo and lipid bilayer experiments [31].

Moraru et al. [20] have presented open time distributions for canine cerebellar IP3R in two different conditions (lipid bilayer experiments, [Ca2+]  = 0.1 and 0.01 M, and [IP3]  = 2 M) (black dashed line in Figures 5 and 6). We simulated the behavior of the selected models in these same experimental conditions (Sim 3 and Sim 4, results in Table 5 and Figure 5) and also with fivefold IP3 concentration (Sim 5 and Sim 6, results in Table 5 and Figure 6). The distributions in the wet-lab experiments are of exponential shape [18]–[20], [22] and simulation results also show exponential shape for all the models. The only distributions that are also otherwise similar to the ones obtained in wet-lab experiments by Moraru et al. [20] are the distributions of the model of Fraiman and Dawson [57] (Figures 5B, 5H, 6B, and 6H). All the simulation conditions used are summarized in Table 6.

Figure 5

Distributions of IP3R open and closed times for all the selected models obtained in simulation conditions Sim 3 (A–F) and Sim 4 (G–L).

(A) Open time distributions of all the models in conditions Sim 3, (B) Enlarged from A, (C) Enlarged from B, (D) Closed time distributions of all the models in conditions Sim 3, (E) Enlarged from D, (F) Enlarged from E, (G) Open time distributions of all the models conditions Sim 4, (H) Enlarged from G, (I) Enlarged from H, (J) Closed time distributions of all the models conditions Sim 4, (K) Enlarged from J, (L) Enlarged from K. Experimental data is from [20]. In simulation conditions Sim 3 [Ca2+] = 0.1 M, [IP3] = 2 M and Sim 4 [Ca2+]  =  0.1 M, [IP3]  = 10 M (as shown in Table 6).

Figure 6

Distribution of IP3R open and closed times for all the selected models obtained in simulation conditions Sim 5 (A–F) and Sim 6 (G–L).

(A) Open time distributions of all the models conditions Sim 5, (B) Enlarged from A, (C) Enlarged from B, (D) Closed time distributions of all the models conditions Sim 5, (E) Enlarged from D, (F) Enlarged from E, (G) Open time distributions of all the models conditions Sim 6, (H) Enlarged from G, (I) Enlarged from H, (J) Closed time distributions of all the models conditions Sim 6, (K) Enlarged from J, (L) Enlarged from K. Experimental data is from [20]. In simulation conditions Sim 5 [Ca2+]  = 0.01 M, [IP3]  = 2 M and Sim 6 [Ca2+]  = 0.01 M, [IP3]  = 10 M (as shown in Table 6).

The Ca2+ concentrations used in the experiments by Moraru et al. [20] are unfortunately at the border or smaller than those observed in a neuron at resting level (i.e., Ca2+ used is 0.1 M or less). As IP3R is, however, known to have functional significance only above the resting level concentrations, more emphasis should be put on physiological conditions in experimental work in the future. In other words, experimental work should additionally be performed with Ca2+ concentrations above the known resting level.

Discussion

In this work, four models of IP3R [35], [55]–[57] were selected among many to examine their steady-state and time series behavior and compare them with experimental data available in literature. We implemented and simulated the selected models using stochastic simulation software STEPS in order to obtain similar data as in single-channel patch-clamp recordings. The open probability curves and statistics, such as the mean open time and open and closed time distributions, were compared to experimental ones obtained in the same conditions. To our knowledge, this is the first detailed evaluation of IP3R model kinetics with stochastic methods. Our comparative study shows significant differences in the behavior and kinetics of the studied models.

Based on our results, the statistical properties of the model of Fraiman and Dawson [57] seem to be the most similar to the ones obtained in wet-lab experiments. The properties of the model of Othmer and Tang [55] are very different when compared to the experimental data. All the models except the model of Dawson et al. [56] produce the bell-shaped open probability curve for Ca2+-dependence and the s-shaped open probability curve for IP3-dependence as seen in the electrophysiological experiments (for example [16], [17], [27]). However, none of the models reproduce the experimental finding presented by Kaftan et al. [19], which shows that Ca2+-dependent open probability curve moves to the right and upward when IP3 concentration increases. This kind of behavior is shown in the original article by Fraiman and Dawson [57]. The reason why the simulation of the same model in this study did not produce similar behavior might be the slight modification that we were forced to make to the model (defining a numerical value for the one parameter that was originally defined as ‘detailed balance’ and neglecting three of the six open states). It is also notable that there is an Errata [67] published for the original article [57] and that we used the parameter set in the Errata [67].

The simulated open and closed time distributions of all the models follow the exponential distribution as does the data from experiments [18]–[20], [22]. However, the distributions are not similar apart from the distribution of Fraiman and Dawson [57]. The reason for this may be the relatively simple structure of the models, insufficiency of modeled states to reproduce the kinetics, and parameter values that do not fit the data.

According to our results, the mean open time of model of Doi et al. [35] is not congruent with the experimental findings. However, the shape and peak value of the open probability curve are in accordance with experimental data. As the model of Doi et al. [35] has originally been published as part of a larger signal transduction model for LTD induction, some inaccuracy in the behavior of the model could have been corrected by other parameters, such as the Ca2+ flux rate and thus the small mean open time does not invalidate the results in the original publication.

As our comparative study points out significant differences in the behavior and kinetics of the studied models, it is of interest to consider reasons for it. We identify four major reasons why the selected models behave differently to each other: 1) the structure (i.e. the equations) and parameter values differ between the models, 2) experimental data that was used in the model development vary, 3) different data handling procedures have been used when developing the models, and 4) model developers did not use automated parameter estimation methods. Next, we will discuss each issue in detail.

Firstly, the most obvious reason for differences in the behavior of models is the structure and parameter values of the models. All the models studied here have different number of states, but this does not cause the differences as such. More importantly, different parameter values and thus the affinities of IP3, as well as activating and inactivating Ca2+, vary between the models. Since the models of Othmer and Tang [55] and Doi et al. [35] reproduce the correct shapes for the open probability curves, re-estimation of their parameters might improve the fitting of models to experimental data. As a general conclusion, all studies neither report the values of all parameters used in simulations nor make it evident which parameter set is used to produce specific results. This makes it difficult to reproduce results (see also discussion in [68]).

Secondly, another reason for the differences in the behavior of the models could be related to the variability in the use of experimental data when constructing the original model. Although the statistical properties of channel kinetics, such as the mean open time and the distributions of open times, are known to be important in properly reconstructing receptor-ion channel kinetics, they are relatively rarely used in developing or evaluating models for IP3R. Furthermore, there exists a clear difference on how experimental data is used to construct (i.e., to define the structure, the number of states, and the number of parameter values in the model) and fine-tune the models (estimation of the unknown parameters). We have noticed that it is not always clear which data is used in modeling and, particularly, how it is used. In general, the models presented for IP3R are constructed based on only some of the data or knowledge obtained from various animal species and experiments. Furthermore, data on kinetics of IP3R have been obtained from various sources: native and recombinantly expressed receptors in cell lines and Xenopus oocytes, and from vertebrate cerebellum or hepatocytes.

Doi et al. [35] use the model of Adkins and Taylor [34] as their starting point and construct the model based on data by Marchant and Taylor [33] and use the open probability curve of Bezprozvanny et al. [16] to study the fitness of their model. The model of Othmer and Tang [55] is also shown to fit the data by Bezprozvanny et al. [16] in addition to data by Watras et al. [17] in [43], but this study does not take the difference in IP3 affinity [31] into account as Doi et al. [35] or study the open or closed time distributions of the model. Fraiman and Dawson [57] and Dawson et al. [56] use several experimental observations when constructing their model, but they do not report using any data for actual fitting of the model parameters. The data that Dawson et al. [56] compare their model to is more dealing with temporal aspect of Ca2+ release and accumulation of Ca2+ to cytosol than actual channel kinetics.

Thirdly, the differences between the simulated and experimentally observed open time distributions and mean open times might also be due to differences in data handling procedures. Experimentally observed open time distributions can be biased due to the limitations and established practices regarding the temporal resolution in the patch-clamp recordings, while in the simulations in this study all the events are recorded exactly at the time they happen. Usually the time resolution in patch-clamp recordings is around 1 ms and thus any opening shorter than that would stay unnoticed or be merged with other channel openings.

Fourthly, to our knowledge, automated parameter estimation methods have not been used in the development of the four models here compared. Studies on IP3R models consider, to some extent, the kinetic ion channel data to define the mathematical structure of the models. However, only a few previous studies use automated parameter estimation techniques and statistical data on ion channel kinetics to fine-tune the IP3R models [36], [69]–[72].

One of the major challenges in modeling the IP3Rs is the lack of access to original raw data, for example from electrophysiological measurements, that could be used in quantitative modeling. This data is not currently available in any public database and as the years pass by it becomes extremely hard to acquire the data from its original sources. This problem is not new or limited just to measurements of ion channels but to all neuroscience data [73], [74]. Some suggestions to improve the situation have been made. For instance, De Schutter [75] suggests that data publishing should be distinguished from paper publishing. Furthermore, Ranjan et al. [76] have established an information management framework for ion channel information, which hopefully will make IP3R experimental data more accessible in the future.

Despite several shortcomings in the development and presentation of models, previous models on IP3R, including the present comparative study on four stochastic IP3R models, will give a good setting for constructing a new, realistic model of IP3Rs for compartmental modeling of neuronal functions. It will be a challenge to develop computationally inexpensive models that can produce realistic stochastic behavior of an individual ion channel. A wealth of evidence indicates, however, an important role of randomly opening ion channels on the global behavior of cells. For example, in neurons the stochastic openings of single ion channels shape the integration of local signals in dendrites or spines [77], stochastic openings of voltage-gated ion channels have an important role in adjusting the transmembrane voltage dynamics [78]–[80], and the reliability of action potential propagation along thin axons is affected by the stochastic opening of voltage-gated ion channels [81]. Furthermore, molecular noise of single ion channel is shown to be translated into global cellular processes in astrocytes [82].

In summary, the development of new IP3R models clearly calls for both steady-state and kinetic data. Fitting of the new computational models should be done using automated estimation techniques, possibly using Bayesian approaches [72], [83]–[85]. Data for model construction and fine-tuning would ideally be acquired from the same neuronal type as the model is built for.