Intracellular calcium signals display an avalanche-like behavior over multiple lengthscales

Many natural phenomena display “self-organized criticality” (SOC), (Bak et al., 1987). This refers to spatially extended systems for which patterns of activity characterized by different lengthscales can occur with a probability density that follows a power law with pattern size. Differently from power laws at phase transitions, systems displaying SOC do not need the tuning of an external parameter. Here we analyze intracellular calcium (Ca2+) signals, a key component of the signaling toolkit of almost any cell type. Ca2+ signals can either be spatially restricted (local) or propagate throughout the cell (global). Different models have suggested that the transition from local to global signals is similar to that of directed percolation. Directed percolation has been associated, in turn, to the appearance of SOC. In this paper we discuss these issues within the framework of simple models of Ca2+ signal propagation. We also analyze the size distribution of local signals (“puffs”) observed in immature Xenopus Laevis oocytes. The puff amplitude distribution obtained from observed local signals is not Gaussian with a noticeable fraction of large size events. The experimental distribution of puff areas in the spatio-temporal record of the image has a long tail that is approximately log-normal. The distribution can also be fitted with a power law relationship albeit with a smaller goodness of fit. The power law behavior is encountered within a simple model that includes some coupling among individual signals for a wide range of parameter values. An analysis of the model shows that a global elevation of the Ca2+ concentration plays a major role in determining whether the puff size distribution is long-tailed or not. This suggests that Ca2+-clearing from the cytosol is key to determine whether IP3-mediated Ca2+ signals can display a SOC-like behavior or not.


INTRODUCTION
Ca 2+ signals are involved in a large variety of physiological processes (Berridge et al., 1998). The diversity of spatio-temporal behaviors that the Ca 2+ concentration can display is fundamental for the universality of the signals (Berridge and Dupont, 1994;Berridge et al., 2000). The prolonged permanence of high free Ca 2+ concentrations in the cytosol leads to cell death. For this reason, the signals are built upon local Ca 2+ elevations due to Ca 2+ entry from the extracellular medium or internal stores through specialized channels. Depending on different factors, these local elevations can remain local or induce the opening of other Ca 2+ channels giving rise to signals that propagate over larger distances (Bootman et al., 1997;Callamaras et al., 1998;Sun et al., 1998). Ca 2+ release from the endoplasmic reticulum (ER) through inositol 1,4,5-trisphosphate (IP 3 ) receptors (IP 3 R's) is a key component that participates of both local (e.g., "puffs") and global (waves) Ca 2+ signals (Berridge, 1993;Choe and Ehrlich, 2006). There is an extensive literature on the observation of this type of signals which have led to the picture that IP 3 R's are organized in clusters on the membrane of the ER . Although the spatial distribution of IP 3 R's is currently under debate (see e.g., Taufiq-Ur-Rahman et al., 2009), it is apparent that those that are activatable (i.e., can become open) are indeed distributed inhomogeneously . Therefore, Ca 2+ signals are characterized by very diverse lengthscales: the typical size of the channel pore (∼1 nm); the dimensions of the IP 3 R cluster (∼100 nm); the extent of the cell (10 μm−1 nm). Whether the signal remains local or propagates throughout the cell depends on Ca 2+ itself, since IP 3 R's become open upon Ca 2+ and IP 3 R binding. Ca 2+ , in turn, plays a dual role, inducing channel opening at moderate concentrations and channel closing at higher ones (Foskett et al., 2007). When a channel becomes open, Ca 2+ rapidly flows in and diffuses inside the cytosol. Diffusion is not free given that the cells have several mechanisms (among them, binding to immobile or slowly moving "buffers") to reduce Ca 2+ concentration. In any case, the local elevation of [Ca 2+ ] that originates at an open channel eventually spreads out in space inducing new channel openings if it finds IP 3 -bound IP 3 R's on its way. This Calcium-Induced-Calcium-Release (CICR) is basic for Ca 2+ signal propagation.
The fact that IP 3 R-clusters are not evenly distributed inside cells introduce a new factor that also modulates the ability of a signal to propagate or not: the mean inter-cluster distance, d. In Keizer et al. (1998); Pearson and Ponce-Dawson (1998); Ponce Dawson et al. (1999), we introduced the simple deterministic "fire-diffuse-fire" (fdf) model to analyze the effect of this additional variable. In particular we found different wave front regimes depending on the ratio between the typical time a cluster remains active (releasing Ca 2+ ) and the time it takes for Ca 2+ to diffuse over the distance d. The existence of an inter-cluster separation can also lead to propagation failure even if some clusters are able to release Ca 2+ . In Pearson and Ponce-Dawson (1998), we analyzed propagation failure within the framework of the deterministic fdf model. We found that the dynamics of the wave became irregular with its temporal pattern undergoing a series of period-doubling bifurcations until it eventually disappeared as an attracting solution of the model. Propagation failure was also studied within the framework of stochastic fdf models (Coombes and Timofeeva, 2003;Calabrese et al., 2010). These stochastic models reproduce better the experimental observations in which, after a series of localized signals, there is a burst of activity that propagates over a much larger spatial region. The relevance of the stochastic localized activity to elicit a global signal has been recently pointed out. In particular, it has been argued that apparently periodic global signals are not truly periodic but are the consequence of stochastic channel openings (Skupin et al., 2008(Skupin et al., , 2010Thurley et al., 2011). In fact, it has long been reported that the occurrence of several local signals can eventually elicit a Ca 2+ wave (Marchant and Parker, 2001). This transition from local to global Ca 2+ signal propagation ressembles that of percolation (Stauffer and Aharony, 1992), a concept that, although originally referred to fluid displacement in a disordered medium, is now used to describe the behavior of connected clusters of any sort. In the case of the Ca 2+ signals that we are focusing on in this paper, the latter would correspond to (active) IP 3 R clusters. In systems where the dynamics involves some type of propagation that could evnetually fail, as in the case of Ca 2+ signals, the relevant concept is that of directed percolation. In particular, it has been observed that the stochastic fdf model of (Coombes and Timofeeva, 2003) is in the universality class of directed percolation (Timofeeva and Coombes, 2004). This means that the correlation lengthscale and other properties of the system scale as |p − p o | α with specific values of α for p ≈ p o , where p is a parameter that controls the transition from a macroscopically active (percolating) to an inactive (non-percolating) state (that occur, e.g., for p > p o and p < p o , respectively). The exponent, α, depends on the quantity whose dependence with p is being analyzed. The values that these exponents take on depend on the type of transition that takes place at p = p o independently of the specificities of the system under study (i.e., they are "universal"). The inhomogenous spatial distribution of IP 3 R's also shapes the kinetics of local signals such as puffs (Diambra and Marchant, 2011). In fact, we investigated the "percolation" of Ca 2+ signals inside an IP 3 R cluster rather than between clusters in (Solovey and Ponce Dawson, 2010). In that paper we called percolation the transition between a situation for which most signals involved the opening of all IP 3 R's with IP 3 of the cluster (percolating states) to a situation for which very few IP 3 R's of the cluster opened (non-percolating states). To that end we introduced a very simple model with which we were able to reproduce the experimentally observed puff size distribution reported in . In particular, by fitting the model to the observations, we determined that, for one open IP 3 R, the radius around which CICR is effective is of the order of the cluster radius (∼250 nm). This means a very strong inter-channel coupling due to CICR inside the cluster, but not strong enough so that all IP 3 R's with IP 3 -bound would become open during each puff. Actually, the experimental data seemed to be near the transition between the percolating (IP 3 -dominated) and nonpercolating (for which Ca 2+ is limiting) states. In this paper we investigate this aspect further by analyzing both new experimental data and a simple model that describes inter-cluster coupling. Again in this case the experiments seem to be close to a transition point.
One could argue that finding a natural system at a transition point is non-generic. However, the existence of systems that behave as if they are naturally near a critical point has been widely observed (Bak, 1996). This is the main idea behind "Self-Organized Criticality" (SOC) (Bak et al., 1987(Bak et al., , 1988, a property of spatially extended dynamical systems that have a critical point as an attractor. In systems displaying SOC, the state is usually identified to be critical when the lengthscale distribution is scale invariant or follows a power law. The paradigmatic example of a system that displays SOC is the sand pile in which avalanches of all sizes can occur as sand grains are added slowly compared to the typical time of energy dissipation. The existence of avalanches of all possible sizes is similar to the occurrence of Ca 2+ signals that propagate over different regions and eventually fail to advance. In fact, noise-induced spiral Ca 2+ waves in astrocytes (Jung et al., 1998) and arrhythmogenic Ca 2+ waves in cardiac myocytes  have been shown to display features of SOC. SOC has been argued to occur in several other biological systems (Bak and Sneppen, 1993;Fraiman et al., 2009;Tetzlaff et al., 2010). Directed percolation and SOC have been associated in the literature (Jovanović et al., 1994). For example, a modified version of a model that displays SOC was shown to evolve to a directed percolation critical point (Sornette and Dornic, 1996). In particular, a SOC-like behavior (characterized by a scale invariant distribution) can be reached in systems with a phase transition onto an absorbing (inactive) state that are driven slowly by infinitesimal additions of the quantity that brings the system close to the transition and are allowed to relax upon activation by decreasing this quantity Vespignani et al., 2000). Although this "SOC recipe" is still debatable Pruessner and Peters (2006) it could be the mechanism that underlies the scaleinvariant behaviors observed in Ca 2+ signals as we discuss in this paper.

OOCYTE PREPARATION
Experiments were performed on immature oocytes from Xenopus laevis previously treated with collagenase in order to defolliculate them. Oocytes were loaded by intracellular microinjection 1 h before recording with Fluo-4 (fluorescent calcium indicator), caged InsP 3 (D-Myo-Inositol 1,4,5-Triphosphate, P4(5)-(1-(2-Nitrophenyl)ethyl) Ester) and the exogenous calcium buffer EGTA (Ethylene glycol-bis(2aminoethylether)-N,N,N ,N -tetraacetic acid). Final intracellular concentrations of Fluo-4 and InsP 3 were 36 μM and 9 μM respectively, assuming 1 μl cytosolic volume. The final concentration of EGTA was varied from 0 μl to 90 μM. Recordings were made at room temperature. All recordings were obtained with the scan line focused at the depth of the pigment granules in the animal hemisphere of the oocyte. Fluo-4 and InsP 3 were from Molecular Probes Inc.; EGTA was from Sigma Aldrich.

CONFOCAL MICROSCOPY
Confocal linescan imaging was performed using a spectral confocal scanning microscope Olympus FluoView1000 that has a spectral scan unit connected to an inverted microscope IX81. Flash photolysis of the caged compound was made with a mercury lamp that comes with the microscope using the modification introduced in Barella et al., 2011. The ultraviolet illumination from the lamp (350-400 nm) enters through an optical fiber that is coupled to a non-conventional port of the microscope in order to perform the photolysis while simultaneously acquiring confocal fluorescence images. The fluorescent indicator was excited using the 488 nm line of a multiline Argon laser that was focused on the oocyte with a ×60 oil immersion objective (NA 1.35). The emitted fluoresecence was detected in the 500-600 nm range. The UV flash was delivered approximately 2 s after having started the image acquisition.

IMAGE ANALYSIS
All images were analyzed using routines written in MATLAB. Images were smoothed out with a filter prior to this analysis. To determine the amplitude, A, of a puff, the mean basal fluorescence at each spatial point, F 0 , was computed before delivery of the UV flash. The region occupied by a puff was identified manually. The puff amplitude, A, was then defined as the maximum value of F/F 0 (F − F 0 /F 0 ) in the region, with F the fluorescence value at each point. We show in Figure 1A an example of an experimental linescan image with time and space in the horizontal and vertical axes, respectively, and using a color code to display F/F 0 that goes from colder to warmer colors with increasing fluorescence. In this example two puffs appear at the same release site but at different times, the first one being more intense than the second one. In order to evaluate the area occupied by a puff in the spatio-temporal record provided by the image, a threshold criterion is applied. The mean fluorescence F m and its standard deviation σ are calculated for the whole ( F/F 0 ) image. Then all pixels with F/F 0 lower than 2F m + 10 σ are set equal to 0 and all pixels with fluorescence above that value are set equal to 1. After appliying the threshold, a binary image is obtained for wich the area of the "1 pixels" is calculated. We show in Figure 1B the thresholded image obtained from Figure 1A using this criterion with the "1 pixels" displayed in black and the rest in white. In this way we determine the area occupied by a puff in the spatio-temporal record provided by the image.
In order to verify that photobleaching was not affecting our data we grouped together puffs that occurred within the same 1 slong interval along the experimental record and computed the mean and standard deviation of their amplitude. We show the results in Figure 2 where we can observe that there are no significant changes as time progresses. In this way, we are confident that photobleaching is not affecting the statistics that we are drawing from the experiments.

DATA FITTING
Some of the data analysed were fitted with different distributions.
To evaluate the goodness of those fits, we used the coefficient of determination (R 2 ) defined by 1 − R ss /T ss where R ss is the residual sum of squares and T ss is the total sum of squares.

Intracluster dynamics
In order to obtain the distribution of puff sizes that occur at a given cluster we follow the model introduced in Solovey and Ponce Dawson (2010)

Intercluster dynamics without calcium accumulation
In order to obtain the distribution of event sizes that involve the opening of IP 3 R's in many different clusters inside the cell we use a modified version of the model just described. This modified version assumes that the separation between two channels in a cluster may be neglected when computing their contribution to the free Ca 2+ concentration at the location of any other cluster and that all the channels of those other clusters sense the same Ca 2+ concentration at any given time. This assumption relies on the expected difference between lateral cluster size (∼200 nm) and inter-cluster distance (∼2 μm). Under this assumption, taking into account that inter-channel coupling is due to the Ca 2+ ions that flow through open IP 3 R's and that Ca 2+ diffuses with an effective coefficient between neighboring cluster sites we assume that the contribution of a cluster with N oi open channels located at position r i to the Ca 2+ concentration at another position, r j , is proportional to N oi /|r i − r j |. Then, since at the intra-cluster level we assume that all available IP 3 R's within a distance r inf of one open channel become open, at the level of how clusters activate one another we assume that a cluster with N o open channels induces the "opening" of all other clusters with available channels that are within a distance r inf N o of it. We also assume that, when a cluster is "open" (or active), all its available channels are open (i.e., releasing Ca 2+ ). In this scheme we do not add up the Ca 2+ concentration coming from several open clusters to decide the opening of a new cluster. This simplification is valid as long as the influence of the closest open cluster is more important than that of the more distant ones. The simulations that we show in this paper were done considering a circular total area of radius R = 5 μm and r inf = 0.250 μm as before. For each realization we selected the number of clusters from a Poisson distribution with mean 25 to guarantee a mean inter-cluster distance of ∼2 μm and then determined the number of available channels for each cluster from a Poisson distribution with mean 8. In this way, the fraction of clusters with available channels is similar to the fraction of available to total channels in the case of the intra-cluster simulations.

Intercluster dynamics with calcium accumulation
The previous model of the intercluster dynamics treats every signal as independent of one another. There is evidence, however, that several local signals usually precede a Ca 2+ (global) wave (Marchant and Parker, 2001). The previous model is then modified to take into account this "coupling" among signals.
To that end a larger region (of radius R = 20 μm) is considered where a certain number of clusters was randomly distributed with uniform probability over the area. The number of clusters is drawn from a Poisson distribution with mean 400 in order to guarantee a mean inter-cluster distance of ∼2 μm and the number of IP 3 -bound IP 3 R's per cluster is drawn from a Poisson distribution of mean 8, as in the previous model. These numbers do not change during the simulation. In order to include signal coupling, we do a time-dependent simulation in which signals are assumed to occur instantaneously and the clusters that participate of it are determined as in the previous intercluster model with some minor modifications. In between signals, the number of active IP 3 -bound IP 3 R's and the Ca 2+ concentration change with time. Thus, there is an implicit time-scale separation with a fast dynamics during a burst of activity and a slower recovery between such bursts. On the other hand, in this new model we distinguish between active IP 3 -bound and inactive IP 3 -bound IP 3 R's. In particular, we assume that only active IP 3 -bound IP 3 R's can become open during a signal and that all those that participate of an event become inactive immediately after it. Each inactive IP 3 R remains inactive for a certain (random) time. Following (Fraiman et al., 2006) we draw the inhibition time of every inactive IP 3 R from an exponential distribution with mean 2.5 s. ] is close to its basal value at a distance of the order of 1 μm. In this new model the initiation of a signal not only depends on the number of active IP 3 -bound IP 3 R's but also on [Ca 2+ ]. In particular, we assume that the probability per unit time that a cluster with N active IP 3 R's become open starting a new signal is 0.225s −1 * N active * [Ca 2+ ]/Ca basal . In this way, this probability takes on the value estimated in (Fraiman et al., 2006) for puff initiation in the presence of basal Ca 2+ . Once the first cluster becomes open, the opening of the other clusters that participate of the signal is decided as before but taking into account the background [Ca 2+ ] that is already accumulated. In particular, an open cluster, k, with N k active channels induces the opening of all other clusters, i, with N i > 0 active channels that satisfy N k 0.414 μMμm/d ik + [Ca 2+ ] ≥ 0.414 μMμm/r inf + [Ca basal ] where d ik is the distance between the i-th and the k-th cluster, [Ca 2+ ] is the (accumulated) background Ca 2+ at the time of the signal and r inf = 0.250 μm.

INTRACLUSTER AND INTERCLUSTER PERCOLATION USING A SIMPLE MODEL
In this section we show the results of using the intra and intercluster Models without Ca 2+ accumulation described in the "Materials and Methods" section to generate distributions of Ca 2+ signal sizes. As explained in that section, one of these models is used for local signals (puffs) that involve the opening of IP 3 R's within a cluster of ∼400 nm sides. The other one models signals that can propagate between clusters that are ∼2 μm apart from one another. Both models include a simplified description of the two main dynamical processes that shape the signals: IP 3 and Ca 2+ binding, treating each signal as independent of one another.
In another section we describe the results obtained when Ca 2+mediated signal coupling is included in the intercluster model. Signal types strongly depend on IP 3 and Ca 2+ since IP 3 R's need to bind these two species to become open. In the experiments with which we compare the models the cells are loaded with caged IP 3 that is eventually released by means of a UV flash. In this way, the experimentalist has control over the IP 3 distribution. Therefore, the models assume that [IP 3 ] is fixed. The dynamics of Ca 2+ , on the other hand, is entirely controlled by the cells. Thus, its concentration is very low at the beginning of the experiment and only after a latency a Ca 2+ ion encounters an IP 3 -bound IP 3 R inducing its opening. Once that happens, the released Ca 2+ diffuses away from the open channel and can induce the opening of nearby IP 3 R's with IP 3 bound. The change in [Ca 2+ ] due to one open channel is noticeable only in its vicinity. Therefore, the Ca 2+ -mediated inter-channel coupling depends on the distance. This feature is taken into account in both models in a simplified way. In the local model, we assume that each channel opens a nearby IP 3 -bound IP 3 R if it is within a fixed distance, r inf . In the model of global signals we assume that each cluster that is releasing Ca 2+ through N o open IP 3 R's can induce Ca 2+ release from a nearby cluster if the latter is within a distance N o r inf of the former. We show in Figure 3A the puff size distribution obtained with 1000 realizations of the local model using the simulation parameters that allowed us to fit the experimental observations of . The distribution is relatively flat for small size events and approaches a Poisson distribution for larger events. This behavior is intermediate between two limiting cases: an IP 3 dominated situation in which all the channels of the cluster that have IP 3 bound when the event starts participate of the signal. This occurs when IP 3 -bound IP 3 R's are densely packed in the cluster, in which case the size distribution is Poisson. The other limiting case corresponds to events in which IP 3 -bound IP 3 R's are so far away from one another that the opening of one them is not likely to induce the opening of any other available IP 3 R in the cluster. In this case Ca 2+ -binding is the limiting factor for the propagation of the signal. Although there is not a clear phase transition, in Solovey and Ponce Dawson (2010) we referred to the change from the Ca 2+ -binding limited cases to the IP 3 -bound limited cases as percolation since it implies a transition from a situation in which very few channels open to another in which all IP 3 -bound IP 3 R's become open. It is interesting to note that the parameters that best fit the experimental observations are more or less at the border between both limiting cases.
We show in Figure 3B the event size distribution obtained with 1000 realizations of the intercluster coupling model without Ca 2+ accumulation. The rationale for the parameter values chosen for these simulations is as follows. Given that the simulation covers a ∼100 μm 2 area we chose a mean number of clusters equal to 25 to guarantee a mean inter-cluster separation of around 2 μm. Then, we chose the mean number of IP 3 -bound IP 3 R's in a cluster equal to 8. This is slightly larger than the mean value derived from the observations of Smith and Parker (2009) (6.02 ± 0.48) which correspond to experiments that focused on localized signals. Thus, by choosing a slightly larger value we expect to generate signals that involve the activation of several clusters. Since we assume that all available (i.e., IP 3 -bound) channels of a cluster become open upon activation of the cluster, the mean numbers of IP 3 -bound IP 3 R's and of open channels during a release event, N o , in a cluster coincide. Having N o = 8 gives N o r inf = 2 μm for r inf = 0.250 μm as in the intracluster model. For each open cluster, N o r inf is to be compared with the distance to all the other clusters to decide on their activation. Therefore, the parameters are such that N o r inf is, on average, equal to the mean inter-cluster distance. This is similar to the situation encountered for the intracluster dynamics model when using the parameters that gave the best fits to the experiments of ). Namely, the mean distance between available channels in the intra-cluster case (∼0.12 nm), and the one between clusters, in the inter-cluster case (∼2 μm), are of the same order as the numbers they need to be compared with to decide channel opening or cluster activation (r inf = 0.250 μm and N o r inf , respectively). In fact, we observe that the event size distributions of Figures 3A,B have similar shapes. As in the case of the intracluster dynamics, these parameters of the intercluster model roughly correspond to a mid situation between a Ca 2+ -binding dominated and an IP 3 -binding dominated distribution. This is illustrated in Figure 4 where we display the event size distributions obtained with the same parameters as in Figure 3B but with different values for the mean number of IP 3 -bound IP 3 R's in a cluster. In Figure 4A this value is 4 and intercluster coupling rarely occurs. In Figure 4B the value is 20 and almost all clusters with available IP 3 R's participate of each event.

RESULTS FROM EXPERIMENTS
In this section we present the experimental results obtained after analizing puffs from 13 different cells as described in the "Materials and Methods" section. In all the experiments 3 compounds (Fluo4, caged IP 3 , EGTA) were injected in the oocytes. The photorelease of the caged compound was then induced with a UV flash and the resulting fluorescence was measured. All images were obtained in the linescan mode which gives information on the fluorescence along a single line inside the oocyte. A variety of events were elicited in the experiments ranging from small puffs to waves. We only analyzed events that were localized in space and time which limits our experimental distributions. Puff sizes were characterized by their amplitude and (spatio-temporal) area which were computed as described in the "Materials and Methods" section. We present in this section the distributions of these two quantities derived from the experiments.
We show in Figure 5A the puff amplitude distribution computed over 202 puffs obtained from experiments performed in 13 cells. Puff amplitudes varied from 0.54 to 3.95. The distribution shows a single maximum at A = 1.5. We show in Figure 5B the same distribution but on a logarithmic scale. A power law fitting of these data is not good.
We show in Figure 6A the corresponding distribution of puff areas (spatio-temporal spread). Reliable area values could only be obtained for 179 puffs observed in 12 oocytes. Puff areas were estimated as explained in the "Materials and Methods" section. The same distribution is shown on a logarithmic scale in Figure 6B. In this case, the distribution can be approximated by a power law of exponent between 1.7 and 2.0 depending on the range of areas included in the fitting. The goodness of the fit was estimated by the coefficient of determination R 2 , which varied from 0.77 to 0.81 depending on the range of areas included. We compared this value with the ones obtained by fitting the data with other distributions. The goodness was worse in all cases with the exception of the log-normal distribution where we obtained R 2 = 0.90. It is always difficult to distinguish power law and lognormal distributions (Clauset et al., 2009) and in this case it is even more difficult given the small range of events spanned by the distribution. In either case, a key feature of both types of distributions is the existence of a long tail. As discussed later, in the case of calcium signals, this is a relevant property for its possible physiological implications.

CALCIUM ACCUMULATION AND IP 3 R INHIBITION CAN LEAD TO A SIZE DISTRIBUTION THAT CAN BE FITTED BY A POWER LAW
In the previous section experimental signal sizes of mainly local signals were characterized by two quantities: their amplitude and their spatio-temporal spread. These quantities not only depend on the number of IP 3 R's that open during the signal but also on the times of the openings and their spatial distribution (Solovey et al., 2011;Diambra and Marchant, 2011). In any case, they can be compared with the distributions of channels that participate of the events obtained from simulations of the different models. In particular, neither the models nor the experimental puff amplitude distributions can be well described by a power law. The experimental distribution of puff areas, however, could be fitted by a log-normal distribution or a power law, albeit over a relatively short range of lengthscales. The limited range of lengthscales could be attributed to the fact that we have mainly analyzed   By fitting a power law relationship to some of the data points we obtain exponents between 1.8 and 2.0 with R 2 between 0.77 and 0.81 (if we include data between 800 and 4500 px or between 800 and 4000 px, respectively).
local signals. Had we included events that involved the recruitment of IP 3 R's from several clusters we might have obtained a power law over a larger range. Signal areas can readily be associated with the size of "avalanches" in typical examples of systems that display SOC. The ability to fit the distribution of Figure 6B with a power law could be an indication that IP 3 Rmediated Ca 2+ signaling in oocytes displays SOC as Ca 2+ waves in astrocytes (Jung et al., 1998) and in cardiac myocytes Nivala et al. (2012). This power law behavior, however, is not captured by the model distributions of Figure 3. The (spatio-temporal) area covered by a signal is strongly dependent on Ca 2+ -mediated channel and cluster coupling. There is experimental evidence, on the other hand, that several local signals usually precede a Ca 2+ (global) wave when cells are continuously subjected to IP 3 www.frontiersin.org September 2012 | Volume 3 | Article 350 | 7 release (Marchant and Parker, 2001). Although cells have very efficient mechanisms to return Ca 2+ back to basal once a signal has finished this does not occur immediately. Thus, there can be a slow Ca 2+ accumulation in between signals which has in fact been observed (Marchant and Parker, 2001). This remaining Ca 2+ effectively couples local signals, an effect that is not taken into account in the intra and intercluster models probed in section 3. Furthermore, IP 3 R's become inactive after having opened. The typical inhibition time has been estimated to be of the order of 2s (Fraiman et al., 2006). This introduces an additional inter-signal coupling. In order to take these two effects into account we modified the intercluster model simulated in section 3. As explained in the "Materials and Methods" section, in order to mimic signal coupling, we do a time-dependent simulation in which we distinguish between active and inactive IP 3 -bound IP 3 R's and where signals are assumed to occur instantaneously. In between signals, the number of active IP 3bound IP 3 R's (the only ones that can become open during a signal) and the Ca 2+ concentration change with time. Thus, the model has an implicit time-scale separation with a fast dynamics during a burst of activity and a slower recovery between such bursts. IP 3 R's become inactive after having opened and go back to being active with a timescale of the order of 2 s. [Ca 2+ ], on the other hand, increases after a signal depending on how many channels were open and then decays back to its basal value with a rate of the order of 200s −1 . For more details on this new model we refer the reader to the "Materials and Methods" section. We show in Figure 7 the results obtained with the intercluster model that takes signal coupling into account. In this case we decided to characterize the signal size by both the number of clusters and the number of IP 3 R's (receptors) that participated of each event. We show the distributions of these two quantities in Figures 7A,B on a logarithmic scale. For this model these two distributions can be approximated by a power law over a relative large range of sizes. Furthermore, the exponent for the distribution of the number of channels is within the range of the one determined from experiments for the puff area distribution. We also fitted these distributions according to the methods described in Clauset et al. (2009) and obtained similar results. The power law approximation fails at large sizes probably due to a finite size effect. Both distributions were not well fitted by other distributions such as log-normal or loglogistic.

MEAN CYTOSOLIC CALCIUM PROVIDES A GLOBAL COUPLING MECHANISM THAT ALTERS THE SIGNAL SIZE DISTRIBUTION
In order to analyze the robustness of the power law behavior encountered in the previous section we simulated the model for other values of its various parameters. More specifically, we studied how the distributions of the number of event-participating clusters and IP 3 R's changed with the mean inter-cluster distance, dm, the mean time an IP 3 R remains inhibited, t inh , the rate at which cytosolic Ca 2+ decays to its basal value, δ Ca and with the radius of influence of an open channel, r inf which were dm = 2 μm, t inh = 2.5s, δ Ca = 200/s and r inf = 0.250 μm in all the simulations of the previous section. We found that in almost all cases the distributions could be fitted by a power law (with exponents between 1 and 1.8) but over a range of event sizes that could be different depending on the parameter values. This is reasonable since the largest event size is strongly dependent on some of these parameters. For example, for large values of the intercluster distance, dm, or of the inhibitory time, t inh and for small values of the mean number of IP 3 -bound IP 3 R's per cluster, λ, the largest event size decreased significantly. Varying r inf between 0.01 μm and 0.1 μm did not produce any noticeable changes. While varying the cytosolic Ca 2+ clearing rate, δ Ca , between 50/s and 100/s did not change the power law behavior of the signal size distribution, using smaller values (we probed 20/s ≤ δ Ca ≤ 30/s) produced a striking change. Namely, the event size distribution ceased to be long-tailed and looked more like Gaussian. We illustrate the results obtained in Figures 8 and 9. We show in By fitting a power law relationship to some of the data points we obtain exponents between 1 and 1.7 with R 2 between 0.76 and 0.92 (if we include data between 20 and 80 IP 3 R's involved, respectively).  (open) IP 3 R's, N o , varies with dm, λ and t inh . We also observe that, when varying any of these 3 parameters the standard deviation of N o decreases with decreasing N o . We have also plotted N o and σ No normalized to the initial number of active IP 3 R's in the simulation obtaining a similar behavior (data not shown). The behavior encountered for small vallues of δ Ca is quite different. Namely, the standard deviation is much smaller than for all the other sets of parameter values explored. Furthermore, the www.frontiersin.org September 2012 | Volume 3 | Article 350 | 9 deviation, σ No , decreases with increasing N o . This is the consequence of the distribution no longer being long-tailed. This can be observed in Figure 9 where we have plotted the distribution of the number of IP 3 R's that participate of the signals when all the parameters are kept at the values used in the previous section except for the cytosolic Ca 2+ clearing rate for which we used δ Ca = 20/s (A), δ Ca = 25/s (B), δ Ca = 30/s (C) and δ Ca = 200/s (D). There we observe that for δ Ca large enough the distribution is well approximated by a power law. As δ Ca decreases, this power law behavior still holds for small enough N o but a new "Gaussianlike component" seems to build in the large N o region. Eventually, when δ Ca is small enough, the power law region disappears and we obtain a Gaussian distribution.

Frontiers in Physiology
In order to interpret this result, we compared the number of receptors that participate of each event, N o and the number of IP 3 R-bound IP 3 R's that are active, N act at any given time. We show a plot of these two numbers as a function of time for δ Ca = 200/s Figure 10A and for δ Ca = 20/s in Figure 10B. We observe that while, in the former N o < N act for most events, for δ Ca = 20/s all active receptors participate of each event (N o = N act ). In this way, the distribution of event sizes coincides with that of active channels which, after a short transient, has a well defined mean value with a relatively small deviation around it. The event size distribution then looks Gaussian.

DISCUSSION AND CONCLUSIONS
Intracellular Ca 2+ signals are ubiquitous across cell types. Among them, those that involve Ca 2+ release from the ER through IP 3 receptors (which are Ca 2+ channels) play a relevant role. Apparently, activatable IP 3 R's are distributed in clusters on the surface of the ER. This uneven distribution together with the kinetic properties of the IP 3 R according to which the channel needs to bind Ca 2+ and IP 3 to become open give rise to a variety of signals that go from very localized ones (puffs) to global waves that propagate throughout the cell. Ca 2+ -mediated channel coupling due to CICR thus plays a key role to determine the spatial size over which the signal spreads. The sequence of channel openings during a signal is then similar to the processes involved in percolation. In particular, we used these ideas in Solovey and Ponce Dawson (2010) to analyze a local signal generation model (i.e., limited to one cluster of IP 3 R's). In Solovey and Ponce Dawson (2010) we called percolation the transition between a situation for which most signals involved the opening of all IP 3 -bound IP 3 R's of the cluster (percolating states) to a situation for which very few IP 3 R's of the cluster opened (non-percolating states). By fitting the model to the experimental results of  (which were obtained in intact mammalian cells of the human neuroblastoma SY5Y cell line) we determined in Solovey and Ponce Dawson (2010) that the experiments seemed to be at a transition point between percolating and non-percolating states. In the present paper we obtained a similar result when looking at the distribution of puff amplitudes determined from experiments performed in Xenopus Laevis oocytes. The experimental distribution of puff (spatio-temporal) areas in this cell type, on the other hand, was slightly different and could be fitted by a power law (goodness of fit) or a log-normal distribution (goodnes of fit). It is always difficult to distinguish power law and log-normal distributions (Clauset et al., 2009). In the present case, it is especially difficult given the small range of event sizes spanned by the experiments which correspond mainly to local signals. However, since puff area is a better reporter than puff amplitude of Ca 2+ -mediated channel coupling when channels are relatively apart from one another we expect a similar type of behavior to be encountered if signals that involve channel recruitment from several clusters are also included in the statistics. In this way, the range of sizes would be larger and a stronger indication on the type of distribution could be drawn. In fact, the power law distributions typical of SOC have been observed for Ca 2+ waves in astrocytes (Jung et al., 1998) and in cardiac myocytes .
In order to investigate the size distribution of signals that involve the opening of IP 3 R's in many clusters we extended the model of Solovey and Ponce Dawson (2010) to include intercluster coupling. The first modification that we tried treated each signal independently of one another. It led to a similar distribution to that of local (intracluster) signals which does not follow a power law. This first extension of the model, however, is not suitable for situations in which cells are constantly challenged with a given stimulus as done in Marchant and Parker (2001). Assuming that the signals act independently of one another is not very realistic in such a case. It has been reported that several local signals usually precede Ca 2+ waves, particularly when the cells are subjected to a constant release of IP 3 as done in Marchant and Parker (2001). This is an indication that, upon a continuous stimulation,

Frontiers in Physiology | Fractal Physiology
September 2012 | Volume 3 | Article 350 | 10 signals can affect one another, especially when they involve the opening of many channels in various different clusters. In order to take this inter-signal coupling into account we further extended the intercluster model including both the dynamics of the global cytosolic [Ca 2+ ] and of the activation and inactivation of IP 3bound IP 3 R's upon channel opening. Using realistic parameter values we obtained distributions for the number of clusters and for the number of channels which participate of the signals that can be fitted by a power law with exponents that are within the values obtained for the experimental spatial puff size distribution. We probed the robustness of this behavior by changing the values of various parameters of the model (the mean inter-cluster distance, dm, the mean time an IP 3 R remains inhibited, t inh , the rate at which cytosolic Ca 2+ decays to its basal value, δ Ca and the radius of influence of an open channel, r inf ). In particular, we determined that the distributions could still be fitted by a power law if we changed any of these parameters around the values used in the original simulations. This robust behavior could lead to the interpretation that our model and thus, the Ca 2+ signals it describes, display SOC. Such a conclusion would not be surprising given that the system under study has the main "ingredients" with which a SOC state can be reached Vespignani et al., 2000). Namely, it is a system with a phase transition onto an absorbing (inactive) state (the state with all channels closed) that is driven slowly by infinitesimal additions of the quantity that brings the system close to the transition to the active state (in our case there are two factors that affect this transition: global cytosolic [Ca 2+ ] and the number of IP 3 -bound IP 3 R's that are active at any given time) and which is allowed to relax upon activation by decreasing this quantity (in our case the number of active IP 3 -IP 3 R's which drops immediately after a signal is elicited). In any case, as already mentioned, it is always difficult to establish with certainty that a given distribution follows a power law or not. In the case of the experimental data that we analyzed the log-normal distribution gave a better fit. However, this could be due to the limited amount of data that we have. It is in fact quite difficult to distinguish a power law from a log-normal distribution (Clauset et al., 2009). It has been shown, for example, that multiplicative processes which naturally lead to log-normal distributions give rise to power law distributions if there is a bounded minimum that acts as a lower barrier to the multiplicative model (Mitzenmacher, 2004). In the case of Ca 2+ signals, there is a lower limit to event sizes (which corresponds to the size of a single channel opening) and one could argue that this could favor a power law over a log-normal distribution. From the point of view of the mechanisms that could lead to one or the other situation, our modeling results seem to favor the power law distribution as well. However, given that our model is only an approximation to the real problem, we cannot tell for sure that IP 3 -mediated Ca 2+ signals are an example of SOC. We are certain, however, that the signal size distribution has a long tail and this could have physiological implications as we discuss later. The event size distribution was not long-tailed for all the parameter values that we probed with our model. In particular, when the cytosolic Ca 2+ clearing time, δ Ca , was decreased by an order of magnitude we obtained a Gaussian-like distribution around a mean for both the number of IP 3 R's ( Figure 9) and of clusters (data not shown) that participate of each signal. We analyzed how the mean and standard deviation of the number of participating IP 3 R's varied with each of the above mentioned parameters. We found that, in most cases probed, the deviation increased with the mean. The exception was the situation encountered for small δ Ca for which we obtained that the deviation, σ No , increased slightly while the mean, N o , decreased when δ Ca increased up to a point at which a large increase in σ No was observed (Figure 8). This signals a transition from a situation with a Gaussian-like distribution to a situation for which the distribution was long-tailed and could be fitted by a power-law (at least over some range of signal sizes). In order to understand why the different behaviors are obtained we compared the number of active IP 3 R's, N act , and that of IP 3 R's that participate of a signal, N o , at any given time. We observed that, while for those cases in which the distribution is long-tailed, N o was less than N act for most events, the number of participating IP 3 R's was equal to that of active ones for almost all events in the case of the Gaussian-like distribution. In particular, we achieve this situation decreasing δ Ca because, in that way, the cytosolic Ca 2+ concentration remains high enough so that no matter how small the addition of the Ca 2+ released by an open IP 3 R is, all IP 3 R's become open as soon as they become active. Thus, in this situation inter-channel coupling plays almost no role. As discussed in Cui et al. (2009), for a model of coupled excitable units that are spatially distributed and that, when "firing" (become open) can induce the "firing" of their active neighbors and then become inactive for a certain refractory time, random openings turn into a random distribution of active units in space. In this way, after a signal that involves the opening of almost all available units, it is very rare that a newly open unit has another active unit in its neighborhood to induce its opening and the signal will fail to propagate. Thus, most likely small signals will follow until a large enough number of IP 3 R's become active again. The random and spatially patchy distribution of active excitable units is a basic ingredient for the existence of a long-tailed distribution of size events as the one observed in our model and in Nivala et al. (2012). This is also key for the large variability of the time interval in between large "spikes" as has been observed in Marchant and Parker (2001); Skupin et al. (2008). Forcing such a system with an external periodic stimulus (like what happens in cardiac myocytes) can lead to a period doubling instability Cui et al.
which can explain Ca 2+ alternans (Jia et al., 2012). When all units are globally coupled (or, rather, they are continously subjected to an environment that make them "fire" as soon as they become active) the random sparse spatial distribution of excitability disappears, the distribution of event sizes ceases to be long-tailed as observed in our model and the period-doubling instability described before disappears. We have observed that the cytosolic Ca 2+ concentration provides a very effective way to switch from a sparsely connected system in which the signal size distribution is long-tailed and may be approximated by a power law to a globally coupled system with a Gaussian-like event size distribution. This last behavior resembles what may observed in the experiments of Marchant and Parker (2001): when a large concentration of IP 3 is continuously released, basal Ca 2+ accumulates and the variability of interspike intervals becomes negligible. Situations with a more efficient coupling mechanism but that had to rely on signal propagation to produce a large event (e.g., setting dm = r inf but keeping λ = 1) failed to give a purely Gaussian distribution giving instead a mixture of a power law for small events and a Gaussian for large ones. Our experiments and those of Marchant and Parker (2001) were performed in immature oocytes. It is known that, with maturation, the ER is reconfigured and the IP 3 R spatial distribution changes which, in turn, leads to Ca 2+ signals that tend to spread over a larger spatial region (Machaca, 2004(Machaca, , 2007. It has been argued that an increase of IP 3 R sensitivity could underlie this behavior of the signals (Ullah et al., 2007). Most likely it is a combination of causes which final effect is to couple IP 3 R's more efficiently preventing the random patchiness of excitability observed in immature oocytes to occur. In this way the mature egg becomes ready to be fertilized and a Ca 2+ wave can propagate without failure (Whitaker, 2006).
Although further studies are necessary to verify our findings, the results obtained in this paper show under what circumstances intracellular Ca 2+ signals can have a long-tailed (SOC-like) size distribution and in which ways this can be changed.