{'fraction_non_alphanumeric': 0.028866821425500816, 'fraction_numerical': 0.04313902501934486, 'mean_word_length': 4.878822340156685, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulation of calcium signaling in fine astrocytic processes: effect of spatial properties on spontaneous activity

Astrocytes, a glial cell type of the central nervous system, have emerged as detectors and regulators of neuronal information processing. Astrocyte excitability resides in transient variations of free cytosolic calcium concentration over a range of temporal and spatial scales, from sub-microdomains to waves propagating throughout the cell. Despite extensive experimental approaches, it is not clear how these signals are transmitted to and integrated within an astrocyte. The localization of the main molecular actors and the geometry of the system, including calcium channels IP3R spatial organization, are deemed essential. However, as most calcium signals occur in astrocytic ramifications that are too fine to be resolved by conventional light microscopy, most of those spatial data are unknown and computational modeling remains the only methodology to study this issue. Here, we propose an IP3R-mediated calcium signaling model for dynamics in such small sub-cellular volumes. To account for the expected stochasticity and low copy numbers, our model is both spatially explicit and particle-based. Extensive simulations show that spontaneous calcium signals arise in the model via the interplay between excitability and stochasticity. The model reproduces the main forms of calcium signals and indicates that their frequency crucially depends on the spatial organization of the IP3R channels. Importantly, we show that two processes expressing exactly the same calcium channels can display different types of calcium signals depending on channels spatial organization. Our model with realistic process volume and calcium concentrations successfully reproduces spontaneous calcium signals that we measured in calcium micro-domains with confocal microscopy. To our knowledge, this model is the first model suited to investigate calcium dynamics in fine astrocytic processes and to propose plausible mechanisms responsible for their variability.

Astrocytes were first characterized as non-excitable cells of the central nervous 2 system since, although they express voltage-gated channels [1], they do not exhibit 3 electrical excitability [2]. Astrocytes excitability instead results from variations of 4 cytosolic calcium concentration [3]. At the cellular level, those calcium signals emerge in 5 astrocytes in response to synaptic activity and may cause the release of molecules called 6 Astrocytic calcium signals are considered to rely mainly on the IP3R calcium 23 channel pathway. Indeed, type-2 IP3R calcium channel is enriched in astrocytes [34] 24 and knocking-out IP3R2 channels abolishes all calcium signals in astrocytic soma and 25 roughly half of them in the cell processes [30]. The molecular origin of the 26 IP3R2-independent signals in processes remains a matter of debate, and could involve 27 calcium fluxes through the plasma membrane [30] and/or other IP3R channel 28 subtypes [35]. At any rate, astrocytes respond to G-protein-coupled receptor (GPCR) 29 agonists with calcium transients [36,37]. The binding of agonists to G q/11 -GPCRs 30 activates IP3 synthesis. In turn, the binding of both IP3 and calcium ions to IP3R 31 channels on the membrane of the endoplasmic reticulum (ER) triggers a calcium influx 32 from the ER to the cytosol [38]. The initiation and propagation of calcium signals 33 within astrocytes then relies on the so-called calcium-induced-calcium release (CICR) 34 mechanism: an increase, even small, of the local calcium concentration increases IP3R 35 opening probability thus increasing the probability for local calcium concentration to 36 rise further. 37 80% of the astrocyte calcium activity in vivo take place in the gliapil, which is 38 mostly formed by astrocytic ramifications that cannot be spatially resolved by 39 conventional light microscopy [39], yet account for 75% of the astrocytic volume [40]. 40 According to electron microscopy studies, the perisynaptic astrocyte projections (PAPs) 41 that belong to the gliapil could be as thin as 30-50nm in diameter [41,42]. At this parameter [48][49][50]. In those models, calcium transients emerge as limit-cycle oscillations 59 from a Hopf bifurcation (or a saddle-node on an invariant circle) beyond a critical value 60 of the IP3 concentration. The first astrocyte-specific calcium signaling models arose a 61 decade later. In those models, the IP3 concentration is usually a dynamical variable 62 coupled to calcium but calcium transients still emerge through the Hopf-bifurcation 63 scenario. Notably, those models focused on intercellular IP3 transport within astrocyte 64 networks via gap junctions [51,52]. Stochastic models of IP3R-mediated calcium 65 signaling have also been proposed, that take into account the stochasticity associated 66 with molecular interactions [53][54][55]. Yet, none of those models take into account the 67 stochastic effects associated with the small volumes and the low copy number of 68 molecules involved in fine processes. Recently, individual-based modeling has been 69 introduced to evaluate the impact of diffusive noise on IP3R opening dynamics [56], but 70 this simplified model disregarded IP3 dynamics and restricted stochasticity to the 71 vicinity of the IP3Rs. 72 Here, we propose an IP3R-mediated calcium signaling model adapted to the 73 dynamics of CICR in small spatial volumes corresponding to thin PAPs. To account for 74 the stochasticity inherent to small sub-cellular volumes and low copy numbers expected 75 in fine processes, our model is both spatially explicit and particle-based: each molecule 76 is described individually, diffuses in space through a random walk and reacts 77 stochastically upon collision with reaction partners. The kinetics of IP3R channels is 78 accounted for with a simplified version of the 8-state Markov model on which most of 79 the previous CICR models are based. In order to explore the range of dynamical 80 behaviors that the model can display, we first focus on a 2D version of our model, that 81 is less compute-intensive than the 3D one. Extensive simulations of the 2D model show 82 that spontaneous calcium signals arise in the model via the interplay between the 83 excitability of the system and its stochasticity. The model accounts for various forms of 84 calcium signals ("blips" and "puffs") and their frequency depends on the spatial 85 organization of the IP3R channels. In particular, we demonstrate that the 86 co-localization of sources of calcium influx plays a crucial role in triggering an effect of 87 IP3R clustering on calcium signaling. Finally, as solute concentrations can hardly be 88 defined in 2D, we use a 3D version of the model in order to compare it to experimental 89 data. We show that the spontaneous calcium signals generated by the 3D model with 90 realistic process volume and astrocytic calcium concentrations successfully reproduce 91 the spontaneous calcium transients measured in calcium micro-domains with confocal 92 microscopy in organotypic culture of hippocampal astrocytes. Our model therefore 93 represents the first validated tool to investigate calcium signals in realistic small 94 sub-cellular volumes such as in PAPs, where astrocytes and synapses communicate. 95 This provides a crucial step towards a better understanding of the spatiotemporal 96 response patterns of astrocytes to neuronal activity and beyond, towards 97 astrocyte-neuron communication. 98 Materials and methods 99 Modeling methods 100 Reaction scheme 101 The model considers cytosolic calcium and IP3 dynamics in the framework of 102 calcium-induced calcium release (CICR) signaling. The reaction scheme considered is 103 shown in figure 1 A. In short, we consider calcium fluxes between the cytosol and the 104 extracellular space or the endoplasmic reticulum (ER), including via IP3R channels. We 105 also take into account the effect of phospholipase C δ (PLCδ ), that, when activated by 106 calcium, synthesizes IP3. To derive simple models for this scheme, we made the 107 following assumptions: 108 • We considered that the extracellular and ER calcium concentrations are constant 109 during the simulation, as well as the electrical potentials across the plasma and 110 ER membranes. In this case, calcium outflow from the cytosol to the ER or to the 111 extracellular medium can be lumped into a single first-order rate α. Likewise, 112 calcium entry from the extracellular medium or any IP3R-independent Ca 2+ 113 influx from the ER can be considered constants, too. We lumped them into a 114 single overall constant flux γ. 115 • PLCδ enzymes remain located in the cytosol (no translocation) and the amount of 116 their substrate PIP2 is present everywhere in large excess. Under this condition, 117 activated PLCδ produces IP3 with constant rate δ.

118
IP3R channels are gated both by calcium and IP3, with a bell-shaped dependence of the 119 open probability to calcium concentration [57]. To model their dynamics, we used the 120 classical 8-state Markov model proposed in [48,57], with two calcium binding sites and 121 one IP3 binding site for each IP3R (see figure 1B). However we used the following We first modeled the kinetic scheme described in figure 1 with a lattice-free 142 spatially-explicit stochastic particle-based model, referred to as "Particle-based" model 143 below, in two spatial dimensions, with reflective boundary conditions. Each molecule of 144 the system was explicitly modeled with its associated position in space. To determine the positions of the N IP3R IP3R molecules, we first chose the centers 150 of N c = N IP3R /η IP3R clusters uniformly at random in the reaction space, where η is 151 the number of IP3R per cluster (as illustrated in figure 1 C). For each cluster, we 152 positioned η IP3R molecules uniformly at random within a distance R c of the cluster 153 center, with R c = d IP3R η/0.91, where d IP3R is the interaction distance of the IP3R, 154 i.e. the maximal distance between IP3R center and a Ca 2+ or IP3 molecule below which 155 binding can occur. According to this algorithm, η = 1 corresponds to randomly 156 distributed independent IP3R molecules (no clustering) whereas IP3R molecules become 157 increasingly clustered when η increases, with constant IP3R density within the clusters 158 and constant total IP3R number in the reaction space.

159
Each MC stimulation step (of duration ∆t) consists in iterating the following steps: 160 1. Diffusion. The position of each mobile molecule (Ca 2+ and IP3) is updated 161 independently according to Brownian motion: Gaussian-distributed random numbers with zero mean and unit variance. In a 164 subset of simulations, the new position of each mobile molecule was chosen at 165 random in the reaction volume, i.e. r(t + ∆t) = ζ were ζ is a vector of i.i.d. 166 random numbers uniformly distributed in [0, L], with L the length of the spatial 167 domain. We refer to this setting as "infinite" diffusion coefficients, D = ∞.  Likewise, each Ca 2+ or IP3 molecule close enough to an IP3R molecule (i.e. 172 within its interaction radius) can bind it depending on its occupancy state. If the 173 IP3 binding site is free, an IP3 molecule binds with probability a 2 ∆t. If one of the 174 Ca sites is free, a Ca 2+ ion binds the free site with probability a 1 ∆t (first Ca site) 175 or a 3 ∆t (second Ca site). If both Ca sites are free, binding occurs to the first site 176 with probability a 1 ∆t and to the second one with probability (1 − a 1 ∆t)a 3 ∆t.   influx from IP3R-independent sources in the ER membrane or through the plasma 189 membrane. Note that the position of this new calcium is not uniform over space 190 but depends on parameter R γ and works as follows: an IP3R molecule is chosen 191 (uniformly) at random and the new Ca 2+ ion is positioned uniformly at random 192 within distance R γ of the chosen IP3R. Therefore low values of R γ emulate 193 co-localization between IP3R-dependent and IP3R-independent Ca 2+ influx 194 sources, whereas the location of IP3R-independent Ca 2+ influx is uniform over the 195 reaction volume when R γ becomes as large as the volume side length.
196 Table 1 gives the parameter values used in our 2D simulations, including the initial 197 numbers of Ca 2+ , PLCδ , IP3 and IP3R molecules.
where [P LCδ] = N plc /V . Finally, the mean-field dynamics of the free Ca 2+ is obtained 209 with: [110] (3) For comparison with the output of the other models, the concentrations were 211 transformed into numbers of molecules by multiplication by the reaction volume V .

212
Perfectly-stirred stochastic temporal dynamics (SSA) 213 For comparison, we also modeled the reaction scheme depicted in Figure 1 using 214 Gillespie's exact Stochastic Simulation Algorithm (SSA) that accounts for stochasticity 215 due to low copy numbers and assumes perfect mixing of the reactants [58,59]. Here, the 216 dynamic variables are the number of Ca 2+ and IP3 molecules in the system, N Ca and sampled from an exponential distribution with mean 1/R T , where R T is the sum of the 222 reaction rates of all reactions. The next reaction to occur at time t + τ is chosen as an 223 integer random variable with point probability given by the ratio of its rate to R T . For 224 instance, for the reaction illustrated above, the probability that this reaction is the one 225 state ijk, 0 else. If the illustration reaction described above on receptor l is chosen, this 229 means N l to as the reaction-diffusion master equation (RDME). In RDME, space is partitioned 238 into voxels inside which perfect mixing is assumed, while diffusion between adjacent 239 voxels is modeled as first order reactions [60,61]. STEPS uses a derivative of the SSA in 240 tetrahedral voxels that allows for a better resolution than the cubic voxels mostly used 241 in voxel-based models [62].

242
Geometry. The main advantage of STEPS in the context of the present study is its 243 automatic handling of external and internal membranes [63]. Moreover, STEPS 244 simulations can easily be parallelized [64], a crucial property given the computational 245 burden of such compartmentalized 3d simulations. This allowed us to explicitly describe 246 the presence of the ER membrane inside the 3d cell cytoplasm and the fact that IP3R 247 channels are located in the ER membrane (see figure 6B). The geometry of the reaction 248 volume consisted in a cylinder of length L astro =1 µm and radius R astro =0.1 µm. The

253
Reactions. In this spatial configuration, we modeled the IP3R-mediated calcium 254 signaling kinetic scheme of figure 1 with IP3R channels positioned on the intracellular 255 ER membrane and according to two model variants:

256
• A first variant, referred to as the "No-GCaMP" model, did not include fluorescent 257 calcium indicators. In this 3d model, parameter values were taken, whenever 258 possible, from the literature (Supplementary Table S1). γ and α values were 259 adjusted to yield basal calcium concentration around 120 nM [65,66]. Likewise, β 260 and µ were adjusted for a basal IP3 concentration of 120nM [67]. Note that this 261 value is based on recent, precise measurements of IP3 concentration and differs by 262 an order of magnitude from IP3 concentration values routinely used in 263 IP3R-mediated calcium models [68][69][70]. IP3R density on the ER surface has been 264 measured from TIRF-microscopy analysis in cell cultures [71] Table 2.

281
Simulations code

282
The code of our ODE, Gillespie, Particle-based and STEPS models is available on 283 ModelDB at https://senselab.med.yale.edu/modeldb/enterCode.cshtml?model=247694. 284 The code will be in public access once the paper will be published.

285
Peak detection and analysis 286 Automated peak detection from the model simulations was based on the statistics of 287 baseline calcium trace. A histogram of Ca 2+ trace was built with a bin size of 0.25 ions 288 and the mode of this histogram was used to define baseline calcium. A peak initiation 289 corresponded to the time step where calcium trace overcame a peak threshold defined as 290 baseline + nσ Ca where σ Ca is the standard deviation of the above histogram. The value 291 of n varied between 2 and 4 and was set by hand for each simulation, depending on its 292 signal/noise ratio. The peak was considered terminated when the calcium trace 293 decreased again below peak threshold. This implies that in case of a second calcium 294 peak starting before the first one terminated, both events were considered as being part 295 of the same peak. Peak duration was defined as the time between peak initiation and 296 termination. Peak amplitude was defined as the maximum number of calcium ions 297 reached during the peak duration. In the 3D model, the peak amplitude A was rescaled 298 to facilitate comparison with experimental data, using   Supplementary Table S1. Parameter values in the 3d model have been adjusted to optimize the match with experimental data as described in the Methods section and shown in Fig 6. Note that the values for calcium and IP3 binding or unbinding to IP3R, i.e. the a i 's and b j 's parameters below, are smaller in our model than in the literature, probably because our model is not cooperative. For GCaMP6s, we used the diffusion coefficient of calmodulin. The initial number of Ca 2+ ions was adjusted so that the measured basal GCaMP6s-Ca concentration was around 300nM [65,66].
Parameter  All experiments were performed as described in [80]. We give below the main  For Ca 2+ imaging, we used a custom-built setup based on an inverted microscope 322 body (Leica DMI6000), as previously described in [83]. We used a 1.3 NA glycerol 323 immersion objective equipped with a correction collar to reduce spherical aberrations 324 and thereby allow imaging deeper inside brain tissue [84]. The excitation light was  Spontaneous calcium events were detected and analyzed automatically by ImageJ 335 plugin LC Pro [85] and then manually verified using Igor Pro (Wavemetrics) [35]. For stochastic models, we generated 20 simulations (with different random numbers) 338 and quantified these simulations as mean ± standard deviation over those 20 simulations. 339 One-way ANOVA was performed to investigate the effect of a given parameter on 340 calcium dynamics. Comparison between two simulation conditions were performed with 341 unpaired Student T test if values followed a Gaussian distribution. Otherwise, a 342 Mann-Whitney test was performed. The same method was applied to compare 343 simulation to experimental results. * is for p ≤ 0.05, ** for p ≤ 0.01, *** for p ≤ 0.001. 344

345
Spontaneous oscillations in the model 346 We first analyzed our particle-based model for the CICR signaling system of figure 1. 347 To that end, we compared Monte-Carlo simulations of the particle-based model in two 348 dimensions with the corresponding Mean-Field and Gillespie's SSA models (see Cytosolic calcium can exit the cytosol to the extracellular space or the endoplasmic reticulum (ER) at a (total) rate α, lumping together the effects of ER and plasma membrane pumps. Likewise, Ca 2+ can enter the cytosol from the extracellular space or from the ER via IP3R-independent flow, with (total) rate γ, emulating calcium channels from the plasma membrane. When an IP3R channel opens, calcium enters the cytosol through the channel at rate µ. Phospholipase C δ (PLCδ), once activated by calcium binding, produces IP3 at rate δ. Like Ca 2+ , IP3 can bind the IP3R channel and is removed with rate β. (B ) Our model of the kinetics of the IP3R channel is an 8-state Markov model adapted from [48,57]. Each IP3R channel monomer is associated with 3 binding sites, two calcium binding sites and one IP3 binding site. Occupancy states are designated by a triplet {i, j, k} where i stands for the occupation of the first Ca binding site (i = 1 if bound, 0 else), j for that of the IP3 binding site and k for the second Ca site. The first calcium binding site has higher affinity than the second. The open state is state {110}, where the first Ca and the IP3 sites are bound but not the second Ca site. (C ) Spatial parameters for the particle-based model. The N IP3R IP3R molecules are positioned within uniformly distributed clusters, with η IP3R in each cluster. Hence η = 1 corresponds to uniformly distributed IP3R (no clustering), while the degree of clustering increases with η (for constant total IP3R number). To account for potential co-localization between IP3R-dependent and IP3R-independent calcium sources, the influx of IP3R-independent calcium (at rate γ) occurs within distance R γ of an IP3R. Thus, low values of R γ emulate co-localization between IP3R-dependent and IP3R-independent Ca 2+ influx sources.
model keeps the perfect mixing hypothesis but assumes stochastic kinetics while the with fluctuations of roughly ± 20 ions, large and fast peaks arise spontaneously with a 360 total amplitude between 20 and 120 ions above the baseline. In strong opposition, the 361 (deterministic) mean-field model does not show these oscillations: one gets a stationary 362 trace, that systematically coincides with the baseline level of the stochastic traces 363 (Fig 2B). Comparing the two stochastic models (SSA and particle-based) indicates that 364 both display the same basal calcium level (Fig 2B) and the same frequency and mean 365 peak amplitude (Fig 2C). A Piecewise-Deterministic Markov Process (PDMP) model of 366 the system of Figure 1, in which Ca and IP3 dynamics obeyed to the mean-field model 367 while IP3R state dynamics was expressed with SSA, also exhibited robust spontaneous 368 Ca 2+ peaks (not shown). Altogether, this suggests that stochasticity is necessary for 369 spontaneous calcium signals to occur in this model. 370 We next searched for the dynamical mechanism that gives rise to those spontaneous 371 peaks. A thorough numerical parameter exploration of the mean-field model failed to 372 evidence the existence of Hopf bifurcations or of any other bifurcation that would 373 generate limit-cycle oscillations in the model (not shown). This is a distinctive feature 374 of our model, since spontaneous oscillations in the vast majority of IP3R-mediated 375 calcium signaling models arise from limit-cycle generating bifurcations [48][49][50]. This is 376 however not unexpected since the simplifications made to derive our model significantly 377 reduced its nonlinearity compared to these models, and the emergence of limit-cycle obtained, before relaxation to the stationary state ( Fig 2D). Figure 2D2 shows how the 385 resulting transient amplitude depends on the amount of injected IP3. For low values of 386 IP3R calcium binding rate (first site), a 1 , the calcium response is basically linear with 387 the number of injected IP3: doubling the amount of IP3 injected only doubles the 388 amplitude of the calcium response. However, as a 1 increases, peak amplitude becomes a 389 strongly nonlinear function of the number of IP3 injected. With a 1 = 5 a.u. for 390 instance, doubling the number of injected IP3 from 50 to 100, results in an almost 391 threefold increase of the calcium response. Therefore the mean-field model with large 392 values of a 1 is an excitable system that amplifies the fluctuations of IP3 in its calcium 393 responses. We conclude that spontaneous calcium transients occur in the system of 394 figure 1 through the interplay of the stochasticity of the SSA or particle-based models 395 and the underlying excitability of the system.

396
Transitions between calcium activity regimes 397 The experimental and modeling literature on intracellular calcium signals 398 distinguishes two classes of localized calcium peaks: puffs and blips [86]. Blips refer to 399 brief and weak peaks that correspond to the opening of a single IP3R channel (or a 400 single IP3R channel tetramer), whereas puffs are longer and higher peaks resulting from 401 the concerted opening of a group of nearby IP3R channels (or tetramers thereof), via 402 the calcium-induced calcium-release principle. We next examined whether our model (C ) Quantification of calcium transients in the stochastic models (calcium peak frequency and mean peak amplitude). No significant difference between the two models was observed. (D) Excitability of the Mean-Field model: increasing quantities of exogenous IP3 molecules were injected at time t = 20∆t, after model equilibration. The amplitude of the resulting calcium response (D1 ) was quantified depending on the amount of IP3 injected and the value of the binding rate constant to the first calcium IP3R site, a 1 (D2 ). Parameter values for the particle-based model: D Ca = D IP3 = ∞ (perfect mixing) and η = 1, R γ =200, i.e. no IP3R channels clustering, and no co-localization of IP3R with IP3R-independent Ca 2+ sources. For SSA and particle-based models, the figure shows the average ± standard deviation over 20 simulations.
was able to reproduce these observations. 404 We carried out parameter exploration of the particle-based model in conditions of 405 perfect mixing for mobile molecules (Ca and IP 3) and uniform spatial distribution of 406 the immobile ones (PLCδ , IP3R). As expected, we found that calcium peaks frequency 407 depends on parameter values (Fig 3A). When the rate of calcium influx through open 408 IP3R channels µ or the binding rate constant to the first Ca IP3R site a 1 are too small, 409 the model does not exhibit calcium peaks at all, only fluctuations around a stationary 410 state (Fig 3C ). This is in agreement with our analysis of the systems excitability 411 above, that evidenced excitability only for large enough values of a 1 ( Figure 2D2). Note 412 however that in the model, IP3R openings do not necessarily lead to a calcium peak, 413 especially for low values of both µ and a 1 (Fig 3C ). Spontaneous calcium transients 414 are obtained in the particle-based model beyond threshold of (µ, a 1 ) values, with a peak 415 frequency that increases with parameters values (Fig 3A). Inspection of the maximal 416 number of open IP3R per peak reveals that not only the frequency, but also the type of 417 these transient signals changes with parameters values (Fig 3A): the less frequent 418 signals are generally associated with a single open IP3R per peak (Fig 3C ), 419 corresponding to blips, whereas the high-frequency spontaneous signals rely on the 420 opening of 2 − 12 IP3R in a peak (Fig 3C•), corresponding to puffs. In agreement with 421 experimental observations [68,87], calcium puffs in the particle-based model are  Impact of calcium diffusion coefficient on calcium signals

428
A modeling study has demonstrated the necessity to account for the stochasticity 429 inherent to calcium diffusion when modeling calcium signaling in small volumes [88]. 430 We next investigated the impact of calcium diffusion on calcium dynamics in the 431 particle-based model. In neurons or astrocytes, the amount of endogenous calcium 432 buffers is large so that the diffusion distance of free calcium is believed to be very small. 433 Many of the endogenous buffers are however mobile. Therefore, our diffusion coefficient 434 for calcium is to be interpreted as an effective diffusion coefficient lumping together 435 calcium buffering by mobile endogenous buffer and diffusion of these buffers.

436
Moreover, several plasma membrane proteins, in particular the Na + -Ca 2+ exchanger 437 (NCX) have been observed to co-localize with ER proteins in neurons and 438 astrocytes [89]. Such a co-localization of calcium signaling molecules might imply 439 spatial organizations including raft-like micro-domains. This organization seems 440 essential for calcium wave propagation in astrocytes [90]. Moreover, mGluR5-ER 441 proteins co-clusters mediated by an interaction with Homer1 scaffold protein have been 442 observed in astrocytic processes [91]. Homer1 is also known for increasing calcium 443 activity in neurons by increasing IP3R-mGluR5 proximity [92]. Those experimental 444 studies suggest that several calcium sources are co-localized with ER proteins in 445 astrocytes and that it might alter calcium dynamics. Such a co-localization could be 446 crucial for calcium signaling, in particular in small volumes. We thus placed our study 447 of the influence of calcium mobility on calcium signaling in a framework where calcium 448 sources (IP3R-dependent and IP3R-independent) can co-localize. 449 To this end, the IP3R-independent calcium influx in the cytosol (from e.g. plasma 450 membrane transporters or channels) was made dependent on parameter R γ , that sets 451 the distance from IP3R receptors within which new calcium ions are injected in the 452 cytosol when they originate from IP3R-independent fluxes (see Methods section). When 453 R γ =0, the initial location of the new calcium ion is shared with an IP3R channel 454 whereas when R γ increases, the injection positions of new calcium ions are increasingly 455 uncorrelated from those of the IP3R channels. When R γ becomes as large as the size of 456 the reaction surface (i.e. for R γ =100), the injection position of the new calcium ion is 457 effectively independent of the positions of the IP3R channels. co-localized with IP3R channels (i.e. for R γ ≥ 5), both the peak frequency and the type 471 of signal were found not to depend on the calcium diffusion coefficient anymore. Those 472 results suggest that calcium diffusivity could control the frequency and type of calcium 473 signals within astrocytes when IP3R channels are co-localized with IP3R-independent 474 calcium sources. IP3R clustering controls calcium signals when co-localized 476 Experimental data demonstrate that IP3R in SH-SY5Y and COS7 cells are not 477 uniformly distributed on the ER membrane but form clusters [68,87]. We next 478 investigated the impact of IP3R clustering on calcium signal dynamics in our 479 particle-based model. Simulations were performed with D ca =0.1 and various amounts 480 of co-localization between IP3R channels and other calcium sources (parameter R γ ).

484
R γ = 10). These traces indicate that the frequency and type of calcium signal in this 485 case is heavily dependent on the spatial distribution of IP3R channels: clustered IP3R 486 seem to exhibit much larger peak frequency and slightly more frequent puffs. However, 487 here again this effect is quite mitigated by the amount of co-localization between IP3R 488 channels and the IP3R-independent calcium sources. In particular, the dynamical range 489 of the modulation by IP3R cluster size η (i.e. the ratio between the frequency at η = 50 490 and η = 1) is maximal for intermediate co-localizations (2 ≤R γ ≤ 10) but the calcium 491 peak frequency is hardly dependent on η when co-localization is either very strong 492 (R γ < 2) or very weak (R γ ≥ 50). Increasing clustering also tends to improve the 493 emergence of puffs, although the effect is significant only for strong co-localization 494 (R γ ≤ 2, Fig 5D). We emphasize that in such cases of strong co-localization, the regime 495 of calcium activity (puffs vs blips) changes by simply rearranging the spatial 496 distribution of the IP3Rs, without changing any of the kinetics parameters of the model. 497 Taken together, these simulation results pinpoint the interplay between calcium 498 source co-localization and the degree of IP3R clustering as a crucial modulator of 499 temporal characteristics of the calcium signals and of the signaling regime. In particular, 500 they suggest that in the presence of certain amount of co-localization between IP3R 501 channels and other sources of calcium influx in the cytosol the spontaneous calcium 502 peak frequency can have a large amplitude variation. Within this range of parameters, 503 calcium peak frequency can be finely tuned by the geometry of the colocalization.

504
Simulations in a compartmentalized 3d geometry reproduce 505 spontaneous calcium microdomains signals 506 The above 2d simulations of the particle-based model have the advantage of a good 507 computational efficiency, which makes them suitable for parametric studies with 508 averaging over a number of Monte-Carlo simulations. However, the 2d setting does not 509 facilitate the comparison of the copy number of molecules in the simulations with 510 species concentrations as measured experimentally. Moreover, it is difficult to 511 investigate with a 2d setting the impact of the fact that IP3R channels are specifically 512 localized at the surface of the ER membrane and not freely diffusing in the cytosol bulk. 513 To tackle those questions, we carried out simulations of our model in a more refined 514 three-dimensional setting (Fig 6), in which we could adjust more precisely molecule 515 concentrations, reaction volume and cytosol compartmentalization to what is expected 516 in fine astrocytic processes. We then compared our simulations to experimental 517 measurements of calcium dynamics in microdomains of comparable dimensions in mice 518 hippocampal organotypic culture (Fig 6A).

519
As 80% of calcium activity occurs in astrocytic ramifications that cannot be resolved 520 by optical microscopy [40], astrocytic calcium signaling models must take into account 521 small volumes associated to it. For that purpose, we created the 3d structure mimicking 522 one process geometry shown in Figure 6B. The reaction volume was chosen to match 523 the range of sizes that are within reach of current imaging methods: a 1 µm-long 524 cylinder of 100 nm radius (i.e. a volume around 0.03 fL), inside which we position a 0.75 525 µm-long cylindrical ER with a radius of 30 nm. In this 3d implementation, calcium and 526 IP3 molecules diffuse in the bulk 3D space located between the external (plasma) 527 membrane and that of the ER, while IP3R molecules are distributed uniformly at 528 random over ER membrane surface.  6A3 and E). Notice that these experimental traces correspond to spontaneous 535 signals to the extent that they were measured in the absence of any neuronal nor 536 astrocytic stimulation. In particular, TTX application in this preparation did not alter 537 peak frequency [80].

538
Our first noticeable result is that our model is able to reproduce the emergence of 539 spontaneous calcium peaks of comparable frequency and signal-to-noise ratio ( fig 6C). 540 This result therefore indicates that spontaneous calcium signals can emerge in the fine 541 processes even with a realistic basal calcium concentration of 83 ± 29 nM, which 542 corresponds to only one to two (1-2) calcium ions in the whole cylinder. However, 543 quantification of the signal properties ( fig 6E) shows that the agreement is only 544 qualitative, since the peak amplitude and frequency are slightly too large while peak  6E3). Note that our experimental statistics are tightly 550 associated with the temporal sampling frequency used in the experiments (2 Hz) since 551 very fast calcium events may be accessible only to higher sampling frequencies [40]. In 552 particular, the experimental peak frequency measured might have been higher with 553 better temporal resolution. In any case, our results show that genetically encoded 554 calcium indicators (GECIs), such as GCaMP6s, may change local calcium concentration, 555 in particular close to open IP3R channels, leading to an increased peak duration.

556
Interestingly, the well-mixed implementation of this model in the cylinder geometry, 557 with the same parameter values, displays a higher peak frequency (data not shown).

558
This reflects the importance of stochasticity and local variations of calcium 559 concentration for tuning IP3R activity and thus the overall calcium signal dynamics.

560
Together these results demonstrate that our model, without any endogenous buffers, is 561 enough to reproduce calcium signals within fine astrocytic processes in a quantitative 562 way, making it a powerful tool to investigate calcium dynamics in the small volumes 563 associated with the PAPs. The 3d geometry used for the 3D model is a cylinder of length L astro =1 µm and radius R astro =0.1 µm, with ER as a thinner cylinder inside. The interior volume is roughly 0.03 fL. (C ) Representative simulations of calcium dynamics within the above cylinder with the "No-GCaMP" simulations. The raw signal corresponds to cytosolic free calcium concentration. (D) Representative simulations of calcium dynamics within the above cylinder with the "GCaMP" simulations. The raw signal corresponds to calcium-bound GCaMP concentration. For both simulation types, parameter values were partly taken from the literature and partly adjusted for fitting calcium traces shown in A (reported in Table 2). (E ) Quantitative comparisons of the spontaneous calcium signals measured experimentally (black bars) or simulated with the "No-GCaMP" or "GCaMP" models (white bars). The compared quantities are peaks amplitude in terms of ∆F/F ratio (E1), their frequency (measured in min −1 for each µm 2 area, E2) and duration E3.

565
Recent experimental reports suggest that the complete dependence of cytosolic 566 calcium transients on IP3R2 is only observed in the astrocyte cell body whereas calcium 567 signals measured within astrocytic processes are a mix of IP3R2-dependent and 568 non-IP3R2-dependent calcium signals [21,30]. The identity, subtype and localization of 569 the receptors responsible for non-IP3R2-dependent calcium signals in astrocytes, in 570 particular their processes, are still to be uncovered. However, our study sheds light on 571 the importance of the localization of these various calcium sources. Our simulation 572 results indeed indicate that whenever IP3R channels are (even moderately) co-localized 573 with IP3R-independent calcium sources, e.g. plasma membrane calcium channels, the 574 degree of IP3R clustering and/or the mobility of the calcium buffers will have a strong 575 impact on the frequency and amplitude of the spontaneous calcium signals. In 576 particular, our simulations predict that two astrocyte processes expressing exactly the 577 same repertoire of channels, pumps and receptors but in a different spatial organization 578 (for instance various degrees of clustering or co-localization), can exhibit very different 579 types and properties of spontaneous calcium events. This could result in significant 580 variability of the calcium response of different processes, even from the same cell.

582
During the past few years, fine astrocytic processes have been regarded as devoid of 583 ER [93,94]. This questions the validity of our model, in which the presence of 584 ER-attached IP3R in the process is crucial for spontaneous activity. We however note 585 that a recent EM study has observed that ER dynamically ramified in astrocyte 586 perivascular processes in vivo and detected contact sites between ER processes and 587 plasma membrane, often positioned in apposition to neuronal synapses [95]. Such 588 contiguous membranous juxtapositions would definitely validate the presence of ER in 589 PAPs. Although dynamical ER remodeling has been reported in dissociated astrocyte 590 culture [96], technical limitations have prevented direct investigation of ER localization 591 within PAPs in vivo or in slices. To our knowledge, it is not even clear whether 592 astrocytic ER is continuous or consists in several independent reservoirs.

593
Super-resolution microscopy of cellular ER and mitochondrial dynamics and structure 594 (resolution ≈ 100nm) has recently been developed and could help solve the controversy 595 regarding the presence of ER in fine processes [97,98]. Correlative super-resolution 596 fluorescence imaging and electron microscopy approaches can yield a resolution of less 597 than 50 nm (down to 10nm) [99], which is very promising avenue to PAPs 598 ultrastructure investigation. At any rate, since the IP3R pathway is involved in calcium 599 dynamics, further investigations regarding ER sub-cellular localization, 600 sub-compartmentalization and dynamics are crucial for better understanding astrocyte 601 information processing. Meanwhile, a straightforward extension of our computational 602 model would be to simulate neuronal stimulation-triggered calcium dynamics.

604
IP3Rs are thought to assemble as tetramers, and a recent experimental study 605 suggested that the four subunits of the tetramer must be simultaneously bound to IP3 606 for the tetramer to allow calcium influx, independently of cytosolic calcium or ATP 607 concentrations [100]. Actually, the original IP3R models predicted that subunit 608 cooperativity for calcium binding is also necessary to fit experimental data of IP3R 609 dynamics [48,50] EPSCs, because this inhibitor of V-ATPases inhibits miniature EPSCs by blocking the 630 refill of presynaptic vesicles. However, the impact of bafilomycin bath application on 631 the frequency of spontaneous calcium signals in astrocytes has proven variable (compare 632 e.g. [105] and [106]). In our preparation, bath-application of bafilomycin strongly 633 decreased peak frequency and amplitude [80]. As bafilomycin has a wide range of effects 634 on calcium signaling that is independent of its effect on the refill of presynaptic 635 neurotransmitter vesicles [107,108], we cannot conclude whether those signals are 636 triggered by EPSCs and further investigation is needed to decipher whether the 637 "spontaneous" calcium signals reported in astrocyte processes are due to spontaneous 638 release of presynaptic vesicles or rely on a synapse-independent mechanism inherent to 639 the CICR system.

641
For simplicity, IP3R clustering in our model was considered static during simulation 642 time. Experimentally, though, IP3R clustering might be highly dynamic [109,110].

643
Several molecules can trigger IP3R clustering, including IP3 and calcium 644 themselves [109,110], through a mechanism that may include the lateral diffusion of 645 IP3R on the ER surface [110] or be independent from it [72]. Beyond this IP3R 646 classification into clustered and un-clustered populations, another approach is to 647 quantify single IP3R channels based on their mobility. A recent study on HeLa cells [43] 648 indicates that calcium signals emerge most of the time from immobile IP3R, which are 649 found in apposition to ER-plasma membrane junctions, whereas the mobile IP3R 650 fraction would not be involved in calcium influx. Our simulation results, in agreement 651 with previous IP3R-mediated calcium models [111,112], indicate that IP3R clustering 652 can lead to an increase of the frequency and amplitude of their calcium signals. This 653 result is in contradiction with a previous modeling study that concluded in favor of a 654 reduction of IP3R channel activity upon IP3R clustering [113]. This discrepancy might 655 rely on the different modeling choices. In particular, the model in [113] incorporates a 656 5-state IP3R model derived from ref [114,115]. All of those modeling studies however 657 agree that dynamical IP3R clustering could be a mechanism used by astrocyte processes 658 to modulate their calcium signals. This could provide astrocyte processes with a 659 capacity for information processing plasticity.

661
In our model, the value of the rate constant for calcium binding on IP3R changes the 662 type of spontaneous dynamics (e.g. blips vs puffs) in addition to its characteristics 663 (frequency, amplitude). Experimentally, several post-transcriptional mechanisms can 664 modulate IP3R affinity. For instance, phosphorylation of type-1 and -2 IP3R by 665 cAMP-activated PKA increases the affinity of IP3R to calcium and IP3 [116]. At a 666 larger time scale, the sensitivity of IP3R to calcium is encoded in a sequence of calcium 667 sensor (Cas) region that differs depending on the IP3R isoform [114,117,118]. Since 668 multiple IP3R isoforms seem to be involved in calcium signaling within astrocytic 669 processes [35], they could assemble into a variety of homo-or hetero-IP3R tetramers 670 that would exhibit a range of calcium and IP3 affinity. In addition, immobile or weakly 671 mobile endogenous calcium buffers are responsible for an effective intracellular calcium 672 diffusion that is an order of magnitude slower than free calcium ions [74]. The 673 variability of these endogenous calcium buffers, with various kinetics and various 674 diffusion coefficients, is large [119]. Some of them are overexpressed in hippocampal and 675 striatal astrocytes, possibly in a region-specific pattern [120]. Our simulation results 676 indicate that the value of the effective Ca 2+ mobility also participates in the 677 determination of the type and characteristics of the spontaneous events, thus confirming 678 previous modeling approaches [88]. Such a regional differential expression of the genes 679 coding for endogenous calcium buffers could therefore be involved in the regional 680 variability of astrocytic calcium signalling [121]. Therefore, the IP3R repertoire, the 681 post-transcriptional regulation of IP3R affinity and the differential expression of 682 endogenous calcium buffers could also be potential determinants allowing a range of 683 responsiveness and spatio-temporal characteristics of calcium signals in astrocyte 684 processes.

686
To conclude, we have presented a spatially-explicit stochastic model to investigate 687 intracellular calcium signaling based on CICR in small sub-cellular volumes. Two recent 688 studies proposed models for the simulation of astrocytic calcium signals in 3d with 689 deterministic differential equation models that correspond to cellular volumes large 690 enough to validate a law of large numbers [94,122]. To our knowledge, our model is the 691 first model suited to reproduce spontaneous calcium signals in the finest astrocyte 692 processes, where low copy number and spatial localization effects are expected to be 693 more prominent than in larger volumes. Since these fine processes are thought to be the 694 place of initiation of neuron-astrocyte interactions, we believe this model might be 695 useful to investigate the initiation and spatiotemporal integration of calcium signals in 696 the sponge-like network of the astrocyte process, a prerequisite to understand 697 neuron-astrocyte communication.