Modeling of Astrocyte Networks: Toward Realistic Topology and Dynamics

Neuronal firing and neuron-to-neuron synaptic wiring are currently widely described as orchestrated by astrocytes—elaborately ramified glial cells tiling the cortical and hippocampal space into non-overlapping domains, each covering hundreds of individual dendrites and hundreds thousands synapses. A key component to astrocytic signaling is the dynamics of cytosolic Ca2+ which displays multiscale spatiotemporal patterns from short confined elemental Ca2+ events (puffs) to Ca2+ waves expanding through many cells. Here, we synthesize the current understanding of astrocyte morphology, coupling local synaptic activity to astrocytic Ca2+ in perisynaptic astrocytic processes and morphology-defined mechanisms of Ca2+ regulation in a distributed model. To this end, we build simplified realistic data-driven spatial network templates and compile model equations as defined by local cell morphology. The input to the model is spatially uncorrelated stochastic synaptic activity. The proposed modeling approach is validated by statistics of simulated Ca2+ transients at a single cell level. In multicellular templates we observe regular sequences of cell entrainment in Ca2+ waves, as a result of interplay between stochastic input and morphology variability between individual astrocytes. Our approach adds spatial dimension to the existing astrocyte models by employment of realistic morphology while retaining enough flexibility and scalability to be embedded in multiscale heterocellular models of neural tissue. We conclude that the proposed approach provides a useful description of neuron-driven Ca2+-activity in the astrocyte syncytium.

Neuronal firing and neuron-to-neuron synaptic wiring are currently widely described as orchestrated by astrocytes-elaborately ramified glial cells tiling the cortical and hippocampal space into non-overlapping domains, each covering hundreds of individual dendrites and hundreds thousands synapses. A key component to astrocytic signaling is the dynamics of cytosolic Ca 2+ which displays multiscale spatiotemporal patterns from short confined elemental Ca 2+ events (puffs) to Ca 2+ waves expanding through many cells. Here, we synthesize the current understanding of astrocyte morphology, coupling local synaptic activity to astrocytic Ca 2+ in perisynaptic astrocytic processes and morphology-defined mechanisms of Ca 2+ regulation in a distributed model. To this end, we build simplified realistic data-driven spatial network templates and compile model equations as defined by local cell morphology. The input to the model is spatially uncorrelated stochastic synaptic activity. The proposed modeling approach is validated by statistics of simulated Ca 2+ transients at a single cell level. In multicellular templates we observe regular sequences of cell entrainment in Ca 2+ waves, as a result of interplay between stochastic input and morphology variability between individual astrocytes. Our approach adds spatial dimension to the existing astrocyte models by employment of realistic morphology while retaining enough flexibility and scalability to be embedded in multiscale heterocellular models of neural tissue. We conclude that the proposed approach provides a useful description of neuron-driven Ca 2+ -activity in the astrocyte syncytium.

INTRODUCTION
Astrocytes of the cortical and hippocampal gray matter are important actors in a number of information processing processes, including synaptic plasticity, long-term potentiation, and synchronization of neuronal firing (Haydon, 2001;Lee et al., 2014;De Pitta et al., 2016;Poskanzer and Yuste, 2016) as well as in coupling neuronal activity to blood flow changes (Otsu et al., 2015). Recent evidence converges on a close connection of these functions with whole-brain processes and systemic regulation pathways. Thus, astrocytes respond to and are able to regulate systemic blood pressure (Marina et al., 2020); they significantly (up to 60%) change their volume during sleep or under anesthesia (Xie et al., 2013); astrocytes play an important role in the clearance of beta-amyloids, a process with mechanisms that are now being actively discussed (Iliff et al., 2012;Abbott et al., 2018;Semyachkina-Glushkovskaya et al., 2018;Mestre et al., 2020); both intracellular and network-level activity of astrocytes are significantly different in sleep and during wakefulness, and activates with locomotion (Bojarskaite et al., 2020;Ingiosi et al., 2020;McCauley et al., 2020). It is important to note that many of the mentioned astrocyte functions are not directly related to neural activity, but are governed by their own regulatory pathways (O'Donnell et al., 2015). Some of these functions are tightly linked to dynamic regulation of astrocyte morphology and volume and depend, for example, on the circadian rhythm of aquaporin expression (Hablitz et al., 2020).
In summary, this frames a new mindset for understanding the function of astrocytes and at the same time poses a challenge for modeling studies. Namely, the morphological features should now be considered as a specific control parameter that significantly contribute to the both single-cell dynamics and network activity patterns. This problem breaks down into three specific tasks: (i) to provide tractable, but still biologically reasonable mathematical account for contribution of subcellular morphological features to intracellular calcium dynamics; (ii) to further develop approaches to modeling of Ca 2+ dynamics on data-driven irregular structures, both for an individual cell and for a network; (iii) to reveal how realistic morphological features are manifested in the spatiotemporal patterns of the calcium dynamics.
We address these tasks in more detail in the rest of the Introduction.

Calcium Signaling in Astrocytes
With plasma membranes enriched in a variety of potassium channels and lacking voltage-gated sodium channels, astrocytes are not electrically excitable (Verkhratsky and Nedergaard, 2018). On the other hand, they display a rich repertoire of Ca 2+activity at multiple spatial and temporal scales (Lind et al., 2013;Volterra et al., 2014;Wu et al., 2014;Bindocci et al., 2017). Although astrocytic Ca 2+ transients can occur spontaneously, their frequency is modulated by neuronal activity (Stobart et al., 2018), changes in local tissue oxygenation Marina et al., 2020), and other factors (Semyanov et al., 2020). As outputs, Ca 2+ -activity in astrocytes leads to release of signaling molecules: gliotransmitters, such as GABA, D-serine, and glutamate, as well as vasoactive metabolites (Serrano et al., 2006;Henneberger et al., 2010;Bazargani and Attwell, 2016). This has been summarized in a concept of "tripartite synapse, " i.e., sensing of synaptic neurotransmitter release by perisynaptic astrocyte processes, encoding this information in Ca 2+ signals and response with secretion of neuroactive molecules (Araque et al., 2014). There is still however an ongoing debate on the mechanisms involved in generating Ca 2+ transients in astrocytes and the extent of effect of astrocyte-derived molecules on synaptic plasticity, e.g., on LTP (Fiacco and McCarthy, 2018;Savtchouk and Volterra, 2018).
Recent experimental evidence obtained with genetically encoded or pipette-loaded Ca 2+ indicators (Tong et al., 2013;Rungta et al., 2016) heralds functional segregation between the less frequent global internal store-operated Ca 2+ transients at the level of cell soma and primary branches, and the more frequent spatially limited microdomain Ca 2+ transients in the thin mesh of astrocytic leaflets-ramified nanoscopic processes, also known as perisynaptic processes (PAPs) due to their proximity to synaptic connections between neurons. The transients located in the leaflets primarily rely on influx of Ca 2+ through plasma membrane, in part because of the high surface-to-volume ratio in this region and in part because the leaflets are often devoid of organelles including ER (Patrushev et al., 2013) and thus can not support exchange with intracellular stores.
The coupling from synaptic activity to local Ca 2+ transients in PAPs and from the latter to global Ca 2+ events is an area of active research. As reviewed in Savtchouk and Volterra (2018), early works attributed this to activation of G-protein coupled receptors to glutamate, but later this pathway has been put to question due to apparent lack of mGluR5 receptor expression in adult astrocytes. Alternative sources of microdomain transients have been proposed, such as via TRP channels (Shigetomi et al., 2011), from mitochondria (Agarwal et al., 2017), etc. One plausible alternative causal pathway can be formulated as follows (Rojas et al., 2007;Verkhratsky et al., 2012;Kirischuk et al., 2016;Parpura et al., 2016): neurotransmitters, released from the presynaptic membranes, primarily glutamate, but also GABA, are cleared from the extracellular space by astrocytic transporters utilizing Na + gradient to drive the neurotransmitters into the cell. This leads to build up of Na + ions in the cytosol, which can lead to temporary reversal of Na + /Ca 2+ -exchanger allowing for Ca 2+ entry via this transporter. Conceivably, if this local Ca 2+ influx happens near the ER and coincides with an increase in inositol trisphosphate (IP 3 ) production by phospholipase C, it can trigger Ca 2+ -induced release of Ca 2+ from intracellular stores via IP 3 receptors (IP 3 Rs) of the ER.
The release of calcium from ER is spatially inhomogeneous due to the non-uniform, clustered, distribution of IP 3 receptors (Smith et al., 2009;Taufiq-Ur-Rahman et al., 2009;Ross, 2012), with clusters spaced at about 0.5-5 µm apart. At a detailed level, calcium release from the receptor clusters has a stochastic character. The effect of the stochastic activation of IP 3 R clusters on the calcium dynamics has been investigated by Shuai and Jung both in point and distributed models Jung, 2002, 2003). In the case of a large enough number of clusters, Ca 2+ release events can be averaged to a lumped deterministic description. Particularly, the increase in IP 3 level transforms stochastic calcium increases into regular waves.
Recapitulating, calcium signaling mechanisms are inhomogeneous across the cell and depend on local morphological parameters, which has to be taken into account in modeling. It seems practical to introduce a metaparameter to describe the relative inputs of store-related and plasma membrane-related Ca 2+ pathways. This metaparameter can reflect local surface-to-volume ratio or the dominant size of processes and can empirically be linked to the astrocyte cytoplasm volume fraction parameter, which can be estimated directly from fluorescent images.

Cell Morphology and Network Connectivity
Astrocytes have intricate and highly complex morphology, which raises computational issues and demands an elaborate approach to modeling. The contribution of the astrocytic spatial segregation and coupling to brain physiology and functions is still not sufficiently understood, especially taking into account that astrocyte-to-neuron and astrocyte-to-astrocyte interaction mechanisms are diverse and depend on brain region. The existence of intercellular Ca 2+ waves traveling across the network of astrocytes suggests a distinct mechanism for long-distance signaling (Cornell-Bell et al., 1990) and plasticity, which operates in parallel to and at much slower time scales than neuronal synaptic transmission (Pirttimaki and Parri, 2013;Sims et al., 2015).
The size of cliques of cortical astrocytes coupled within a local network is estimated around 60-80 cells (Haas et al., 2006;Houades et al., 2006Houades et al., , 2008, but several networks can also connect via a limited number of "hub" astrocytes (Carmignoto, 2000). The implications of inter-astrocyte connectivity have been analyzed in a modeling study by Lallouette et al. (2014) with the main conclusion that sparse short-range connections can promote Ca 2+ wave propagation along the network. This allows to conjecture that once initiated, a wave of excitation can propagate over long distances in the brain cortex and affect (activate or inhibit) postsynaptic neurons at distant synaptic terminals, although most Ca 2+ events are confined to a single astrocyte spatial domain. Propagating calcium waves can travel distances of more than 100 µm with speed from 7 to 27 µm/s in culture and brain slices (Dani et al., 1992). However, the waves observed in vivo rarely spread more than 80 µm (Hoogland et al., 2009;Brazhe et al., 2013), although this observation can be influenced by imaging protocol, as Kuga and colleagues reported large-scale Ca 2+ glissandi in vivo that were only observable under low laser intensity .
It follows that for meso-scale problems related to brain tissue physiology, it is computationally cumbersome to build a groundup model starting from individual processes. We propose a more pragmatic approach based on texture-like volume segmentation to classes such as "soma, " "large branches, " and "gliapil" or a mesh of unresolved thin processes. This rasterization radically simplifies model implementation and scales to large networks. At the same time, by defining morphology-based spatial distribution of a metaparameter, one can study the effects of spatial heterogeneity at different scales. Indeed, the spatial distributions used for simulations are ideally data-driven. Because it is not always possible to infer the astrocytic network structure or even individual domain boundaries from experimental data, and because the networks can be variable anyway, it seems inviting to generate variable astrocytic tilings from images of individual cells.

Modeling Studies
Models of IP 3 -mediated Ca 2+ oscillations have been extensively reviewed both in general (Dupont et al., 2011) and in application to astrocytes (Riera et al., 2011;Manninen and Havela, 2017;Oschmann et al., 2017a), which included both point-and spatially extended models. In particular, the De Young-Keizer model stemmed several currently popular models of Ca 2+ dynamics in literature. This model allows to simulate IP 3sensitive calcium dynamics in cytoplasm and ER occurring at the constant level of IP 3 including also a variant of the model with the positive-feedback mechanism of Ca 2+ on IP 3 production (De Young and Keizer, 1992). Li and Rinzel (1994) reduced De Young-Keizer model to a two-variable system introducing the experimentally observed time scale difference between fast and slow inactivation of IP 3 receptor by Ca 2+ . Adding the dynamics for [IP 3 ] with synthesis dependent on activation of metabotropic glutamate receptors and [Ca 2+ ] degradation leads to a three variable model (Ullah et al., 2006). Also building on Li-Rinzel model and providing a more detailed description of IP 3 degradation, De Pitta and coauthors proposed a three-variable model for glutamate-induced intracellular calcium dynamics caused by the synaptic activity in astrocytes (De Pitta et al., 2009).
One of the first models for intercellular propagation of calcium waves has been described in (Sneyd et al., 1994) by adding diffusion of IP 3 and cytosolic Ca 2+ to the two-pool Ca 2+ model. The effect of Ca 2+ diffusion rate on spatiotemporal patterns of Ca 2+ signaling was studied by Shuai and Jung (2003) in a lattice-based model. Later, Kang and Othmer (2009) regarded networked astroglial Ca 2+ signaling in a 2D model using spatial patterns in form of sparsely connected irregularly branching cells with simplified morphology. Both intracellular diffusion of IP 3 via gap-junctions and extracellular purinergic signaling was regarded as a mechanism of intercellular communication in Kang and Othmer (2009), Edwards and Gibson (2010); intercellular Ca 2+ diffusion was however disregarded in most modeling studies, e.g., Ullah et al., 2006;Kang and Othmer, 2009;Edwards and Gibson, 2010, primarily based on the notion of a much faster diffusion of IP 3 than Ca 2+ , see Allbritton et al., 1992, and small permeability of gap junctions to Ca 2+ . More recently Savtchenko et al. (2018) suggested an advanced NEURON-based modeling environment for detailed spatially extended models of astrocytes. However, they did not address full calcium dynamics models or morphology-defined variations of mechanisms. Specifically, the relative weights of plasma membrane-dependent mechanisms (IP 3 synthesis and Ca 2+ influx) and store-dependent mechanisms scale with astrocytic process morphology, as defined by surface to volume ratio, cytoplasm volume fraction and the physical presence of ER in the process. This has been studied in point-models by Oschmann et al. (2017b) and in 1D extended model by Wu et al. (2018). Recently, Brazhe et al. (2018) studied the implications of the spatial segregation between IP 3 synthesis and plasma membrane exchange and the IP 3 -mediated ER exchange in discrete spatial templates of variable complexity.
The tripartite synapse concept and the computational role of astrocytes in neural network activity has early attracted the attention of modeling studies, pioneered by papers by Jung (2004, 2007). Understanding of the tripartite synapse from the viewpoint of non-linear dynamics and functional models have been developed in works of Postnov et al. (2007Postnov et al. ( , 2008Postnov et al. ( , 2011, and Tewari and Majumdar (2012). Later, tripartite synapses have been adopted in more formal neural network models (Alvarellos-González et al., 2012;Sajedinia and Hélie, 2018;Lenk et al., 2020). Because there is still no consensus based on experimental evidence on mechanisms of Ca 2+ transients in PAPs and gliotransmission effects (Savtchouk and Volterra, 2018), it is hard to formulate a comprehensive model that would include all conceivable pathways and still remain tractable. While at this stage refraining from closing the loop from astrocytes to neurons, we believe it is important to understand the spatiotemporal patterning of astrocytic Ca 2+ signaling at levels from microdomains to networks.

The Proposed Modeling Approach
Our main motivation in this study is to learn if uncorrelated background synaptic activity, when sensed by astrocytes, will be shaped into morphology-defined patterns of Ca 2+ signaling. We present a model of multi-cellular network of astrocytes based on realistic spatial templates. We start from a single-cell model, which is considerably simpler than in Savtchenko et al. (2018), allowing for smaller computational costs, and move on to connect separate cells together to obtain a network model.
We focus on the implications of the morphology-dependent spatial segregation of the Ca 2+ signaling mechanisms between astrocytic leaflets and branches. We follow the lines set out in Brazhe et al. (2018) toward more realistic and larger scale spatial templates, ranging from single astrocytes to networks. In contrast to astrocytes in culture or retina in vivo, cortex astrocytes are non-flat and occupy some volume in 3D space. Nevertheless, we chose to reduce dimensions to 2D and flatten astrocyte images used as spatial templates. One reason for this was to reduce computational cost, especially when addressing network models. Another reason that most existing Ca 2+ imaging data are obtained as time series in single focal plane, and the experimentally obtained dynamics is confined to flat 2D anyway. The work of Bindocci et al. (2017) demonstrated the richness of Ca 2+ dynamics within the whole astrocytic domain in 3D, but volumetric imaging is not yet widely used in the context of astrocytes. We therefore contemplated that using 2D templates for simulations would not restrict us from observing diverse and physiologically relevant Ca 2+ signaling patterns, event if real astrocytes have more degrees of freedom. The rest of the paper is organized as follows: we start from a description of the proposed model in a top-down order: the general concept is followed by proposed algorithm of creating spatial templates for modeling and then continues with description of the differential equations for dynamics of intracellular and ER Ca 2+ , intracellular IP 3 , and extracellular glutamate concentrations. Having defined the model, we test its plausibility on single-astrocyte templates and after quantification of Ca 2+ event statistics we proceed to behavior of astrocyte networks, where we observe noise-driven regular activation patterns.

Model Design and Overview
In this work we aim to conceptualize our current understanding of spatial organization of the astrocytic Ca 2+ dynamics in a form of a spatially detailed model of individual and networked astrocytes excited by stochastic background neuronal activity. In the light of the striking differences between Ca 2+ signaling in astrocytic leaflets and thin processes on the one hand and global somatic signaling on the other, we start with segregation of the modeling space into three major classes as shown in Figure 1: astrocyte soma with thick branches (I), a mesh of astrocytic thin processes (II) and extracellular space (III). The continuum between the two extreme classes I and II is defined as local fraction of astrocytic cytoplasm volume (AVF) and a related parameter-local surface-to-volume ratio (SVR) of the astrocytic processes. In extreme class I regions, such as soma, Ca 2+ dynamics are dominated by exchange with intracellular stores, and a unit of modeled space (template pixel/voxel) contains only astrocyte, while in extreme class II regions (leaflets), Ca 2+ dynamics is dominated by exchange with plasma membrane and each modeled pixel contains a mesh of extremely thin astrocyte processes tangled with neuropil. We thus define a mapping of each pixel in the spatial model template to either class III (no astrocyte) or to a continuous variable between the extreme cases of class I and II with implications in local calcium dynamics and diffusion.
With regard to the local calcium dynamics, the extreme complexity and sheer number of cellular pathways involved, makes the detailed and comprehensive modeling of every Ca 2+related mechanism extremely challenging. Not surprisingly, there is a substantial body of published models that aim to account for the essential features of calcium dynamics in astrocytes, which do not completely agree with each other (Manninen and Havela, 2017). To choose the best model we build upon a model proposed by Ullah et al. (2006) as a prototype, while other models could fit in the proposed approach as well, for example the "ChI" model by De Pitta et al. (2009), which is similar to that of Ullah et al.

Spatial Structure
To represent astrocyte networks with realistic geometry of the regions I-III, one needs to create such templates algorithmically or, alternatively, obtain them from experimental data. Each of the two variants has its benefits and drawbacks. To provide just two examples, the experiment-based approach was employed in Wallach et al. (2014) and an algorithmic creation of network templates was employed in Postnov et al. (2009). Here, we draw advantages from both approaches by suggesting a simple stochastic data-driven algorithm to create realistic surrogate spatial templates of astrocyte networks. Specifically, we use experimental images of astrocytes obtained from a public database, and arrange network structure using Voronoi partition and simple geometrical transformations (see section 2.2.1).

Neuronal Activity
We assume that astrocytic Ca 2+ response to local neuronal activity is primarily driven by the transporter-mediated uptake of neurotransmitters released from presynaptic membranes. One of the possible coupling mechanisms is the reversal of the Na + /Ca 2+ -exchanger transport due to an increase in [Na + ] allowing for a Ca 2+ influx. Here, we sacrifice biophysical details in favor of model simplicity and assume that astrocyte calcium dynamics is excited directly by glutamate released from the presynaptic terminals, causing transient fluxes of Ca 2+ through the plasma membrane. A typical cortex astrocyte is associated with 300-400 individual dendrites and is in contact with about 10 4 -10 5 synapses (Bushong et al., 2002;Halassa et al., 2007). Judging by these numbers and taking into account sparsity of neuronal signaling in the cortex it seems safe to treat each pixel in the distributed model template as associated with a single or just a few individual synapses. For as long as we are not focused on information processing in the cortex, we can assume independent stochastic nature of spiking activity in any of the presynaptic units and describe local activity only statistically, neglecting any complex spike timing patterns. Consequently, we describe the synaptic glutamate drive to the model in each pixel as triggered by presynaptic spike trains drawn from independent homogeneous Poisson process ξ p (t) with intensity p Hz.

Data-Driven Network Generation
Astrocytes, like neurons, have complex morphology. Ideally, an algorithm to create spatial templates should provide means to "grow" realistic branching 3D shapes of astrocytes from a set of randomly placed "seed" locations. Indeed, there are many experimental and modeling studies of the branching patterns for various types of neurons (Ascoli et al., 2007;Donohue and Ascoli, 2008;Cuntz et al., 2010;Polavaram et al., 2014), providing means for creation of realistic surrogate shapes of as many neurons as needed. However, unlike neurons, there is less data available on the statistics of astrocyte branching, which makes it harder to create surrogate spatial templates of astrocyte networks. This hindrance can be circumvented by using a public database of microscopic images of cortical and hippocampal astrocytes (Martone et al., , 2008. To create a library of realistic spatial templates for individual cells, we downloaded a set of 27 fluorescent confocal 3D stacks of hippocampal astrocytes (4-week old rats, microinjection loaded with lucifer yellow in acute slices) (Bushong et al., 2004). The stacks have average lateral resolution of ≈ 0.07 µm/px and vertical (Z-axis) resolution of 0.2 µm; there are 45-60 Zplanes in each stack, thus encompassing the thickness of about 10 µm along the Z-axis. Because our model is set in twodimensional space, the stacks were flattened along the Z-axis by max-projection, Figure 2. The projections were downsampled 4× before simulations, resulting in lateral resolution of ≈ 0.28 µm/px. Each of the experimental astrocyte images then serves as a progenitor of randomized offsprings obtained by applying 250 random rotations (from 0 • to 360 • , shearings and stretchings (within ±20% of original size, uniform distribution), which results in a collection of 6,750 randomized pseudoexperimental astrocyte templates, used to tile the model space. Such data set expansion from a limited number of "real-world" objects is a popular approach in machine learning (Simard et al., 2003;Krizhevsky et al., 2012) helping to prevent overfitting and providing for transformation-invariant feature learning.
Inspired by the fact that astrocytes establish distinct nonoverlapping territories, we employ an algorithm based on Voronoi partitioning and active contours to tile the model space with astrocytes. First, we create a lattice of "seed points" regularly spaced at some intervals corresponding to average cell density, typically around 50 µm, shown in light gray in FIGURE 2 | Algorithm to create surrogate templates of astrocyte network. First, a set of seeding points on a regular grid (light-gray) is perturbed with random shifts (dark-gray points). Voronoi diagram is then drawn for these points (blue lines). Each patch in the Voronoi partitioning is then filled with the best shape-matching template from an augmented collection of astrocyte images. The lookup collection is created from a set of experimental images taken from CCDB (Martone et al., , 2008 by applying multiple different random rotations and shears to each experimental image.

Figure 2.
The resulting regular grid is then deformed by jittering x and y coordinates of every point by a random displacement value drawn from Gaussian distribution with σ = 10 µm (dark gray points in Figure 2). Different values of spacing and jitter can be used, the ones used here tended to give the most realistic tiling results. Next, a Voronoi diagram, which for each seed point delineates territories closer to it than to any other seed point, is drawn for the jittered points. We then iteratively pick a polygonal area patch from the Voronoi partitioning, look up a template astrocyte from the randomized collection, with a convex hull best matching the shape of the given Voronoi patch, and place this template into the model space. Repeated for all patches in the Voronoi partition, this creates a preliminary tiling with partially overlapping domains of neighboring astrocytes and occasional empty spaces. Next, this draft tiling is optimized with an active deformable model: the perimeter of each cell template is treated as an elastic two-dimensional curve, which optimizes an energy functional designed to promote repulsion between overlapping regions and adhesion between neighboring cells, with a penalization of the major cell shape deformation. After all domain boundaries are settled, the spatial templates are interpolated into the deformed contours. The described process of the network template creation is visualized in Supplementary Video 1.

Computational Design
Our simulations are based on compiling an encoded raster image representation (a template) of the model space to regionspecific equations. For the sake of computational simplicity, we use two-dimensional spatial layout-each pixel of the spatial template can be interpreted as a thin slab, occupied either exclusively (e.g., in the soma) or partly by astrocyte cytosol; or as belonging to extracellular space. As follows from this approach, each pixel in the model space has to be assigned to either astrocyte-free space (class III) or astrocyte-occupied space, ranging from class I, astrocyte soma and thick branches, to class II, i.e., elements of volume containing a tangle of thin astrocyte processes and unresolved neuronal structures, e.g., synaptic boutons. We account for a graded transition from thick branches to thin processes to leaflets by introducing a local astrocyte volume fraction (AVF) parameter, which defines the landscape of how much of each pixel volume is occupied by astrocyte in the 3D prototype. AVF here is defined as a ratio between local fluorescence intensity of the template and the intensity at the soma AVF = max(I/I max , AVF min ). An example of the described mapping from image intensity to AVF is shown in Figures 3A,B, where the colormap in Figure 3B is such that the non-zero blue channel delineates the presence of astrocyte cytoplasm (non-zero AVF), and the intensity of the red channel encodes the AVF value. To describe the relative input of the storeoperated calcium flux and plasma membrane flux, we introduce a surface-volume ratio (SVR) parameter, which inversely depends on AVF ( Figure 3C). The SVR value is maximal at the edges of the leaflets and minimal in the soma. Accordingly, a simple raster RGB image serves as a spatial template to encode the model space. Specifically, non-zero values in the blue channel define astrocyteoccupied pixels, while intensity in the red channel encodes AVF and ranges from minimal value AVF min (class II) to 1 (class I). Thus, one can set up computation for a specific spatial template by simply drawing it algorithmically or with an indexed palette using a graphical editor.
At each integration step the master program module optionally compiles the provided image into a set of equations by mapping each pixel color to equation set following the colorcoded dictionary. For each pixel first the point dynamics are applied, i.e., right-hand terms are evaluated. Next, diffusion of ions and molecules and any other short-range interactions is taken into account based on the class of the neighboring pixels. This approach is flexible, but has an overhead of compiling the color-to-equation mapping. To improve the computational performance we employ NVIDIA CUDA, a parallel GPU-based computing technology.

Intracellular Calcium Dynamics: Principal Quantities and Flows
The model for local Ca 2+ dynamics is based on that of Ullah et al. (2006) with a few modifications previously introduced in Brazhe et al. (2018), which sum up to treating ER calcium as a dynamic variable, adding neurotransmitter-dependent calcium influx via plasma membrane, and segregation between thick and thin processes. Below we describe the proposed model, focusing on the differences with the Ullah et al. (2006)  To account for the morphology-based spatial heterogeneity of astrocytes, we introduce a parameter r ∈ (0 < r min . . . 1]-a scalar quantity, roughly representing local AVF. This parameter also defines a linked parameter s representing local SVR of the astrocyte processes; SVR is inversely related to AVF: s = 1/(1 − exp[0.1(r − 0.5)]). SVR scales relative inputs of Ca 2+ exchange through plasma membrane and with ER as described below, while AVF scales effective diffusion coefficients for Ca 2+ and IP 3 .
The equation set for the principal variables reads: where J ER is the total flow of calcium ions in exchange between the cytosol and endoplasmic reticulum; J pm is the total flow of calcium ions through the plasma membrane in exchange between the cytosol and extracellular space; I Glu and I Ca stand for glutamate-and calcium-dependent inositol trisphosphate production mechanisms, mediated by phospholipases β and δ; I eq is a simplified first-order equilibration of inositol trisphosphate concentration to the basal level [IP 3 ] 0 ; [Glu] amb is the ambient concentration of extracellular glutamate, and τ Glu is the timescale of its clearance and return to the baseline level; ξ p (t) is stochastic source of glutamate from nearby located neuronal synapses triggered by Poisson spike trains in each pixel. Additionally, J diff , I diff , and G diff describe the finite-element approximation of diffusion of cytosolic Ca 2+ , IP 3 and extracellular glutamate, respectively and are described below. The weighting coefficient s accounts for the stratification of intracellular dynamics according to Figure 3: in the leaflets r = r min ≪1 and s ≈ 1, while for deep cytosol locations r = 1 and s ≈ 0. With this we (i) describe that input from all plasma membrane calcium currents is maximal in the leaflets and (ii) assume that endoplasmic reticulum does not invade leaflets much, thus ER exchange is large only in thicker branches and soma. We also assume that all IP 3 is produced in the plasma membrane by means of G-protein coupled phospholipase β or Ca 2+ -dependent phospholipase δ.

Calcium Exchange Between the Cytosol and ER
Total calcium flow across the endoplasmic membrane is composed of IP 3 R-mediated current J IP 3 , leak of Ca 2+ from endoplasmic reticulum J leak , and contribution of ER membrane Ca 2+ pump J pump : Ca 2+ current via IP 3 receptors J IP 3 is modeled in the same way as in Ullah et al. (2006) (Equations 2, 4, 5, 9-12). Two other terms in Equation (5) stand for the leak of calcium from ER, and for the pumping it back, respectively, following Equations (3,6) in Ullah et al. (2006).

Transmembrane Calcium Flows
Transmembrane calcium exchange J pm consists of three flows: where J in describes the sum of background constant Ca 2+ influx and agonist-dependent IP 3 -stimulated Ca 2+ influx across plasma membrane from the extracellular space and J out is an extrusion current (eqs. 7-8 in Ullah et al.); J Glu = γ [Glu] describes the direct effect of extracellular glutamate on additional calcium influx.

Inositol Trisphosphate Turnover
The dynamics for IP 3 concentration from Equation (3) second, we account for the Ca 2+ -stimulated IP 3 production in the same way as in Ullah et al. (2006), (eq. 14), and third, we account for glutamate-driven IP 3 production I Glu following Ullah et al. (2006), (eq. 15).

Synaptic Glutamate Drive
Stochastic glutamate source ξ p (t) in each pixel is modeled as quantal release triggered by a spike train drawn from a homogeneous Poisson process with intensity p, which agrees with statistics of neuronal firing (Softky and Koch, 1993). Accordingly, the ξ p (t) term in Equation (4) is given by where A is the instantaneous increase in glutamate release rate associated with each presynaptic event and t k are times of presynaptic spikes in the given pixel following Poisson process with intensity p syn .

Intracellular Diffusion
Elevated cytoplasmic Ca 2+ can remain confined to the spatial domain of a single astrocyte, but can also spread to the neighboring astrocytes (Nedergaard, 1994;Carmignoto, 2000;Falcke, 2004) in a wavelike manner. The involvement of a large number of cells into a wave is still not fully understood though it may be an important aspect of information processing in the brain (Haas et al., 2006). At least two main mechanisms can account for the intercellular wave propagation: (i) secretion and diffusion of extracellular ATP and its action on P2Y receptors on astrocytic membranes and (ii) diffusion of intracellular IP 3 and Ca 2+ between contacting astrocytic leaflets, via gap junctions. Relative input of the two mechanisms differs across brain regions and for the cortical astrocytes the one mediated by the gap junctions has been reported to prevail (Haas et al., 2006). The current work is therefore focused on the latter mechanism. Accordingly, astrocytes in the model are networked by an analog of gap junctions dispersed over the parts of the cell perimeter and simulated as the connection of areas with low AVF.
Here we employ a rather simplified description of diffusion in the cytoplasm where the region occupied by astrocyte is considered as a continuous space. As corroborated by evidence for autologous gap junctions between the processes of the same astrocyte (Wolff et al., 1998;Nagy and Rash, 2003;Genoud et al., 2015), astrocyte cytosolic volume can be described as a porous sponge-like medium rather than a branched structure or acyclic graph. Thus, possible hindrance to IP 3 or Ca 2+ diffusion through the intricate mesh of astrocytic processes due to tortuosity and porosity of the astrocytic volume can be accounted for by a simple scaling of the apparent diffusion coefficient. Though an interesting issue, a detailed account for intracellular diffusion and connectivity between neighboring points in an astrocyte is outside the scope of the current study and here we resort to a rather minimalistic description.
Following the study of diffusion coefficients of IP 3 and Ca 2+ in Xenopus laevis oocytes (Allbritton et al., 1992), most modeling studies assume a much faster diffusion of IP 3 (≈ 300 µm 2 /s)  Ullah et al. (2006).
than Ca 2+ (≈ 10 µm 2 /s) due to Ca 2+ buffering and often disregard Ca 2+ diffusion altogether. However, the effective rate of IP 3 diffusion can occur on a much slower timescale (Dickinson et al., 2016), equalizing the signaling range and speed of the two signaling factors. Moreover, gap junctions formed by connexin43, characteristic for astrocytes, are permeable to Ca 2+ (De Bock et al., 2012). We thus account for intra-and intercellular diffusion of both IP 3 and Ca 2+ in the model with the same effective diffusion coefficients. This includes exchange at borders between neighboring astrocytes to imitate the function of gap junctions, which is supported by evidence that IP 3 can diffuse through the gap junctions along with Ca 2+ (Yule et al., 1996). Specifically, the diffusive term in Equation (1), e.g., for Ca 2+ reads: where D * Ca is the diffusion rate defined as the diffusion coefficient for Ca 2+ scaled by spatial grid step δx and local AVF value r: D * Ca = r 2 D Ca /δx 2 , and i enumerates the N nearest neighboring astrocyte-containing units. Here, we regard diffusion in porous media and assume that larger AVF is associated with larger cross-sectional area open for diffusion, as the astrocyte process diameter increases, and simultaneously with less tortuous paths taken up by diffusing molecules as the processes become less entangled. This leads to approximately quadratic scaling of D * with r. The diffusive term for IP 3 is defined in a similar way.
Finally, neighboring pixels with different AVF values obviously contain unequal volumes of astrocyte cytoplasm; hence, small concentration changes in areas with high AVF should cause larger diffusion-mediated concentration changes in the neighboring pixels with low AVF. This was accounted for by scaling the concentration rates of change due to diffusive exchange by the ratio of the AVF values of the two neighboring pixels.

Model Parameters and Numerical Details
The basic set of model parameters is given in Table 1. Only new parameters and values different from that in Ullah et al. (2006) are shown. For convenience, the full set of parameters is provided in Supplementary Table 1. The few values that are different were adjusted in order to provide the reasonable dynamics with the introduced treatment of [Ca 2+ ] ER as a variable in our model and the spatially extended layout. Diffusion coefficient for Ca 2+ is taken as a lower-bound estimate in Allbritton et al. (1992). Slow diffusion coefficient of IP 3 is based on Dickinson et al. (2016). We also added new parameters, specifically, A, τ Glu , and D Glu . The latter was chosen as a small value to describe only minimal spillover from a release site and buffering by binding to transporters. The pair of parameters describing instantaneous glutamate release rate and slower decay could be varied, because it is hard to assess the actual transmitter concentration and decay time as sensed by astrocyte leaflets. Extracellular glutamate transients occurring due do quantal synaptic release as estimated by fluorescent glutamate sensor have decay timescale in close to 100 ms (Jensen et al., 2019), and this value was used for the simulations shown below. This led to local glutamate transients peaking at 1.2 µM and decaying within 200 ms. We note that qualitatively similar Ca 2+ signaling dynamics could be obtained with a shorter τ Glu value, compensated by higher release rate A.
Numerical integration of the model differential equations is done in an explicit scheme (4th order Runge-Kutta method adopted for stochastic differential equations with a fixed timestep dt = 0.002 s) implemented in AGEOM-CUDA software (Postnov et al., 2012). Spatial grid step was δx = 0.275 µm/pixel for single-cell templates and δx = 0.55 µm/pixel for network templates. For reproducibilty, a reference implementation of spatial template generation and model simulation is available at https://zenodo.org/record/ 4552726#.YDAz1nUzZQ8 in form of Jupyter notebooks, Python and C code.
To quantify spatiotemporal properties of the simulated Ca 2+ dynamics, we examined complementary cumulative distribution functions (CCDF) of areas and lifetimes of individual Ca 2+ transients in all single-cell spatial templates to avoid selection bias. Ca 2+ transients were thresholded at 25% change from the local baseline level. The resulting contiguous TXY volumes of suprathreshold Ca 2+ concentration were treated as discrete events. CCDFF X (x) of some random variable X is defined as the probability P that X is greater than some x value:F X (x) = P[X > x]. We present the CCDF curves in double logarithmic coordinates to test if they can be approximated by a straight line, implying a power-law behavior. If a given variable is distributed according to a power law with probability density function (PDF) P X (x) ∝ x −α , then the CCDF also has a power-law behavior, but with a smaller exponentF X (x) ∝ x −(α−1) .

RESULTS
The proposed model, including the modifications to the local calcium dynamics and spatial mapping, was tested in a number of simulation experiments with different parameter settings and different spatial templates (which we call "cells" for shorthand below). To test for agreement between model behavior and the experimentally observed dynamics, first, we looked at the effect of the level of mean neuronal firing rate (local rate of the Poisson point process in terms of our model) on spatio-temporal dynamics of astrocytic calcium, and second, we tested whether the artificial spatial templates could provide for realistic intercellular calcium waves or other collective variants of astrocytic calcium dynamics.  Figure 2. At low excitation (p syn = 0.005 Hz) most Ca 2+ events were spatially confined and tended to start at a small number of sites, as shown with max-span contours and red labeling in Figure 4A (left) for 25 largest events. At higher excitation (p syn = 0.01 Hz) many Ca 2+ events spread to occupy the whole cell domain and again tended to initiate at the same sites. Synaptic signaling events were integrated into a spatial glutamate profile as shown in Figure 4A (bottom): local surges of extracellular glutamate are sparse at low excitation, while their instantaneous spatial density increases at high excitation, with a tendency of nearby sparks to blend. The tendency, exemplified by a single template in Figure 4A, was supported by the majority of single-cell templates ( Figure 4B): the number of events (during 2,500 s simulation time) covering more than 25% of the cell area increased with excitation for nearly all cells except six, which were incapable of generating whole-cell transients at p syn = 0.01 Hz. There were no obvious differences between these cells and all the others in overall morphology or AVF statistics. All cells did generate whole-cell transients at a higher excitation of p syn = 0.02 Hz. Most cells demonstrated a decline in the number of small events, covering less than 25% of the cell area with excitation, as a larger proportion of events was enabled to spread over larger areas, while the rate of event initiation could remain stable. The area of large (> 25%) events increased for all cells which generate large events under low drive conditions except the three, which did not generate large events at all, apart from other cells. The average area of the small (< 25%) events increased with excitation strength for all cells.

Wave Patterns in Single-Cell Templates
Stochastic local glutamate surges initiate two parallel processes: fast localized Ca 2+ transients and slower IP 3 production. Both integrate over time to steady-state levels of the model variables. Because the relative input of plasma membrane transport is defined by AVF in our model, we expected that the steady state levels and the probability of Ca 2+ event initiation should depend on AVF as well. Steady-state values of [Ca 2+ ] and [IP 3 ] decreased with growing AVF (Figures 4C,D), forming an uneven spatial profile. Calcium, as well as IP 3 levels, were higher in the periphery and lower near the soma, which is due to Ca 2+ entry during the synaptic excitation and due to higher IP 3 production in the regions with lower AVF. Different cells varied in steady-state levels of the variables, while intensified stimulation lead on average to a slight elevation of steadystate [Ca 2+ ] i , due to increased Ca 2+ entry and did not affect [IP 3 ] i , as the increase in [Ca 2+ ] i was insufficient for activation of PLC δ .
An example of a single large Ca 2+ -event is shown in Figure 4E. The expanding wave of elevated [Ca 2+ ] i initiates in a small location and spreads over the whole spatial domain (top row). Due to the large range of steady-state [IP 3 ] i concentrations, the event is unclear in absolute [IP 3 ] i values (middle row), but is obvious in the relative scale (bottom row).
Despite the stochastic nature of excitation, Ca 2+ activity in most cells is self-organized in a repeated pattern of Ca 2+ transient initiation and spreading ( Figure 5); event initiation sites were tightly clustered. Interestingly, activation in some clusters lead to spatially confined events, unlinked to activity in the rest of the cell, while transients originating in other sites tended to spread over the whole spatial domain in a repeated fashion. This is illustrated in Figures 5A-D for an example spatial template (see also Supplementary Video 2). This cell is markedly anisotropic, which defines the dominant wave spreading properties. The two active initiation sites labeled as #1 and #2 display different properties: events, starting in the site #1 often spread over large portions of the cell, as shown for a line-scan path in Figure 5B, while line-scan along the path starting in #2 was either activated by a wave coming from #1 or-very locally and with a higher frequency-by confined transients initiating in #2. Averaging small temporal windows around Ca 2+ spikes at the origin of path #1 shows that activation along this path is time-locked to activation of the initiation site. On the other hand, averaging similar temporal windows triggered by Ca 2+ spike at the origin of path #2 did not reveal any structured activation patterns. Figure 5D shows max-span contours of 25 largest events with their initiation sites mapped in red, as well as 5 peak-delay maps for a repeated pattern of activation. In these maps color indicates delay in seconds between the Ca 2+ peak at the initiation site and the Ca 2+ peak at each point of the cell. The precise wave initiation site could vary within 3-5µm, but the overall spatiotemporal pattern remained similar, with activation spreading mainly along the long thick processes toward the soma. Additional examples  of repeated patterns for other cells are provided in Figure 5E. In some cells preferred wave initiation sites could alternate between two polar positions or a few neighboring regions. There was also some scatter in the maximum delay between the initiation of the wave and its full expansion.
Though localized Ca 2+ events could be initiated in the low-AVF regions due to direct influx through the plasma membrane, these event seeds needed to reach a tip of a thicker branch with higher AVF to be amplified by IP 3 -mediated ER exchange. Thus, initialization of a global Ca 2+ wave critically depends on a coincidence of exactly the right spot in the AVF profileallowing both for a high enough [IP 3 ] i baseline and sufficient ER exchange-and a wide and long enough cluster of glutamate release due to local increase in synaptic activity. Areas containing FIGURE 6 | Complementary cumulative distribution functions for areas (left) and durations (right) of Ca 2+ events in all single-cell templates. Red lines-low excitation (p syn = 0.005 Hz), blue lines-high excitation (p syn = 0.01 Hz). Event area CCDFs: slopes of the fits for the low and high excitations are 3.0 and 2.3, correspondingly; the bend at the large areas corresponds to transition to whole-cell activation. Event duration CCDFs: slopes of the fits are 4.1 and 3.5 for the low and high excitations, correspondingly.
only thin processes with low AVF will display only frequent local Ca 2+ sparks, unable to invade the neighboring regions and thus will primarily set the baseline levels of [Ca 2+ ] i and [IP 3 ] due to diffusion. Indeed, astrocytes in hippocampal slices display frequent localized Ca 2+ events in the cell periphery, often termed "microdomains, " with a characteristic size of high-Ca 2+ spots much smaller than the cell size and originating in the thin processes region (Rungta et al., 2016). At the same time, the "lifespan" of the Ca 2+ wavefront increases with the already invaded area of the thick branch region due to regenerative Ca 2+ release, which in turn relies on background [IP 3 ] level. Thus, the specific topology of the cell template predicts the "hot spots" for the probability of Ca 2+ wave seeding.
We next describe statistics of areas and lifetimes of individual Ca 2+ transients (thresholded at 25% deviation from baseline). The left column of Figure 6 provides CCDF of areas, covered by individual events, while the right column describes durations of calcium events. Red curves correspond to the case of low neuron activity, the blue ones correspond to high neuron activity. In both cases the CCDF curves are presented in double logarithmic coordinates and can be approximated by a straight line within some ranges of event areas or durations, suggesting a power-law behavior in agreement with experimental data (Wu et al., 2014). After recalculating from CCDF to PDF exponents, the resulting parameters were different from that reported in Wu et al. (2014): for areas, α ≈ 3.3 . . . 4.0 in the model vs. α ≈ 2.1 . . . 2.4 in cultured astrocytes, and for durations α ≈ 4.5 . . . 5.1 in the model vs. α ≈ 1.97 . . . 2.16. These discrepancies can be explained by imperfection of the model and the 2D spatial embedding in the model. Defining model parameters that govern the shape of the event size and duration distributions is a potentially interesting outlook for further studies. Increase in synaptic excitation favored larger areas occupied by waves and longer durations of each event. The kinks in the CCDFs for event areas correspond to transition to whole-cell waves, i.e., events covering more than 30-50 µm 2 were likely to expand further and cover the whole cell.
We so far examined Ca 2+ transients in 2D spatial templates, which we use to generate astrocytic network in the next section to reduce computational load. We note however that the presented modeling approach is extendable to 3D with minimal modifications. An example of simulation for a 3D single-cell spatial template is presented in Supplementary Video 3. The features reported above for 2D templates remain in 3D: there are Ca 2+ transients at different spatial and temporal scales, waves that expand toward soma and vice versa, there are visible repeated patterns of activation.

Collective Dynamics of Astrocytes
After testing model behavior at microdomain and single-cell level, we turned to ensembles of connected cells. Collective dynamics of astrocytes was simulated using spatial templates containing about 40 cells as shown in Figure 7A; see also Supplementary Video 4. The simulated network is smaller than the size of cliques or networks of connected astrocytes in the neocortex (Houades et al., 2008), but still at the same order of magnitude. Simulations of larger spatial network templates are also possible but are more computationally demanding. Ca 2+ activity was not uniform across the spatial template: there were active "hot" brims in the periphery of most of the cell domains and some cells were less active than others, as shown in Figure 7B. We observed Ca 2+ transients originating in some astrocytes and expanding to their neighbors in a wavelike manner, an example of such a wave is shown in Figure 7C as a sequence of contours of elevated [Ca 2+ ] i , separated by 1 s.
To simplify the description of emerging spatiotemporal patterns, each cell can be considered as an element that "fires, " i.e., producing a global whole-cell spike, or remains silent. This cell-wise activity is shown for [Ca 2+ ] i and [IP 3 ] i in Figure 7D, where line colors correspond to domain colors in Figure 7A. We observed individual cell spikes as well as packs of tightly grouped spikes corresponding to multi-cellular waves. A considerable scatter is clear in the baseline levels of IP 3 , which reflects individual properties of each cell. IP 3 concentration peaks are wider and occur with a small delay in comparison to Ca 2+ peaks, which reflects the slower kinetics of IP 3 production timescale and a lag due to diffusion of IP 3 from the periphery to the somatic region.
A large Ca 2+ transient in one astrocyte can spread to other cells. In single cell templates we observed self-organized repeated patterns of activation. We were curious, if there will also emerge repeated patterns of intercellular dynamics at a network level. As a simple test for repeatable patterns, we calculated averages of the domain-averaged Ca 2+ rasters in a short time window, triggered by Ca 2+ spikes in different domains. We collected Ca 2+ traces, where each spatial domain served as a large region of interest (ROI), selected one cell as a seed, created time windows around Ca 2+ spikes in this cell and averaged Ca 2+ dynamics snapshots in all other ROIs across the time windows. In the case of stable activation sequence, the spikes appear with the same time-lag relative to the seeding astrocyte, and are thus visible in the raster plot, while if spiking in other astrocytes is not time-locked to the seeding astrocyte, the average Ca 2+ signal will be faint. Examples of such spike-triggered averages are shown in Figure 7E for three cells labeled as "#, " "*, " and "&" in Figure 7A used to center the time windows. Here, Ca 2+ spikes in cell "#" was time-locked to activation of a single other cell, activation of cell "*" lead to repeated activation of several other cells with a stable delay, while activation of cell "&" did not repeatedly lead to activation of other cells.
Inspired by this cell-wide "repeatability" measure based on the contrast of spike-triggered averages, we defined a similar score for more a detailed mapping of whether activation in some region repeatedly lead to activation in other areas with stable time lags. To this end, we split the spatial domain into overlapping square windows of size 5×5 µm and extracted Ca 2+ dynamics from these patches. We then selected the patches where there were more than 10 Ca 2+ spikes reaching at least 0.5 µM [Ca 2+ ] i and created spike-triggered 50 s-long averages from the raster of Ca 2+ signals in all patches. Percentage of all points with [Ca 2+ ] i > 0.2 µM in such spike-triggered windows was used as the repeatability score for a given patch. The patches were then projected back onto the spatial template, with averaging of the values in overlapping areas between patches. This resulted in an automated mapping of areas leading to repeated downstream activation at the network level, revealing cells, serving as hubs in spreading multicellular events ( Figure 7F). As suggested in Brazhe et al. (2018) for a simpler model, some spatial configurations of thick branches and leaflets can trigger persistent pacemaker-like activity, taking over the control of dynamics at the network level. We thus tested if alterations to the spatial template could lead to self-organization in a different spatial pattern of the centers of high repeatability (Figure 8). The cell labeled as "*" in Figure 7 was often activated from its direct neighbor to the right. Unlinking this cell from the network by setting to zero all contacts at domain boundaries ( Figure 7A) lead to a change in the repeatability map, decreasing the score for the cell "*" and increasing it for the three cells in the center. Unlinking the cell "*" from the network effectively silenced it, leaving only three cells with relatively high repeatability scores in the template.
A subtler alteration of the spatial template can be directed at the sites of frequent Ca 2+ transient initiation sites shown in red in Figure 8. Substituting the natural AVF profile at these sites with the average AVF value (Figure 8C) or ablating cell content from these areas (Figure 8D) lead to a dramatic reorganization of the Ca 2+ initiation sites and the centers of high repeatability. In both cases the number of initiation sites was reduced, with some of the former sites being silenced and some new sites formed in the neighborhood of the previous sites. Also in both cases the inequality of repeatability score was increased, with a few cells showing very high values. We attribute this to the reduced number of initiation sites, leading to a more repeated activation sequences.

DISCUSSION
We proposed a spatially detailed model of astrocytic calcium activity, which reflects current understanding of the two distinct mechanisms of Ca 2+ dynamics: excitable IP 3 -mediated exchange with ER in astrocyte soma and branches and plasma membrane exchange in the fine astrocytic processes and leaflets, sensitive to external conditions. Specifically, we suggest (i) an algorithm for data-based generation of 2D spatial templates matching realistic astrocyte morphology, and (ii) morphology-dependent spatially non-uniform parameter landscape for the calcium dynamics. To this end, we introduce the AVF parameter, which sets locally the relative input of the plasma membrane and ER based pathways and scales effective intracellular diffusion coefficients. The central idea underlying this separation is that astrocytes "sense" synaptic activity with fine processes, and it is where Ca 2+ transients are relying on extracellular Ca 2+ rather than intracellular stores, and where the bulk of IP 3 can be produced, while thicker branches and somata provide the positive feedback gain mechanism for IP 3 -mediated Ca 2+ -induced Ca 2+ release from ER. This mechanism separation is directly mapped to cell morphology in our approach.
We tested the suggested framework both at individual cell level and for algorithmically created multicellular astrocytic network templates. Our results show that the model is able to reproduce characteristic spatiotemporal patterns of Ca 2+ dynamics driven by synaptic activity, represented by spatially uncorrelated point-sources of glutamate release coupled to focal Ca 2+ entry, triggered by independent stochastic Poisson spike trains. At the single-cell level, the statistics of Ca 2+ event durations and expansion areas turned out to have a powerlaw distribution resembling experimental data (Wu et al., 2014). Power-law statistics of the Ca 2+ transients does not directly FIGURE 8 | Modification of spatial templates changes activation patterns. Color-coding: repeatability score, red: wave initiation sites. (A,B) Unlinking single cells from the neighbors; Blue contours-isolated cells, red regions denote wave initiation sites as in Figure 4A. Isolation of the cell often activated the cell denoted "*" in follow from the model equations and is an emergent property of the interplay between astrocyte morphology and Ca 2+ dynamics.
The presented model is a rather simplified representation of native astrocyte morphology and Ca 2+ dynamics. It is less detailed, but also less computationally demanding than the framework proposed by Savtchenko et al. (2018), allowing for simulations on a GPU-equipped laptop rather than on a supercomputer or a cloud. A major simplification of our model, dictating its limitations, is the reduction of real 3D astrocytic morphology to flat 2D spatial templates. The flattening was primarily done for the sake of computational tractability, but also conceptually matches single-plane imaging regime. The emergence of repeated activation patterns should not depend on the embedding space dimensionality, although the repeated propagation patterns can become more complex and elaborate in 3D; one can also expect different expansion rates and initiation probabilities for Ca 2+ events in 3D. Notwithstanding, we argue that using 2D patterns can be a useful approximation. First, some astrocyte network systems can be regarded as effectively twodimensional, e.g., astrocyte cultures or retinal astrocyte networks. Second, "true" astrocyte morphology can be regarded as less than three-dimensional, with the degrees of freedom limited by branching connectivity of the astrocytic processes, creating more or less independent astrocyte lobes and sub-domains.
A 2D morphology-based modeling approach has earlier been employed by Kang and Othmer (2009) to study calcium waves in retinal astrocytes. Our interpretation of astrocyte morphology is however different from the one in Kang and Othmer (2009) in several key aspects. First, Kang and Othmer used a GFAP-positive immunostained micrograph of retinal astrocytes, but GFAP stains only a small fraction of astrocyte cytosol volume, while the rest, especially the mesh of thin processes, got excluded from consideration of calcium dynamics, leading for the reconstructed morphology to resemble that of cultured astrocytes (see e.g., Wu et al., 2014 for comparison of morphology). On the other hand, an account for spatial segregation of calcium activity in fine astrocytic processes and thick branches was key to this study. Second, we account for a morphology-based balance between Ca 2+ entry via plasma membrane and from ER, which is absent in the work of Kang and Othmer. Third, we use a larger spatial template and describe Ca 2+ dynamics on a larger spatial scale. Our model is also more focused on intercellular communication through gap junctions, as we do not account for extracellular purinergic signaling to explain Ca 2+ wave propagation, which was done in earlier works focused on retinal astrocytes (Kang and Othmer, 2009;Edwards and Gibson, 2010), keeping in mind that the role of the extracellular pathway can be more pronounced in the retina than in neocortex.
Another limitation comes from our approximation of the 3D mesh of fine astropil sponge by a continuous active medium, parameterized by astrocytic cytosol volume fraction. Effectively, we "glue" together the individual branches and leaflets, assuming that they are at least partly interconnected by autologous gap junctions and branch-to-branch loops. The idea of "loopy" or sponge-like organization of the astropil has experimental support (Wolff et al., 1998;Genoud et al., 2015;Arizono et al., 2020), and we therefore adopt the AVF framework to represent unresolved astrocyte processes, also accounting for the tortuosity of the sponge by AVF-dependent scaling diffusion coefficients.
Indeed, the mentioned simplifications are expected to limit the predictive power of the model with regard to event frequencies, scaling characteristics and propagation speed. On the other hand, the simulated patterns of single-cell calcium transients are qualitatively similar to that observed in a single focal plane, which suggests that this reduction seems to preserve main features of astrocyte dynamics, while it is worth to be investigated further at sub-cellular spatial scales in future work. A proof-of-concept possibility to use the proposed modeling framework to simulate Ca 2+ signaling in 3D is presented in Supplementary Video 3.
A notable feature of simulated Ca 2+ dynamics in this system is spontaneously emerging stable patterns in both initiation and propagation of calcium transients, in tune to the co-active neuronal and astrocytic cells or repeating sequences of neuronal activation reported in slices (Sasaki et al., , 2014. In agreement with Brazhe et al. (2018) we observed morphologydependent emergence of hotspots with persistent pacemaker-like activity, taking over control of the dynamics at larger scales. In single-cell templates these preferred initiation sites could lead to activation of either spatially confined microdomains or larger expanding Ca 2+ areas, covering up to the whole cell. In multicellular systems we observed self-organized patterns of repeated calcium activity involving multiple cells. All cells in the template sharing the same equations and parameters, local differences in morphology of single astrocytes and geometry of astrocyte-to-astrocyte contacts favored initiation of multicellular Ca 2+ waves in some cells, followed by repeated sequences of cell activation as the Ca 2+ wave sweeped across the network. Excluding some cells from the original template caused dramatic reshaping of repeating activation patterns; removing or mashing up cell content in event initiation hotspots effectively reduced the number of active initiation sites and led to more stereotyped network activity. Our simulations showed that different cells in the network templates had different levels of activity and could develop patterns of Ca 2+ events with preferred directionality. Recently, Wang et al. (2019) have demonstrated heterogeneity of individual neocortical astrocytes with respect to the properties of their Ca 2+ activity. They also report a tendency for anatomical directionality during network-wide bursts of Ca 2+ activity, accompanying locomotion and presumed to result from adrenergic input from locus coeruleus. In the light of the present study, it is interesting to ask, whether the experimentally observed heterogeneity of astrocytic signaling is defined by different patterning and levels of local synaptic activity or by individual astrocyte morphology, or the interplay of both? At a larger spatial scale it is interesting whether the observed population-wide directionality is a product of the afferentation delays from locus coeruleus or some previously unstudied anatomic directionality of astrocyte morphology? Going even further, it seems an exciting possibility that-if Ca 2+ activity indeed affects LTP and Ca 2+ activity patterns are shaped by astrocytic morphology template-some part of memory in the neural network can be engraved in the fingerprint of astrocyte morphology. We conclude that with the same parameter set, the specific dynamical regime and role of an individual cell within astrocytic network is to a large extent defined by its morphology and this may have implications for computational performance of the underlying neuronal network.
The concept of spatially segregated oscillator regarded in (Brazhe et al., 2018) attracts growing interest as a new class of bio-inspired dynamic models. For example, Vanag (2020) explores a model of spatially segregated Belousov-Zhabotinsky oscillating reaction, where the catalyst (analogous to calcium induced calcium release in Ca 2+ models) is confined to small beads, whereas the rest chemical components of the reaction can diffuse freely, which leads to complex dynamical regimes and interaction between neighboring beads with immobilized catalyst. Another similar concept is that of volume transmission and intercellular communication via diffusing signaling molecules secreted by excitable cells (Sykova and Nicholson, 2008). There also appears to be an interesting connection between the self-organizing Ca 2+ signaling in coupled astrocytes and the concept of "reaction-diffusion computing" or "chemical computing" based on oscillatory chemical (e.g., Belousov-Zhabotinsky) systems in coupled reactor volumes, e.g., microemulsions (Adamatzky et al., 2005;Epstein et al., 2012;Showalter and Epstein, 2015;Torbensen et al., 2017;Vanag, 2019;Proskurkin et al., 2020). It is conceivable that astrocytes can provide a substrate for such Ca 2+ -based reaction-diffusion computing in parallel to, and organizing sparse network-based neuronal connectivity. Interestingly, it is common for self-organizing spatially distributed chemical systems to rely on inhibitory diffusive coupling (Li et al., 2015), while only excitatory coupling of neighboring astrocytes was modeled. Finding a mechanism for negative diffusive coupling, competing with the excitatory one can be an interesting outlook stemming from this analogy. Finally, it is inviting to believe that a combination of "analog" chemical computing approaches and "digital" neuronal ones, as well as a combination of graph-based and neighbor-based connectivity will inspire new algorithms for machine learning.

DATA AVAILABILITY STATEMENT
The model code and data, used for network spatial template generation and simulations are available for download at https://zenodo.org/record/4552726#. YDAz1nUzZQ8 (doi: 10.5281/zenodo.4552726).

AUTHOR CONTRIBUTIONS
DV: computer simulations, analysis of simulation results, figure preparation. AV: computer simulations, analysis of simulation results, figure preparation, writing the manuscript. DP: program design and architecture for simulations, writing the manuscript. AB: conceptualization of the model, network generation algorithm, figure preparation, writing the manuscript. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
DP and AB acknowledge the support from the Russian Foundation for Basic Research, grant 19-515-55016, which covered the development of the modeling approach for spatial segmentation of inter-astrocyte functions. AV and AB acknowledge support from Russian Science Foundation, grant 17-74-20089, which covered development of the numeric algorithm for creation of realistic spatial templates, biophysical interpretation of the model, and the collective astrocyte signaling analysis. AB also acknowledges support from Interdisciplinary Scientific and Educational School of Moscow University Molecular Technologies of the Living Systems and Synthetic Biology, which covered 3D model simulations. DV acknowledges support from the Russian Foundation for Basic Research, grant 16-32-50221 for mobility of young scientists, which covered model interpretation in terms of non-linear dynamics, CUDA-based simulations of the model, and analysis of single-cell dynamics.