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 (Ca²⁺) signals, a key component of the signaling toolkit of almost any cell type. Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ concentration plays a major role in determining whether the puff size distribution is long-tailed or not. This suggests that Ca²⁺-clearing from the cytosol is key to determine whether IP₃-mediated Ca²⁺ signals can display a SOC-like behavior or not.

1. Introduction

Ca²⁺ signals are involved in a large variety of physiological processes (Berridge et al., 1998). The diversity of spatio-temporal behaviors that the Ca²⁺ 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²⁺ concentrations in the cytosol leads to cell death. For this reason, the signals are built upon local Ca²⁺ elevations due to Ca²⁺ 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²⁺ channels giving rise to signals that propagate over larger distances (Bootman et al., 1997; Callamaras et al., 1998; Sun et al., 1998). Ca²⁺ release from the endoplasmic reticulum (ER) through inositol 1,4,5-trisphosphate (IP₃) receptors (IP₃R's) is a key component that participates of both local (e.g., “puffs”) and global (waves) Ca²⁺ 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₃R's are organized in clusters on the membrane of the ER (Smith and Parker, 2009). Although the spatial distribution of IP₃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 (Smith et al., 2009). Therefore, Ca²⁺ signals are characterized by very diverse lengthscales: the typical size of the channel pore (~1 nm); the dimensions of the IP₃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²⁺ itself, since IP₃R's become open upon Ca²⁺ and IP₃R binding. Ca²⁺, 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²⁺ 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²⁺ concentration. In any case, the local elevation of [Ca²⁺] that originates at an open channel eventually spreads out in space inducing new channel openings if it finds IP₃-bound IP₃R's on its way. This Calcium-Induced-Calcium-Release (CICR) is basic for Ca²⁺ signal propagation.

The fact that IP₃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²⁺) and the time it takes for Ca²⁺ 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²⁺. 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, 2010; Thurley et al., 2011). In fact, it has long been reported that the occurrence of several local signals can eventually elicit a Ca²⁺ wave (Marchant and Parker, 2001). This transition from local to global Ca²⁺ 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²⁺ signals that we are focusing on in this paper, the latter would correspond to (active) IP₃R clusters. In systems where the dynamics involves some type of propagation that could evnetually fail, as in the case of Ca²⁺ 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₃R's also shapes the kinetics of local signals such as puffs (Diambra and Marchant, 2011). In fact, we investigated the “percolation” of Ca²⁺ signals inside an IP₃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₃R's with IP₃ of the cluster (percolating states) to a situation for which very few IP₃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 Smith and Parker (2009). In particular, by fitting the model to the observations, we determined that, for one open IP₃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₃R's with IP₃-bound would become open during each puff. Actually, the experimental data seemed to be near the transition between the percolating (IP₃-dominated) and non-percolating (for which Ca²⁺ 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, 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²⁺ signals that propagate over different regions and eventually fail to advance. In fact, noise-induced spiral Ca²⁺ waves in astrocytes (Jung et al., 1998) and arrhythmogenic Ca²⁺ waves in cardiac myocytes (Ko et al., 2012) 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 (Dickman et al., 2000; Vespignani et al., 2000). Although this “SOC recipe” is still debatable Pruessner and Peters (2006) it could be the mechanism that underlies the scale-invariant behaviors observed in Ca²⁺ signals as we discuss in this paper.

2. Materials and Methods

2.1. 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₃ (D-Myo-Inositol 1,4,5-Triphosphate, P4(5)-(1-(2-Nitrophenyl)ethyl) Ester) and the exogenous calcium buffer EGTA (Ethylene glycol-bis(2-aminoethylether)-N,N,N′,N′-tetraacetic acid). Final intracellular concentrations of Fluo-4 and InsP₃ 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₃ were from Molecular Probes Inc.; EGTA was from Sigma Aldrich.

2.2. 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.

2.3. 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₀, 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₀ (F − F₀/F₀) 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₀ 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 and its standard deviation σ are calculated for the whole (ΔF/F₀) image. Then all pixels with ΔF/F₀ 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.

FIGURE 1

Figure 1. Linescan image with two puffs in the same site. (A) ΔF/F₀ fluorescence image with the corresponding color code. (B) Δ F/F₀ thresholded binary image.

In order to verify that photobleaching was not affecting our data we grouped together puffs that occurred within the same 1 s-long 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.

FIGURE 2

Figure 2. Puff amplitude as funcion of puff occurrence. Mean amplitude (symbols) and corresponding standard deviation (vertical lines) of puffs that occurred within the same 1 s interval along the experimental record as a function of the time elapsed since the delivery of the UV pulse.

2.4. 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 (ℛ²) 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.

2.5. Mathematical Models

2.5.1. 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) using the parameter values that reproduced the observations of Smith and Parker (2009). Briefly, we consider a circular cluster of fixed radius, R = 250 nm, over which there are N = 18 randomly distributed IP₃R's. Each IP₃R has a probability p = 5/18 of having IP₃ bound (i.e., of being available). The dynamical phase starts by picking an available channel and assuming it is open. We then open all other available IP₃R's within a distance r_inf = 0.250 μm of the open IP₃R. We iterate the scheme for every new open channel until there are no available IP₃R's within r_inf of an open channel that is still closed. We characterize the size of the event by the number of IP₃R's that became open. Repeating the simulation of our local model for many realizations we produce puff size histograms as shown in other sections. This modeled “size” distribution can be compared with the puff amplitude distribution obtained in our experiments provided that all open IP₃R contribute equally to the [Ca²⁺] at the cluster. Care must be taken in doing this since this approximation does not always hold (Solovey et al., 2011).

2.5.2. Intercluster dynamics without calcium accumulation

In order to obtain the distribution of event sizes that involve the opening of IP₃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²⁺ concentration at the location of any other cluster and that all the channels of those other clusters sense the same Ca²⁺ 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²⁺ ions that flow through open IP₃R's and that Ca²⁺ 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²⁺ 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₃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²⁺). In this scheme we do not add up the Ca²⁺ 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.

2.5.3. 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²⁺ (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₃-bound IP₃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₃-bound IP₃R's and the Ca²⁺ 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₃-bound and inactive IP₃-bound IP₃R's. In particular, we assume that only active IP₃-bound IP₃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₃R remains inactive for a certain (random) time. Following (Fraiman et al., 2006) we draw the inhibition time of every inactive IP₃R from an exponential distribution with mean 2.5 s. At t = 0 we assume that all IP₃-bound IP₃R's are active and that [Ca²⁺] is equal to [Ca_basal] = 0.1 μM. It then varies by adding a fixed (spatially homogeneous) amount immediately after the occurrence of a signal per each IP₃R that opened during the event. [Ca²⁺] subsequently decays exponentially with time towards [Ca_basal] with rate 200/s. We add a spatially uniform amount to [Ca²⁺] for each open channel under the assumption that Ca²⁺ diffusion occurs fast enough so that [Ca²⁺] becomes approximately uniform in the vicinity of the clusters on a timescale that is much shorter than the mean separation time between events (Fraiman et al., 2006). In order to determine the fixed contribution of each open channel to the total [Ca²⁺] we reason as follows. First, we assume that, while an IP₃R is open, its Ca²⁺ current is I = 0.1pA and that Ca²⁺ diffuses away from it with D = 200 μm²/s [approximately the free Ca²⁺ diffusion coefficient (Allbritton et al., 1992)]. Thus, when the channel is open it adds a non-uniform contribution Δ[Ca²⁺] to the background [Ca²⁺] that we approximate by the stationary solution of the diffusion equation in the presence of a 2I point source. Namely, Δ[Ca²⁺] ≈ β I μ M μ m³/(4πDr) where β is a unit conversion factor and r is the distance to the open channel. Given the values of I and D considered, it is Δ[Ca²⁺] ≈ 0.414 μ M μ m/r. Then, we assume that the total amount of ΔCa²⁺ contained in a cylinder of height h and radius R which encloses the ER's area of the simulation (π R²) can be approximated by $h{\displaystyle {\int }_0^R\Delta }\left[{\text{Ca}}^{2+}\right]2\pi dr\text{ = 2}\pi hR0.414\mu \text{M}\mu \text{m}$ while the IP₃R is open and gets uniformly distributed over the same volume immediately after the channel becomes closed. In this way, we estimate the fixed spatially uniform contribution to [Ca²⁺] of each open channel immediately after it becomes closed as 2 π h R 0.414 μ M μ m/(hπ R²) = 0.828 μ M μ m/R = 0.0414 μ M. Therefore, after a signal with N_o open channels we add N_o0.0414 μ M to [Ca²⁺]. As mentioned before, this concentration subsequently decays towards its basal value with a rate 200/s that has been estimated as D/1 μ m² given that, for each open channel, Δ[Ca²⁺] 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₃-bound IP₃R's but also on [Ca²⁺]. In particular, we assume that the probability per unit time that a cluster with N_active IP₃R's become open starting a new signal is 0.225s⁻¹*N_active*[Ca²⁺]/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²⁺. 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²⁺] 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_k0.414 μ M μ m/d_ik + [Ca²⁺] ≥ 0.414 μ M μ m/r_inf + [Ca_basal] where d_ik is the distance between the i-th and the k-th cluster, [Ca²⁺] is the (accumulated) background Ca²⁺ at the time of the signal and r_inf = 0.250 μm.

3. Intracluster and Intercluster Percolation Using a Simple Model

In this section we show the results of using the intra and intercluster Models without Ca²⁺ accumulation described in the “Materials and Methods” section to generate distributions of Ca²⁺ signal sizes. As explained in that section, one of these models is used for local signals (puffs) that involve the opening of IP₃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₃ and Ca²⁺ binding, treating each signal as independent of one another. In another section we describe the results obtained when Ca²⁺-mediated signal coupling is included in the intercluster model. Signal types strongly depend on IP₃ and Ca²⁺ since IP₃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₃ that is eventually released by means of a UV flash. In this way, the experimentalist has control over the IP₃ distribution. Therefore, the models assume that [IP₃] is fixed. The dynamics of Ca²⁺, 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²⁺ ion encounters an IP₃-bound IP₃R inducing its opening. Once that happens, the released Ca²⁺ diffuses away from the open channel and can induce the opening of nearby IP₃R's with IP₃ bound. The change in [Ca²⁺] due to one open channel is noticeable only in its vicinity. Therefore, the Ca²⁺-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₃-bound IP₃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²⁺ through N_o open IP₃R's can induce Ca²⁺ release from a nearby cluster if the latter is within a distance N_or_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 Smith and Parker (2009). 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₃ dominated situation in which all the channels of the cluster that have IP₃ bound when the event starts participate of the signal. This occurs when IP₃-bound IP₃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₃-bound IP₃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₃R in the cluster. In this case Ca²⁺-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²⁺-binding limited cases to the IP₃-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₃-bound IP₃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.

FIGURE 3

Figure 3. Event size distributions obtained with simulations of the intra and intercluster models without Ca²⁺ accumulation. (A) Distribution obtained with 1000 realizations of the intracluster (local) model using the parameters described in the “Materials and Methods” section (B) Similar to (A), but for the intercluster model without Ca²⁺ accumulation. See main text for the simulation parameter values.

We show in Figure 3B the event size distribution obtained with 1000 realizations of the intercluster coupling model without Ca²⁺ accumulation. The rationale for the parameter values chosen for these simulations is as follows. Given that the simulation covers a ~100 μm² 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₃-bound IP₃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₃-bound) channels of a cluster become open upon activation of the cluster, the mean numbers of IP₃-bound IP₃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_or_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_or_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 (Smith and Parker, 2009). 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_or_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²⁺-binding dominated and an IP₃-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₃-bound IP₃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₃R's participate of each event.

FIGURE 4

Figure 4. Transition between Ca²⁺-binding limited and IP₃-binding limited cases in the event size distribution of the inter-cluster model. (A) Similar to Figure 3B but using a mean number of IP₃-bound IP₃R's equal to 4. (B) Similar to (A) but using a mean number of IP₃-bound IP₃R's equal to 20.

4. 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₃, 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.

FIGURE 5

Figure 5. Distribution of puff amplitudes obtained experimentally. (A) Distribution obtained by pooling together data from 202 puffs observed in 13 cells. Experiments were performed as explained in the “Materials and Methods” sections. (B) Same as (A) but on a logarithmic scale.

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 ℛ², 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 ℛ² = 0.90. It is always difficult to distinguish power law and log-normal 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.

FIGURE 6

Figure 6. Distribution of puff areas obtained experimentally. (A) Similar to Figure 5 but for the area covered by the puffs. Areas were determined as explained in the “Materials and Methods” section pooling data from 179 puffs observed in 12 cells. (B) Same as (A) but on a logarithmic scale. By fitting a power law relationship to some of the data points we obtain exponents between 1.8 and 2.0 with ℛ² between 0.77 and 0.81 (if we include data between 800 and 4500 px or between 800 and 4000 px, respectively).

5. Calcium Accumulation and IP₃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₃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 local signals. Had we included events that involved the recruitment of IP₃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₃R-mediated Ca²⁺ signaling in oocytes displays SOC as Ca²⁺ 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²⁺-mediated channel and cluster coupling. There is experimental evidence, on the other hand, that several local signals usually precede a Ca²⁺ (global) wave when cells are continuously subjected to IP₃ release (Marchant and Parker, 2001). Although cells have very efficient mechanisms to return Ca²⁺ back to basal once a signal has finished this does not occur immediately. Thus, there can be a slow Ca²⁺ accumulation in between signals which has in fact been observed (Marchant and Parker, 2001). This remaining Ca²⁺ effectively couples local signals, an effect that is not taken into account in the intra and intercluster models probed in section 3. Furthermore, IP₃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₃-bound IP₃R's and where signals are assumed to occur instantaneously. In between signals, the number of active IP₃-bound IP₃R's (the only ones that can become open during a signal) and the Ca²⁺ 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₃R's become inactive after having opened and go back to being active with a timescale of the order of 2 s. [Ca²⁺], 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⁻¹. 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₃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 log-logistic.

FIGURE 7

Figure 7. Distributions of clusters and of receptors involved obtained with simulations of intercluster dynamics with calcium accumulation. (A) Distribution of the number of clusters involved in each signal displayed on a logarithmic scale. By fitting a power law relationship to some of the data points we obtain exponents between 1 and 1.4 with ℛ² between 0.72 and 0.85 (if we include data between 30 and 80 clusters involved, respectively). (B) Distribution of the number of IP₃R's involved on a logarithmic scale. By fitting a power law relationship to some of the data points we obtain exponents between 1 and 1.7 with ℛ² between 0.76 and 0.92 (if we include data between 20 and 80 IP₃R's involved, respectively).

6. 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₃R's changed with the mean inter-cluster distance, dm, the mean time an IP₃R remains inhibited, t_inh, the rate at which cytosolic Ca²⁺ 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₃-bound IP₃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²⁺ 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 Figure 8 the mean and standard deviation of the number of IP₃R's that participate of the signals, N_o, as a function of a varying parameter when all the others are kept fixed at the values used in the previous section. We can observe in this Figure the trends we have just described with respect to how the mean number of participating (open) IP₃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₃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 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₃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²⁺ 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 “Gaussian-like 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.

FIGURE 8

Figure 8. Dependence of the number of IP₃R's that participate of a signal on various parameters. Mean (symbols) and standard deviations (vertical lines around mean) of the number of IP₃R's that participate of a signal, N_o, as a function of the inter-cluster distance, dm (A), of the mean number of IP₃-bound IP₃R's per cluster, λ (B), of the mean IP₃R inhibition time, t_inh (C), and of the cytosolic Ca²⁺ clearing rate, δ_Ca (D). In all the subfigures, the parameters that are not varied are fixed at the values used in section 5.

FIGURE 9

Figure 9. Event size distributions obtained for different values of the cytosolic Ca²⁺ clearing rate. Distribution of the number of IP₃R's that participate of the signals when all the parameters are kept at the values used in section 5 except for the cytosolic Ca²⁺ clearing rate for which we used: (A) δ_Ca = 20/s, (B) δ_Ca = 25/s, (C) δ_Ca = 30/s, (D) δ_Ca = 200/s.

In order to interpret this result, we compared the number of receptors that participate of each event, N_o and the number of IP₃R-bound IP₃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.

FIGURE 10

Figure 10. Number of active IP₃R's and of IP₃R's that are participating of a signal as a function of time. Number of active IP₃R's, N_act (with solid lines) and of open IP₃R's, N_o (with open circles) as a function of time for (A) δ_Ca = 200/s and for (B) δ_Ca = 20/s.

7. Discussion and Conclusions

Intracellular Ca²⁺ signals are ubiquitous across cell types. Among them, those that involve Ca²⁺ release from the ER through IP₃ receptors (which are Ca²⁺ channels) play a relevant role. Apparently, activatable IP₃R's are distributed in clusters on the surface of the ER. This uneven distribution together with the kinetic properties of the IP₃R according to which the channel needs to bind Ca²⁺ and IP₃ 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²⁺-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₃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₃-bound IP₃R's of the cluster (percolating states) to a situation for which very few IP₃R's of the cluster opened (non-percolating states). By fitting the model to the experimental results of Smith and Parker (2009) (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²⁺-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²⁺ waves in astrocytes (Jung et al., 1998) and in cardiac myocytes (Nivala et al., 2012).

In order to investigate the size distribution of signals that involve the opening of IP₃R's in many clusters we extended the model of Solovey and Ponce Dawson (2010) to include inter-cluster 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²⁺ waves, particularly when the cells are subjected to a constant release of IP₃ as done in Marchant and Parker (2001). This is an indication that, upon a continuous stimulation, 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²⁺] and of the activation and inactivation of IP₃-bound IP₃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₃R remains inhibited, t_inh, the rate at which cytosolic Ca²⁺ 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²⁺ 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 (Dickman et al., 2000; 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²⁺] and the number of IP₃-bound IP₃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₃-IP₃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²⁺ 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₃-mediated Ca²⁺ 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²⁺ 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₃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₃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₃R's, N_act, and that of IP₃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₃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²⁺ concentration remains high enough so that no matter how small the addition of the Ca²⁺ released by an open IP₃R is, all IP₃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₃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. (2009) which can explain Ca²⁺ 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²⁺ 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₃ is continuously released, basal Ca²⁺ 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₃R spatial distribution changes which, in turn, leads to Ca²⁺ signals that tend to spread over a larger spatial region (Machaca, 2004, 2007). It has been argued that an increase of IP₃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₃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²⁺ 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²⁺ signals can have a long-tailed (SOC-like) size distribution and in which ways this can be changed.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Guillermo Solovey for providing the code to simulate the intracluster dynamics. This research has been supported by the following grants: UBACyT 20020100100064, PICT 2010-1481, PICT 2010-2767, and CONICET (PIP). Silvina Ponce Dawson is a member of Carrera del Investigador Científico.

References

Allbritton, N., Meyer, T., and Stryer, L. (1992). Range of messenger action of calcium ion and inositol 1,4,5-trisphosphate. Science 258, 1812–1815.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bak, P. (1996). How Nature Works. New York, NY: Copernicus.

Bak, P., and Sneppen, K. (1993). Punctuated equilibrium and criticality in a simple model of evolution. Phys. Rev. Lett. 71, 4083–4086.

Pubmed Abstract | Pubmed Full Text

Bak, P., Tang, C., and Wiesenfeld, K. (1987). Self-organized criticality: an explanation of the 1/f noise. Phys. Rev. Lett. 59, 381–384.

Pubmed Abstract | Pubmed Full Text

Bak, P., Tang, C., and Wiesenfeld, K. (1988). Self-organized criticality. Phys. Rev. A 38, 364–374.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Barella, M., Espada, R., Ponce, M. L., Sigaut, L., and Ponce Dawson, S. (2011). Custom-modification of a commercial confocal microscope to photolyze caged compounds using the conventional illumination module and its application to the observation of Inositol 1,4,5-trisphosphate-mediated calcium signals. J. Biomed Opt. 16, 066013.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Berridge, M. J. (1993). Inositol trisphosphate and calcium signalling. Nature 361, 315–325.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Berridge, M. J., Bootman, M. D., and Lipp, P. (1998). Calcium – a life and death signal. Nature 395, 645–648.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Berridge, M. J., and Dupont, G. (1994). Spatial and temporal signalling by calcium. Curr. Opin. Cell Biol. 6, 267–274.

Pubmed Abstract | Pubmed Full Text

Berridge, M. J., Lipp, P., and Bootman, M. D. (2000). The versatility and universality of calcium signalling. Nat. Rev. Mol. Cell Biol. 1, 11–21.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bootman, M., Niggli, E., Berridge, M., and Lipp, P. (1997). Imaging the hierarchical Ca2+ signalling system in HeLa cells. J. Physiol. (Lond.) 499, 307–314.

Pubmed Abstract | Pubmed Full Text

Choe, C. U., and Ehrlich, B. E. (2006). The inositol 1,4,5-trisphosphate Receptor (IP3R) and its regulators: sometimes good and sometimes bad teamwork. Sci. STKE 2006, re15.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Calabrese, A., Fraiman, D., Zysman, D., and Ponce Dawson, S. (2010). Stochastic fire-diffuse-fire model with realistic cluster dynamics. Phys. Rev. E 82, 1–12.

Pubmed Abstract | Pubmed Full Text

Callamaras, N., Marchant, J. S., Sun, X.-P., and Parker, I. (1998). Activation and co-ordination of InsP3-mediated elementary Ca2+ events during global Ca2+ signals in Xenopus oocytes. J. Physiol. (Lond.) 509, 81–91.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Clauset, A., Shalizi, C. R., and Newman, M. (2009). Power-law distributions in empirical data. SIAM Rev. 51, 661–703.

Coombes, S., and Timofeeva, Y. (2003). Sparks and waves in a stochastic fire-diffuse-fire model of ca²⁺ release. Phys. Rev. E 68, 021915.

Pubmed Abstract | Pubmed Full Text

Cui, X., Rovetti, R. J., Yang, L., Garfinkel, A., Weiss, J. N., and Qu, Z. (2009). Period-doubling bifurcation in an array of coupled stochastically excitable elements subjected to global periodic forcing. Phys. Rev. Lett. 103, 044102.

Pubmed Abstract | Pubmed Full Text

Diambra, L., and Marchant, J. S. (2011). Inositol (1,4,5)-trisphosphate receptor microarchitecture shapes ca(2+) puff kinetics. Biophys. J. 100, 822–831.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Dickman, R., Muñoz, M., Vespignani, A., and Zaperi, S. (2000). Paths to self-organized criticality. Braz. J. Phys. 30, 27–41.

Foskett, J. K., White, C., Cheung, K. H., and Mak, D.-O. D. (2007). Inositol triphosphate receptor Ca2+ release channels. Physiol. Rev. 87, 593–658.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fraiman, D., Balenzuela, P., Foss, J., and Chialvo, D. R. (2009). Ising-like dynamics in large-scale functional brain networks. Phys. Rev. E 79, 061922.

Pubmed Abstract | Pubmed Full Text

Fraiman, D., Pando, B., Dargan, S., Parker, I., and Dawson, S. P. (2006). Analysis of puff dynamics in oocytes: interdependence of puff amplitude and interpuff interval. Biophys. J. 90, 3897–3907.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jia, Z., Bien, H., Shiferaw, Y., and Entcheva, E. (2012). Cardiac cellular coupling and the spread of early instabilities in intracellular Ca2+. Biophys. J. 102, 1294–1302.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jovanović, B. C. V., Buldyrev, S. V., Havlin, S., and Stanley, H. E. (1994). Punctuated equilibrium and “history-dependent” percolation. Phys. Rev. E 50, R2403–R2406.

Pubmed Abstract | Pubmed Full Text

Jung, P., Cornell-Bell, A., Madden, K. S., and Moss, F. (1998). Noise-induced spiral waves in astrocyte syncytia show evidence of self-organized criticality. J. Neurophysiol. 79, 1098–1101.

Pubmed Abstract | Pubmed Full Text

Keizer, J., Smith, G. D., Ponce-Dawson, S., and Pearson, J. E. (1998). Saltatory Propagation of Ca2+ Waves by Ca2+ sparks. Biophys. J. 75, 595–600.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ko, C. Y., Nivala, M., Nivala, M., Weiss, J. N., and Qu, Z. (2012). Self-organized criticality underlies arrhythmogenic calcium waves in cardiac myocytes. Biophys. J. 102, 100a.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Machaca, K. (2004). Increased sensitivity and clustering of elementary Ca2+ release events during oocyte maturation. Dev. Biol. 275, 170–182.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Machaca, K. (2007). Ca 2+ signaling differentiation during oocyte maturation. J. Cell. Physiol. 213, 331–340.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Marchant, J. S., and Parker, I. (2001). Role of elementary Ca2+ puffs in generating repetitive Ca2+ oscillations. EMBO J. 20, 65–76.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mitzenmacher, M. (2004). A brief history of generative models for power law and lognormal distributions. Internet Math. 1, 226–251.

Nivala, M., Ko, C. Y., Nivala, M., Weiss, J. N., and Qu, Z. (2012). Criticality in intracellular calcium signaling in cardiac myocytes. Biophys. J. 102, 2433–2442.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pearson, J. E., and Ponce-Dawson, S. (1998). Crisis on skid row. Phys. A Stat. Mech. Appl. 257, 141–148.

Ponce Dawson, S., Keizer, J., and Pearson, J. E. (1999). Fire-diffuse-fire model of dynamics of intracellular calcium waves. PNAS 96, 6060–6063.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pruessner, G., and Peters, O. (2006). Self-organized criticality and absorbing states: lessons from the ising model. Phys. Rev. E 73, 025106.

Pubmed Abstract | Pubmed Full Text

Skupin, A., Kettenmann, H., and Falcke, M. (2010). Calcium signals driven by single channel noise. PLoS Comput. Biol. 6:e1000870. doi:10.1371/journal.pcbi.1000870

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Skupin, A., Kettenmann, H., Winkler, U., Wartenberg, M., Sauer, H., Tovey, S. C., Taylor, C., and Falcke, M. (2008). How does intracellular Ca2+ oscillate: by chance or by the clock? Biophys. J. 94, 2404–2411.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Smith, I. F., and Parker, I. (2009). Imaging the quantal substructure of single IP3R channel activity during Ca2+ puffs in intact mammalian cells. Proc. Natl. Acad. Sci. U.S.A. 106, 6404–6409.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Smith, I. F., Wiltgen, S. M., Shuai, J., and Parker, I. (2009). Ca(2+) puffs originate from preestablished stable clusters of inositol trisphosphate receptors. Sci. Signal. 2, ra77.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Solovey, G., Fraiman, D., and Dawson, S. P. (2011). Mean field strategies induce unrealistic nonlinearities in calcium puffs. Front. Physiol. 2:46. doi: 10.3389/fphys.2011.00046

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Solovey, G., and Ponce Dawson, S. (2010). Intra-cluster percolation of calcium signals. PLoS ONE 5:e8997. doi: 10.1371/journal.pone.0008997

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sornette, D., and Dornic, I. (1996). Parallel bak-sneppen model and directed percolation. Phys. Rev. E 54, 3334–3338.

Pubmed Abstract | Pubmed Full Text

Stauffer, D., and Aharony, A. (1992). Introduction to Percolation Theory, 2nd Edn. London: Taylor & Francis.

Sun, X.-P., Callamaras, N., Marchant, J. S., and Parker, I. (1998). A continuum of InsP3-mediated elementary Ca2+ signalling events in Xenopus oocytes. J. Physiol. (Lond.) 509, 67–80.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Taufiq-Ur-Rahman Skupin, A., Falcke, M., and Taylor, C. W. (2009). Clustering of InsP3 receptors by InsP3 retunes their regulation by InsP3 and Ca2+. Nature 458, 655–659.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tetzlaff, C., Okujeni, S., Egert, U., Wrgtter, F., and Butz, M. (2010). Self-organized criticality in developing neuronal networks. PLoS Comput. Biol. 6:e1001013. doi: 10.1371/journal.pcbi.1001013

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Thurley, K., Smith, I. F., Tovey, S. C., Taylor, C. W., Parker, I., and Falcke, M. (2011). Timescales of IP3-evoked Ca2+ spikes emerge from Ca2+ puffs only at the cellular level. Biophys. J. 101, 2638–2644.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Timofeeva, Y., and Coombes, S. (2004). Directed percolation in a two-dimensional stochastic fire-diffuse-fire model. Phys. Rev. E 70, 062901.

Pubmed Abstract | Pubmed Full Text

Ullah, G., Jung, P., and Machaca, K. (2007). Modeling Ca2+ signaling differentiation during oocyte maturation. Cell Calcium 42, 556–564.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Vespignani, A., Dickman, R., Muñoz, M. A., and Zapperi, S. (2000). Absorbing-state phase transitions in fixed-energy sandpiles. Phys. Rev. E 62, 4564–4582.

Pubmed Abstract | Pubmed Full Text

Whitaker, M. (2006). Calcium at fertilization and in early development. Physiol. Rev. 86, 25–88.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text