Comparison of Models for IP3 Receptor Kinetics Using Stochastic Simulations

Inositol 1,4,5-trisphosphate receptor (IP3R) is a ubiquitous intracellular calcium (Ca2+) channel which has a major role in controlling Ca2+ 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 Ca2+ transients in dendrites or overall Ca2+ 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 Ca2+ and IP3 should be more carefully addressed when new models for IP3R are developed.


Introduction
Inositol 1,4,5-trisphosphate receptor (IP 3 R) is a ligand-gated calcium (Ca 2z ) release channel typically expressed on the endoplasmic reticulum (ER) in neurons and many other cell types. It has a major role in intracellular Ca 2z 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 Ca 2z concentration ([Ca 2z ]) are involved, among others, in neuroplasticity and development (see recent reviews [1,2]), and in neurodegeneration (see [3,4]). Transient, repetitive changes in cytosolic Ca 2z concentration are crucial for synapse modification and plasticity, including long-term potentiation (LTP) and long-term depression (LTD) [5][6][7][8]. These phenomena constitute the biological basis for learning and memory formation in the brain [8,9]. Particularly in the cerebellum, IP 3 Rs are relatively highly expressed in Purkinje cells [10]. Ca 2z 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 IP 3 R and their splice variants from which homo-or heterotetramers can form [12]. IP 3 R is activated and opened by both IP 3 and Ca 2z . Ca 2z can also act as the inhibitor of IP 3 R in higher concentrations. IP 3 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 Ca 2z store, and while open, IP 3 R can release Ca 2z from ER lumen to the cytosol. Transient rises or oscillations in Ca 2z concentration can then activate various enzymes and even induce changes in the transcriptional level. IP 3 Rs are known to be responsible for the phenomenon called Ca 2z -induced Ca 2z 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 IP 3 R has been studied with electrophysiological techniques. However, since IP 3 Rs are prevalently located on the endoplasmic reticulum of a cell, performing the recordings is not straightforward. The first recordings performed on IP 3 Rs involved isolated microsomes from smooth muscle cells incorporated into artificial lipid bilayer [15]. Later, the same technique has been used, for example, for IP 3 R in canine cerebellum [16][17][18][19][20], in mouse cerebellum [21], and in HEK cells [22] (IP 3 R recombinantly expressed). IP 3 Rs have also been recorded from the plasma membranes of DT40 cells [23] (IP 3 R endogenously expressed (native)) and DT40-3KO cells [24,25] (stably expressed IP 3 R construct, native IP 3 R ablated). Since the nuclear membrane is a continuation of the ER, IP 3 Rs have also been recorded from isolated nuclei of Xenopus oocytes (for example [26] (recombinantly expressed and native IP 3 Rs), Purkinje neurons and granule cells [27,28] (IP 3 R 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 IP 3 R 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 Ca 2z or IP 3 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, Ca 2z imaging and radioactive assays have also been used to study the behavior of IP 3 R in vitro. For example, Fujiwara et al. [31] analyzed the kinetics of Ca 2z release via IP 3 R in controlled cytoplasmic environment in permeabilized cerebellar Purkinje cells. In addition, superfusion and 45 Ca 2z release assay (radioactive assay) have been used for studying the Ca 2z release and inhibition of IP 3 R by Ca 2z in hepatic microsomes [32][33][34]. These kind of studies give more detailed information on the IP 3 R regulation by IP 3 and Ca 2z and their affinities than electrophysiological studies. In some cases, the data obtained from Ca 2z 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 IP 3 R, as well as its involvement in various cellular processes, it is of interest to build models of IP 3 R. 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 IP 3 R (for a comprehensive review, see, for example [39]). There are models presented for different types of IP 3 R (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 IP 3 R have been compared either analytically or by means of simulation [36,[43][44][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][48][49][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 IP 3 R 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 IP 3 R models the selection of a specific model for describing IP 3 R 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 IP 3 R 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 IP 3 R because it is most commonly expressed in neurons [10]. After a preliminary study, we chose four models [35,[55][56][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 IP 3 R in canine cerebellum or type 1 IP 3 R. 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 IP 3 R 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 IP 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 IP 3 R (open and closed times and the open probability) with different concentrations of both Ca 2z and IP 3 should be more carefully addressed when new models for IP 3 R are developed.

Materials and Methods
In our present work, after a preliminary review on existing IP 3 R models, we selected four models [35,[55][56][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 IP 3 R, 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 Ca 2z 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 Ca 2z or IP 3 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 IP 3 has to bind to its binding site before Ca 2z can bind and the channel can open, as well as the activating Ca 2z has to bind to its site before the inhibition by Ca 2z 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].  [56] is applicable to type 1 and 2 IP 3 Rs and with some modification to type 3. Dawson et al. [56] assume that IP 3 R has two conformations, R and P. The conformation R can bind four IP 3 molecules rapidly, but with low affinity, to reach an open state. The conformation P, on the other hand, slowly binds four IP 3 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 Ca 2z (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 Ca 2z flux and diffusion. Moreover, we used constant Ca 2z concentration and the simulated reactions happened in well-mixed system and in the present work only the kinetics of the IP 3 R, not Ca 2z dynamics was studied. We used the the model presented in Figure 1D and the parameter values given in Table 2.
The model of Fraiman and Dawson. The IP 3 R model of Fraiman and Dawson [57] was originally built to study the effects of different Ca 2z concentrations inside the ER to the kinetics of IP 3 R. It is the only model included in the present study that has a Ca 2z 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.
Originally, six states, O a , O b , O c , P a , P b , and P c , were considered open. However, it has been experimentally shown that  [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). doi:10.1371/journal.pone.0059618.g001 Table 1. Rate constants for IP 3 R model of Othmer and Tang [55].  [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

Simulations and data analysis
In the present study, the simulations were designed to reproduce the data produced in experimental electrophysiological measurements from neuronal IP 3 Rs. 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 IP 3 Rs 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    [35]. 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 IP 3 R was placed on the surface and the cytosolic concentrations of Ca 2z and IP 3 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

Open probability
It has been experimentally shown that the open probability (P o ) of IP 3 R is dependent on the cytosolic Ca 2z 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 P o (P o = 0.33) at cytosolic Ca 2z 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)  The open probability of IP 3 R is also dependent on cytosolic IP 3 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 IP 3 Rs on Purkinje cell nuclear membrane, Marchenko et al. [27] have shown that the P o stays close to 0 until IP 3 concentration reaches 0.3 mM and keeps rising until IP 3 concentration is 3 mM ([Ca 2z ] = .25 mM). Watras et al. [17] have shown that the rise starts when IP 3 concentration is 0.03 mM and settles after 1 mM. The P o 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 IP 3 concentration as P o 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), P o 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 IP 3 R that the bell-shaped Ca 2z -dependence curve moves upward and to the right when IP 3 concentration is increased. They used IP 3 concentration values of 0.02, 0.2, 2, and 180 mM. We used the same concentrations, in addition to their fivefold values, except 180 mM in our simulation for all the models (results in Figure 3  [66] and Tang et al. [43]. None of the models reproduced the results presented by Kaftan et al. [19].

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 IP 3 R is 2.9 + 0.2 ms and Kaznacheyeva et al. [22] that the mean open time of wild-type rat cerebellar IP 3 R is 4.2+0.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, [IP 3 ] = 2 mM, [Ca 2z ] = 0.2 mM). 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 IP 3 concentration (Sim 2, results in Table 5 and Figure 4G [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] 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 IP 3 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][19][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.
The Ca 2z 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., Ca 2z used is 0.1 mM or less). As IP 3 R 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 Ca 2z concentrations above the known resting level.

Discussion
In this work, four models of IP 3 R [35, [55][56][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 IP 3 R 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 Ca 2zdependence and the s-shaped open probability curve for IP 3dependence 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 Ca 2z -dependent open probability curve moves to the right and upward when IP 3 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][19][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 Ca 2z 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 The different simulation conditions (Sim 1 -Sim 6) are presented in Table 6. doi:10.1371/journal.pone.0059618.t005 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 IP 3 , as well as activating and inactivating Ca 2z , vary between the models. Since the models  Table 6). doi:10.1371/journal.pone.0059618.g006 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 IP 3 R. 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 IP 3 R are constructed based on only some of the data or knowledge obtained from various animal species and experiments. Furthermore, data on kinetics of IP 3 R 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 IP 3 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 Ca 2z release and accumulation of Ca 2z 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 IP 3 R 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 IP 3 R models [36,[69][70][71][72].
One of the major challenges in modeling the IP 3 Rs 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 IP 3 R experimental data more accessible in the future.
Despite several shortcomings in the development and presentation of models, previous models on IP 3 R, including the present comparative study on four stochastic IP 3 R models, will give a good setting for constructing a new, realistic model of IP 3 Rs 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][79][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 IP 3 R 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][84][85]. Data for model construction and fine-tuning would ideally be acquired from the same neuronal type as the model is built for. The simulations were done in the same conditions as wet-lab experiments [20,22] and with five times greater IP 3 concentration in order to take into account the affinity difference between in vivo and lipid bilayer experiments [31]. doi:10.1371/journal.pone.0059618.t006