Mean Field Strategies Induce Unrealistic Non-Linearities in Calcium Puffs

Mean field models are often useful approximations to biological systems, but sometimes, they can yield misleading results. In this work, we compare mean field approaches with stochastic models of intracellular calcium release. In particular, we concentrate on calcium signals generated by the concerted opening of several clustered channels (calcium puffs). To this end we simulate calcium puffs numerically and then try to reproduce features of the resulting calcium distribution using mean field models were all the channels open and close simultaneously. We show that an unrealistic non-linear relationship between the current and the number of open channels is needed to reproduce the simulated puffs. Furthermore, a single channel current which is five times smaller than the one of the stochastic simulations is also needed. Our study sheds light on the importance of the stochastic kinetics of the calcium release channel activity to estimate the release fluxes.

. This dependence of the open probability of IP 3 Rs on cytosolic Ca 2+ creates communication between channels. As a result of the spatial organization of channels and the process of Ca 2+ -induced Ca 2+ -release (CICR), cytosolic Ca 2+ signals display a hierarchical spatiotemporal organization which in some cells (e.g., oocytes) have length scales that span over six orders of magnitude (Callamaras and Parker, 1994;Mak and Foskett, 1997;Shuai and Jung, 2003a). The hierarchy of Ca 2+ signals go from local ones such as "blips" (Ca 2+ release through a single IP 3 R) and "puffs" (Ca 2+ release through several IP 3 Rs in a cluster; Sun et al., 1998;Bootman et al., 1997;Yao and Parker, 1991;Parker et al., 1996) to waves that propagate globally across the cell by successive cycles of CICR and Ca 2+ diffusion between clusters (Yao and Parker, 1991;Parker et al., 1996;Callamaras et al., 1998;Sun et al., 1998). These signals have been observed using fluorescence microscopy and Ca 2+ sensitive dyes (Bootman et al., 1997;Callamaras et al., 1998;Sun et al., 1998). However, fluorescence microscopy experiments cannot follow the kinetics of individual IP 3 Rs or the detailed Ca 2+ dynamics within a cluster of channels in most cell types. The spatial resolution of conventional light microscopy is limited by diffraction. This determines that features below the 250-nm length scales are unobservable in typical confocal images. This limit can be beaten using super-resolution techniques. Techniques of this type have been applied to Ca 2+ signals using total internal reflection fluorescence (TIRF) and a fast CCD Parker and Smith, 2010;Wiltgen et al., 2010). The need of using TIRF to guarantee a good axial resolution restricts the application of this approach to cells, as those of the human neuroblastoma SH-SY5Y cell line, where IP 3 Rs are close enough to the plasma membrane. In particular, it does not work in oocytes. A vast body

IntroductIon
An important problem in different areas of science is to relate phenomena that occur at different scales. In biology, an observed behavior at a certain level of organization is often the result of interactions at a lower level, with very different physical and temporal scales. Calcium (Ca 2+ ) signals are an extraordinary example of a multiple scale biological system suitable for being studied with different experimental techniques, each of which focus on a particular temporal and spatial scale. Therefore, it is important to understand the interactions of processes occurring at different scales and how to infer properties at unobserved scales from experimental observations. In this respect, mathematics is invaluable (Cohen, 2004).
Ca 2+ ions are important signaling elements that regulate a large variety of cellular functions such as fertilization, proliferation, development, learning and memory, contraction and secretion (Berridge et al., 2000). The versatility of Ca 2+ as a intracellular cell messenger is possible because cells have mechanisms to control efficiently and precisely the distribution of Ca 2+ in space and time. Within cells, Ca 2+ is kept at high concentrations in stores such as the endoplasmic reticulum (ER). Release of Ca 2+ from the ER occurs through specialized channels among which inositol triphosphate (IP 3 ) receptors play a most relevant role . IP 3 Rs are spatially organized in clusters located on the membrane of the ER and separated by a few microns (Yao et al., 1995). The open probability of IP 3 Rs depends on both the IP 3 and cytosolic Ca 2+ concentration ([Ca 2+ ]; Taylor, 1998;Patel et al., 1999). A key feature is the well-established biphasic action of Ca 2+ in both facilitating and inhibiting the opening of IP 3 Rs. For relatively low [Ca 2+ ], the Ca 2+ released by one channel increases the open probability of neighboring channels, whereas at high [Ca 2+ ] it favors a closed state of the channel (Iino, 1990; Bezprozvanny of literature on Ca 2+ signalsin this cell type exists which is still awaiting a clarification of the underlying Ca 2+ and IP 3 R dynamics. Oocytes, on the other hand, are intact cells and can be used to explore changes in the morphology of the signals with maturation or with depth (Callamaras and Parker, 1999;Machaca, 2007). Thus, their study is of great interest. Unfortunately, fluorescence microscopy experiments in this cell type only give, so far, average (mean field) information of the underlying process. Mathematical tools have been widely used to link experimental observations at the level of puffs with the single channel activity.
As with experiments, there are different modeling approaches that differ on the spatial or time scale they try to resolve. Some models rely on the hypothesis that channels are in such close contact that the concentration of Ca 2+ can be considered homogeneous throughout the cluster Jung, 2002, 2003b;Thul and Falcke, 2004;DeRemigio and Smith, 2005;Diambra and Guisoni, 2005;Ullah and Jung, 2006). This assumption allows for relatively fast simulations of cytosolic Ca 2+ dynamics since Ca 2+ diffusion within the cluster can be neglected. Other models inferred basic puff properties using a deterministic dynamics of spatially localized channels within a cluster combined with an analysis of fluorescence experiments Baran, 2007;Baran and Popescu, 2009;Bruno et al., 2010). These approaches have been useful as a first step for more elaborate analyses that incorporate the stochastic transitions between open and closed states of the channels. Some useful puff properties inferred with deterministic models include estimates of the current of underlying puffs, duration, and spatial width of the release sites.
We have recently analyzed experimental calcium puffs observed in Xenopus laevis oocytes using linescan experiments (Bruno et al., 2010). Using a previously developed backward algorithm (Ventura et al., 2005), we estimated the calcium flux that underlies fluorescent images. Since the algorithm is model independent, the estimates of Bruno et al. (2010) are independent of the dye or the exogenous buffer that are used in the experiments. The results we presented in Bruno et al. (2010) and those of Rose et al. (2006) could be explained by a model in which the total calcium current released through a cluster of IP 3 Rs during a puff is a non-linear function of the number of open channels: linear when the number is small (between 1 and ∼20 channels) and proportional to the square root of the number of channels for larger puffs. This type of scaling could be attributed to local Ca 2+ depletion inside the ER. In fact, in Thul and Falcke (2004), the authors presented simulations of the current and concentration profiles generated by Ca 2+ release from the ER and compared them with signal mass measurements in X. laevis oocytes. Taking into account both the luminal and cytosolic Ca 2+ dynamics, they found that the release current was approximately proportional to the square root of the number of open channels. In view of these results, the change of scaling could be explained as follows (Bruno et al., 2010). Puffs of small amplitude correspond to very few open channels and relatively large distances from one another in which case depletion does not affect the single channel current. Puffs of large amplitude correspond to much smaller inter-channel distances in which case local depletion affects the luminal concentration in the vicinity of any other open channel and as a result, the total Ca 2+ current scales as the square root of the number of channels. Thul and Falcke (2004). However, these models do not describe the stochastic nature of the transitions between open and closed states of individual channels, so they can be viewed as mean field descriptions (Thul and Falcke, 2004;Bruno et al., 2010). Recently, the need of using stochastic and spatially resolved models to unveil the full repertoire of Ca 2+ dynamics has been emphasized (Solovey et al., 2008;Thul et al., 2009;Rüdiger et al., 2010;Solovey and Ponce Dawson, 2010a,b;Diambra and Marchant, 2011;Thurley and Falcke, 2011).
In this work we compare mean field and stochastic puff models, with the latter including the description of individual channel openings and closings. In particular we show that the stochastic model reproduces the experimental observations of Bruno et al. (2010) with a fixed single channel current of ∼0.1 pA regardless of the number of open channels. This implies that the instantaneous total current is proportional to the number of open IP 3 Rs at any given time. It is the average Ca 2+ current released during the whole puff duration that scales non-linearly with the number of open or available channels, as does the (instantaneous) total current in the mean field models of Bruno et al. (2010), Thul and Falcke (2004). The studies of the present paper also show that in order to reproduce the experimental observations under the assumption that channels open and close simultaneously a smaller single IP 3 R current is needed than the actual one. Therefore, using mean field models to derive the underlying channel current from an image can lead to an underestimation of its value. The present study shows the limitations of mean field models and the importance of including a detailed description of the intra-cluster dynamics in order to infer realistic single channel properties from collective signals such as puffs.
The organization of the paper is as follows. In Section 2 we describe the models and the procedure to simulate puffs. Section 3 is divided in four parts. First we compare the time course of the current and fluorescence obtained with simulations of the mean field and stochastic puffs (see Section 3.1). In Section 3.2 we show a non-linear property of the system at the level of puffs, both in experimental data and simulations. In Section 3.3 we present the main result of this work which is the artifactual non-linearities that appear in mean field puff models and in mean field based analyses of experimental puffs. A simple semi-analytical model (see Section 3.4) helps to understand the results. Discussions and conclusions are presented in Section 4.

Puff Models based on a cluster of sPatIally localIzed IP 3 rs
We consider two kinds of models that differ mainly on the way we deal with the channel temporal dynamics. In mean field models, all available channels open and close simultaneously at times chosen a priori while in stochastic models the transitions between the open and closed conformations of the channels is stochastic and follows a channel kinetic model. The number and spatial distribution of the channels is based on the following evidence (Bruno et al., 2010): (a) puffs result from the release of Ca 2+ through a few tens of IP 3 Rs within a cluster; (b) the size of the Ca 2+ release region is ∼500 nm 2 , independent of the number of channels. As we describe in Section 2.1.2, we do not include the dynamics of IP 3 binding and unbinding. Thus, by available channels we mean IP 3 Rs with IP 3 bound.
Clusters consist of N p available channels (N p ≤ 60) randomly distributed within a circle of radius R = 230 nm (stochastic models) or a square of side ∼250 nm (mean field models) at the beginning of the simulation. We only use cluster configurations that gave nearest neighbor channel distances larger than 40 nm in the stochastic simulations or 45 nm in the deterministic case.
The stochastic transitions between open and closed states of the channels occur at time scales much faster than the time scale that the experiments can resolve. For this reason, inferred values of the puff calcium current should be taken as averaged values. In contrast, in the simulations we know exactly how many channels are open with a resolution of ∼0.1 ms, much better than the resolution with which we could infer the Ca 2+ currents from the experimental records in Bruno et al. (2010). Therefore, in order to compare the maximum released current inferred from the experiments with the current used in the simulations we have to average the time course of the simulated current. This was done using a square moving average window.

Mean field models
The number of open channels, N p , is constant during a puff and each open channel releases a constant Ca 2+ current I ch . However, the value of I ch may depend on N p . We analyze two situations: the mean field linear model where I ch is independent of the number of open channels, in which case the total Ca 2+ current is with I ch = 0.1 pA and the mean field non-linear model (Bruno et al., 2010) where I ch depends non-linearly on the number of open channels and the total Ca 2+ current is given by: with I 01 = 0.017 pA, I 02 = 0.08 pA, and N pt = 22. This expression corresponds to the one derived from the analysis of experimental puffs presented in Bruno et al. (2010).

Stochastic models
The single channel current is 0.1 pA, regardless of the number of open channels. We use this value for the current because it agrees with previous estimates of the type 1 IP 3 R current in the ER membrane of X. laevis oocytes . It is also within the range of currents found recently for single recombinant type 3 IP 3 Rs (Vais et al., 2010). The transition between the open and closed conformations of the channels is stochastic. We use two simple IP 3 R kinetic models which differ solely on the behavior of the mean open time as a function of [Ca 2+ ]. We call these models CI (for which the mean open time is independent of [Ca 2+ ]) and CD (for which the mean open time decreases with [Ca 2+ ]). We use these models to test whether a Ca 2+ dependent or a Ca 2+ independent open to closed transition leaves a detectable signature on puffs or not. Both models assume that the channel can be in three states, closed (C), open (O), or inhibited (I). The stochastic transitions between states occur according to the following schemes: These models do not consider the kinetics of IP 3 binding or unbinding to the IP 3 R. According to most IP 3 R kinetic models, these processes occur so fast that the fraction of IP 3 Rs with IP 3 bound can be assumed to be fixed and given by an equilibrium relationship (Young and Keizer, 1992). In that sense, the kinetic scheme actually describes the transitions between states of the IP 3 R with IP 3 bound. In addition, the kinetic schemes do not consider either that the channels can reopen once they enter the inhibited states. This simplification is reasonable to investigate basic puff properties. Namely, the time an IP 3 R remains inhibited under the conditions for which puffs are observed has been estimated as ∼2 s (Fraiman et al., 2006), a time much larger than the duration of a puff. This long inhibition time is compatible with the IP 3 R open probability obtained in reconstituted bilayer experiments (Bezprozvanny et al., 1991) but seems contradictory with the results of patch clamp experiments performed in isolated Xenopus oocyte nuclei using K + as the channel carrier (Mak et al., 1998). The latter type of experiments seems less disruptive than the former one since IP 3 Rs remain in a native membrane. However, the Ca 2+ concentrations that are used on the luminal side of the channel in the patch clamp experiments are much lower than the ones encountered under physiological conditions. As discussed in Fraiman and Ponce Dawson (2004), this could explain the apparent contradiction between the long inhibition time estimated in Fraiman et al. (2006) and the observations of Mak et al. (1998). In fact, the existence of a broad variety of experimental results derived from single channel recordings calls for the definition of an optimal system where the source of the observed discrepancies can be unambiguously established .
The values of the parameters of the kinetic models we consider in this paper are listed in Table 1. The open transition rate, k CO , is taken from the De Young-Keizer IP 3 R model (Young and Keizer, 1992) and the transition rate from the open to the inhibited state is defined so that the mean open time of a single IP 3 R is 10 ms (a typical IP 3 R open duration; Mak et al., 1998 (Solovey et al., 2008). The transitions between states of the channels are computed using a discrete time Markov Chain procedure with a time step 1 μs.
We calculate the duration of the Ca 2+ release during a puff as: CaB is an averaged version of CaB due to the microscope point spread function (Bruno et al., 2010).

seMI-analytIcal Puff Model
In order to gain a more intuitive understanding of the behaviors encountered with the numerical simulations, we also analyze a cartoon-like model that is amenable to analytic calculations. We consider the IP 3 R kinetic models introduced in Section 2.1.2 but the Ca 2+ dynamics is simplified in two ways. First, at time zero all N p channels are open. All they can do is make a transition to the inhibited state. The puff ends when all channels are inhibited. This simplifies enormously the complexity of the Ca 2+ dynamics, specially in the CI model, in which the transition from the open to the inhibited state is independent of Ca 2+ . To further simplify the dynamics of the CD model, we make a second approximation: all channels of the cluster are subject to the same [Ca 2+ ] which is proportional to the number of open channels n, so that the (Solovey et al., 2008) and n = N p , N p − 1, …, 1, 0.
The two versions of the semi analytic model can be considered as extreme cases regarding the channels "coupling strength" (DeRemigio and Nguyen et al., 2005;Groff and Smith, 2008). The one with the CI kinetics corresponds to a fully uncoupled case where channels close stochastically and independently after a exponentially distributed open time with mean value 1/k OI . The one with the CD kinetics, channels are fully coupled by cytosolic Ca 2+ .
The mean puff duration can be computed analytically for both models using Equation 4 and taking t M such that n(t M ) = 0 for the first time (see details of the derivation in the Appendix). The result is that the mean duration of puffs that start with N p open channels is: for models CI and CD, respectively. The average released current, I ave , is given by Equation 5. We do not have analytic expressions for the mean average current. Therefore, in order to compute it, we perform stochastic simulations of the semi analytic models (see Appendix for further details) that show that the mean average current can be approximated by the ratio of mean values: It corresponds to the current of a puff that releases the same total amount of Ca 2+ but at a constant rate during the whole duration of the Ca 2+ release. It is equivalent to a case where all channels of the cluster open and close simultaneously and release the same total amount of Ca 2+ during the same time as the original puff.

Simulations
We solve numerically the set of coupled reaction-diffusion equations that result from the diffusion of cytosolic Ca 2+ in the presence of localized Ca 2+ channels. The details of the methods to simulate puffs are described in Solovey et al. (2008). Therefore, here we only present a brief description. Besides Ca 2+ , we consider the following species: an immobile endogenous buffer, a cytosolic fluorescent Ca 2+ indicator (B) and an exogenous mobile buffer. We use the same parameters we used in Solovey et al. (2008).
In the case of the mean field models, we use a fine spatial grid where individual channels are located. We use clusters of 1, 3, 5, 10, and 15 channels for the linear mean field model and 1, 3,5,10,15,20,25,30,40, and 50 channels for the simulations with the nonlinear mean field model. The release duration was set to 18 ms, the mean duration of Ca 2+ release inferred from experiments in Bruno et al. (2010). In order to do a qualitative description of the differences between stochastic and mean field models (Section 3.1), we used long release durations (∼50 ms), as illustrated in Figure 1.
To generate puffs with the stochastic model we use a simplified method also introduced in Solovey et al. (2008). The method is based on a scale separation, a similar approach as the one taken more recently by Skupin et al. (2010). Namely, we use a fine grid to determine the [Ca 2+ ] distribution within the cluster but describe its dynamics with a quasi-stationary approximation. In order to determine the contribution of each open channel to the [Ca 2+ ] distribution within the cluster, we first analyze this distribution when there is only one channel in the same cytosolic environment. The total [Ca 2+ ] within the cluster when there are several open channels is then approximated by a linear combination of the contributions due to each individual open channel. This detailed description of [Ca 2+ ] within the cluster is used to determine the openings and closings of the various channels. In this way we obtain the total Ca 2+ current that is released from the cluster as a function of time. The spatiotemporal [Ca 2+ ] distribution outside the cluster is determined using a coarser grid in which the cluster is represented by a point source whose current is proportional to the number of open channels determined before. A reaction-diffusion system is solved on this coarser grid. We simulated 10 puffs for a cluster of a given number of channels that we vary between 10 and 50.
The fluorescence F is related to the concentration of Ca 2+ ions bound to the fluorescent dye (CaB): simulated with the mean field model the current is steady and the fluorescence increases during the total duration of Ca 2+ release (Figures 1C,D).
The examples of Figures 1A,B are representative examples of two basic features observed in stochastic puffs: (a) the time to peak current shortens if the number of available channels increases and (b) the difference between puffs simulated with the CI and CD IP 3 R model is not significant. There is a rapid recruitment of IP 3 Rs when Ca 2+ release starts. This initial phase is followed by a usually larger falling phase in which channels close one by one, in agreement with . All subsequent stochastic simulations will be done using model CI unless otherwise mentioned because the fluorescence signals of puffs simulated with the CI and CD models are statistically equivalent.
The amplitude of the puffs simulated with the non-linear mean field model is similar to the amplitudes of the puffs obtained with the stochastic model. The number of channels is also similar. However, for models CI and CD, respectively.

tIMe course of sIMulated Puffs
Simulations can display the time course of puffs at different scales: from the single channel activity to the fluorescence signal. Useful parameters that characterize puffs are the amplitude (A, peak fluorescence) and the rising time (t f , time until the amplitude is reached). Qualitative differences between the stochastic and the mean field models are evident in Figure 1. In the stochastic case (Figures 1A,B), the current varies stepwise and the peak fluorescence is achieved before the end of Ca 2+ release. However, in puffs that dye saturation occurs at much larger amplitudes (A ∼ F max / F 0 = 48) than those of the puffs analyzed in Bruno et al. (2010). In the stochastic simulations, the Ca 2+ -bound dye concentration is never larger than 35 μM (for N p = 50) while the total dye concentration is 40 μM. The reason for the observed non-linearity is the combination of the various processes that modulate the Ca 2+ -bound dye concentration within the volume corresponding to the pixel of the image. This processes depend on the kinetic rates of the reaction, the Ca 2+ current, and the rate at which the Ca 2+ bound dye leaves the volume due to diffusion.
Although puff amplitudes depend non-linearly on the underlying maximum currents, both for the simulations and for the experiments, the mean values of the puff current can be approximated by a function of the mean puff amplitudes with a relatively small dispersion. This indicates that, even if the interaction with the dye acts as a non-linear filter, the puff amplitude could eventually be used to estimate the underlying Ca 2+ current.

the fInIte tIMe resolutIon of the exPerIMents fIlters the observatIons dIfferently dePendIng on the nuMber of oPen channels
We have analyzed in the previous Section the relationship between puff amplitude and the Ca 2+ current that underlies the signal. We have shown, in particular, that these quantities depend non-linearly on one another in mean field models (linear and non-linear), stochastic models (independently of the IP 3 R kinetic scheme), and experiments. Even though the amplitude-current relationship obtained in all these cases is similar, we have observed that the number of channels involved in the linear and the non-linear mean field models are very different. In order to analyze this aspect in more detail, we now explore the dependence between number of channels and puff amplitude.
We show in Figure 3A a plot of puff amplitude as a function of the number of available channels, N p , obtained with simulations of the stochastic puff model. We see that, as N p grows, the mean puff amplitude saturates. This indicates that the mean puff amplitude observed at a given cluster is not a good estimator of the number of channels of the cluster that can become open. We can distinguish an approximately linear regime in which the fluorescence signals contributed by each open channel summate linearly, in agreement with the recent results of . Interestingly, these results show that, besides the non-linear filter that the Ca 2+ -dye reaction introduces between Ca 2+ current and fluorescence, there is an additional source of non-linearity when going from the number of channels to puff amplitude. This additional source of non-linearity is related to the temporal coarse-graining imposed by the experimental temporal resolution acting on data that is characterized by different time courses depending on the size of the puffs. As shown in Figure 1, the time it takes for the current and the fluorescence to achieve their maximum values is shorter the larger is the puff (and the number of open channels that are involved). This means that the fluorescence that is measured experimentally is integrated over a shorter time interval. The dye response time and the finite time resolution of the experiments act as filters of this rapid increment producing the observed non-linear dependence. The fact that the time course is different depending on puff size is in turn related to the spatial organization of IP 3 Rs within the cluster. All models we the linear mean field model reproduces the same amplitude as the stochastic puffs but with significant fewer channels. The examples displayed in Figure 1 illustrate this difference.
Interestingly, the ratio of puff amplitudes is similar in all models when the number of channels is duplicated, as it is shown in Figure 1, comparing puffs from panels Figures 1A-D. This suggest that the ratio of amplitudes can give ambiguous results when trying to determine the "correct" puff underlying model.
As illustrated in Figures 1A,B, the rising time is shorter as more channels participate of the puff. This is a consequence of the fact that the coordination of channels due to CICR evoked by local [Ca 2+ ] elevations is more effective as the distance between IP 3 Rs decreases. Approximately 80% of the puffs simulated using the stochastic model are such that all available channels of the cluster participate although only approximately half of them were simultaneously open (data not shown). In the remaining ∼20% of the puffs, only ∼2 channels did not open during the 150-ms of the simulation time. This is more likely to occur in small clusters (〈N p ∼ 13〉) where the spacing between IP 3 Rs is relatively large and the coordination between channels due to CICR is less effective.
It is important to note that the average current and fluorescence time course of several stochastic simulations does not converge to the mean field counterpart. The mean field model assumes that all channels open and close simultaneously, therefore the release current is constant and the fluorescence increases during the puff.

the aMPlItude and current of Puffs are non-lInearly related
Observed puff amplitudes do not grow linearly with the maximum underlying current estimated from the observations (Bruno et al., 2010). Simulated puffs using the stochastic models show the same behavior, as shown in Figure 2. For large enough currents, the amplitude increases as the current increases, but at a lower rate than for small currents. One possible explanation to these observations is that the fluorescent dye is saturated at the position of the cluster due to the large amount of Ca 2+ ions that are released during a puff. Nevertheless, this is not the case. Calibration experiments showed FIgURe 2 | Non-linear relationship between the amplitude and the release current. Puff amplitude as function of the maximum released current for experimental puffs (adapted from Bruno et al., 2010) and puffs simulated with the stochastic, linear and non-linear models. In the simulations of the stochastic puff model, IP 3 Rs follow the CD kinetic mode. use in the simulations assume that the spatial cluster size is constant regardless of the number of channels that they contain. Thus, the simulated puffs that involve more open channels and have larger amplitudes, occur in clusters where the mean distance between available channels is shorter. This implies that CICR occurs on a faster time scale and the peak current is achieved within a shorter time. This non-linearity is absent in the linear mean field model.
In order to compare the new non-linearity that arises in the simulations of the stochastic puff model with the one that is intrinsic to the non-linear mean field model we compute the average current (I ave , defined in Equation 5) released during a puff simulated using the stochastic model. This average current corresponds to a puff that releases a constant current during the whole release duration (as if all N p channels of the cluster opened and closed simultaneously and released the same total amount of Ca 2+ during the same time as the original puff). Therefore, it is equivalent to replacing the stochastic puff by one that could have been obtained with a mean field model. We found that the derivative of I ave (N p ) decreases as N p grows ( Figure 3B). The deviation from linearity for large N p occurs although the simulations were performed with the stochastic model using a constant single channel current of 0.1 pA, regardless of the number of channels. The exact scaling between the average current and N p is not so relevant. However, it is remarkable that we can fit the data with the non-linear function I ave = I 01 Np for N p ≤ 22 and I I N ave p = 02 1 2 / for N p ≤ 22 (black line in Figure 3B). The square root scaling for large N p agrees with the value predicted in the simulations of Thul and Falcke (2004). From the fitting we find an effective single channel current for small N p , I 01 = 0.019 pA, which is five times smaller than the actual single channel current used in the simulations. Interestingly, it agrees with the value that the non-linear mean field model needed to reproduce the experimental observations (see Equation 2).

understandIng the results usIng the seMI-analytIcal Puff Model
The fact that the effective single channel current is smaller than the actual one can be understood by considering two extreme (hypothetical) puffs with N p available channels and total release duration, t r . In the first extreme case all N p channels open and close simultaneously. Thus, they all remain open during a time t r . The corresponding average current is I ave = I ch N p t r /t r I ch N p , where I ch is the single channel current. In this case, I ave is a linear function of the number of channels that participate of the puff. The other extreme case is a puff in which the channels open and close sequentially so that there is one open channel at any given time. In this case, the average current is I ave = I ch , independent of N p . The simulations using the stochastic puff model give results in between these two extreme behaviors so that I eff = I ave /N p ≤ I ch .
In order to get an intuitive understanding of why the slope between I ave and N p decreases with N p we simulated 500 puffs using the semi-analytical model (Section 2.2). We illustrate the nonlinear dependence between I ave and N p for this model in Figure  4. In this case, we use both kinetic models of IP 3 Rs because they lead to different results. Although it is clearer for CD channels, in both cases the current grows non-linearly as N p increases. For example, in the CI case, 〈I ave 〉 ≅ 0.4 pA for N p = 10 and 〈I ave 〉 ≅ 1 pA for N p = 40. This non-linear relationship is equivalent to that of the non-linear mean field model. Results from simulations of the stochastic puff model that use the CD kinetic model lie in between the semi-analytical results for the two IP 3 R models for N p > 10 (red squares and error bars in Figure 4A). The effective single channel current, I eff = I ave /N p of the semi-analytical models is plotted in Figure 4B as a function of N p . It is clear that as the number of channel increases, I eff approaches the value ∼0.02 pA, similarly to the puffs of the stochastic model.
The semi-analytical model allows us to go one step further in the analysis. We can use the approximation 8 to estimate I ave . According to this expression, if the mean time with n simultaneously open channels, 〈t n 〉, is independent of n, then the average current is linearly related to the number of open channels, 〈I ave 〉 ≅ I ch (N p + 1)/2. However, in most cases 〈t n 〉 does depend on n. In particular, for the two IP 3 R kinetic models that we considered, CI and CD, we obtain 〈t n 〉 ∼ 1/n and 〈t n 〉 ∼ 1/n 2 respectively. Approximating p , we conclude that 〈 〉 ∼ I ave I N lnN p p and 〈I ave (N p )〉 D ∼ ln N p for models CI and CD, respectively. Therefore, we conclude that I ave is a non-linear function of N p with a slope that decreases with N p . Within this view, the ability to approximate I ave by a linear function  In (A,B), black circles and squares correspond to the semi-analytical model using the CD and the CI kinetic model, respectively. Small red squares and error bars represent mean and SD for puffs simulated with the stochastic model. for small N p and a relationship of the form I N ave p ∼ 1 2 / instead of another non-linear function is only due to the limited interval of N p values over which we can observe the function. These estimates also show why I ave is smaller for the CD than for the CI model. In fact, the reason for this is that the transition from the open to the inhibited state is faster for the CD model and that the difference of time-scales increases with the number of channels that are simultaneously open. Thus, the state of the cluster with a large number of simultaneously open channels is shorter for the CD than for the CI model which results in a smaller average current.

dIscussIon
Intracellular signals that involve Ca 2+ release from the ER through clustered IP 3 R-Ca 2+ channels are crucial in a variety of cell types. A complete understanding of the processes involved in such signals requires a combination of experiments and mathematical modeling. Given the multiple time and spatial scales involved, mean field approaches are frequently used in models and data analysis. In a sense, data from fluorescence microscopy experiments is intrinsically of mean field type since it is averaged over the finite space and time resolution of the experimental method. In this work we have compared mean field approaches to stochastic models of Ca 2+ puffs, a type of localized signal that involves Ca 2+ release from several IP 3 Rs in a cluster and which has been observed in intact cells using optical techniques. We found that mean field models only match more accurate stochastic models under two unrealistic assumption: (a) a non-linear relationship between the current and the number of open channels and (b) a single channel current five times lower.
We have used two mean field models that we presented in a previous work (Bruno et al., 2010) where we showed that they approximate the inferred distribution of puff currents and the relationship between the puff amplitude and the underlying current. Both models assume that all available IP 3 Rs of the cluster open and close simultaneously. They differ on how the current depended on the number of channels: linear (Equation 1) and non-linear (Equation 2). We compared simulations of mean field models with simulations of a stochastic puff model that includes the stochastic transitions between open and close states of the channels (Solovey et al., 2008). In order to compare both models, we defined a mean field analog of the stochastic puffs as one with a constant current given by the average current released during the whole puff duration.
The first observation is that the amplitude of puffs is nonlinearly related to the released current, as it was shown in Bruno et al. (2010). This occurred although the fluorescent dye was not saturated. This result was also obtained with the simulations of the stochastic model (Figure 2). Only the non-linear mean field model was able to reproduce this non-linear relationship for the full range of the number of available channels considered here. The linear mean field model described the same behavior but using a much smaller number of channels than the number inferred from the experiments (Bruno et al., 2010; Figure 2).
Recent experimental results concluded that the total calcium current released through a cluster of IP 3 Rs during a puff is a nonlinear function of the number of open channels (Bruno et al., 2010): linear when the number is small and proportional to the square root of the number of channels for larger puffs. Although this type of scaling could be attributed to local Ca 2+ depletion inside the ER (Thul and Falcke, 2004), we argue that this non-linear relationship occurs due to the assumed simultaneously action of the release sites intrinsic to mean field approaches. In fact, even though our stochastic simulations were done with a constant single channel current of 0.1 pA, we found that the average released current of stochastic puffs depended non-linearly on the number of available channels with a slope that decreased with the number of available channels (3B). Furthermore, we showed that the average current could be approximated by Equation 2 with I 01 approximately equal to the single channel current of the non-linear mean field model (∼0.02 pA). This results implies that the mean field assumption also underestimates the single channel current by a factor of 5 (Figure 3).
We investigated the differences between stochastic and mean field approaches analytically as well. For this purpose, we used a simplified stochastic model to describe the dynamics of clustered IP 3 Rs (semi-analytical model). Results using this model confirmed that the assumption that channels open and close simultaneously leads to an underestimation of the single IP 3 R channel current and that the relationship between the average current and the number of available channels is non-linear (Figure 4).
In a previous work (Bruno et al., 2010), we have shown that the non-linear mean field model reproduces the puff to trigger amplitude distribution reported in Rose et al. (2006) better than the linear model. Our results shed light on why this is could be the case. Namely, the puff to trigger amplitude distribution does not take dynamics into account. However, the non-linear mean field model captures the additional non-linear filtering that the experiments impose. Indeed, it is remarkable that the relationship between puff amplitude and current obtained with the non-linear mean field model and the stochastic model conceded using the same range of number of available channels (Figure 2).
There is some controversy on how IP 3 Rs behave (Bezprozvanny et al., 1991;Mak et al., 1998). Consequently, there is a variety of IP 3 R kinetic models (Young and Keizer, 1992;Fraiman and Ponce Dawson, 2004;Shuai et al., 2007) and it is not completely clear which one describes the behavior of the channel under the conditions encountered in intact cells. In particular, it is not clear if the mean open time depends on Ca 2+ or not. Therefore, in the case of the stochastic models, we did simulations using two simple kinetic models of the IP 3 R that only differed on whether the transition from the open to the closed state of the channel was Ca 2+ dependent or not (CI and CD model respectively). The differences between both models are significant in the context of the semi-analytical model (Figure 4). This is easy to understand because the semi-analytical model maximized the differences between them. When using the CD model, the semi-analytical model is such that the [Ca 2+ ] sensed by each channel is additive and proportional to the number of open channels, as if all IP 3 Rs were in the same spatial location. In contrast, when using the CI model, the underlying assumption is that the Ca 2+ released by one channel has no influence on the other open channels. The real situation is somewhere in between both, as it is shown in Figure 4. In general, we did not found any significant difference between the puffs simulated with the stochastic model using any of the two kinetic IP 3 R models. We therefore conclude that this single channel property (a Ca 2+ dependent or a Ca 2+ independent open to closed transition) does not leave a detectable signature on puffs.
Ca 2+ puffs are the result of a complicated underlying spatiotemporal dynamics. The relationship between puff amplitude and the underlying Ca 2+ current can only be understood in terms of the interplay between the intensity and the dynamics of the release. The dynamics of the release, on the other hand, is strongly dependent on the spatial organization of the channels that become open. The combination of all these factors and the experimental limitations can lead to (unexpected) non-linear relationships. Modeling efforts that do not include the variability of the individual channel openings and closings provide averaged descriptions that are still valid. In particular, they can be used as building blocks of mathematical models of more global signals such as waves. However, if one is willing to understand the intra-cluster dynamics or to extract information on single channel kinetics from puff observations, a stochastic description is unavoidable. = and therefore m d (n) = k OI n. The assumption underlying the expression for m d is that the channels are at the same spatial location and [Ca 2+ ] is additive.
The variance of the puff duration can be obtained using Equation A.1 and the fact that t n are independent exponentially distributed random variables. Therefore: It is also possible to calculate the distribution function of t r , F(t). Since t r is the sum of N p exponentially distributed random variables (t n , n = 1, …, N p ) with different mean values [〈t n 〉 = (nm) −1 ], following (Cox, 1962) the distribution function is: where t r is the total release time given by Equation A.1. Therefore, the puff is a random process that is fully defined by t n (n = 1, …, N p ) and each t n is the minimum of n exponential random variables, Equation A.2. So we numerically generate n exponentially distributed random numbers with mean m −1 (in the case of the CI model, m = m i =k OI and in the case of the CD model, m = m d =k OI n) and take the minimum between them which we define as t n . Repeating this operation for each n completes the random process. With these values, we calculate the total release time given by Equation A.1 and the average current given by Equation A.6. Repeating this random process 500 times for each N p we obtain the mean I ave as a function of N p that we plot in Figure 4A. The approximate solutions in Equation 9b are obtained calculating the ratio of mean values in Equation A.6. Although this is not true in general, the results of the simulations show that the relative error of those approximations is always below 15% and decreases with N p (for the CI model) and between 20 and 25% (for the CD model) for the values of N p considered.

aPPendIx seMI-analytIcal Puff Model
We describe here the derivation of the formulas for the mean duration of the puff in the semi-analytical model, Equation 7b, the variance of the duration which we use to calculate effective current in Figure 4B, the simulation procedure to obtain the average current, Equation eq: stoch.metodos.semianalitico.Iave, and the approximate solutions, Equations eq: stoch.metodos.semianalitico. Iavemedio. eq: stoch.metodos.semianalitico.Iavemedio In this model we consider that the puff starts with N p open channels that may inhibit at any time with a rate m that, in principle, may depend on the number of open channels n. After t = 0, the number of channels decreases stepwise, one by one, so the puff duration is the sum of the time elapsed with N p open channels ( ), t N p plus the time elapsed with N p − 1 open channels ( ) t N p −1 and so on until the last channel closes. Therefore, the duration of the Ca 2+ release, Equation 4, can be expressed as: The minimum between n exponential random variables is also exponentially distributed but with mean value 〈t n 〉 = nm −1 . Finally, the mean Ca 2+ release duration of a cluster of N p channels is: Replacing m by the corresponding open to inhibited transition rate, it reduces to Equations 7b. We consider two models as in Section 2.1.2. If the channels of the cluster follow the CI model, they behave independently and close stochastically after a mean dwell time ( ), k OI −1 so m = m i =k OI . On the other hand, model CD considers that IP 3 Rs are coupled by Ca 2+ , so that Ca 2+ release from one channel affects the inhibiting probability of the rest opened channels. Different levels of coupling were considered in Groff and Smith (2008), DeRemigio and Smith (2005)