Bistability and Chaos Emergence in Spontaneous Dynamics of Astrocytic Calcium Concentration

: In this work, we consider a mathematical model describing spontaneous calcium signaling in astrocytes. Based on biologically relevant principles, this model simulates experimentally observed calcium oscillations and can predict the emergence of complicated dynamics. Using analytical and numerical analysis, various attracting sets were found and investigated. Employing bifurcation theory analysis, we examined steady state solutions, bistability, simple and complicated periodic limit cycles and also chaotic attractors. We found that astrocytes possess a variety of complex dynamical modes, including chaos and multistability, that can further provide different modulations of neuronal circuits, enhancing their plasticity and ﬂexibility.


Introduction
Recent experimental and theoretical findings demonstrated that glial cells, particularly astrocytes, can participate in neuronal signaling and brain information processing. Astrocytes operate at the time scale of seconds and can release different gliotransmitters, e.g., neuroactive chemicals, modulating neuronal activity from the level of single synapses up to coordination of different networks. Chemical calcium excitability represents the characteristic feature of astrocytic signaling [1][2][3][4][5][6][7][8]. Ca 2+ elevations in astrocytes are vital for the optimal functioning of the CNS [9]. Astrocytic Ca 2+ signals also control K + uptake [10], contribute to the regulation of local blood flow [11] and morphological plasticity of these cells [12] and induce the release of gliotransmitters [13]. It is believed that this release which is crucial for neuronal regulation is directly connected with the generation of intrinsic calcium pulses. That is, the calcium signals and their characteristics further define synaptic transmission efficiency and, hence, the dynamics of the accompanying neuronal networks.
Astrocytes display both triggered and spontaneous calcium signals. Spontaneous Ca 2+ events are generated intrinsically without any external stimuli, whereas triggered calcium events occur in response to changes in the astrocytic environment, such as synaptic or neuronal activity and physiologically relevant internal and external triggers. The mechanisms underlying the generation of spontaneous Ca 2+ elevations in astrocytes are still poorly understood [14]. In experiments, it was shown that astrocytes can generate spontaneous calcium signals during blocking the neuronal activity [15], blocking vesicular release from neurons and astrocytes [16,17] and, in the case of inositol 1,4,5-trisphosphate receptor type 2 (IP3R2), genetic deletion in astrocytes [18,19]. However, genetic knockout of IP3R2 leads to a considerable reduction of spontaneous calcium events supporting the importance of the endoplasmic reticulum (ER) and the IP3R. key bifurcation scenarios of transitions between different dynamical modes of calcium signaling including the appearance of multistability and chaotic calcium oscillations. The outline of the paper is as follows. In Section 2, we describe the considered mathematical model. In Section 3, we focus on the role of an input calcium flow and examine the properties of the regimes observed for the spontaneous calcium concentration. The coexistence of different regimes is discussed in Section 3.2. The role of the output calcium flow in the emergence of complicated regular and irregular mixed-mode oscillations of calcium concentration is studied in Section 4. In Section 5, compliance with experimental data and astrocytic chemical activity contribution to information processing are discussed. Finally, Section 6 summarizes the conclusion.

Description of the Mathematical Model
The Lavrentovich-Hemkin [31] model describes the dynamics of the intracellular calcium concentration and Ca 2+ -dependent concentration of inositol-1,4,5-triphosphate (IP 3 ), by the following three-dimensional dynamical system: where the expressions for J serca , J CICR and J PLC are: .
Dynamical variables include the intracellular calcium concentration [Ca 2+ ] cyt , the calcium concentration in the internal storage endoplasmic reticulum (ER), [Ca 2+ ] ER and the concentration of the secondary messenger inositol-1,4,5-triphosphate [IP 3 ] cyt , which controls the ionic channels and removes calcium into the cytosol. Most of the parameters of the system (1) were chosen according to the data of [31]: [31], the authors provided the computational simulation of conditions carried out in real experimental studies and illustrated the dynamics that were qualitatively similar to the experimental data. With this aim, the authors demonstrated the impact of the influx of the Ca 2+ ions from the extracellular matrix on both (i) the period change of the regular spontaneous calcium oscillations in astrocytes and (ii) the complicated dynamics appearance. In [30], the authors expanded this system of nonlinear differential equations by combining it with different types of voltage-gated calcium channels. At the same time, the mechanisms responsible for the change of the regimes still remained incomprehensible. To our knowledge, refs. [37][38][39] are the first works with the strict analysis of the Lavrentovich-Hemkin model. Independently, the authors of both groups were interested in the bifurcation mechanisms of the spontaneous calcium oscillatory dynamics emergence and studied the peculiarities of Andronov-Hopf bifurcations in the model. Furthermore, for both scientific groups, the chaotic dynamics in spontaneous calcium oscillations had become the subject of close attention. As a result, in [40], the authors studied ways to control such complicated behavior as well as in [41], the authors demonstrated that even small changes in the parameters of the system can significantly modify the bifurcation diagram, revealing the absence of complicated oscillations.
In this work, particular attention was paid to the influence of the parameters controlling transmembrane calcium transport. In particular, we focused on the role of extracellular calcium flow J in and the rate of the calcium escape k out on the change of calcium dynamics. Additionally, the following parameters were also varied: the rate of calcium flow through serca to the ER from the cytosol, v M2 , the rate of calcium flow through IP 3 from the ER to the cytosol, v M3 , and the amount of feedback between calcium in the cytosol and IP 3 and v p , Figure 1. In contrast to previous works, in this study, we focused on the emergence of bistable behavior in spontaneous calcium dynamics and on the ways for the emergence of various types of chaotic dynamics. In our calculations, the system of nonlinear differential equations (1) was numerically solved by the fourth order Runge-Kutta integration scheme with a time step of 0.005 s. To study the peculiarities of transition from the steady state solution to the oscillatory mode, different one-and two-parametric diagrams were calculated. In all calculations to obtain more accurate results; the transients were discarded. To study various chaotic regimes, we calculated the largest Lyapunov exponent by the use of the numerical method introduced by Wolf et al. [42].

Input Calcium Flow
Let us start from the parameter set considered in the original paper by Lavrentovich and Hemkin in [31]. For this case, weak flow of an extracellular calcium through the astrocytic membrane leads to a steady state regime at which the concentration of calcium in the ER can exceed the concentration of calcium in the cytosol by tens of times. As illustrated in Figure 2a, for J in = 0.02 µM/s, for instance, at the time t = 1000 s the concentration of [Ca 2+ ] cyt ≈ 0.0397 µM and [Ca 2+ ] ER ≈ 3.6175 µM. Meanwhile, IP 3 concentration in the cell approaches [IP 3 ] cyt ≈ 0.01 µM. The increase of J in leads to the appearance of oscillations in the system dynamics. Particularly, from the data obtained for J in = 0.05 µM/s, Figure 2b, it follows that decreasing the phase of [Ca 2+ ] ER during the oscillatory mode provokes a sharp increase in both [Ca 2+ ] cyt and [IP 3 ] cyt up to the values of 0.65 µM and 0.33 µM, respectively. Following the increase of J in leads to another steady state regime at which the concentrations of calcium in the ER and in the cytosol weakly differ. For J in = 0.07 µM/s, Figure 2c, at the time t = 1000 s these concentrations differ by four times only, 0.14 µM and 0.57 µM, respectively.

Steady States and Oscillatory Modes
To illustrate the peculiarities of the transition from steady state mode to oscillatory regime, Figure 3a-i show the two-parametric bifurcation diagrams. To present the whole dynamical picture, all three variables of the considered system were studied. The diagrams shown in the first column of Figure 3 were obtained for [Ca 2+ ] cyt , the results for [Ca 2+ ] ER are presented in the middle column, and the right column shows the diagrams obtained for [IP 3 ] cyt . The calculations were carried out for the parameters J in , v M2 , v M3 , and v p . Different shades of blue show the values of the system variables within the quiescent modes. The darker tone of blue corresponds to larger values of the corresponding concentrations. As seen from the diagrams, the concentration of calcium in the cytosol [Ca 2+ ] cyt is always increased with the growth of extracellular calcium flow through the astrocytic membrane. In Figure 3a,d,g, the darker tone of blue is observed for larger values of J in . The law for this growth can be obtained analytically by equating the right parts of the system (1) to zero (see Appendix A for details): where As follows from (2), the steady state concentration of calcium in the cytosol depends on only two parameters, J in and k out . Thus, the increase of J in leads to the linear growth of Ca 2+ * cyt as seen in Figure 3a,d,g. A similar change in the dynamics of [IP 3 ] cyt concentration can also be observed in Figure 3c,f, where the increase of J in provides the monotonous growth of the [IP 3 ] steady state concentration for any values of v M2 and v M3 , respectively. Indeed, it follows from (5), that the nonlinear law for this growth is: in ) 2 is positive for any J in > 0, equality (5) is the monotonously increasing function. Similar analysis of the Formula (5) allows us to explain the diagram shown in Figure 3i. Here, the increase of v p leads to the corresponding scaling only: larger values of v p give larger values of [IP 3 ] * cyt . From the diagrams shown in Figure 3b,e,h, it follows that, for [Ca 2+ ] ER , a more complicated behavior is observed. As expected, the larger values of v M2 (the rate of calcium flow through serca to the ER from the cytosol) lead to larger concentrations of [Ca 2+ ] ER . In contrast, for a larger value of v M3 (the rate of calcium flow through IP 3 from the ER to the cytosol), the smaller steady state value of [Ca 2+ ] ER is observed. Indeed, from equality (3), Note that for small input calcium flow J in , particularly, for J in < 0.005 µM/s, it seems that the level of the calcium concentration in the endoplasmic reticulum does not depend on J in , especially for the diagrams shown in Figure 3e,h. However, that is not the case. It can easily be shown that for small J in , k * is small enough, and equality (5), therefore, yields revealing that (i) [Ca 2+ ] ER (J in )-dependence increases nonmonotonously (its derivative increases with the increase of J in ) and (ii) both parameters, v M3 and v p do not change the value of calcium concentration in the endoplasmic reticulum. Thus, the analysis of the steady state domains shows that two ranges in J in can be contingently highlighted: One range is for J in before the oscillatory mode (for small J in ), and the other is after it (for large J in ). Within these ranges, the intracellular calcium [Ca 2+ ] cyt and IP 3 concentrations are different. Moreover, within these ranges, the ER calcium concentration [Ca 2+ ] ER even changes differently.
In the oscillatory regime, the difference between the minimal and the maximal value of the corresponding variable is shown by shades of red. Thus, the domains with a redto-yellow gradient present the parameter-dependent evolution of ] min cyt , respectively. In this case, the darker tone of red corresponds to larger values of the corresponding difference. From Figure 3a-c, it follows that the increase of v M2 leads to an increase of all the considered differences ∆ cyt , ∆ ER , and ∆ IP 3 . Moreover, an analysis of the oscillatory solutions calculated for various v M2 showed that the increase of this parameter also leads to the increase of the maximal values of all variables and to the significant growth of the oscillations period. For the considered range of the parameter v M3 , the difference ∆ cyt (together with the maximal value of the corresponding variable) monotonously grows with the increase of v M3 . Meanwhile, ∆ ER and ∆ IP 3 demonstrate nonmonotonous behavior with the increase of v M3 , Figure 3e,f. Analysis of the solutions calculated for various v p showed that the increase of this parameter leads to the increase of ∆ IP 3 (as well as the maximal values of [IP 3 ] cyt ), while ∆ cyt and ∆ ER are decreased. Moreover, with the increase of v p , the period of oscillations is also decreased.
Note that for small values of the considered parameters v M2 , v M3 , and v p , for any value of the external calcium concentration J in , only a steady state regime can be observed. Thus, only for the parameters exceeding certain threshold values v M2 > v th M2 , v M3 > v th M3 , and v p > v th p , the emergence of an oscillatory mode is possible. Particular attention should be paid to the analysis of the system behavior near the boundaries between the regions with steady states and oscillatory modes, i.e., between the domains shown in Figure 3 by the gradients of blue and red(yellow). It was shown that the nature of the equilibrium point stability change can be different near the different parts of this boundary. To further show the peculiarities, we consider Figure 3a in more detail.

Bistability Emergence: Coexistence of Steady State and Oscillatory Mode
Examining the diagram shown in Figure 4a Figure 4d,e, the coexisting solutions near the lower boundary (for J in = 0.0238 µM/s) are demonstrated. Figure 4d was obtained for initial conditions 3.96, 0.01). Thus, for spontaneous calcium dynamics, the multistability being the crucial in sudden switchings of the dynamical regimes [43][44][45][46][47], is also possible.
To show the change of amplitude in oscillatory mode, we calculated the differences To obtain the width of the ranges with bistable dynamics, we focused on ∆ cyt only. Near both boundaries of the oscillations-to-quiescence transition, the calculations were performed twice. Namely, for the lower boundary, the calculations of ∆ cyt with the increase of the parameter J in were carried out using the initials in a vicinity of the steady state (blue solid line in Figure 4f). These initials can be obtained automatically because for small J in , the steady state is globally stable (due to its uniqueness). Therefore, in numerical calculations, small changes in J in with taking the initials in the steady state observed for the previous value of J in , provide the initials within the small vicinity of the shifted steady state. The calculations of ∆ cyt with the decrease of the parameter J in were carried out by using the initials in a vicinity of the limit cycle (orange dashed line in Figure 4f). Here, for large enough values of J in (but still close to the boundary), the limit cycle is a unique attracting set in the phase space of the system. Thus, the initials near the limit cycle with the small change of J in can be obtained automatically.
Similarly, for the upper boundary, the calculations of ∆ cyt with the increase of the parameter J in were carried out using the initials in a vicinity of the limit cycle (blue solid line in Figure 4h). In contrast, the calculations of ∆ cyt with the decrease of the parameter J in were carried out using the initials in the vicinity of the steady state (orange dashed line in Figure 4h). Thus, the width of the ranges with the bistable type of behavior is ≈0.0001 µM/s and ≈0.002 µM/s for the lower and upper boundaries, respectively.
Obviously, the equilibrium point stability is lost via subcritical Andronov-Hopf bifurcations leading to the birth of unstable limit cycles near the red points depicted in Figure 4f,h. Note that this result is consistent with theoretical studies reported in [38]. For the lower boundary, the increase of the unstable cycle amplitude with the decrease of J in occurs, while for the upper boundary, the unstable cycle amplitude is increased with the increase of J in .
Approaching the stable limit cycles, both unstable cycles collide with the latter and disappear via fold limit cycle bifurcations at J in ≈ 0.02374 µM/s and J in ≈ 0.0615 µM/s, respectively. Taking into account the numerical results presented above, we can summarize the following proposition.

Proposition 1. For the Lavrentovich-Hemkin model with the parameters taken as in [31]:
(a) The change of extracellular calcium flow reveals the existence of two J in -values where the Andronov-Hopf bifurcation occurs; (b) Both the Andronov-Hopf bifurcations are subcritical, i.e., the change of the equilibrium stability occurs with an unstable limit cycle emergence; (c) These unstable limit cycles coalesce with the stable cycles and disappear via the fold limit cycle bifurcations with further change of extracellular calcium flow. The J in -parameter ranges between the Andronov-Hopf and the fold limit cycle bifurcations define the ranges of bistability, where the coexistence of an oscillatory and steady state modes is observed.
Note that for small values of v M2 , the transition from the steady state mode to the oscillatory regime occurs via the supercritical Andronov-Hopf bifurcation [48]. Here, the emergence of small amplitude oscillations is observed, and the boundary crossing corresponds to a transition from one monostable regime to another.

Output Calcium Flow
Similarly to the analysis given in Section 3, the peculiarities of transition from steady state mode to oscillatory regime were also examined for the parameters k out , and v M2 , v M3 and v p .  (Figure 5a,d,g), (Figure 5b,e,h), and (Figure 5c,f,i). Here, for ∆ cyt , the darker tone of color is observed for smaller values of k out because from equality (2), it follows that [Ca 2+ ] cyt ∼ 1/k out .  Figure 5c,f, the increase of k out provides the monotonous decrease of the [IP 3 ] steady state concentration for any values of v M2 and v M3 , respectively. Indeed, it follows from (5) that the nonlinear law for this decrease is [IP 3 ] * cyt (k out ) ∼ α/(β + k 2 out ), where α, β > 0. This means that for any values of parameters α and β, equality (5) defines the monotonously decreasing function because its derivative −2αk out /(β + k 2 out ) 2 is negative for any k out > 0. A similar analysis of Formula (5) allows us to explain the diagram shown in Figure 5i. Here, the increase of v p leads to the corresponding scaling only; larger values of v p give larger values of [IP 3 ] * cyt . From the diagrams shown in Figure 5b, In an oscillatory regime, as in Section 3, the difference between the maximal and the minimal value of the corresponding variable is shown by the gradient of red-to-yellow. The darker tone of the red in the two-parametric diagrams Figure 5a

Chaotic Spontaneous Calcium Dynamics
From the experimental data, it follows that a peculiarity of astricytic chemical activity is in the presence of peaks with different amplitudes of calcium concentration. In accordance with terminology given in [49], blips are short and weak peaks that correspond to the opening of one IP3R channel (or one tetramer in an IP3R channel), while the puffs are longer and higher peaks resulting from the coordinated opening of a group of neighboring IP3R channels (or their tetramers) through the calcium-induced calcium release principle (CICR). The emergence of complicated alternations of such peaks is possible within the framework of the Lavrentovich-Hemkin mathematical model [31]. Furthermore, we present three types of chaotic spontaneous elevation of astrocytic calcium.

Burst-Type Dynamics: Small-Peaks Irregularity
Note that in [31], the authors presented an example of complicated calcium dynamics that can emerge with the change of an extracellular calcium level. To obtain such complicated chemical activity, few parameters were changed in the model. Namely, the authors considered the case when k p = 0.164 µM (instead of k p = 0.3 µM) and k CaA = k CaI = 0.27 µM (instead of k CaA = k CaI = 0.15 µM). Physically, these changes of the parameters correspond to the situation where the PLCδ1 dynamics have a faster response to Ca 2+ in astrocytic cytoplasm, while the rate of calcium release via the IP3R (through the CICR) occurs at a higher concentration of cytosolic calcium, and this rate would not drop off as quickly. Using the same assumptions, we consider the parameters as in [31] and study the role of the output calcium flow in the emergence of chaotic chemical activity. Figure 6a,b show an example of irregular time series numerically obtained for the concentration of the cytoplasmic calcium [Ca 2+ ] cyt and 3D pictures of the corresponding attractor in the phase space of the system (1), respectively. As seen from the time series, being on this attractor, the phase point returns on a large amplitude loop after different time intervals defined by the time spent in the region of low-amplitude oscillations. Such activity looks like bursts with various durations of small-amplitude, during which the amplitude is chaotically varied. It is remarkable that such types of bursting or mixedmode oscillations are quite widespread in nonlinear neurodynamical systems [50][51][52][53][54][55][56][57][58]. Note that for the parameters taken in Figure 6a  To study the dynamical mechanism leading to emergence of the attractor shown in Figure 6a, the one-parametric bifurcation diagram was obtained. As a control parameter, the rate of the output calcium flow k out was considered. For each value of k out taken from the interval k out ∈ [0.42; 0.8] s −1 , the local maximums of [Ca 2+ ] cyt (t)-dependence were plotted on the diagram, Figure 6c. It should be noted that, simple limit cycle with local maximum located near ≈ 0.2 µM (for small k out ), disappears with the following decrease of k out (for k out < 0.439 s −1 the steady state can be obtained).
As seen from the enlarged part of the diagram, Figure 6d, with the increase of k out , the period-doubling cascade is observed. As a result of such a bifurcation scenario, a chaotic small-amplitude attractor emerges in the phase space of the system (not shown, but similar oscillations were observed, for instance, with the change of v p in [37] or v M2 in [41]). A further increase of k out leads to the appearance of a large-amplitude loop, and the attractor becomes similar to that shown in Figure 6a. On the bifurcation diagram, this corresponds to the appearance of an upper line for k out > 0.49 s −1 . As seen from the enlarged part shown in Figure 6d, the change of k out within the interval k out ∈ [0.4967; 0.4992] s −1 leads to the appearance of several alternating periodic windows and chaotic bands. The observed transitions are well characterized by the maximal Lyapunov exponent shown in Figure 6e.
With a further increase of k out , only periodic mixed-mode oscillations take place. The evolution of such oscillations can be described by the use of the α β notation with α and β being the integers indicating the numbers of large (puffs) and small (blips) maximum values of oscillations in one period, respectively. Thus, for k out > 0.5 s −1 , the following sequence can be written: 1 6 → 1 5 → 1 4 → 1 3 → 1 2 → 1 1 . Note that contrasting similar transitions can be obtained with the decrease of the input calcium flow [31].

Small-Amplitude Chaotic Chemical Activity
Considering the smaller value of the parameter k p , e.g., k p = 0.133 µM, here we examine the peculiarities of astrocytic chemical activity when the PLCδ1 dynamics have an even faster response to Ca 2+ in astrocytic cytoplasm than it was assumed in the previous case. For k out = 0.62 s −1 , Figure 7a  changed values for all its maximums. For the enlarged part of the diagram, three relatively wide periodic windows alternating with chaotic bands were obtained. In the widest window k out ∈ [0.629; 0.62967] s −1 , a sequence of reverse period-doubling bifurcations yields a period-7 attractor which is later destroyed at k out ≈ 0.62967 s −1 in a crisis event. A period-5 attractor observed for the large k out in Figure 7d, is also destroyed in a crisis event that occurs for k out ≈ 0.631 s −1 . For k out > 0.635 s −1 , similar to the notation used before, the following sequence can be written: 1 4 → 1 3 → 1 2 .  Figure 8e confirm the observed transitions. For k out > 0.485 s −1 , similar to the notation used before, the following sequence can be written: 1 2 → 1 1 . Taking into account the numerical results presented above, we can summarize the following proposition.

Proposition 2.
Within the framework of the Lavrentovich-Hemkin mathematical model, various scenarios of chaos emergence can be realized, and therefore, various types of chaotic spontaneous calcium oscillatory dynamics in astrocytes can be simulated.

Discussion
Small variations in the Ca 2+ entry through plasma membrane and regimes of IP 3 production in astrocytes can display a surprisingly rich dynamical repertoire of spontaneous Ca 2+ signaling. Keeping in mind that the astrocytic calcium induces modulations of synaptic transmission and neuronal activity, such variety of different calcium dynamics, from simple periodic oscillations to complex chaotic bursts, opens possibilities to control neuronal signaling indirectly through the neuron-astrocyte interaction.
Our study suggested several biophysical mechanisms underlying the emergence of spontaneous Ca 2+ oscillations with various amplitudes, frequencies, shapes, and other crucial features. Two important predictions follow from the model analysis hitherto discussed. The first one concerns the spontaneous Ca 2+ signal emergence in the certain range of transmembrane Ca 2+ flux values. This result is consistent with experimental findings that discovered that most spontaneous Ca 2+ events start in an optimal range of thin distal processes [27]. It was shown that the mechanism underlying such subcellular distribution of the Ca 2+ events is that the level of Ca 2+ entry through plasma membrane in astrocytic branchlets depends on their surface-to-volume ratio. Surface-to-volume ratio is highest in the distal branchlets, where Ca 2+ entry into the cytosol therefore produces the largest Ca 2+ level elevations, which then can induce CICR by activating IP3Rs. Second is the observation that it is sufficient to slightly vary IP 3 production rates by PLCδ to induce dramatic changes in the subsequent Ca 2+ dynamics. This could result in the emergence of regular, self-sustained stable oscillations, bursting or chaotic. Note that chaotic oscillations of different types combining irregular pulse sequences with variable amplitudes further provide different levels of modulations of neuronal activity, both in time and in space in the context of neuron-astrocyte interaction. Consequently, they facilitate or depress particular signal transmission pathways in neuronal nets, and chaos can be viewed as a tool to enhance degrees of freedom in astrocytic modulation of neuronal signal transmission. We also note that our results are consistent with other recent theoretical studies reporting on chaotic astrocyte signaling [38,39]. In this context, the dependencies of Ca 2+ and IP 3 oscillatory dynamics on the regime of IP 3 production suggested different modes of stimuli encoding by astrocytes. Periodic Ca 2+ oscillations could represent a mechanism of frequency encoding. In turn, chaotic Ca 2+ signaling could perform more complex encoding, employing frequency, phase, and amplitude encoding characteristics. It follows from our analysis that chaotic oscillations are more likely to appear for large rates of IP 3 production by PLCδ. Experimental studies show that the proteins PLCδ tether to the plasma membrane and various intracellular structures and mainly associate with the cell periphery [59]. Thus, there is overlapping in subcellular distributions of PLCδ and values of surface-to-volume ratio parameter, which determine the optimal level of Ca 2+ entry through the plasma membrane for the emergence of the spontaneous Ca 2+ activity in astrocytes. This evidence supports the hypothesis of the key role of PLCδ in the generation of the spontaneous Ca 2+ signals mediated by the Ca 2+ flux through the plasma membrane. On the other hand, our results confirm that Ca 2+ dynamics in individual astrocytes may be significantly dependent on the astrocytic morphology and spatial distribution of the cellular and subcellular components.
Spontaneous Ca 2+ events can be modulated by different external stimuli corresponding, for example, to neuronal activity and changes in the cell environment [14]. Experimental studies reveal that Ca 2+ dynamics do not simply replicate synaptic activity but are actually much more complex [60]. This may indicate that the properties of spatiotemporal Ca 2+ dynamics both spontaneous and triggered by neuronal inputs are likely to be gov-erned by the interplay of intrinsic astrocytic cellular properties, characteristics of neuronal inputs, and environmental changes. Understanding the complex dynamic mechanisms of intracellular Ca 2+ activity has remained a major challenge and is required due to recently identified roles of astrocytic signaling in synaptic, neural network, and memory functions [13,61]. Astrocytic Ca 2+ activity contribution to information processing remains unclear. Computational models can help with this issue. It was shown that astrocytes can induce neuronal firing synchronicity and synaptic coordination [62][63][64][65], can enhance generation of information in neuronal ensembles [66][67][68][69][70], and can contribute to memory formation [71][72][73][74][75].
Abnormal astrocytic signaling can indicate pathological conditions and cause synaptic and network imbalances, leading to cognitive impairment [61]. Understanding the mechanisms underlying calcium activity in astrocytes and their role in physiological and pathological neuronal activity will open up a perspective of new therapeutic opportunities [76][77][78].

Conclusions
In this work, we present a variety of dynamical regimes available within the framework of the Lavrentovich-Hemkin model in a wide and meaningful region of its parameter space. We have shown that the system has a unique equilibrium point and, for this steady state, the parameter dependent laws of the variables changes were analytically obtained. We numerically determined the domains where the dynamics fall to the unique equilibrium point of the system, where the simple limit cycles are observed, and where the bistability can emerge. We have shown the variety of the observed bursting regimes of calcium oscillations and demonstrated the attractors corresponding to various types of chaotic chemical astrocytic activity. For complicated periodic bursting activity, the parameter regions with different oscillatory dynamical behaviors were classified using the α β notation with α and β being the integers indicating the numbers of large and small maximum values of chemical activity in each period, respectively. Such classification according to the periodicity of the observed mixed-mode oscillations, allows illustrating possible dynamical transitions. These transformations in chemical activity might be crucial elements used in astrocytic coding and, therefore, implied the mechanisms modifying the neuronal activity. Providing useful background information about the mechanisms of possible astrocytic activity-based coding, the results obtained in this study can be used in further simulation of neuron-astrocytic networks and might be helpful in understanding their complicated interplay.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.