Particle-based Multiscale Modeling of Calcium Puff Dynamics

Intracellular calcium is regulated in part by the release of Ca$^{2+}$ ions from the endoplasmic reticulum via inositol-4,5-triphosphate receptor (IP$_3$R) channels (among other possibilities such as RyR and L-type calcium channels). The resulting dynamics are highly diverse, lead to local calcium"puffs"as well as global waves propagating through cells, as observed in {\it Xenopus} oocytes, neurons, and other cell types. Local fluctuations in the number of calcium ions play a crucial role in the onset of these features. Previous modeling studies of calcium puff dynamics stemming from IP$_3$R channels have predominantly focused on stochastic channel models coupled to deterministic diffusion of ions, thereby neglecting local fluctuations of the ion number. Tracking of individual ions is computationally difficult due to the scale separation in the Ca$^{2+}$ concentration when channels are in the open or closed states. In this paper, a spatial multiscale model for investigating of the dynamics of puffs is presented. It couples Brownian motion (diffusion) of ions with a stochastic channel gating model. The model is used to analyze calcium puff statistics. Concentration time traces as well as channel state information are studied. We identify the regime in which puffs can be found and develop a mean-field theory to extract the boundary of this regime. Puffs are only possible when the time scale of channel inhibition is sufficiently large. Implications for the understanding of puff generation and termination are discussed.

1. Introduction. Intracellular calcium plays a major role in many signaling pathways and regulates enzymatic activity, gene expression [7], and neural activity [60,36]. In order to control a variety of cell functions, cells control the local cytoplasmic calcium concentration via the exchange of Ca 2+ ions with the extracellular space and the release of Ca 2+ ions through channels situated on the membrane of reservoirs, such as the endoplasmic and the sarcoplasmic reticulum. Modulation of the Ca 2+ release regulates muscle contraction, pathway cross-talk, and mitochondrial activity, and disruption of these processes is associated with various diseases such as early-onset Alzheimer's [5,57], heart failure [3,2], and such psychological conditions as bipolar disorder and schizophrenia [6]. Hence, detailed knowledge about the underlying processes governing this signaling mechanism is necessary for progress in our understanding of these diseases.
In this paper, we focus on the release of calcium ions from the endoplasmic reticulum (ER) via inositol-4,5-triphosphate receptor (IP 3 R) channels. Upon the binding of Ca 2+ to binding sites on its cytosolic part, a channel opens, and calcium ions flow from the ER into the cytoplasm. Channels occur in clusters of 10 to 20 channels [12,47]. The opening of a single channel usually triggers the release of Ca 2+ from other channels in the same cluster due to the increased Ca 2+ concentration in their vicinity [46]. This mechanism results in a highly localized increase of the cytosolic calcium concentration. These "puffs" of calcium ions have been detected and analyzed in experiments by using fluorescent calcium buffers [58,23].
Traditional approaches to modeling intracellular calcium dynamics are based on deterministic macroscopic rate equations [27,11]; however, the intrinsically random, erratic nature of calcium signals in many cells or cell domains necessitates an approach going beyond the deterministic regime [59,24,31,44]. Progress has been made by recognizing the importance of number fluctuations in the binding to the channels [48,18] and by using hybrid models, where the deterministic calcium concentration is coupled to stochastic channel binding models [42]. Recently, it was found that local fluctuations stemming from diffusive noise (i.e., noise originating from the random movement of ions) have a crucial influence on calcium dynamics for clusters of intracellular channels, particularly on the interpuff waiting time [22] but also on single-channel equilibrium behavior [56]. Stochastic effects were also investigated for L-type calcium channels and RyR channels in the dyadic cleft [34,49,28]. Other studies consider a Langevin equation governing the fraction of open channels in a cluster [54]. However, tracking the exact diffusion of individual ions in the complete computational domain is computationally intensive.
Here, for the first time, we apply spatial stochastic multiscale methods to model the dynamics of calcium puffs, including the release of Ca 2+ from IP 3 R channels, and track individual ion positions in order to accurately incorporate diffusive noise. We take into account both activating and inhibitory channel properties. We study the dynamics of the calcium concentration as a function of the ion binding affinities and explore the regime where there exist puffs, as well as the regime in which channels do not close after their first opening. The paper is organized as follows: In section 2 we discuss the spatial stochastic model for diffusion and channel gating used throughout this study. Section 3 discusses the multiscale approach that we employ in order to reduce the computational effort required to track single ions and hence make this study feasible. In section 4 we present results on the statistics of puffs extracted from simulated time series data, and study the transition between perpetually open channel clusters and the parameter regime in which puffs can be observed. In addition, we develop a mean-field model for channel dynamics and use it to extract the boundary between the regimes. Finally, we summarize our findings in section 5.
2. Spatial stochastic model for intracellular Ca 2+ release. The spatial extent of our computational model consists of the three-dimensional domain Ω which models a part of the intracellular space. Ca 2+ ions are able to undergo free diffusion in Ω. They bind to and dissociate from binding sites on the channels, which are positioned on a small area of the domain boundary, corresponding to the membrane of the ER. The domain geometry and boundary conditions are specified in section 2.2. In the following subsections we discuss the components of this model in detail. The parameter values used throughout this study can be found in Table 1.

Diffusion-Brownian Dynamics.
A versatile method for the simulation of particles in a solution is given by Brownian Dynamics (BD). Collisions of particles with solution molecules lead to overdamped dynamics and random forcing on a sufficiently long time scale [15]. Assuming that there are Q(t) free ions in the simulation where D is the diffusion coefficient of ions in solution, X j (t) ∈ Ω describes the trajectory of the jth ion, and W j is a three-dimensional vector of independent Wiener processes. This approach dramatically reduces the dimensionality of the problem compared to molecular dynamics approaches wherein the degrees of freedom of every participating molecule need to be taken into account. Nevertheless, the computational load is still high compared to deterministic PDE-based approaches to diffusion. There are a number of approaches for simulating (2.1) in the literature, ranging from discretization with a fixed time step [4] to event-based methods [52,37]. In this paper, we discretize time with the time step ∆t and use the Euler-Maruyama discretization of (2.1); i.e., the position of the jth ion is updated according to where ξ j is a vector of three independent normally distributed random numbers with zero mean and unit variance.
2.2. Domain geometry and boundary conditions. In our simulations, the computational domain is given as cube Ω = [0, L] 3 , where the value of L is specified together with other parameters in Table 1. Ca 2+ ions are able to undergo free diffusion in Ω, which we simulate using (2.2). A cluster of nine IP 3 R channels is positioned in c 2016 SIAM. Published by SIAM under the terms of the Creative Commons 4.0 license a 3 × 3 grid with channel spacing = 0.15 µm, centered in the z = 0 plane; i.e., the positions of nine channels in the cluster are given as We found no significant effects from varying the channel spacing ; hence we chose this particular cluster configuration for ease of implementation. Ions bind to and dissociate from binding sites on the channels. The boundaries of the computational domain at x = 0, x = L, y = 0, y = L, and z = L are constant-concentration boundaries, and hence they absorb and introduce ions such that in the absence of open channels the concentration of Ca 2+ in the computational domain is held at its equilibrium value, c 0 , on average. Hence, ions can enter and leave the domain via the boundaries, in addition to being introduced through open channels. The boundary at z = 0 is reflective, corresponding to the membrane of the ER. If the channels are closed, then the average number of free ions in the computational domain, Q(t) , is equal to c 0 |Ω|, where |Ω| = L 3 is the volume of Ω. Using our parameter values (see Table 1), c 0 |Ω| ≈ 1.5 × 10 3 . The average number of ions in the simulation domain, Q(t) , increases when the channels are open (by around two orders of magnitude). In order to simulate the system over sufficiently long time intervals, we will use a multiscale approach, which we describe in detail in section 3.

2.
3. Stochastic channel binding model. The conformational changes between the open and closed states of IP 3 R channels are controlled by the binding of Ca 2+ to activating and inhibitory binding sites. These channels consist of four subunits, each of which can be in an active or a neutral/inhibited state. Each subunit has three different binding sites: An activating binding site for Ca 2+ ions, an inhibitory binding site for Ca 2+ ions, and an IP 3 binding site.
To accurately capture channel gating events, we employ a simplified DeYoung-Keizer model [11]. Here, we disregard IP 3 dynamics and consider the effects of IP 3 only via their influence on the dissociation constant b i of Ca 2+ ions from inhibitory binding sites. Free Ca 2+ ions can bind to the activating and inhibitory sites, while bound ions can dissociate from occupied sites. Therefore, there are two reversible reactions in our model for each subunit of a channel: The rates a a and a i describe the binding affinity of a Ca 2+ ion to an activating or inhibitory binding site, respectively, while the off-rates b a and b i describe the corresponding dissociation reaction rates. The variables S a and S i describe the binding site state and can take values of 0 and 1. Channels consist of four subunits, with each subunit having one activating and one inhibitory Ca 2+ -binding site; see Figure 1. A model subunit can then be in three distinct states: neutral (no binding site occupied), active (only the activating site occupied), and inhibited (if the inhibitory site is occupied regardless of the state of the activating site). A channel then opens when at least three of its four subunits are in the active state.  to their corresponding binding sites. Ions become binding candidates as soon as they enter a half-sphere of radius around a binding site. (The channels are positioned on the ER membrane, and hence only the half-sphere above them is available for ion binding.) They are then allowed to bind to the site with a probability P λ per time step spent in the binding region [35,16]. An ion bound to an activating site (resp., inhibitory site) is allowed to dissociate with a probability 1 − exp(−b a ∆t) (resp., 1 − exp(−b i ∆t)) and placed at a distance of σ (unbinding radius) from the binding site. We use the values of binding and unbinding radii, and σ, as given in Table 1. Their values are chosen to be reasonable, such that no overlaps occur between channels and that they are large enough such that the BD simulation time step can be chosen reasonably large. Then we can precalculate the remaining parameter P λ (binding probability) before the start of simulations using the approach described in [35, section 5].
2.5. Channel opening and calcium flux. Each channel has four subunits. Let S j,k a (t) ∈ {0, 1} (resp., S j,k i (t) ∈ {0, 1}), j = 1, 2, . . . , 9, k = 1, 2, 3, 4, be the state of the activating (resp., inhibiting) site of the kth subunit of the jth channel in the cluster. Then the number of active subunits of the jth channel at time t is   Table 1, that our computational domain may contain on the order of 10 5 ions during the peak of a puff.  implementation of BD binding dynamics, described in section 2.4. Therefore, we split our computational domain into two regions: a cube containing the channel sites (red region in Figure 2), as well as the remaining space Ω \ Ω BD in which we will use a coarser description of ion movement, as described in the next subsection. Here, L BD < L is the length of the edge of the cube Ω BD . We use L BD = L/5 in our simulations (see Table 1). In particular, we use BD simulations in a small fraction of 1/5 3 ≈ 0.8% of the computational domain Ω.
3.1. Compartment-based model for diffusion. We subdivide the region Ω \ Ω BD into compartments (cubes) with size h (illustrated in the bottom left of Figure 2) and employ a compartment-based method to simulate the movement of ions [17]. Ions are allowed to move between adjacent compartments with a rate of d = D/h 2 . The compartment-based algorithm stores and evolves only the number of ions in each compartment (rather than following individual ions). An event-based stochastic simulation algorithm is used to efficiently simulate this system. Several equivalent methods have been developed in the literature, such as the Gillespie algorithm [26], the Next Reaction Method [25], the Next Subvolume Method [14,30], as well as the Optimized Direct Method [9].
In this paper, we employ the next reaction method [25]. For each possible move between two neighboring compartments, a putative time for the next jump of an ion to occur is calculated. It is given by where t is the current simulation time, d = D/h 2 is the jump rate for one ion, A c (t) is the current number of ions in the compartment from which an ion is jumping, and r is a random number uniformly distributed in the interval (0, 1). Clearly, the putative jump time (3.1) is infinity if A c (t) = 0, i.e., if the corresponding compartment is empty. The putative times are smaller (on average) if the corresponding compartment contains more ions.
The putative jump times are inserted into a priority queue (a heap data structure [10]), which enables us to efficiently extract the earliest jump time and thus the next jump. This move is then performed, and the numbers of ions in compartments (and the corresponding putative times and their entries in the priority queue) are updated. We then iterate this process by finding the minimal putative time and performing the corresponding ion jump at each iteration.
At the boundaries of the computational domain Ω (except for the boundary at z = 0, which is reflective), ions are absorbed (i.e., they leave the domain) or can enter with rates consistent with a constant equilibrium concentration of c 0 outside the domain. The jump rate from outside the domain into a compartment just inside the domain boundary is c 0 Dh, which is used, instead of dA C (t), in (3.1) to compute the corresponding putative times.

3.2.
Coupling the BD simulation in Ω BD with the compartment-based approach in Ω \ Ω BD . Several methods exist in order to couple the BD and compartment-based methods across their interface, such as the two-regime method (TRM) [19,20] or the ghost-cell method [21]. Here, we employ the TRM. In order to accurately capture diffusion across the interface, the jump rates from adjacent compartments into Ω BD need to be adjusted from the bulk rate d = D/h 2 to the interface rate [19] (3.2) These jump rates are used, instead of d, in (3.1) to compute the corresponding putative times. If the chosen move in the compartment regime is a jump from a compartment adjacent to the interface into Ω BD , the occupancy number of the compartment is reduced by one, and a new ion is added in Ω BD at the distance x from the interface which is sampled from the probability distribution [19] where erfc is the complementary error function. The rate (3.2) and the distribution (3.3) are used to transfer ions from Ω \ Ω BD into Ω BD . In the opposite direction, the TRM transfers any ion from Ω BD which during the time step interacts with the interface. For a detailed discussion and the derivation of the probabilities and jump rates above, please see references [19,20,21]. We employ the particle-based simulation library package Tyche, which implements the TRM [38]. The TRM has also been recently implemented in Smoldyn [39]. where the elements in S (resp., E) are ordered, i.e., s 1 < s 2 < · · · < s Np (resp., e 1 < e 2 < · · · < e Np ). The interpuff times are then given by where the threshold t p = 0.25 s. This threshold filtering decreases the number of puffs considered by removing the puffs which are not well separated. Indeed, the inset in Figure 4(a) shows that the distribution of interpuff times below the threshold follows an exponential distribution and is due to random channel reopenings stemming from residual calcium. The distribution of interpuff times (histogram of set T ) is shown in Figure 4(a) for a set of 954 puff intervals. Its shape is similar to a Gamma distribution: Indeed, when choosing κ = 2.82 and θ = 0.64 s in (4.4) (such that the resulting distribution has the same mean and variance as our data), the comparison is good. As the Erlang distribution (the special case of a Gamma distribution with κ ∈ N) results from summing κ exponentially distributed random variables, we can speculate that an effective model of calcium channel opening might be described by a threestage binding/unbinding process. The maximum of the interpuff time distribution, τ m ≈ 1.81 s, can be interpreted as the typical time between puffs in this system.
The distribution of puff amplitudes is determined by taking the maximum value in a previously determined puff interval, The amplitude distribution (histogram of set A) is shown in Figure 4(b) and mirrors the broad distribution characterized in experiments [13].
The puff durations are calculated by considering the full duration at half maximum, which is a common criterion for determining the puff duration in the literature [13]. To this end, we calculate the indices of crossings of the half maximum threshold on the rising slope of a puff, where s k ∈ S, e k ∈ E, a k ∈ A, 1 ≤ k ≤ N p , and the falling slope of a puff, where s k ∈ S, e k ∈ E, a k ∈ A, 1 ≤ k ≤ N p .
The puff durations are then calculated via The distribution of puff durations (histogram of set D) is plotted in Figure 4 Figure 5 shows the results of a simulation run with two different values of inhibitory site dissociation rate b i , illustrating two different modes of behavior. We use the same value of a i in all four panels; the rest of the parameter values are given in Table 1. Figures 5(a)-(b) show the results of a simulation run with b i = 5 s −1 (corresponding to a dissociation constant of K D = b i /a i = 50 µM). These simulation parameters are consistent with the modeling of calcium puffs in the literature and are used in hybrid PDE-based models [43]. Figure 5(a) shows the Ca 2+ concentration in the computational domain Ω, while Figure 5(b) displays the fraction of occupied inhibitory sites in the channel cluster as a function of time. Due to the large dissociation rate, the necessary level of bound inhibitory sites can never be sustained long enough for all the channels in the cluster to close at the same time. Therefore, the Ca 2+ concentration is kept at a level where empty activating sites are immediately filled and the channel cluster stays perpetually open. Hence, puff termination in this case can be speculated to be facilitated by a mechanism other than channel inhibition, such as ER Ca 2+ reservoir depletion or a mechanism involving dissociation of IP 3 [40]. We checked the influence of lowering the binding radius to 6 nm and found that it had no effect on the channel cluster closing.

Inhibition transition.
In contrast, a lower inhibitory site dissociation constant yields well-defined calcium puffs in our particle-based simulation scheme.
where the averages are again take over all values of j, j = 1, 2, . . . , N (compare with (4.2)). The puff score (4.5) can take values in [0, 1]. A puff score greater than 0.25 indicates channel excitability and therefore the existence of puffs in the system. In order to visualize the two parameter regimes, we performed simulations for different inhibitory site binding parameters a i ∈   Table 1. Note that the time scale of inhibition decay is long enough for the calcium concentration to decay completely; thus the channel cluster closes and we can identify well-defined puffs.  Table 1. The black line shows the numerically determined phase boundary from the meanfield model (4.6)-(4.8).
color indicates the value of [PS](a i , b i ). The transition between the two regimes is not sharp, but gradual, especially for higher a i . This is due to prolonged channel reopenings becoming more likely due to faster dissociation of bound inhibiting ions when b i approaches the transition. The phase boundary is consistent with a dissociation constant of 4 µM < K D < 10 µM.
To directly compare our results with previously used schemes from the literature, we devised a simplified hybrid scheme: The main difference between our method and hybrid simulation algorithms is the much lower calcium concentration in the vicinity of open channels and the resulting weaker channel inhibition. Therefore, in the simplified hybrid scheme, whenever a channel is open, we do not use the particle-based binding described in section 2.4. Instead we assume a constant high calcium concentration of c B = 150 µM (consistent with results from hybrid simulations [43]) to generate random binding events to inhibitory sites with a rate of a i c B , irrespective of any ions in the vicinity. In all other respects, the simulation proceeds as described previously. Figure 7 shows the resulting map of puff scores.
In order to study the influence of buffers (whose main effect is to slow down ion diffusion [33]) on the boundary of the puff regime, we performed a similar set of simulations with a lower Ca 2+ diffusion constant of D = 20 µm 2 /s. Because ions bound to buffer molecules can be viewed as slowly diffusing ions, lowering the diffusion constant is a way to model the effects of buffers [53]. The result is shown in Figure 8. With a lower diffusion constant, the calcium concentration in the channel cluster nanodomains decays more slowly. For puffs to exist, the time scale of the decay of inhibitory site binding needs to be longer than the time scale of Ca 2+ decay. Hence, the boundary separating the two regimes is pushed to smaller values of the inhibitory dissociation rate b i , corresponding to approximately 1 µM < K D < 2 µM. Hence, the effective diffusion constant plays an important role in determining the boundary between the two regimes. Note that the overall variation of the puff score [PS](a i , b i ) here is smaller compared to the case of D = 220 µm 2 /s; therefore Ca 2+ puffs become less pronounced with slower ion diffusion. Ca 2+ buffers lead to additional extrinsic  Table 1. For comparison, the black line shows the numerically determined phase boundary from the mean-field model (4.6)-(4.8) for Figure 6.  Table 1. Color indicates the puff score (4.5). The black line shows the numerically determined phase boundary from the mean-field model (4.6)- (4.8). Note that the scale of the y-axis here is different from that in Figures 6 and 7. noise due to the binding and unbinding of calcium ions to buffer molecules, which enhances Ca 2+ fluctuations, which might have effects on the puff statistics in addition to what we have presented here [55].
Here, Θ(x) is the Heaviside function with the properties The first term on the right-hand side of (4.6) describes the above-mentioned channel openings: The channels open only if three subunits are active and not inhibited. The influx rate is determined via the channel current ν = (2eV ) −1 I C = 518.28 µM s −1 (where e = 1.602 × 10 −19 C is the electron charge; the in-flowing ions are assumed to be spread over a volume of V = 1 µm 3 ). This value is also consistent with influx rates extracted from the rising flanks of puffs in our simulations. The exponential decay parameter λ was determined by fitting an exponential decay to calcium puff simulation data. The parameters of the mean-field model are summarized in Table 2.  Concentration c(t)  c(t) returns to its equilibrium value c 0 , while for higher values, it oscillates at a high concentration level. These two states correspond to the "puff" regime and the "always open" regime, respectively. Figure 9(b) shows the corresponding trajectories in a(t)b(t) space. When the trajectory, after its initial transient phase, does not touch the line segment {(a, 2) for a ∈ (3, 4)}, highlighted in gray in Figure 9(b), the concentration will decay to c 0 and the system returns to its original state. If, however, the a(t)b(t) trajectory hits this boundary, the system stays excited, and the channels stay perpetually open. This can be translated into a temporal criterion by viewing the time evolution of the system as a two-step process: (1) Find the time τ such that b(τ ) = 2 (i.e., the time until the second Heaviside function in (4.6) becomes zero); (2) find τ 1 > τ such that b(τ 1 ) = 2, and τ 2 > τ such that a(τ 2 ) = 3. If τ 2 < τ 1 , welldefined puffs are possible; otherwise, c(t) stays elevated and channels stay perpetually open.
We now proceed to find the mean-field phase boundary between the "puff" regime and the "perpetually open" regime. Steady-state analysis of (4.6)-(4.8) does not yield the correct results in the latter regime due to the discontinuous nature of the Heaviside functions in (4.6) and the resulting temporal criterion above. The system rapidly switches between the inhibited state with b(t) > 2 and the noninhibited state with b(t) < 2, leading to the observed decaying oscillations of c(t) in Figure 9(a). Hence, we numerically integrate (4.6)-(4.8) to test whether a given parameter set a i , b i falls into one of the two regimes. We let the system first evolve until b > 2 and continue until either a(t) becomes smaller than 3 or b(t) becomes smaller than 2. We then identify the regime the system is in according to the criterion described above. A bisection algorithm is used to find the boundary b i (a i ) between the two regimes.
The black lines in Figures 6, 7, and 8 show the extracted phase boundaries for the parameter values given in Table 2. 5. Conclusions. In this paper, we have reported a novel application of a particlebased spatial algorithm for diffusion to investigate the influence of diffusive noise on the dynamics of intracellular calcium release. Particle number noise in calcium microdomains has attracted interest in recent studies [22,55,56,32], also for L-type and RyR calcium channels [34,49,28], and it is important to clarify whether calcium diffusion as an additional noise source needs to be incorporated to obtain a better understanding of subcellular calcium signals.
In order to make this study feasible, we split the domain into a compartmentbased regime and a Brownian dynamics regime, coupled via the TRM. This allowed us to model the dynamics of full puffs, including release of a realistic number of ions and inhibition dynamics. We extracted concentration time traces and analyzed the resulting puff statistics. The interpuff time distribution as well as the distributions of puff amplitudes and lifetimes agree qualitatively with experimental data in the literature [23].
We then proceeded to analyze the binding parameter regimes under which welldefined Ca 2+ puffs are possible. We found that, surprisingly, an inhibitory binding site dissociation constant K D = 50 µM, consistent with the literature and patch-clamp experiments [45], does not yield puffs in our model. In this parameter regime, channels stay perpetually open. In order to investigate the transition between well-defined puffs and perpetually open channels, we characterized calcium concentration traces for various combinations of the inhibitory site binding parameters. The phase boundary visible in our data is consistent with a dissociation constant in the region 4 µM < K D < 10 µM. Lower values of the Ca 2+ diffusion constant yield a phase boundary at smaller inhibitory site dissociation rates and thus an even smaller dissociation constant. Given the reliable puff generation and termination in previous studies based on fitted gating models with large K D [43,50], this is an unexpected result. Why does our model not show the same robust termination at large K D as the hybrid approaches [41,43,44]? In the latter models the binding and unbinding to the receptors is stochastic, but the calcium distribution is calculated from deterministic reaction-diffusion equations. To test for the differences of puffs in both models we have performed a set of simulations using large local Ca 2+ concentrations for open channels, similar to what is obtained in the hybrid method. These simulations confirm the robust termination in the hybrid scheme even at large K D . In a hybrid approach, a deep inhibitory state is achieved owing to the large local nanodomain around each open channel displayed in the solution of deterministic calcium equations [41]. Thus, in the hybrid model, a channel is not inhibited by shared calcium in the domain but by its own released calcium. This effect has been termed self-inhibition [50].
In our simulations at large K D , however, there is insufficient inhibition of the channels during the early phase of a puff. This means that no or perhaps only one subunit per channel binds inhibitory calcium, while reliable inhibition requires binding of three or four calcium ions. Our results suggest that diffusive noise mixes calcium in the cluster domain and diminishes the localized domains around open channels and self-inhibition. Thus we are led to a model of Ca 2+ puffs that is very different from the previous hybrid model. Our diffusive model allows inhibition only from a much less localized concentration profile and can therefore only be achieved at a much smaller dissociation constant of inhibitory binding sites. This conclusion is also supported by the mean-field ODE model that we have devised and that captures the average c 2016 SIAM. Published by SIAM under the terms of the Creative Commons 4.0 license binding state of the cluster's activating and inhibitory sites as well as the resulting Ca 2+ concentration in a shared and well-mixed microdomain. This model displays a sharp phase boundary between the two regimes, which agrees well with the data from our spatial simulations.
Evidence for an inhomogeneous or homogeneous calcium distribution in the cluster is hard to obtain directly from our simulations because of the short lifetime of nanodomains. In any case, our study highlights the role of the local calcium concentration in the termination of puffs and shows that puffs are very sensitive to fluctuations of residual calcium remaining after channel closing. It has to be noted though, that, apart from the diffusive noise, there are other differences between the current BD setup and the former hybrid approaches. Notably, these include the presence of calcium-binding buffers and the effect of the terms describing the reuptake of Ca 2+ into the ER [43], and it remains to analyze the extent to which these differences affect puff termination.