Intra-Cluster Percolation of Calcium Signals

Calcium signals are involved in a large variety of physiological processes. Their versatility relies on the diversity of spatio-temporal behaviors that the calcium concentration can display. Calcium entry through inositol 1,4,5-trisphosphate (IP) receptors (IPR's) is a key component that participates in both local signals such as “puffs” and in global waves. IPR's are usually organized in clusters on the membrane of the endoplasmic reticulum and their spatial distribution has important effects on the resulting signal. Recent high resolution observations [1] of Ca puffs offer a window to study intra-cluster organization. The experiments give the distribution of the number of IPR's that open during each puff without much processing. Here we present a simple model with which we interpret the experimental distribution in terms of two stochastic processes: IP binding and unbinding and Ca-mediated inter-channel coupling. Depending on the parameters of the system, the distribution may be dominated by one or the other process. The transition between both extreme cases is similar to a percolation process. We show how, from an analysis of the experimental distribution, information can be obtained on the relative weight of the two processes. The largest distance over which Ca-mediated coupling acts and the density of IP-bound IPR's of the cluster can also be estimated. The approach allows us to infer properties of the interactions among the channels of the cluster from statistical information on their emergent collective behavior.


Introduction
The calcium (Ca 2z ) ion is a universal second messenger that is involved in a large number of physiological processes [2]. To this end, cells regulate cytosolic Ca 2z concentration ([Ca 2z ]) very precisely. At basal conditions free cytosolic [Ca 2z ] is very low (*100nM). [Ca 2z ] is several orders of magnitude higher in the extracellular medium and in internal reservoirs, such as the endoplasmic reticulum. Different signals can induce the opening of specific Ca 2z channels located on the plasma membrane or on the membrane of the internal reservoirs leading to local increments of the cytosolic [Ca 2z ] of various durations. This [Ca 2z ] change evokes different end responses depending upon the spatiotemporal distribution of [Ca 2z ]. Thus, it is of interest to measure the latter and how different factors shape it.
One of the Ca 2z channels involved in intracellular Ca 2z signals is the inositol 1,4,5-trisphosphate (IP 3 ) receptor (IP 3 R) which is expressed in many cell types and is located at the surface of intracellular membranes such as the endoplasmic reticulum (ER), the sarcoplasmic reticulum (SR) and the nucleus. The IP 3 R is biphasically regulated by Ca 2z , with a bell-shaped open probability as a function of [Ca 2z ]. Kinetic models of the IP 3 R take this dual effect into account by assuming that the receptor has at least one activating and one inhibitory site such that Ca 2z binding to the first one induces channel opening (provided that IP 3 is also bound to the receptor) and Ca 2z binding to the second one induces channel closing [3][4][5]. Given that the affinity for Ca 2z of the activating site is larger than that of the inhibitory site, a local increase of cytosolic Ca 2z in the vicinity of an IP 3 R with IP 3 bound induces channel opening first. This leads to a phenomenon called Ca 2z -induced Ca 2z -release (CICR) because the Ca 2z ions released by one channel can in turn trigger the opening of other nearby channels with IP 3 bound. Ca 2z channels are not uniformly distributed in the cell. IP 3 R's, in particular, are usually organized in clusters on the membrane of the ER that are separated by a few microns [6]. These clusters have been estimated to be 400nm|400nm in size in oocytes [7,8]. The simulations of [7] showed that previous observations could be reproduced assuming that between 25 and 35 IP 3 R's opened simultaneously during puffs. A similar estimate was obtained in [8] using a mean-field model that assumed that all channels opened and closed simultaneously. Simulations that include a stochastic description of the individual channel openings and closings, however, show that at most half of the channels with IP 3 bound are simultaneously open during a puff [8]. This implies that even in clusters with 50 IP 3 R's with IP 3 bound, the maximum number of simultaneously open channels is around 20. These results are consistent with observations of Ca 2z signals in the human neuroblastoma SY5Y cell line in which puffs of up to 20 simultaneously open channels were observed [1]. Measurements performed using patches of the outer nuclear envelope of the DT40 cell line give smaller numbers of IP 3 R's in each patch [9]. The non-uniform spatial organization of the IP 3 R's together with the channel coupling induced by CICR gives rise to a large variety of intracellular Ca 2z signals that go from very localized ones to waves that propagate throughout the cell [10].
The hierarchy of intracellular Ca 2z signals that includes Ca 2z ''blips'' (Ca 2z release through a single IP 3 R), ''puffs'' (Ca 2z release through several IP 3 R's in a cluster) and waves that propagate globally across cells by successive cycles of CICR has been observed using fluorescence microscopy and Ca 2z sensitive dyes [10][11][12][13]. The Xenopus laevis oocyte has been frequently used for this purpose because of its relatively large size and because the only Ca 2z channels that are present on the surface of the ER are IP 3 R's. Fluorescent images of these signals obtained with confocal microscopy do not resolve the inner-cluster structure. Therefore, different modeling strategies have been presented in order to determine the properties of the dynamics and spatial organization of IP 3 R's within clusters that are compatible with these experimental observations [7,8,14,15]. In particular, in [8,16] we made the very simple assumption that the number of IP 3 R's that open during the first puff that occurs at a site is given by the number of IP 3 R's with IP 3 bound. The underlying assumption was that the Ca 2z released by the first open channel would induce the opening of all the other IP 3 R's of the cluster with IP 3 bound. Therefore, if all the clusters had approximately the same number of IP 3 R's and all IP 3 R's were equally sensitive to IP 3 , the distribution of the number of channels that opened during a puff could be approximated by a binomial or Poisson distribution [8], provided that the probability that the channels become open were the same immediately before the occurrence of each puff. This last condition would not be satisfied in a non-stationary situation, e.g. if the concentrations of the agonists right before the release event differed significantly from puff to puff. It would not hold, in particular, for data containing sequences of puffs that are coupled through CICR or to puffs in which the inhibitory effect of the Ca 2z released in a previous event was noticeable, as described in [16]. In oocytes, the latter is only relevant for very long records containing many puffs at a site, which is usually not the case in most experiments. Calcium induced calcium release is also affected by buffers that can trap Ca 2z ions as they diffuse. This not only reduces the [Ca 2z ] but also alters the rate of Ca 2z transport [17]. The distances that separate IP 3 R's within a cluster are very small (10-20nm) [9]. Thus, only large concentrations of very fast buffers could affect Ca 2z -mediated inter-channel coupling in cases with many active channels [18,19]. The assumption that all the channels with IP 3 bound participate of the first puff of their site is the simplest way of approaching the complex problem of Ca 2zmediated inter-channel communication. Yet, it is applicable as long as the distance between IP 3 -bound channels is not too large. In the present paper we drop this assumption and analyze how Ca 2z -mediated inter-channel coupling affects the distribution of puff sizes. Our approach provides a simple tool to study some of the effects of buffers on the intra-cluster dynamics.
The quantal properties of Ca 2z release during puffs have recently been revealed in [1] using total internal reflection fluorescence (TIRF) microscopy in intact mammalian cells of the human neuroblastoma SY5Y cell line. The proximity of IP 3 R's to the plasma membrane in this cell type allowed the use of TIRF microscopy in which fluorescence can be elicited in a very small (attoliter) volume. This, together with the use of a fast CCD camera, permitted a much better temporal resolution than the one achieved with confocal microscopy. In this way, abrupt step-wise transitions between fluorescence levels were observed during the falling phase of puffs. Furthermore, many puffs could be elicited at each release site due to the use of a membrane-permeable form of IP 3 [20]. The authors then inferred that the step-wise transitions between fluorescence levels occurred in multiples of a basic unit that they identified with the amplitude contribution of each channel at the site [1]. Using this relationship they could readily obtain the distribution of the number of channels that open during a puff. Given that there is a large variability among cluster sites, they analyzed the subset of events that occur in clusters with a similar number of IP 3 R's. The authors did not find any sign of an inhibiting effect of the Ca 2z released in their records. In spite of that and even constraining the data set as mentioned before, they found that a Poisson distribution failed to reproduce the observed histogram of event sizes particularly in the region of small events (i.e., puffs with very few open channels). They could approximately describe the distribution with a model that assumes a weak cooperativity among channels. Inter-channel cooperativity is mediated by the Ca 2z released through an open IP 3 R that subsequently diffuses to a neighboring channel. Thus, the distance between channels is a key factor that regulates the cooperativity strength [21]. The approach of [1], however, does not take space into account.
In the present paper we introduce a simple model that takes into account both the stochasticity due to IP 3 binding and the distancedependent Ca 2z -mediated cooperativity. It can reproduce the event size distribution reported in [1] for events involving any number of open channels. The distribution obtained with our model approaches a binomial or Poisson distribution as the cooperativity strength increases so that the opening of one IP 3 R induces the opening of all other IP 3 R's with IP 3 bound. This transition from Ca 2z -dominated to IP 3 -binding dominated stochasticity is similar to a percolation transition. It also occurs if the number of IP 3 R's with IP 3 bound increases. Therefore, the transition can be reflected on the distribution of the number of IP 3 R's that open at a given release site.
Percolation in connection with Ca 2z signals has been invoked to explain the transition from abortive to propagating waves in cells [22][23][24]. Our paper is the first to identify two limiting regimes of the intra-cluster dynamics that underlies puffs and to characterize the change between them as a percolation transition. Furthermore, we show how information on the transition between both regimes (the IP 3 -binding and the Ca 2z dominated behaviors) can be extracted from the distribution of the number of IP 3 R's that open during a puff. Knowledge on this transition can, in turn, yield information on the largest distance over which Ca 2zmediated cooperativity acts and on the mean density of IP 3 -bound IP 3 R's of the clusters. In this way, we can estimate biophysical parameters that affect the intra-cluster dynamics from statistical information on the emergent collective behavior of the channels of the cluster.
The aim of the simple model that we introduce in this paper is to characterize the basic mechanisms that shape the distribution of the number of channels that open during puffs. In particular, we identify the competition between two stochastic processes as the main determinant of the form of the distribution. Therefore, an analysis of this form may give information on the relative weight of the two competing processes. The model does not include a detailed description of the dynamics that takes place during or between events. For some time most models of intracellular Ca 2z dynamics were deterministic (see e.g. [25]). The observation of local signals such as puffs led to the development of several models that included a stochastic description of Ca 2z release [14,15,26,27] or of the spatial location of the IP 3 R's [28]. It is currently clear that stochastic effects are not only relevant for local release events but are a fundamental aspect of the Ca 2z dynamics for the full range of observed signals, including waves [29][30][31][32]. More information on stochastic models of Ca 2z signals can be found in a recent focus issue on the subject [33]. Simulations of these stochastic dynamic models could be used to probe the main findings of the present paper.

The Model
We introduce here a simple model to describe the distribution of puff sizes that occur at sites with similar numbers of IP 3 R's. The model is simple in the sense that it does not include a detailed description of the dynamics of the individual channel openings and closings or of IP 3 or Ca 2z binding and unbinding. However, it does include the stochasticity associated to IP 3 binding and channel coupling via CICR. Given the estimates of [8], the model assumes that clusters occupy a fixed size region (more specifically, a circle of radius, R) and that N IP 3 R's are randomly distributed over the cluster region with uniform probability. Each IP 3 R of the cluster has a probability p of having IP 3 bound. N p is the random variable that represents the number of available IP 3 R's (i.e., of IP 3 R's with IP 3 bound) of the cluster before a puff starts. The distribution of this variable is a binomial of parameters N and p, that can be approximated by a Poisson distribution of parameter l~pN for N large and p small enough (for example, with N~20 and p~0:05 the absolute difference between both cumulative distributions is lower than 0:01 for each N p ). The model considers that if an IP 3 R with IP 3 bound becomes open and Ca 2z starts to flow through its pore all other IP 3 R's with IP 3 bound that are within a distance, r inf , of the open IP 3 R will also become open. These newly opened IP 3 R's in turn trigger the opening of new IP 3 R's with IP 3 bound that are within the distance, r inf , from an open one. This scheme triggers a cascade of openings that stops when there are no more available IP 3 R's within the radius of influence (i.e., the distance r inf ) of any open IP 3 R. This cascade determines the number, n, of channels that open during a puff. We call P A (N p ) the probability that there are N p available IP 3 R's (with IP 3 bound) in a cluster and P o (n=N p ) the conditional probability that n channels open during an event given that there are N p with IP 3 bound in the cluster. Given that we are interested in the distribution of event sizes, we only consider the situations for which 1ƒnƒN p . Therefore, we renormalize the probabilities so that P N Np~1 P A (N p )~1 and P Np n~1 P o (n=N p )~1. In this way, P A (N p ) is a binomial or a Poisson distribution divided by one minus the probability that there are not IP 3 R's with IP 3 bound in the cluster. Using these renormalized versions of P A (N p ) and P o (n=N p ), the probability, P n , of having a puff with n open channels is given by: is approximated by a Poisson distribution of parameter l~pN when N??. P n can be readily compared with distributions obtained from experimental observations as the one displayed in Fig. 4D of [1].

Factors That Shape the Distribution of Event Sizes
The two stochastic components of the model are evident in the expression of P n . P A (N p ) reflects the stochasticity of IP 3 binding and P o (n=N p ) the one due to inter-channel coupling via CICR. The relative weight of both factors on the resulting P n depends on the relationship between two typical lengthscales of the problem: the radius of influence, r inf (the maximum distance between channels at which CICR still works) and the mean distance between channels with IP 3 bound, D, which is a random variable that can be computed in terms of R and the number of IP 3 R's with IP 3 bound, N p , as: D can take values between 0 and R ffiffiffi p p =2. Closely related to D is the density of available IP 3 R's which is given by: The relationship between D and r inf determines the relative weight of both stochastic components on P n . In particular, if r inf =D is very large, the opening of any channel of the cluster will eventually lead to the opening of all other available channels. If such a situation holds for most events, then P o (n=N p )&d nNp and P n will mainly be determined by the stochastic component due to IP 3 binding, i.e., P n &P A (n). If, for most events, r inf =D is very small, then most of P n will be concentrated near n~1, regardless of how many available IP 3 R's there are in each realization. We will refer to both extreme behaviors as IP 3 or Ca 2z limited. Depending on the parameters of the model (R, r inf , N, and p), one or the other situation is favored. However, in many situations one or the other behavior is favored depending on the value of N p , i.e., on the realization. In those cases, the dominant stochastic component of P n depends on the value of n.
We first illustrate how the distribution, P n , varies with the number of IP 3 R's of the cluster, N, while all other parameters are fixed. As N increases, the most likely values that N p can take on also increase. This means that it is more probable to have more available IP 3 R's at any given instance. On the other hand, since the spatial dimensions of the cluster are unchanged (R is fixed) the mean distance between available IP 3 R's, D, is more likely to be smaller (see Eq. 2). Given that the typical distance for CICR to occur, r inf , is also fixed, it is more probable that r inf =D be larger. Therefore, P n approaches P A as N is increased. This is illustrated in Fig. 1 where we have plotted the distributions P n (solid circles) and P A (N p ) (bars) obtained with 1000 realizations of our model using p~5=18, R~250nm, r inf~2 30nm and three values of N. In A, N~10, the number of available channels is small for most realizations (its mean value is m~pN~2:8) so that P n is dominated by inter-channel Ca 2z -coupling and concentrated around small values of n. In C, N~50, the number of available channels is large for most realizations (its mean value is m~pN~13:9) so that D is typically smaller than r inf (D(13:9)~59:43nm). In this case, P n is dominated by the IP 3binding stochasticity and almost indistinguishable from the distribution of available channels, P A . The example of Fig. 1 B corresponds to a situation in between these two extreme cases with N~20. We can observe how, as the number of available channels is more likely to be larger, P n approaches P A . We also observe that for N~10 and N~20, P n and P A differ mainly in the region of small values of n. This occurs because it is difficult for one open channel to induce the opening of another one if the mean interchannel distance is large. Thus, if N p is small it is very rare that all available channels become open. In this way, the relative frequency of small events becomes larger than the fraction of instances with a small number of available channels.
A transition from Ca 2z -dominated to IP 3 -binding dominated stochasticity also occurs as r inf is increased, while all other parameters are fixed. In this case, P A remains unchanged and so does the mean distance between available IP 3 R's, D. By changing r inf it is possible to go from a situation in which r inf =D is small for most events and P n is Ca 2z -limited to a situation in which r inf =D is large and P n is IP 3 -binding limited. This is illustrated in Fig. 2 where we have plotted the distribution of event sizes that we obtain with our model for three different values of r inf . For r inf~1 50nm, the distribution is Ca 2z -coupling limited and is concentrated around n~1. As r inf is increased, the relative frequency of events with small n decreases. For r inf~5 00nm, the distribution is IP 3 -binding limited. In this example, P A is well approximated by a Poisson distribution of parameter l~5 (data not shown). The situation in between these extreme cases corresponds to r inf~2 50nm and is able to reproduce reasonably well the experimental distribution of Fig. 4D of [1] (superimposed with bars in Fig. 2).
In the Ca 2z limited behavior the number of open channels, n, is small for most events, regardless of the value of N p . This implies nvN p for almost all events. In the IP 3 -binding limited behavior all available IP 3 R's become open (n~N p in most cases). Therefore, in order to analyze the transition between the Ca 2z -dominated to IP 3 -binding dominated stochasticity, we study how often events occur for which all available IP 3 R's become open. This happens trivially for events with N p~1 . Here we are interested in situations with N p w1. To this end, we compute numerically the probability that all available IP 3 R's, N p , become open, P(n~N p =N p ), which is a function of N p and of only one independent parameter, the dimensionless radius of influence, r inf =R, (see Methods). We plot in Fig. 3 A P(n~N p =N p ) as a function of r inf =R, for N p~1 0 (circles), N p~3 0 (squares) and N p~1 00 (triangles). As expected, P(n~N p =N p ) is an increasing function of N p for each value of r inf =R. We also observe that P(n~N p =Np) is an increasing (sigmoidal-like) function of r inf =R that goes from 0 (i.e. nvN p in almost all cases, which corresponds to Ca 2z -dominated stochasticity) to 1 (i.e. n~N p in almost all cases, which corresponds to IP 3binding dominated stochasticity) and that such transition occurs over a smaller interval of r inf =R values the larger N p is.
We can think of the Ca 2z -limited and the IP 3 -binding limited situations as two phases and the transition between them as a phase transition in the limit of very large N p . This percolation-like transition occurs at a well defined value of r inf =R in this limit. For finite values of N p we introduce two quantities, r (1) inf (N p ) and r (2) inf (N p ), that determine the type of regime that we can expect (Ca 2z -limited if r inf vr (1) inf or IP 3 -binding limited if r inf wr (2) inf ) for each value of N p (see Methods). The arrows in Fig. 3 A indicate the values of r (1) inf =R (*) and r (2) inf =R (**) for the N p~1 0 case. We show in Fig. 3 B plots of r (1) inf =R, r (2) inf =R and D=R as functions of r A (Eq. (3)). It is important to note that these curves are the same, regardless of the specific parameter values of the model. We observe that all of them are decreasing functions of r A or, equivalently, of N p . N p is a stochastic variable that changes from realization to realization. Therefore, even for a given cluster (characterized by fixed values of N, p and r inf ) r (1) inf and r (2) inf may take on different values depending on the realization. In this way, depending on r inf and the values that r A may take on, a subset of the events that occur at a cluster may be IP 3 -binding limited (those for which r (2) inf vr inf ) while others are not. An analogous situation may hold regarding the Ca 2z -limited behavior. Furthermore, for some clusters, the Ca 2z -limited condition may hold for some events and the IP 3 -binding limited for others. If the parameters N, p, R and r inf are such that most realizations satisfy r (2) inf (r A )wr inf , then most events will be IP 3 -binding limited. This happens if r inf or pN are large enough, in which case the distribution of event sizes, P n , approaches the distribution of available channels, P A .

Observing Percolation as a Function of Event Size
The results of Fig. 3 B imply that there are clusters that can display different types of behaviors depending on the event. For these clusters, we expect to find, in their distribution, P n , a trace of the transition to the limiting behavior that they can display. Here we are interested in the percolation transition, i.e., the transition to the IP 3 -binding dominated stochasticity. As already discussed, the larger N p the more likely it is that all IP 3 R's become open during the puff (see Fig. 3 A). Thus, the transition to the IP 3 -binding dominated stochasticity should occur as N p and, consequently, n become larger. To study this transition we consider a cluster with fixed parameters R, r inf , p and N (or l~pN in the Poisson limit) and define n Ã as the minimum value of N p such that r inf §r (2) inf (N p ). The definition of r (2) inf is based on the conditional probability, P(n~N p =N p ), which is independent of N. In cases with finite N, n Ã is meaningful provided that it be smaller than N. Since r (2) inf decreases with N p (see Fig. 3 B), taking into account the definitions of n Ã and of r (2) inf (N p ) (see Methods) we conclude that r inf wr (2) inf (N p ) and P o (n~N p =N p ) §0:9 for all N p §n Ã . Thus, we can approximate: with less than 10% error. Inserting this approximation in Eq. (1) we obtain: We then conclude that the n §n Ã tail of P n corresponds to IP 3binding dominated events. Therefore, it should be possible to approximate it by a (renormalized) binomial (provided that n Ã ƒN) or Poisson distribution in the region of large n. The left border of this IP 3 dominated behavior, n Ã , gives information on r inf , i.e. on the maximum distance for which CICR-coupling can work effectively. Therefore, it should be possible to estimate r inf by analyzing P n , i.e., to infer a biophysical parameter that characterizes the intra-cluster dynamics from statistical information on the emergent collective behavior of the channels of the cluster.

Determining Intra-Cluster Properties from Observations of the Cluster as a Whole
We now discuss how we can estimate r inf from an experimental distribution of event sizes, P n . For the sake of simplicity, we assume that P A (n) can be approximated by a renormalized Poisson distribution, P l (n):l n exp ({l)=((1{ exp ({l))n!), of unknown parameter l. The goal of this section is to provide a way to estimate l and n Ã , the value of n at which P n and P A (n)&P l (n) depart from one another (see Eq. (6)). Once n Ã is inferred, we estimate r inf =R as r inf =R&r (2) inf (n Ã )=R using the function displayed in Fig. 3 B. To this end, we focus on the large n tails of P n and P l by computing the complementary cumulative distribution functions: for ' §2. Given that P l is proportional to a Poisson distribution, there is an analytic expression for T l . Namely, T l (')1 {C(½',l)=((1{ exp ({l))(½'{1)!), where C(x,y) is the incomplete C function and ½' is the integer part of '. If the cluster is such that n Ã exists so that r inf is larger than r (2) inf (N p ) for N p §n Ã and it is smaller otherwise, then, according to the calculation of the previous section, P n &P A (n)&P l (n) for n §n Ã . Therefore, the complementary cumulative distribution functions of Eqs. (7)-(8) also satisfy T o (')&T l (') for ' §n Ã .
We now describe how to estimate n Ã and l. The aim is to obtain a (renormalized) Poisson distribution, P l that can approximate P n in the large n region. If we find it, we assume that it is a good approximation of the distribution of available channels, P A . As illustrated in Fig. 1, the mean value, SnT obs that is obtained using the experimental distribution, P n , is smaller than the one that would be obtained if P A (n) was used instead. On the other hand, if P A is a good approximation of P n in the large n region, then the mean value obtained with P A should be smaller than the size of the largest observed event, n max . This implies that  (n~N p =Np), as a function of the dimensionless radius of influence, r inf =R, for N p~1 0 (circles), N p~3 0 (squares) and N p~1 00 (triangles). B: r (1) inf =R (circles), r (2) inf =R (squares) and D=R (triangles) as functions of r A (N p )~N p =pR 2 . The values of r (1) inf =R and r (2) inf =R for the case with N p~1 0 are indicated in A with one and two asterisks, respectively. doi:10.1371/journal.pone.0008997.g003 if P A can be approximated by a renormalized Poisson distribution of parameter l. Therefore, we look for the best P l within a finite set of renormalized Poisson distributions of parameters l i satisfying (9). In order to estimate r inf from the observations the relevant quantity that we need to obtain is n Ã , which is an integer. For this purpose, it is possible to use a rather coarse grid of l values within the interval defined in (9). In particular, we have mainly used integer values of l obtaining good results. Once the values l i are chosen, we compute the complementary cumulative distribution functions, T li (k) given by (8) for each l i and 1ƒkƒn max . We then calculate the error of approximating T o (k) by T li (k) over the interval 'ƒkƒn max as a function of ': We set a threshold for the error, e, and choose n Ã i for each l i as the smallest value of ' for which error(T o ,T li )ƒe. Finally, we choose the best l i as the one with the smallest n Ã .
The procedure is illustrated in Fig. 4 where the ''experimental'' distribution comes from a simulation of our model with R~250nm, p~5=18, r inf~2 30nm and N~18. In this case, n max~1 0. We show in Fig. 4 A the complementary cumulative distribution functions and in Fig. 4 B the errors for the values of l that we have considered: l 1~3 (inverted triangles), l 2~4 (triangles), l 3~5 (squares) and l 4~6 (rhombes). Larger values of l give very bad approximations and are not shown. We show in Fig. 4 C the values, n Ã , obtained for each l using the threshold, e~0:005 (shown with a horizontal line in Fig. 4 B). In this example, the best value is l~4 for which n Ã~8 . We estimate the density of IP 3 -bound IP 3 R's at which the departure between the experimental and the Poisson distribution occurs as r Ã A~n Ã =(pR 2 )&41mm {2 , where we have used R~250nm.
Using the r (2) inf =R vs r A relationship displayed in Fig. 3 B, we estimate r (2) inf (r Ã A )=R&0:9+0:08 from which we get r (2) inf (r Ã A )& 220nm+20nm. This provides an estimate of the radius of influence which compares very well with the value that was used to generate the data, r inf~2 30nm. Using the same procedure, we analyzed the data presented in Fig. 4D of [1] and obtained r inf~2 20nm+20nm assuming R~250nm.

Discussion
Intracellular Ca 2z signals are built from localized release events in which Ca 2z enters the cytosol through one or several channels. Ca 2z release from the endoplasmic reticulum through IP 3 R's is a key component of the Ca 2z signaling toolkit in many cell types. IP 3 R's are Ca 2z channels that need to bind IP 3 and Ca 2z to become open and are usually organized in clusters on the membrane of the endoplasmic reticulum. The intra-cluster organization and the interactions of the channels within it affect the dynamics and extent of the signals. Therefore, their study is a matter of active research.
Recent experiments [1] that use super-resolution optical techniques are providing detailed data on elementary IP 3 Rmediated Ca 2z release events in mammalian cells. In the experiments, the number of IP 3 R-Ca 2z -channels that open during each event can be inferred from the observed puff amplitudes without much processing. The observations of [1] showed that the variability among clusters affected the shape of the event size distribution, P n . In order to get rid of this variability, the distribution coming from sites with similar properties was computed in [1]. The distribution, P n , obtained in this way was not Poisson, as we might have expected if the number of channels that opened during each event was proportional to the number of IP 3 R's with IP 3 bound in the cluster [8]. The authors of [1] could reproduce P n approximately (for events larger than a certain size) assuming a weak cooperativity among channels. Namely, they assumed that the probability that a channel became open scaled as some power of the number of open channels and obtained that the exponent was 1/3 from a fit to the data. The rationale for the cooperativity assumption relied on the fact that the IP 3 R's of a cluster may be coupled via CICR induced by the Ca 2z ions that travel from an open IP 3 R to a neighboring one. The model of [1], however, did not take space into account and did not provide a mechanistic explanation for the obtained scaling.
In this paper we have presented a simple model that includes a description of the intra-cluster spatial organization with which we can reproduce the observed distribution over all event sizes. In the model the distribution, P n , is the result of the competition of two stochastic processes: IP 3 binding and distance-dependent CICR. The model assumes a stationarity condition, namely, that the agonists concentration at the release site is the same immediately before the occurrence of each puff. This condition holds as long as puffs are independent of one another. This is consistent with the observations reported in [1] where cluster coupling was prevented using the slow Ca 2z buffer EGTA and where IP 3 R-Ca 2zinhibition does not play a significant role. In any case, our model is adequate to describe the distribution of first event sizes that occur at each cluster before Ca 2z can exert any inhibiting effect.
There are two limiting cases in which one of the two stochastic processes considered in the model is the main determinant of the distribution shape. If the mean distance between IP 3 R's with IP 3 bound in the cluster is much smaller than the typical distance of inter-channel coupling due to CICR for most events, the distribution is IP 3 -binding limited and it can be approximated by a binomial or Poisson distribution. In the opposite case, CICR dominates and the distribution is peaked around n&1. The Ca 2zlimited and the IP 3 -binding limited situations can be thought of as two phases and the transition between them as a percolation-like transition in the limit in which the number of IP 3 R's with IP 3 bound, N p , is very large. This interpretation of the factors that shape the observed distribution can be tested with simulations of some of the stochastic models of intracellular Ca 2z signals reported in the literature (see e.g. [33] and references therein). They can also be tested experimentally. One possibility is to change the most likely values of N p by changing the amount of IP 3 that is photo-released in the cell. An alternative option is to analyze P n for events coming from clusters that give rise, on average, to larger events. According to the model, the distribution should approach a binomial or Poisson distribution as the mean value of N p becomes larger while other parameters remain the same. Another way to affect the balance between both stochastic components is to disrupt Ca 2z -mediated inter-channel coupling by means of a fast buffer such as BAPTA.
Given that N p is a stochastic variable that varies from event to event, the transition between the Ca 2z -dominated and IP 3binding dominated stochasticity described by the model may be reflected in the way that P n depends on the event size, n. In fact, we have used this property to show how a fingerprint of this transition may be encountered in P n and how information on the inter-channel coupling distance may be extracted from it. This means that a parameter that characterizes the communication between pairs of channels can be estimated from statistical information on the emergent collective behavior of the channels of the cluster. This information could be used to analyze the effect of buffers on the intra-cluster dynamics, a matter that is of active current research [19,34]. Our model provides a simple tool with which this effect can be analyzed in experiments.

Methods
Each term of the sum that defines Eq. (1) is the product of two functions. We have an analytic expression for one of them, P A (N p ), but not for P o (n=N p ). Thus, we compute P n numerically performing realizations of the model with fixed values of N, p, R and r inf . The location of the channels within the cluster and which of them have IP 3 bound vary among realizations and are chosen randomly (see Results). We only keep realizations with N p §1. Once we have the spatial distribution of available IP 3 R's, we start each event by picking at random one of the IP 3 R's with IP 3 bound and assume it is open. If N p~1 , we assume it gives rise to an event with n~1. By changing the values of N, p, R and r inf we analyze how P n varies with them. In this way we can determine the values of the parameters that best reproduce the experimental observations. r inf could be measured in units of the cluster spatial extent, R, in which case we would get rid of one parameter of the problem, R. We keep it to make a connection with the experimental data. However, it is important to note that the number of independent parameters of the model is 3, for finite N and 2 in the limit in which P A can be approximated by a Poisson distribution.
For each value, N p , of available IP 3 R's, we estimate the fraction of events such that the N p IP 3 R's become open. This fraction is one for N p~1 . For N p w1, we compute the probability that all available IP 3 R's become open, P(n~N p =N p ), numerically, performing 500 stochastic realizations of our model for each of which we fix the value of N p a priori. Namely, we fix at the beginning the values of R, r inf and N p and then pick N p locations at random over the circle where we assume there are available IP 3 R's. From there on, the model goes on as before, generating the cascade of openings that determines n. The distribution of events with n open channels for each value of N p gives P(n~N p =N p ). This function of N p depends on only one independent parameter, r inf =R. As expected, it is an increasing function of r inf =R (see Fig. 3 A).
We define two quantities, r (2) inf and r (1) inf , which are values of r inf for which P(n~N p =N p ) is either close to 1 or to 0, respectively. We compute them as follows. We first calculate a lower bound for r (2) inf as the minimum value of r inf such that, if r inf is larger than this lower bound, then P o (n~N p =N p )w0:8. We calculate an upper bound for r (2) inf as the minimum value of r inf for which P(n~N p =N p )~1. Then, we compute r (2) inf as the mean between these two bounds. We assume that the distance between the bounds is the error with which r (2) inf can be determined. We proceed analogously in the case of r (1) inf , but in this case the lower bound is the largest value of r inf for which P(n~N p =N p )~0 and the upper bound is the maximum value of r inf for which P(n~N p =N p )ƒ0:2. We compute r (1) inf and r (2) inf in this way using the numerical estimations of P(n~N p =N p ) for various values of r inf .