Time to blip – stochastic simulation of single channel opening

The stochastic dynamics of the inositol-1,4,5-trisphosphate (IP3) receptor (IP3R) is key to understanding a wide range of observed calcium (Ca2+) signals (Falcke 2004). The stochastic nature results from the constant binding and unbinding of Ca2+ and IP3 to and from their respective binding sites and is especially important in the initiation of a Ca2+ puff, i.e. the release of Ca2+ through a cluster of IP3Rs. Once the first IP3R opens, the Ca2+ concentration rises significantly around the ion channel and hence increases the open probability for neighboring IP3Rs. In turn this may trigger the activation of further receptors giving rise to a Ca2+ puff (Thul et al. 2009; Thurley et al. 2012). In this protocol, we determine the time that it takes for a single IP3R to open from rest. We explicitly take into account the tetrameric structure of the IP3R and the fact that multiple subunits need to be active before the channel opens (Bezprozvanny et al. 1991; Watras et al. 1991). We develop code for a stochastic simulation of the IP3R and simulate it using the software package Matlab (Attaway 2011). This protocol demonstrates the basic form of a stochastic simulation algorithm and may serve as a starting point to investigate more complex gating dynamics.


Materials
• Matlab (http://www.mathworks.com/products/matlab/)Method 1. Derive the transition scheme for the stochastic dynamics.Suppose that each of the four subunits of the IP3R exists in only two states (as e.g. in the Li-Rinzel model (Li and Rinzel 1994)), viz.active and inactive.Then there are five states for the entire IP3R when we assume that all subunits act independently: 4 inactive subunits (C4) to no inactive subunit (C0).When the transition rate between the active and the inactive state is  and between the inactive and active state is  we arrive at the transition scheme depicted in Figure 1.  4. Run the code by either typing stochastic2state at the command line or by pressing the run-button (green triangle) in the Matlab editor.5.The code shown in step 3 simulates 50000 (Nr) channel openings and saves the time it takes for each individual opening in the array openarr.
As suggested by experiments, it is assumed that the channel opens when three subunits are in the active state (while open<4).Every channel opening starts from the configuration when all four subunits are in the closed state (open=0).The if statement represents the stochastic transitions that are illustrated in step 2. The last three lines produce the plot shown in Figure 2. Solution: Check the combinatorial factors, which here are the numbers 1 to 4 in the transition probabilities.Make sure that the time step  is small enough.The sum of all transition probabilities needs to be smaller than 1.

Problem (
Step 3): The produced figure differs from Figure 2. Solution: Since we here study stochastic simulations, no two simulations are the same.To generate reproducible stochastic simulations, we need to initiate the random number generator with the same seed every time we run stochastic2state.m(see the rng command in Matlab).

Problem (
Step 4): The code does not run from the command line.Solution: Make sure that the path in the command line is the same as the one to where stochastic2state was saved.

Discussion
Figure 2 demonstrates that there is a time lag before IP3Rs open from rest.The most likely delay corresponds to the peak of the distribution.The exact shape of the blip distribution depends on the gating scheme used and parameter values.For this protocol, we have chosen the simplest IP3R model.For more sophisticated models, there are further considerations.Firstly, the transition scheme becomes more involved and often cannot be represented for the entire ion channel, but instead for a single subunit.Secondly, we have to use more sophisticated numerical methods such as developed in (Gillespie 1977).The method used here presents an approximation of (Gillespie 1977) and works well for small time steps and a limited number of states.For more detail, I refer the reader to (Thul and Falcke 2006;2007;Higgins et al. 2009)

Figure 1 :
Figure 1: Gating scheme.The state space of an IP3R with four independent subunits and two states per subunit can be represented by the five states C4 to C0, where the subscript indicates the number of inactive subunits.The arrows present transitions between the states along with the corresponding transition rates.

Figure 2 :
Figure 2: The histogram shows the distribution of times for a single IP3R to open.Note the delay and the sharp rise of the blip initiation times.
Write down the stochastic transition step.Suppose e.g. the channel is in the state C3, then there is a probability of 3 of making a transition to C2, a probability of  of making a transition to C4 and a probability of 1 − 3 +   of remaining in the state C3 during a small time interval .When we draw a random number  that is uniformly distributed between 0 and 1, if • 0 ≤  < 3 then make transition to C2, • 3 ≤  < 3 +   then make transition to C4, • 3 +   ≤  stay in C3. 3. Code up the ideas from step 1 and 2 in Matlab as shown below.Save the code as stochastic2state.m.