∎

11institutetext: Thais Appelt Peres Bartiê 22institutetext: 22email:
<EMAIL_ADDRESS>33institutetext: Leonel Teixeira Pinto 44institutetext:
44email<EMAIL_ADDRESS>55institutetext: Department of Chemical
Engineering and Food Engineering, Federal University of Santa Catarina,
Florianópolis, Santa Catarina, Brazil

# Characterization of biophysical determinants of spatio-temporal calcium
dynamics in astrocytes

Thais Appelt Peres Bartiê  Leonel Teixeira Pinto

(Received: date / Accepted: date)

###### Abstract

Most of the functions performed by astrocytes in brain information processing
are related to calcium waves. Experimental studies involving calcium waves
present discrepant results, leading to gaps in the full understanding of the
functions of these cells. The use of mathematical models help to understand
the experimental results, identifying chemical mechanisms involved in calcium
waves and the limits of experimental methods. The model is diffusion-based and
uses receptors and channels as boundary conditions. The computer program
developed was prepared to allow the study of complex geometries, with several
astrocytes, each of them with several branches. The code structure allows easy
adaptation to various experimental situations in which the model can be
compared. The code was deposited in the ModelDB repository, and will be under
number 266795 after publication. A sensitivity analysis showed the relative
significance of the parameters and identifies the ideal range of values for
each one. We showed that several sets of values can lead to the same calcium
signaling dynamics. This encourages the questioning of parameters to model
calcium signaling in astrocytes that are commonly used in the literature, and
it suggests better experimental planning. In the final part of the work, the
effects produced by the endoplasmic reticulum when located close to the
extremities of the branches were evaluated. We conclude that when they are
located close to the region of the glutamatergic stimulus, they favor local
calcium dynamics. By contrast, when they are located at points away from the
stimulated region, they accelerate the global spread of signaling.

###### Keywords:

Glutamatergic stimulus Ion channels dynamics Mathematical modeling Calcium
signaling Astrocyte dynamics

##  1 Introduction

There is experimental evidence to suggest that neurons and astrocytes
communicate with one another. This suggests that astrocytes are involved in
brain information processing. Researchers have proposed several experimental
approaches in order to understand the role of astrocytes in brain
communication (Halassa et al., 2007a; Perea et al., 2009; Pirttimaki et al.,
2017). These approaches are primarily designed to analyze calcium signaling,
because this is the main activity of astrocytes (Bernardinelli et al., 2014;
Bindocci et al., 2017; Howarth, 2014; Navarrete et al., 2012; Sasaki et al.,
2014). The results suggested some information regarding the role of
astrocytes; however, they are quite discrepant. The discrepancies of the
experimental results hinder development of a global view of the contributions
of astrocytes in brain communication (Bazargani and Attwell, 2016; Lallouette
et al., 2018; Volterra et al., 2014). Recently, the very existence of
gliotransmission was discussed in two perspectives papers (Fiacco and
McCarthy, 2018; Savtchouk and Volterra, 2018), with opposite conclusions. The
former do not believe in its existence under physiological conditions, while
the latter consider it to be well established. These discrepancies arise from
difficulties in detecting and interpreting calcium signaling. Elaboration of
sophisticated experimental apparatus and critical analysis of results can
improve understanding of the role of astrocytes (Losi et al., 2017; Poskanzer
and Molofsky, 2018; Rusakov, 2015). Mathematical models can contribute to
experimental data interpretation, because they help to identify the mechanisms
involved in specific phenomena as well as the limits of experimental methods.
Furthermore, validated mathematical models allow the elaboration of new
hypotheses about each phenomenon involved and the prediction of a system’s
dynamic behavior in various situations (Min et al., 2012; Myung and Pitt,
2002; Suffczynski et al., 2006).

Both mathematical descriptions and interpretation of experimental results of
calcium signaling are not trivial. They face problems such as geometric
complexity of astrocytes as well the large number of mechanisms involved in
calcium signaling generation and propagation. Mathematical models also face
another challenge: how to mathematically describe cellular processes such as
the dynamics of the various channel types and their opening probabilities. The
main problems associated with interpretation of calcium signaling and its
mathematical description will be elaborated below.

The astrocytes are considered in this work. They are called ”stars” due to
their high level of branching (Matyash and Kettenmann, 2010). The branches are
important because they allow astrocytes to contact neurons and other
astrocytes. In addition, the branches allow the astrocytes to involve distant
synapses, building communication paths between them (Halassa et al., 2007b).
Lack of reflection on this structural complexity may result in antagonistic
interpretations of astrocyte phenomena (Bazargani and Attwell, 2016; Volterra
et al., 2014). In order to minimize this antagonism, it is important to
consider where the signaling is detected, because cell processes occur at
different amplitudes and periods depending on the region where they start
(Shtrahman et al., 2017; Wu et al., 2014). Another important question is the
way astrocytic connections are spatially distributed. The spatial distribution
can directly affect calcium signaling dynamic (Perez-Alvarez et al., 2014).

Calcium signaling is triggered by a combination of various mechanisms
(McCarthy and Salm, 1991). The endoplasmic reticulum (ER) is the main
organelle contributing to signaling through calcium release by channels
(Berridge et al., 2003; Bootman et al., 2002). In addition to IP3-activated
Ca2+ channels, the endoplasmic reticulum has ryanodine-sensitive receptors.
Both IP3 and ryanodine receptors are modulated by a mechanism called calcium-
induced calcium-release (CICR) as a consequence of the increase of
intracellular Ca2+ concentrations (Simpson et al., 1998; Skupin et al., 2008).
Other biological processes participate in calcium signaling in astrocytes,
including channel opening triggered by activation of glutamate ionotropic
receptors (Koh et al., 1995) and ATP (Neary et al., 1988); activation of
Na+/Ca2+exchanger (Golovina et al., 2003), and Ca2+-ATPase pump (PMCA)
(Verkhratsky et al., 2012), and opening of voltage-dependent calcium channels
(MacVicar, 1984). In addition, passive channels are present in cytoplasmic and
ER membranes (Clapham, 1995). Also, there are channels activated by lack of
Ca2+ in ER linked to transient receptor potential channels (TRPC) (Verkhratsky
and Parpura, 2014).

Based on several mechanisms described above, it is possible to conclude that
calcium signaling is a complex phenomenon. Thus, development of a mathematical
model that describes it requires selection of the most relevant mechanisms.
This is not an easy task; however, good selection is fundamental such that the
mathematical model may permit plausible inferences.

Several mathematical models of calcium signaling in astrocytes have been
proposed in recent years; however, as was the case with experimental results,
these models have not yet led to definitive conclusions regarding the role of
astrocytes in brain information processing. Most of these models were derived
from one of three classic models describing calcium signaling: Keizer and De
Young (1992), Li and Rinzel (1994) or Höfer et al. (2002). In addition to
using one of these models as a basis, recent models used the same values as
their parameters. Keizer and De Young (1992), for example, stated that the set
of parameters they used was not the only one to satisfy the established
conditions. In addition, parametric sensitivity analysis has been rarely
performed in papers proposing mathematical modeling of calcium signaling in
astrocytes (Manninen et al., 2018). Lack of parametric analysis may be
responsible for comparability problems observed by Manninen et al. (2017).
Another relevant question is why the three basis papers do not make direct
comparisons with experimental results. Recent papers rarely do this either. It
is possible that this phenomenon is related to the many gaps remain in
experimental results (Poskanzer and Molofsky, 2018; Durkee and Araque, 2019;
Oheim et al., 2018; Shigetomi et al., 2016; Rusakov, 2015). It is also
important to note that Manninen et al. (2018) and Manninen et al. (2017) point
to difficulties in reproducing literature models, either on account of errors
in model creation itself or by absence of information.

Given the discrepancies and inconsistencies of experimental results and
mathematical models, deep reflection and engagement in the creation of new
approaches are necessary. These approaches must consider peculiarities of
astrocytes and the various challenges they present.

Against this background, the objective of this work was to propose a
biologically coherent mathematical model for prediction of calcium signaling
dynamics triggered by glutamatergic stimuli in astrocytes. This model
considers the most relevant characteristics and processes regarding astrocytes
to be used as a tool for further studies on the roles of these cells. A
parametric sensitivity analysis was performed and the source code of
computational program will be available in ModelDB (McDougal et al., 2017) at
http://modeldb.yale.edu/266795, allowing for its reproduction. The code is
available on the second author’s ResearchGate. Program availability and the
manner in which it was built enables researchers worldwide to contribute in
many ways, both experimental and mathematical. As the model is being enriched,
it will be possible to make effective contributions to understanding of the
role of astrocytes in brain information processing.

##  2 Mathematical model

Calcium signaling depends on mechanisms that occur mainly in three regions:
synaptic cleft, cytosol, and endoplasmic reticulum. In the first region, only
the temporal dynamics of glutamate that reaches the astrocytes proximity is
modeled. In the cytosol, the dynamics of calcium and IP3 are modeled. The
region is considered heterogeneous and the model includes diffusive partial
differential equations. The endoplasmic reticulum acts as a functional calcium
storage site. Due to the proportionally small dimensions of the ER and its
high calcium concentrations, it is possible to disregard its concentration
profiles and to use a homogeneous approach with ordinary differential
equations. Homogeneity is used here, not in the mathematical sense, but in the
chemical sense, as a mixture that has a uniform concentration at all points.
If this consideration is possible, the benefit would be reduced computational
time and memory requirements. The two approaches, homogeneous and
heterogeneous, will be tested and compared.

Fig. 1 shows a schematic drawing of regions included in the mathematical
model, highlighting the mechanisms considered in each. As the objective of
this work was to study calcium signaling triggered by glutamatergic stimulus,
the main mechanisms considered were: activation of ionotropic and metabotropic
glutamate receptors, i-GluR and m-GluR (Verkhratsky et al., 2012); activity of
ATPases and Ca2+ passive channels in both cell and ER membranes (Clapham,
1995); opening of IP3-activated channels inthe ER membrane (Verkhratsky et
al., 2012); TRPC channel operation, and a contributor to mechanism called
store-operated calcium entry (SOCE) (Verkhratsky and Parpura, 2014).

Figure 1: Calcium signaling mechanisms considered in the mathematical model.
The diagram illustrates some mechanisms involved in calcium signaling. In the
cell membrane, there are metabotropic glutamate receptors (m-GluR), which are
activated when neurons release glutamate into the synaptic cleft during
synaptic transmission. Its activation triggers production of IP3, which moves
by diffusion. Upon reaching the endoplasmic reticulum membranes, IP3 activates
receptors that open calcium channels. Calcium concentration decreases in the
endoplasmic reticulum promote activation of STIM, which opens store-operated
calcium channels (SOC) located in the cell membrane. In both cell and
endoplasmic reticulum membranes, passive channels and pump operate to maintain
or restore resting state. In the cell membrane, glutamate released by neurons
also activates ionotropic receptors (i-GluR), opening calcium channels. All
mechanisms mentioned above were included in the mathematical model.

After choosing the relevant mechanisms, mass balance was performed in regions
where ordinary or partial differential equations are used (cytosol and ER).

The cytosol is a colloidal medium that contains organelles and other
structures that make the free movement of ions and molecules difficult (Luby-
Phelps, 1999). In this medium, diffusion is the main mass transport mechanism
(Höfer et al., 2002). Consequently, calcium and IP3 dynamics in cytosol ($C$
and $I$, respectively), can be described by Eq. 1, where $D_{W}$ is the
compound $W$ diffusion coefficient in cytosol (cyt):

$\frac{\partial W_{cyt}}{\partial t}=D_{W}.\nabla^{2}{W_{cyt}},W=C,I$ (1)

The same equation can be used to describe calcium dynamics in ER, when this
organelle is considered heterogeneous. In the equations throughout the text,
the variables are the concentrations of IP3 (I), calcium (C) and glutamate
(G). The main subscript used refer to the regions where the variables are
calculated, cytosol (cyt), endoplasmic reticulum (ER) and synaptic cleft (sc)
and to the membranes where flows are calculated, endoplasmic reticulum (ERM)
and cell body (CM). The description of the parameters is presented in Table 1.

As previously mentioned, the endoplasmic reticulum can be considered
homogeneous. Because of its narrow thickness, the diffusive flow of calcium
quickly smooths out the concentration gradients created by outflow. In this
case, the calcium fluxes occurring in its membrane must be directly computed
in mass balance. Assuming the variation of calcium concentration in the ER
depends on the fluxes produced by calcium pumps ($v^{pp}_{ERM}$) IP3-activated
calcium channels ($v^{lc}_{ERM}$) and calcium passive channels
($v^{pc}_{ERM}$) (Ullah et al., 2006), Eq. 2 can be written to describe the
calcium concentration profile in ER. It is necessary to calculate the integral
of calcium flux over the surface of ER membrane ($S_{ERM}$) to reconcile the
homogeneous ER approach with the heterogeneous cytosol approach.

$V_{ER}\frac{\partial C_{ER}}{\partial
t}=\sum_{g=1}^{3}\int\int_{S_{ERM}}v^{g}_{ERM}$ (2)

where g = 1(pp), 2(pc), 3(lc)

The consideration of the homogeneity of the mixture within the endoplasmic
reticulum is typical in the mathematical models found in the literature. The
review made by (Manninen et al., 2018) listed 106 mathematical models
generated between 1995 and 2018. Only three of these consider diffusion of
calcium in the ER. The question is not whether the flux of calcium through ER
membranes creates concentration profiles, which obviously happens, but whether
they need to be considered in the mathematical model. To help understand this
issue, both approaches were analyzed. The same program configuration was
performed for a duration of two seconds, with both considerations and their
results were compared. The differences were less than 0.0001%, considered
insignificant. The differences in the concentration curves were not
distinguishable, which is why they are not shown in this work. According to
Occam razor, it was determined to maintain the consideration of homogeneity.

The buffering effects of calcium can be considered in two ways: by considering
the free calcium concentration, or by considering the calcium effective
diffusion coefficient in the cytosol. The effective diffusion coefficient was
the form chosen for this work (Höfer et al., 2002).

Mathematical expressions for almost all of fluxes produced by the mechanisms
represented in Fig. 1 can be found in literature. However, they were developed
to be used in homogeneous models. Even in some already proposed heterogeneous
models, such as thosw of Höfer et al. (2002) and Wu et al. (2014), these
fluxes were incorporated in mass balance. This was done to simplify the
mathematical model and to avoid the need to incorporate these fluxes as
boundary conditions. In this paper, we consider the places where fluxs occur,
which implies their use as boundary conditions. They are represented by
equations 3 to 5 and 7 to 9. Therefore, the mathematical expressions proposed
in the literature have been modified to include the density of the channels
and, therefore, the fraction of the area occupied by the channels in the
membranes of the cell and the ER.

IP3 production depends not only on glutamate coupling with metabotropic
receptor, but also on the concentration of intracellular calcium, because Ca2+
keeps the effector protein activated (Pawelczyk and Matecki, 1997). Ullah et
al. (2006) shows an example for calculation of IP3 production summing the
effects of these two compounds. This appears inappropriate, given that the
protein amplifies, but does not initiate the activation (Thore et al., 2005).
Based on this fact, the present model multiplied the two contributions, as
shown in Eq. 3. In relation to the term corresponding to calcium contribution,
the proposal of Höfer et al. (2002) was used. In the equation the gradient
$(G_{sc}-G^{*}_{sc})$ is used because it is assumed that IP3 production is not
relevant at rest concentration. Thus, Eq. 3 represents the IP3 production in
an infinitesimal element of area, dA.

$u^{p}_{CM}=f^{p}_{CM}.dA_{CM}K^{p}\frac{(G_{sc}-G^{*}_{sc})^{n}}{(G_{sc}-G^{*}_{sc})^{n}+k_{p,G}^{n}}\frac{C_{cyt}^{2}}{C_{cyt}^{2}+k_{p,C}^{2}}$
(3)

IP3 degradation, Eq. 4, was described using first-order kinetics, as proposed
by Ullah et al. (2006). Elimination of IP3 occurs until its concentration
returns to the resting value.

$u^{d}_{ERM}=f^{d}_{ERM}.dA_{ERM}.h^{d}.(I_{cyt}-I^{*}_{cyt})$ (4)

Calcium flux through passive channels, both in the cellular or ER membranes
(Di Garbo et al., 2007), was considered proportional to the area fraction of
these channels, to mass transference coefficient, and to concentration
gradient between both sides of the membranes. This relation replaces the
concept of constant current of calcium entry in the cell, which is used in
several articles (Amiri et al., 2011; Höfer et al., 2002; Ullah et al., 2006).
Eq. 5 shows calcium flux through passive channels in the ER membrane.

$v^{pc}_{ERM}=f^{pc}_{ERM}.dA_{ERM}.h^{pc}_{ERM}.(C_{ER}-C_{cyt})$ (5)

We incorporated transient receptor potential channels (TRP channels) activated
from a reduction in ER calcium concentration. It is possible consider these
channels present identical dynamics to those of passive channels after
activation and opening. It can be assumed that there is an effective increase
in passive channel fraction in the cell membrane that depends on Ca2+
concentration in the ER. Mathematically, this idea is expressed in the Eq. 6.

$f^{TRPC}_{CM}=\left\\{\begin{array}[]{rc}\frac{1}{100}.\frac{C^{*}_{ER}-C_{ER}}{C^{*}_{ER}},&\mbox{for}\quad{C^{*}_{ER}-C_{ER}}>0,\\\
0,&\mbox{for}\quad{C^{*}_{ER}-C_{ER}}\leq 0.\end{array}\right.$ (6)

The reduction factor (100) employed in Eq. 6 causes the fraction of these
channels to be at the same order of magnitude as those of other channels.

Eq. 7 gives the adaptation of Eq. 5 to the cytoplasmic membrane, taking into
account the TRP channels.

$v^{pc}_{CM}=(f^{pc}_{CM}+f^{TRPC}_{CM}).dA_{CM}.h^{pc}_{CM}.(C_{sc}-C_{cyt})$
(7)

The equation for describing calcium flux produced by the pumps in the ER
membrane was developed from experimental data fit. The equation considers that
pump operation depends of intracellular calcium concentration, up to a
saturation limit, given by the $k_{pp}$ constant (Ullah et al., 2006). Calcium
influence is of second order, because the pump transports two calcium ions
simultaneously. According to Verkhratsky et al. (2012), these pumps are also
located in cell membrane. Its action was included in the present model. Eq. 8
describes the calcium flux due to pump action.

$v^{pp}_{m}=f^{pp}_{m}.dA_{m}.K^{pp}_{m}.\frac{C^{2}_{cyt}}{C^{2}_{cyt}+k^{2}_{pp}},m=CM,ERM$
(8)

Eq. 9, based on (Ullah et al., 2006), represents calcium flux through ligand-
activated channels present in ER and cell membranes.

$v^{lc}_{m}=p^{o/c}_{m}.f^{lc}_{m}.dA_{m}.h^{lc}_{m}(C_{r}-C_{cyt})\left\\{\begin{array}[]{rc}m=CM,r=sc\\\
m=ERM,r=ER\end{array}\right.$ (9)

The mathematical description of opening and closing of ionic channels, among
all the previously described phenomena, is the most challenging, because the
channels are not simultaneously activated (Falcke, 2003) and each one has its
own particularities (Eisenberg, 1999). There are interesting proposals in
several papers; nevertheless, problems remain, including the number of
parameters and the absence of clear criteria for the choices of their values
(Keizer and De Young, 1992; Means et al., 2006). Although they have been used
in later investigations (Höfer et al., 2002; Taheri et al., 2017; Ullah et
al., 2006; Wu et al., 2014), these values were obtained in studies of
different cell types; therefore, their generalization for study of astrocytes
is questionable.

Figure 2: equations application domain The diagram illustrates some mechanisms
involved in calcium signaling. A. Eq. 12 describes the diffusive movement of
Ca2+ and IP3 in the cytosol of all regions. It is known as Fick’s second law
of diffusion. Eq. 2 is used to calculate the variation in ER Ca2+
concentration. It says that any variation in this concentration is result of
Ca2+ fluxes that occur in the ER membrane, whether produced by calcium pumps,
passive channels or ligand-activated channels. B. IP3 is formed by m-GluR
receptors, as a result of glutamatergic stimulus (Eq. 3). It moves through the
cytosol by diffusion (Eq. 12) and is degraded in the ER membrane (Eq. 4). C.
Calcium flux through passive channels is calculated on ER membranes (Eq. 5)
and on plasmatic membranes at the ends of the branches (Eq. 7). The difference
is that the latter includes Transient Receiver Potential Channels, TRPC. The
lateral membranes of the branches are considered isolated. D. The fluxes
produced by calcium pumps (Eq. 8) and ligand-activated channels (Eq. 9) are
calculated for both membranes. The probability of opening the channels is
calculated by Eq. 10.

The discussion about opening and closing of ionic channels is related to the
dynamics of the parameter $p^{o/c}_{m}$ in Eq. 9, indicating whether the
channel allows, at a specific moment, the passage of calcium ($p^{o/c}_{m}=1$)
or not ($p^{o/c}_{m}=0$). Results obtained by Tu et al. (2005a) and Tu et al.
(2005b) were used for description of dynamics parameters. These researchers
conducted a detailed study of IP3 activated channels and they proposed
equations that describe the channel probability of opening according to
cytosolic concentrations of Ca2+ or IP3, respectively $P_{ca}$ and $P_{I}$.
The authors calculated the probability of opening either maintaining constant
IP3 concentration for the first, or Ca2+ concentration for the last.

They obtained two experimental curves, fitted with specific functions. In this
model, both concentrations vary simultaneously. The total probability of
opening channels must meet the condition $P_{I}\in[0;1]$. Applying the theorem
for independent probabilities, $P(A\cap B)=P(A)$ x $P(B)$, the two expressions
proposed by Tu et al. have been multiplied. The result is Eq. 10, which
represents the total probability of opening channels activated by IP3, for the
situation of simultaneous and independent variation of calcium and IP3
concentrations.

$P_{I}=4.\frac{I^{n_{I}}_{cyt}}{I^{n_{I}}_{cyt}+ks^{n_{I}}_{I}}.\frac{ks^{n_{C}}_{C}.C^{n_{C}}_{cyt}}{(ks^{n_{C}}_{C}+C^{n_{C}}_{cyt}).(Ks^{n_{C}}_{C}+C^{n_{C}}_{cyt})}$
(10)

For glutamate-activated calcium channels (i-GluR), a similar equation was
used. Because there is no information on the participation of Ca2+ in the
activation or inactivation process, a simpler equation was used, solely
dependent on glutamate. The Eq. 11 is similar to the first term of $P_{I}$;
however, the concentration of glutamate is now considered.

$P_{G}=\frac{G^{n_{G}}_{sc}}{G^{n_{G}}_{sc}+ks^{n_{G}}_{G}}$ (11)

Fig. 2 shows in detail the locations where each equation in the mathematical
model is used. It mainly seeks to detail the application of boundary
conditions.

##  3 Aspects related to computer implementation

The geometry to represent the astrocytes should basically contain several
branches and a body to which the branches can be coupled. Shapes, such as
prisms and spheres, could be used to represent cell body geometry,
facilitating the use of macroscopic equations to describe molecular motion.
Prisms or cylinders can represent the branches. The computational
implementation of the mathematical model considered that most of available
data comes from experiments with astrocytes cultures seeded in monolayers. The
use of a two-dimensional model (Höfer et al., 2002) is reasonable in these
conditions, what makes possible to consider the cell body and its branches as
rectangular prisms containing the same depth. The Cartesian coordinate system
is appropriate to the chosen geometry. Thus, considering this coordinate
system and two-dimensional variation, Eq. 1 can be rewritten as the Eq. 12.

$\frac{\partial W_{cyt}}{\partial t}=D_{W}.\frac{\partial^{2}W_{cyt}}{\partial
x^{2}}+D_{W}.\frac{\partial^{2}W_{cyt}}{\partial y^{2}},W=C,I$ (12)

The location and geometry of endoplasmic reticulum influence calcium signaling
(Berridge et al., 2003). According to Pivneva et al. (2008), the endoplasmic
reticulum locates in vicinity of the nucleus and is also dispersed in the
cytosol, near the plasma membrane. In astrocytes in situ, the investigators
have found most of reticulum in vicinity of plasma membrane, that is, in the
branches. Patrushev et al. (2013), By contrast, points out that the volume
occupied by ER increases from the end of branches towards the cell body. In
addition, Oschmann et al. (2018) demonstrated that endfeet processes are
triggered mainly by calcium input from the extracellular environment and not
by release of calcium stored in RE. Literature shows that there is no single
opinion on this subject.

Based on this information this work was carried out in two steps. The first
stage considered the study of effects produced by cytosol scattered ERs
requiring a great computational effort. Not only by increasing the number of
variables, but also by the number of possible alternatives, because the
positions of ERs become almost random.

Given this, were performed calcium signaling analyzes without the existence of
ER in the branches. Each cell has a single ER, positioned around the nucleus.
Nucleus and ER shapes and positions are shown in Fig. 3. In the second stage a
simplified structure was elaborated, in which astrocytes are composed by a
central region with two opposite arms. The ER positioned in the branches are
formed by rectangular prisms. This configuration will be detailed later.

Fig. 3 shows the representation of the computationally implemented geometry,
highlighting the dimensions and the distribution of the branches.

Figure 3: Cell geometry used in mathematical modeling. The diagram shows the
structure used. It consists of four branches connected to the cell body. The
cell body has a nucleus surrounded by the endoplasmic reticulum. The
dimensions, in $\mu m$, are, 8.0 x 8.0 x 1.0 for the cell body, 3.0 x 0.7 x
1.0 for branch elements, 0.7 x 0.7 x 1.0 for connecting elements, and 4.3 x
4.3 x 1.0 for the nucleus. The nucleus is surrounded by the ER, whose external
faces, with dimensions of 5.4 x 1.0, have one hundred of IP3-activated
channels. These channels are numbered clockwise from the center of the bottom
face. The ER is isolated on its internal faces.

It is important to note that the computer program developed is versatile
enough to build several cells, each with several complex branches, of
different sizes and in different positions. This versatility is best described
in a file attached to the software in ModelDB.

For the solution of Eq. 12, a mesh was constructed filling all the cytosol
area, with $\Delta x$ and $\Delta y$ dimensions. In the computer program,
$\Delta x=0.1\mu m$ and $\Delta y=0.05\mu m$ were chosen in all regions. The
fluxes produced by the channels illustrated in Fig. 1 only reach this mesh in
its contours. For this reason, they must be considered as boundary conditions.
These conditions will now be detailed

###  3.1 Boundary conditions

The boundary conditions describe the effect of calcium and IP3 fluxes on the
control volumes that belong to cell membranes. As shown in Fig. 1, these
membranes delimit the branches, the cell body and the ER. The branches lateral
faces are considered insulated. At its extremities, by contrast, there are
calcium fluxes through channels activated by glutamate, passive channels and
pump. In addition, there is production flux of IP3. The boundary condition in
these cases is based on the idea that mass flux crossing membrane moves by
diffusion within the cell. Mathematically, this condition is represented by
Eq. 13 (IP3) and Eq. 14 (Ca2+).

$v^{p}_{CM}=-dA_{CM}.D_{I}.\frac{\partial I_{cyt}}{\partial u}$ (13)

$v^{lc}_{CM}+v^{pc}_{CM}+v^{pp}_{CM}=-dA_{CM}.D_{C}.\frac{\partial
C_{cyt}}{\partial u}$ (14)

In membrane surrounding the cell body, the existing fluxes are calcium through
passive channels and the PMCA pump. Considering again the movement inside the
cell occurring by diffusion, the boundary condition at this location is given
by Eq. 15.

$v^{pc}_{CM}+v^{pp}_{CM}=-dA_{CM}.D_{C}.\frac{\partial C_{cyt}}{\partial u}$
(15)

Finally, in membrane surrounding the endoplasmic reticulum there are calcium
fluxes through channels activated by IP3, and through passive channels and
pumps. In addition, there is degradation flux of IP3. The boundary conditions
at this location, are given by Eq. 16 (IP3) and Eq. 17 (Ca2+), considering
again diffusion in cytosol.

$v^{d}_{ERM}=-dA_{ERM}.D_{I}.\frac{\partial I_{cyt}}{\partial u}$ (16)

$v^{lc}_{ERM}+v^{pc}_{ERM}+v^{pp}_{ERM}=-dA_{ERM}.D_{C}.\frac{\partial
C_{cyt}}{\partial u}$ (17)

The two sets of equations in the model, equations 3-9 and equations 13-17 are
complementary. We will show its application in the model using equations 4 and
16, which represent the degradation of IP3 in the ER membrane. Substituting
the first in the second and simplifying, we obtain:

$\frac{\partial I_{cyt}}{\partial
u}=-\frac{f^{d}_{ERM}.h^{d}}{D_{I}}.(I_{cyt}-I^{*}_{cyt}),u=x,y$ (18)

Eq. 18 means that any IP3 that arrives by diffusion in the ER membrane,
activates calcium channels and is degraded. The discretization of the
equations was performed using the Finite Difference Method. Considering a
point $(i,j)$ belonging to the mesh used in the cytosol to discretize Eq. 12,
and assuming that point $i$ belongs to the ER membrane, and that a regressive
discretization is used at this point, we would have:

$I_{i-1,j}=I_{i,j}+\alpha.(I_{i,j}-I^{*}_{cyt})$ (19)

where the dimensionless $\alpha=f^{d}_{ERM}.h^{d}.\Delta x/D_{I}$

Note that, in the discretization of Eq. 12, applied at point $i$ belonging to
the ER membrane, the concentration of I appears in position $(i-1,j)$. This
position does not belong to the equation domain, which is why it needs to be
replaced by the boundary condition.

###  3.2 Opening probability

The probability of channel opening requires an additional explanation about
how it was code implemented. This is necessary because probability value
calculation is not enough to determine whether the channel will be opened or
not. Another criterion must be established in order to decide the event
occurrence. It was decide the channel would be opened only if a
number,$\gamma\in$ (0;1), randomly chosen from an equiprobable distribution,
was less than the probability calculated for each computational step by Eq. 10
or Eq. 11. The channel, once opened, remains on this condition for a period of
time, $t_{o}$. After that, it is closed for another period, corresponding to
the refractory period, $t_{r}$. Under these specific conditions there is no
further need to calculate the probability, because $p^{o/c}_{m}$ is equal to
one or zero, respectively.

However, the probability of opening is dependent on the time interval used for
each computational step. This is because, if two time intervals were used, the
first being double the second, in the second the number of tests doubles and
the probability of opening also doubles. To avoid this, all simulations were
carried out with the same time interval, $\Delta t=10^{-6}s$

###  3.3 Model parameters values

It was assumed that the rest concentrations of calcium and IP3 in cytosol and
calcium in the ER correspond to 0.1 $\mu$M, 0.16 $\mu$M and 10.3 $\mu$M,
respectively (Ullah et al., 2006). In addition, calcium concentration in the
extracellular fluid corresponds to 2 mM (Bear et al., 2007). Calcium
diffusivity is 25 $\mu m^{2}.s^{-1}$ and the IP3 diffusivity, 280 $\mu
m^{2}.s^{-1}$ (Höfer et al., 2002).

The volumetric fraction occupied by reticulum was estimated based on data
provided by Alberts et al. (2017). The authors declared that cytosol occupies
54% of cell total volume. Thus, all organelles occupy 46% of cell body. As the
cell fraction occupied by ER is about 17%, the organelles fraction
corresponding to ER is 37%.

The calculation of the area fraction occupied by metabotropic receptors was
based on the estimation of receptors density presented by Nadkarni and Jung
(2005): 300 $\mu m^{-2}$. According to this information and because this type
of receptor has one pore of approximately 5 nm of diameter (Traynelis et al.,
2010), the $f^{p}_{CM}$ parameter was calculated (0.006). The same value was
used to $f^{lc}_{CM}$, $f^{lc}_{ERM}$ and $f^{d}_{ERM}$ for lack of better
information. The $k_{p,C}$ value (0.3 $\mu$M) was extracted from Höfer et al.
(2002) and $n$ (0.3) and $k_{p,G}$ (0.78 $\mu$M) were from Ullah et al.
(2006).

The constant associated to ER passive channels, $f^{pc}_{ERM}$.$h^{pc}_{ERM}$,
was estimated by equaling the flux through channels used by Ullah et al.
(2006) to the rest flux used in this work. The obtained value was 0.12 $\mu
M.s^{-1}$. Same value was used for cytoplasmic membrane channels,
i.e.$f^{pc}_{CM}$.$h^{pc}_{CM}$=0.12 $\mu M.s^{-1}$. This flux is ten times
lower for cell body. Otherwise, the exit would be massive and the restitution
period for ER would be too long, which is not biologically feasible.

The pump saturation constant, $k^{pp}$ = 0.1 $\mu$M, was also used by Ullah et
al. (2006). The constants $f^{pp}_{ERM}$. $K^{pp}_{ERM}$ and
$f^{pp}_{CM}$.$K^{pp}_{CM}$ were calculated by equaling the flux of passive
channels and the flux of pumps at resting. Similarly, to passive channels,
this flux is ten thousand times lower in cell body membrane.

The opening probability constants for IP3-activated channels were extracted
from Tu et al. (2005a) and Tu et al. (2005b), considering that the astrocytes
present, mainly type 2 IP3 receptors (Holtzclaw et al., 2002). Thus, $ks_{I}$
= 0.1 $\mu M$; $n_{I}=2.2$; $ks_{C}$ = 0.16 $\mu M$; $n_{C}=2.05$ and $Ks_{C}$
= 0.16 $\mu M$. For glutamate-activated channels, initially, $n_{G}$ equal to
$n_{I}$ was adopted. The value $ks_{G}$ = 1000 $\mu$M was determined based on
knowledge of interval in which glutamate concentration varies in the synaptic
cleft. The concentration of glutamate at rest in the cleft is 3 $\mu$M. Given
that, on average, 3000 glutamate molecules are released in cleft during
synaptic transmission and cleft has a volume of 2.0 x $10^{6}$ $nm^{3}$ (Allam
et al., 2012), the average concentration of this compound reaches 3003 $\mu$M.
Considering the glutamate release through the vesicle in a synapses central
point, this concentration in the areas surrounding the astrocyte has an
amplitude of $\delta$=1000 $\mu$M (Allam et al., 2012). By plotting a curve
with 1000, 2000, 3000 $\mu$M for $ks_{G}$ value, the most coherent behavior
was found at 1000 $\mu$M. Refractory period (10.8 ms) and opening maintenance
time (7.6 ms) of each channel were extracted from Tu et al. (2005b).

The constants associated with production and decay rates of IP3, $K^{p}$ = 6.0
x $10^{-12}$ $\mu mol.\mu m^{2}/s$ and $h^{d}$ = 6.0 x $10^{3}$ $\mu
m.s^{-1}$, were adjusted to make calcium signaling occurs in about 1.0 s as
indicated in Araque et al. (2014). The $h^{lc}_{ERM}$ (150 $\mu m/s$) and
$h^{lc}_{CM}$ (0.3 $\mu m/s$) coefficients were adjusted assuming to be
relevant the reticulum contribution to increase of intracellular calcium
concentration and considering the comparative study with a calcium signaling
presented by Araque et al. (2014), this concentration doubles. It is
noteworthy that after calculation of relation between $v^{lc}_{CM}$ and
$v^{lc}_{ERM}$ flows obtained by Ullah et al. (2006), $h^{lc}_{CM}$
corresponded to 0.002 x $h^{lc}_{ERM}$.

Now, after established the global parameters, the results obtained with the
model and its potentiality will be exposed.

Table 1: Parameters description and adopted values Symbol |  Description |  Range |  Value |  Unit
---|---|---|---|---
Rest concentrations
$C_{cyt}^{*}$ |  Calcium concentration at rest in cytosol |  |  0.1 |  $\mu M$
$C_{ER}^{*}$ |  Calcium concentration at rest in ER |  |  10.3 |  $\mu M$
$C_{sc}$ |  Calcium concentration in synaptic cleft |  |  2.0 |  $mM$
$I_{cyt}^{*}$ |  IP3 concentration at rest in cytosol |  |  0.16 |  $\mu M$
$G_{sc}^{*}$ |  Glutamate concentration at rest in synaptic cleft |  |  3.0 |  $\mu M$
Diffusivities
$D_{C}$ |  Calcium diffusivity in cytosol |  |  25 |  $\mu m^{2}s^{-1}$
$D_{I}$ |  IP3 diffusivity in cytosol |  |  280 |  $\mu m^{2}s^{-1}$
IP3 production
$f^{p}_{CM}$ |  Cell membrane area fraction where occurs IP3 production |  |  0.006 |  -
$K^{p}$ |  Maximum rate of IP3 production |  1 to $20x10^{-12}$ |  $6.0x10^{-12}$ |  $\mu mol\mu m^{-2}s^{-1}$
$k_{p,C}$ |  Dissociation constant for calcium stimulation of IP3 production |  |  0.3 |  $\mu M$
$k_{p,G}$ |  Dissociation constant for glutamate stimulation of IP3 production |  |  0.78 |  $\mu M$
$n$ |  Hill coefficient for IP3 production |  |  0.3 |  -
IP3 degradation
$f^{d}_{ERM}$ |  Area fraction of ER membrane where occurs IP3 degradation |  |  0.006 |  -
$h^{d}$ |  IP3 degradation coefficient |  1 to 20 x $10^{3}$ |  6.0 x $10^{3}$ |  $\mu ms^{-1}$
Calcium fluxes through ligand channels
$f^{lc}_{CM}$ |  Cell membrane area fraction occupied by ligand channels |  |  0.006 |  -
$f^{lc}_{ERM}$ |  ER membrane area fraction occupied by ligand channels |  |  0.006 |  -
$h^{lc}_{CM}$ |  Mass transfer coefficient for cell membrane ligand channels |  |  0.3 |  $\mu ms^{-1}$
$h^{lc}_{ERM}$ |  Mass transfer coefficient for ER membrane ligand channels |  10 to 200 |  150 |  $\mu ms^{-1}$
Calcium fluxes through passive channels
$f^{pc}_{CM}$ |  Cell membrane area fraction occupied by passive channels |  |  0.001 |  -
$f^{pc}_{ERM}$ |  ER membrane area fraction occupied by passive channels |  |  0.001 |  -
$h^{pc}_{CM}$ |  Mass transfer coefficient for cell passive channels |  |  120 |  $\mu ms^{-1}$
$h^{pc}_{ERM}$ |  Mass transfer coefficient for ER passive channels |  |  120 |  $\mu ms^{-1}$
Calcium fluxes through pumps
$f^{pp}_{CM}$ |  Cell membrane area fraction occupied by calcium pumps |  |  0.001 |  -
$f^{pp}_{ERM}$ |  ER membrane area fraction occupied by calcium pumps |  |  0.001 |  -
$K^{pp}_{CM}$ |  Pump flux maximum rate at cell membrane |  |  0.3 |  $\mu mol\mu m^{-2}s^{-1}$
$K^{pp}_{ERM}$ |  Pump flux maximum rate at ER membrane |  |  0.3 |  $\mu mol\mu m^{-2}s^{-1}$
$k_{pp}$ |  Pump saturation constant |  |  0.1 |  $\mu M$
Channel open probability
$ks_{C}$ |  Apparent affinity for calcium inhibitory site |  |  0.16 |  $\mu M$
$Ks_{C}$ |  Apparent affinity for calcium activating site |  |  0.16 |  $\mu M$
$ks_{I}$ |  Apparent affinity for IP3 activating site |  |  0.1 |  $\mu M$
$ks_{G}$ |  Apparent affinity for glutamate activating site |  |  1000 |  $\mu M$
$n_{C}$ |  Hill coefficient to calcium |  |  2.05 |  -
$n_{I}$ |  Hill coefficient to IP3 |  |  2.2 |  -
$n_{G}$ |  Hill coefficient to glutamate |  |  2.2 |  -
Channel opening and closing
$t_{r}$ |  Ionic channels refractory time |  |  10.8 |  $ms$
$t_{o}$ |  Ionic channels opening time |  |  7.6 |  $ms$
$t^{p}_{test}$ |  Computing time for test probability |  |  10.0 |  $ms$
Stimulus characteristics
$\sigma$ |  Glutamate oscillation amplitude in synaptic cleft during a stimulus |  |  10.0 |  $ms$
$\delta$ |  Glutamate oscillation amplitude during a stimulus |  |  1000 |  $\mu M$

##  4 Results and discussion

The calcium and IP3 dynamics in the cytosol, as well as calcium dynamic in the
ER, were obtained from numerical simulation of the proposed model and they are
shown in this section. For this simulation, it was developed a computational
interface using FORTRAN® language.

The free parameters of the model were adjusted aiming to reach medium values
indicated by Araque et al. (2014) for calcium signaling amplitude and
duration. It was considered an astrocyte branch extremity stimulated due to
passage of a signaling from a synapse which it surrounds. Here, the stimulus
was considered glutamatergic, making necessary to consider glutamate
concentration in synaptic cleft near the astrocyte. According to Budisantoso
et al. (2013) the concentration of Glutamate in the synaptic cleft grows
exponentially close to the opening of a vesicle. The article also shows that,
due to the spread of the neurotransmitter caused by diffusion, this growth
softens as we move away from this point. Next to an astrocyte connected to the
synapse, at a distance greater than 100 nm from the starting point, the growth
of the glutamate concentration curve approaches a normal curve. The decay of
the curve depends on diffusion and on the density of m-GluR and i-GluR
channels. Given the difficulty of obtaining data on the dynamic behavior of
glutamate close to the astrocytic termination, it was proposed a sine wave
function shown in Eq. 20. This function has a behavior similar to that
suggested by Budisantoso et al. (2013) close to the astrocyte and guarantees a
stimulus with a duration of 10 ms ($\sigma$) and a concentration amplitude of
the order of $\delta$. Another advantage is to easily allow the generation of
stimulus sequences, or a burst. The sine function represents an increase
followed by a decrease in glutamate concentration. For negative values of this
function, the synaptic cleft glutamate concentration assumes the rest
concentration value.

$G_{sc}=G^{*}_{sc}+\delta.sin(\frac{2\pi}{\sigma}.t),\sigma=10ms$ (20)

For the graphs shown in this section, it was used the following dimensionless
variables:

$\overline{y}^{C}_{cyt}=\frac{\overline{C}_{cyt}}{C^{*}_{cyt}},\overline{y}^{C}_{ER}=\frac{C_{ER}}{C^{*}_{ER}},\overline{y}^{I}_{cyt}=\frac{\overline{I}_{cyt}}{I^{*}_{cyt}}$
(21)

###  4.1 Cell calcium signalling

The curves representing average concentration profiles of Ca2+ and IP3 in the
cytosol and Ca2+ in the ER (Fig. 4) show coherent qualitative behavior. The
glutamatergic stimulation activates metabotropic receptors, triggering IP3
production, which moves through diffusion to ER, where opens calcium channels.
If the metabotropic receptors deactivate and IP3 production reduces, the
concentration increases more slowly until the outflux from the cell and the
return flux to the ER become more relevant, reducing intracellular calcium
concentration. These phenomena can be observed on Fig. 4, which shows
concentration profiles after stimulus.

Figure 4: Calcium and IP3 dynamics. The graph shows calcium and IP3 mean
concentrations in cytosol and calcium concentration in endoplasmic reticulum.
IP3 concentration increases due to its production triggered by activation of
metabotropic receptors in cell membrane. When it reaches endoplasmic
reticulum, IP3 starts to be degraded and its concentration decreases until it
returns to resting value. Cytosol calcium concentration increases mainly
because IP3 reaches ER and opens calcium channels IP3-activated, allowing its
release into cytosol. Calcium concentration in ER decreases because calcium is
released into cytosol. Both calcium concentration in cell body and in ER
return to resting value due to calcium activity pumps and to passive channels.

The results obtained can be compared to the ones published in several
experimental papers like Venance et al. (1997) and Ullah et al. (2006).
However, the experimental conditions used, such as astrocyte cultures and
long-term stimulations, in opposition to the situation simulated in this work,
from one isolated astrocyte, stimulated only at one of its extremities, can
only allow qualitative comparisons.

The diversity of data in the literature, obtained from different cells,
regions and experimental methods, makes it impossible to establish precise
values for model parameters. In addition, most of them are related to specific
local situations. An example is the action of a specific channel that depends
not only on its dynamics, but also on its density and, most of times, on its
availability in the cell membrane. In this context, it is important to
evaluate the relative influence of each parameter of the model on calcium
signaling through a parametric sensibility analysis.

###  4.2 Parameter sensitivity

After an analysis of Araque et al. (2014) paper, it is possible to infer that
calcium cytosolic average concentration doubles in response to an elementary
stimulation. It is still possible to notice that, for this case, signaling
time is about 1.0 s. Thus, the choice of model free parameters values aimed
that calcium cytosolic average concentration at rest doubled and signalization
time was, approximately, 1.0 s.

In the model proposed here there are three free parameters, the constants
associated to IP3 production and degradation, $K^{p}$ and $h^{d}$,
respectively, and the mass transfer coefficient through IP3 activated channels
in ER membrane, $h^{lc}_{ERM}$. Initially is presented a parameter analysis
considering the desired values for signaling time and calcium cytosolic
average concentration. The deviation, here designed, $E$, was calculated by
the least squares method:
$E=(\overline{y}^{C}_{cyt,max}-2)^{2}+(t_{C}-1)^{2}$. In the equation,
$\overline{y}^{C}_{cyt,max}$ is the maximum value reached by dimensionless
average concentration of Ca2+ cytosolic, and $t_{C}$ is the time for return of
this concentration to the rest value.

Fig. 5 shows the results obtained for $10^{12}$ x $K^{p}\in[1;10]$, $10^{-3}$
x $h^{d}\in[1;15]$ and $h^{lc}_{ERM}$ = 100, 150 and 200 $\mu m.s^{-1}$. Fig.
5a and Fig. 5b show $E$ for two values of $h^{lc}_{ERM}$, clearly
demonstrating its dependency of the free parameters of the model. In Fig. 5c
the indicated sets represent lowers values for $E$. For $h^{lc}_{ERM}=150$
$\mu m.s^{-1}$, Fig. 5d shows the maximum values obtained for IP3 cytosolic
average concentration. The continuous line links points to lowers values of
$E$.

Figure 5: Effect of three free parameters in values of
$\overline{y}^{C}_{cyt,max}$ and ${t}_{C}$. (a) Error for $h^{lc}_{ERM}=200$
$\mu m.s^{-1}$. (b) Error for $h^{lc}_{ERM}=150$ $\mu m.s^{-1}$. It may be
noted several parameter sets make possible to obtain error values close to
zero. Those sets are shown in part (c). The four highlighted points correspond
to lowest error values obtained. The set chosen was $h^{lc}_{ERM}=150$ $\mu
m.s^{-1}$, $h^{d}$ = 6 x $10^{3}\mu m.s^{-1}$ and $K^{p}$ = 6 x $10^{-12}$
$\mu mol.\mu m.s^{-1}$. (d) Maximum values of IP3 mean concentration in
cytosol for $h^{lc}_{ERM}=150$ $\mu m.s^{-1}$. Each dotted line joins points
for same value of $h^{d}$, considering just odd values varying from 1 to 15 x
$10^{3}$. Continuous line joins points with lowest E for $h^{d}$ ranging from
1 to 15 x $10^{3}$.

It is important to highlight that, for sets containing lower values of $E$,
the IP3 cytosolic average maximum concentration is restricted to a narrow
range. This characteristic is interesting when assessed together with the
linear behavior of the relation between the parameters that allow fewer
errors.

Considering the results presented at Fig. 5, the values chosen for the free
parameters were: $K^{p}$ = 6 x $10^{-12}$ $\mu mol.\mu m.s^{-1}$, $h^{d}$ = 6
x $10^{3}$ $\mu m.s^{-1}$ and $h^{lc}_{ERM}=150$ $\mu m.s^{-1}$.

Now, after the analysis of the parameter set, it will be exposed how each one
affects the system behavior.

####  4.2.1 IP3 production

The global coefficient of IP3 production rate affects directly the maximum
values of calcium and IP3 average concentrations in cytosol. Fig. 6a and Fig.
6b show temporal profiles of those concentrations. It can be observed that
curves pattern and restitution time does not change. Fig. 6c shows that, as
the value of $K^{p}$ increases, the cytosolic calcium maximum concentration
value increases until it reaches a limit, from which any increase in the
analyzed parameter does not cause relevant changes in calcium concentration
peak.

Figure 6: Dynamics of Ca2+ and IP3 mean concentrations in cytosol as function
of global coefficient of IP3 production rate. a Dynamics of dimensionless mean
concentrations of Ca2+ in the cytosol after the glutamatergic stimulation. The
dotted lines show the $\overline{y}_{cyt}^{C}$ dynamics for several values of
$K^{p}$ (1 to 19 x 10-12, considering just the odd values). The continuous
line shows the calcium dynamics for $K^{p}$ = 6 x $10^{-12}\mu mol.\mu
m.s^{-1}$. b Dynamics of dimensionless mean concentrations of IP3 in the
cytosol. The dotted lines show the IP3 dynamics for odd values of $K^{p}$,
from 1 to 19 x 10-12. The continuous line shows the IP3 dynamics for $K^{p}$ =
6 x $10^{-12}$ $\mu mol.\mu m.s^{-1}$. For both concentrations,
$\overline{y}_{cyt}^{C}$ and $\overline{y}_{cyt}^{i}$, the maximum values are
directly proportional to $K^{p}$. c Maximum values of dimensionless Ca2+
concentrations in cytosol for values of $K^{p}$ from 1 to 20 x 10-12. For
$K^{p}$ = 6 x $10^{-12}$ $\mu mol.\mu m.s^{-1}$, the maximum concentration
reaches the double of resting concentration. d Derivative of curve of the
maximum values of dimensionless Ca2+ concentrations in cytosol in relation to
$K^{p}$.

The observed behavior in Fig. 6c is expected, because the increase in IP3
production results in a higher concentration of this compound in cytosol (Fig.
6b), which allows the activation of a greater number of channels in ER,
releasing more calcium. This release is limited by quantity of channels
present in ER membrane, which reduces this compound influence in calcium
concentration peak, after specific value. Fig. 6d shows the derivative graph
of maximum calcium curve (Fig. 6c) in relation to $K^{p}$. From this last
curve, it is possible to conclude the model is more sensitive to parameter
$K^{p}$ in interval of 1 to 10 x 10-12 $\mu mol.\mu m.s^{-1}$. Sensitivity is
significantly lower above the last value. The value chosen, $K^{p}$ = 6 x
$10^{-12}$ $\mu mol.\mu m.s^{-1}$, is within this region. It is important to
show that IP3 production requires energy expenditure (O’Hayre et al., 2013),
reason to choose values resulting in lower producing for same desired limits
of cytosolic calcium.

####  4.2.2 IP3 degradation

The IP3 degradation, as its production, affects cytosolic calcium
concentration peak, as shown in Fig. 7a and Fig. 7c. How much higher is the
parameter associated to degradation, lower is the value reached by calcium
concentration. Fig. 7d shows restitution time for cytosolic Ca2+ average
concentration. Observed restitution time varies significantly for low values
of $h^{d}$, reaching a region with small variations for higher values. The
beginning of this region can be defined in $h^{d}=6000$ $\mu m.s^{-1}$. As
shown in Fig. 7c, at this point, maximum average concentration reaches twice
the rest value. That is the chosen value, because it reduces restitution time
and maximizes cytosolic Ca2+ maximum concentration.

Figure 7: Dynamics of Ca2+ and IP3 mean concentrations as function of IP3
degradation coefficient. Dynamics of a Ca2+ and b IP3 dimensionless mean
concentrations in cytosol for odd values of $10^{-3}$ x $h^{d}$ between 0 and
20. For both concentrations, $\overline{y}_{cyt}^{C}$ and
$\overline{y}_{cyt}^{i}$, the maximum values are inversely proportional to
$h^{d}$. c Maximum values and d restitution time for cytosolic Ca2+ mean
concentration for values of $10^{-3}$ x $h^{d}$ from 1 to 20.

####  4.2.3 Ca2+ release

Calcium release from endoplasmic reticulum depends on mass transfer
coefficient through ligand dependent channels, $h^{lc}_{ERM}$. Fig. 8a shows
its influence in cytosolic Ca2+ average concentration evolution. Fig. 8c shows
maximum values of this variable for different values of $h^{lc}_{ERM}$,
evidencing calcium release is directly proportional to this parameter. The
constant growth observed at Fig. 8c is justified by a higher mass transfer
coefficient, producing more calcium exit from ER to cytosol. As the process is
limited by ER concentration, which is significantly high, no saturation were
observed in analyzed interval. Also IP3 cytosolic average concentration is not
affected by this parameter (Fig. 8b). Fig. 8d shows little influence of it in
relation to ER Ca2+ minimum concentration. It is important to mention that
calcium concentration in ER reduces about 4.0% when in cytosol increases 100%.
The value chosen for the parameter, $150$ $\mu m.s^{-1}$, is, again,
associated to desired calcium maximum cytosolic concentration ($0.2$ $\mu M$).

Figure 8: Dynamics of dimensionless Ca2+ mean concentrations in cytosol and in
ER as a function of mass transfer coefficient of ligand channels of ER
membrane. (a) Dynamics of Ca2+ mean concentration in cytosol for $10^{-3}$ x
$h^{lc}_{ERM}\in[10;200]$. The continuous line correspond to $10^{-3}$ x
$h^{lc}_{ERM}=150\mu m.s^{-1}$. (b) Dynamics of IP3 mean concentration in
cytosol. This dynamics does not depends on the values of $h^{lc}_{ERM}$. (c)
Maximum values in cytosol and (d) minimum values in ER for dimensionless Ca2+
concentration.

###  4.3 Implication of TRPC incorporation

Incorporation of the TRPC, while not altering cellular dynamics, helps Ca2+
concentration in ER to return to resting value in a reasonable time interval
(Malli et al., 2007; Verkhratsky and Parpura, 2014). However, as shown in Fig.
9, the effect is not significant in case of a single glutamatergic stimulation
lasting 0.01 s. By contrast, in case of burst, with 10 stimulations in
sequence, TRPC incorporation allows calcium baseline concentration returning
to rest value in ER. It means this incorporation is important for an astrocyte
immersed in cell medium, with several connections and submitted to several
stimuli.

Figure 9: Comparison between dynamics of calcium concentration in ER, obtained
with and without the activation of TRPC channels. In the figure, the white
circles represent the dynamic without the activation of the TRPC channels. The
black circles represent the activation of this mechanism. a Curves obtained in
the case of one single stimulation, with 0.01 s duration. b Curves obtained in
the case of a burst, with a series of 10 stimulations of 0.01 s each one.

###  4.4 Opening probability

The use of $P_{I}$ probability (Eq. 10) for opening IP3-activated channels in
endoplasmic reticulum membranes requires their individualization and
numbering. There are one hundred equally spaced channel positions distributed
along the ER membrane. Fig. 3 shows numbering starting in opposite face and
going around membrane clockwise. When occur stimulus in a synapse connected to
an astrocyte arm, the ER side facing that arm has channels with positions 38
to 62, 50 being the central position. On the opposite side, the channel in
central position and farthest from stimulation has number 100. The channels do
not open in sequence but obey the probability test. It is possible to see at
Fig. 10a opening times for all channels, when parameter values are those
described on Section 4.2. An analysis of parameters influence at Eq. 8,
$n_{I}$, $n_{C}$, $ks_{C}$, and $Ks_{C}$, shows that $P_{I}$ is mainly
affected by the first one. At Fig. 10b it is possible see this effect. The
values used for exponent $n_{I}$ were 1.0, 2.2, 5.0 and 10.0. The parameter
reduces probability and increases opening time.

Figure 10: IP3 activated channel opening time in different points in the
contour of ER membrane. a Opening time for all ER channels, with parameter
values indicated at Sect. 4.2. b Opening time for the same channels, for
several values of $n_{I}$. The values of this parameter, from the lighter to
darker points are 1.0, 2.2, 5.0 and 10.0, respectively. As the parameter value
increases, the opening probability reduces and the channels take more time to
open.

###  4.5 Influence of ER in branches

ER is found in astrocyte branches (Golovina and Blaustein, 2000; Golovina,
2005; Holtzclaw et al., 2002). However, its contribution to calcium dynamics
is still unclear (Oheim et al., 2018). In order to study this influence, the
computational program used in this work was modified to allow the allocation
of reticulum in the branches.

Initially, calcium dynamics was simulated in only one branch to verify how
reticulum would influence it. The branch had 6 $\mu m$ length and 1 $\mu m$
width. The ER therein was placed at a distance of 0.4 $\mu m$ from the
astrocyte membrane around the synapse and had 1 $\mu m$ length. Other lengths
were tested as well: 1.5 $\mu m$, 2.0 $\mu m$, 2.5 $\mu m$, and 3.0 $\mu m$.

First we analyzed the ER effect in astrocyte branch in IP3 dynamics, because
this is mainly responsible for propagation of calcium signaling. It was
observed that ER decreases the IP3 mass that leaves the branch, as can be seen
in Fig. 11. All ER sizes cause this same effect; however, the larger the
endoplasmic reticulum size, the smaller the amount of IP3 leaving the
branching. This is due to the IP3 degradation phenomenon that occurs in ER
contours. The larger the endoplasmic reticulum, the greater the degradation.
Consequently, the mass of IP3 available to go ahead in the cell and to release
calcium from stocks is lower.

Figure 11: Effect of ER length in the mass of IP3 leaving the astrocyte
branch. The curves show the mass sum over time that leaves the astrocyte
branching for each reticulum size evaluated.

As showed in Fig. 11, endoplasmic reticulum allocation to astrocyte branch may
not be good for global Ca2+ signaling, because the ER decreases the IP3
available amount for opening calcium channels in cell body ER. Thus, it is
necessary to determine why cells would have ER in astrocyte branchs.

Calcium has local signals not spread to whole cell (Di Castro et al., 2011;
Kanemaru et al., 2014). These signals are short and fast. The endoplasmic
reticulum may contribute to them. Fig. 12 shows the dynamics of dimensionless
calcium concentration in branches. From this figure, it is possible to
conclude that ER promotes local calcium signaling, increasing the levels of
Ca2+ in a smaller time interval.

Figure 12: Effect of ER length in the dimensionless Ca+2 concentration in
astrocyte branch. The curves show the dimensionless Ca+2 concentration over
time to each reticulum size evaluated.

Calcium dynamics observed in Fig. 12 can be compared to experimental results
of Di Castro et al. (2011) with respect to local calcium events. Studying
astrocyte branches, researchers concluded that local signaling extends at
distances of 3 to 8 $\mu m$, having 0.14 s rise time and 0.7 s duration. In
this work, the branch has 6 $\mu m$, calcium concentration rise time was 0.2 s
and duration was 0.8 s. These results can validate the parameters set defined
in the previous simulation.

In addition to contributing to local calcium signaling, ER could contribute to
calcium wave advancement, because the released calcium go ahead in the cell
when reaches the astrocyte branch end. As can be seen in Fig. 11 and Fig. 13,
the amount of calcium leaving the branching is smaller and slower than IP3;
therefore, this calcium possibly does not contribute significantly to calcium
wave advancement.

Figure 13: Effect of ER length in the mass of Ca+2 leaving the astrocyte
branch. The curves show the mass of Ca+2 that leaves the astrocyte branch over
time to each reticulum size evaluated.

The production of IP3 depends on calcium concentration surrounding
metabotropic receptors. Another effect of ER on astrocyte branch endfeet could
be enhancing IP3 production. In the cases discussed above, this effect is
imperceptible, because the time for IP3 production is shorter than time IP3
takes to reach endoplasmic reticulum to release calcium and for calcium to
arrive in the membrane to strengthen IP3 production. It was necessary to
evaluate a burst corresponding to successive synaptic transmissions. Three
stimuli were performed in order to determine whether the calcium released in
the first stimulus was able to influence IP3 production in second and third
stimuli. An insignificant effect was found. IP3 production increases by about
0.001%.

In addition to reticulum size, endoplasmic reticulum position may influence
IP3 production. To evaluate the ER position influence on IP3 production, the
smallest reticulum tested (1.0 $\mu m$) was allocated at 0.4 $\mu m$, 0.8 $\mu
m$ and 1.4 $\mu m$ away from cell membrane where IP3 is produced. Again, the
effect on IP3 production was less than 0.001%.

With this considerations it is possible to verify that ER allocated in
astrocyte branches mainly contributes to local calcium signaling. In relation
to global signaling, its effect is negative because it reduces amounts of IP3
available for channel opening in cell body ER. It was also possible to note
that increasing calcium local concentration does not significantly influence
IP3 production.

To analyze parameters influence in this new simulation, the same parameters
from previous simulation were evaluated: production ($K^{p}$) coefficient,
degradation ($h^{d}$) coefficient and calcium release ($h^{lc}_{ERM}$)
coefficient. The results can be seen in Fig. 14.

Figure 14: Effect of the parameters in the calcium flux through the channels
coupled to ionotropic receptors. Influence of the parameters related to the
production of IP3, degradation of IP3 and release of calcium in the calcium
flux through the channels coupled to ionotropic IP3 receptors present in the
endoplasmic reticulum. The original parameters values are $K^{p}$ = 6.0 x
$10^{-12}\mu mol$, $\mu m^{-2}$$s^{-1}$, $h^{d}$ = 6.0 x $10^{3}$ $\mu
m$$s^{-1}$, $h^{lc}_{ERM}$ = 150 $\mu m$$s^{-1}$.

As previously discussed, parameters related to IP3 have opposite effects: the
higher the production coefficient, the greater the calcium amount released by
ER. The higher the degradation coefficient, the less calcium is released from
ER. Fig. 14 shows that the influence of these parametrs is not significant
when evaluating the flux through channels coupled to ionotropic IP3 receptors.
Regarding the calcium release parameter, the model sensitivity is much higher.
This is because it is directly related to calcium output flux. If the
parameter doubles, calcium output flux doubles. If the parameter is halved,
the output flux is reduced by half as well.

Another interesting aspect to emphasize is channel opening process. In the
cell body ER, IP3 arrives slowly, so channels are opening sequentially (Fig.
10). In branching, by contrast, given the high IP3 diffusion coefficient and
the small space for movement, channels open almost simultaneously. This
phenomenon is evidenced in Fig. 15. For all ER sizes, all channels are opened
simultaneously when the concentration of IP3 is high. When IP3 concentration
decreases around 0.2 s, it is possible to observe more clearly the opening
probability effect in the channels.

Figure 15: Calcium flux through the calcium channels coupled to ionotropic
receptors of IP3. The channels coupled to the ionotropic receptors of IP3 open
simultaneously in the endoplasmic reticulum. The larger the reticulum, the
greater the number of channels and, consequently, the greater the outflux of
calcium. With the passage of time, the number of channels that opens is
smaller and the effect of the probability of opening becomes evident.

It is also possible to observe in Fig. 15 small decreases in output flux over
time despite the fact that all channels open simultaneously. This can be best
understood by observing flux over time through a single channel (Fig. 16). As
the gradients diminish (increasing calcium concentration in the cell and
decreasing calcium concentration in the ER), the flux also decreases.

Figure 16: Calcium flux through a single channel coupled to an ionotropic IP3
receptor. The flux of calcium through a channel, for multiple stimuli,
decreases over time due to decreased gradients. The blank intervals represent
the periods of time that the channel remains closed.

Computational simulation of an astrocyte branch allowed analysis of ER effects
on local calcium signaling. Analysis of its effects on global calcium
signaling requires a more complex geometry. Therefore, we simulateed the
geometry shown in Fig. 17, containing a cell body and two branches connected
to it. The two branches contain ER. This configuration was named L1R1 (one ER
at left branch and one ER at right branch) to facilitate results presentation.
For better result analysis, we tested configurations with the ER in only one
of the branches (L1R0 and L0R1). Finally, for comparison purposes, we also
evaluate the situation with branches without reticulum (L0R0).

Figure 17: Geometric configurations implemented computationally in way to
analyze endoplasmic reticulum effects on astrocytes branches. The figure shows
the four configurations used. From top to bottom, L1R1 has two endoplasmic
reticulum, in left and right branches; L1R0 only has the reticulum in the left
branch; L0R1 only in the right branch; L0R0 has no endoplasmic reticulum in
the branches. The results generated by these four configurations were analyzed
and compared. The length of the endoplasmic reticulum is 1 $\mu m$. Figure 18:
Instantaneous IP3 outflux at the end of the right branch. The graph shows the
outflux of IP3 on the branch opposite the branch where the stimulus is
applied. The flux is smaller for the situation where there are ER in the two
branches, intermediate for the situations with ER in only one of the branches
and greater for the case in which the branches do not have ER.

First we evaluated branches endoplasmic reticulum effect on IP3 production. It
has been observed that configurations L1R1 and L1R0 have 0.000086% greater IP3
production than L0R1 and L0R0. This shows that reticulum located near stimulus
increases IP3 production. Although the increase is not significant, it does
occur. This IP3 production increase occurs because calcium released from
reticulum by IP3 travels to cell membrane and therein can potentiate IP3
production (Pawelczyk and Matecki, 1997). This can be evidenced when the first
stimulus does not cause distinction among fluxes produced in the four cases.
This distinction can be noticed in subsequent stimuli (results not shown).
What makes so small IP3 production increase is the slow temporal calcium
dynamics. Thus, it is possible to state that reticulum at branch ends does not
significantly contribute to increase IP3 production.

With respect to IP3 output flux, that is, to IP3 that would advance in the
cell or even in astrocytes network, it is possible to conclude that the most
favorable configuration would be the one in which there is no endoplasmic
reticulum in the branches (Fig. 18). This is because the ER promotes this
compound degradation. The IP3 degradation will be proportional to endoplasmic
reticulum amount in cell.

The L1R0 and L0R1 configurations have different IP3 flux profiles at the right
branch end. The reason is that, at the stimulated end, L1R0 configuration, a
large concentration of IP3 reaches the ER and is degraded. The degradation of
IP3 is proportional to its concentration in the cytosol, close to the ER
membrane. In the L0R1 configuration, the IP3 flow that reaches the last ER is
diluted by diffusion and because it has already been partially degraded in the
central ER. The lower concentration in the proximity of the last ER membrane
leads to less degradation of IP3 and a greater amount reaches the outlet at
the end of the right branch.

Regarding calcium, it is interesting to note calcium flux in cell membrane in
which the stimulus occurs, showed in Fig. 19. In this region, there is a long
time outflux, because increase of intracellular calcium concentration promotes
pump action attempting to restore the resting state. In cases when ER is in
the left branch, outflux occurs faster and with great oscillation. This
oscillation is result of probabilistic dynamics of channel opening in the
reticulum. As the ER is close to where flux is being measured, this
probabilistic effect is more clearly perceived. In cases where there is no ER
in vicinity of the membrane where stimulus occurs, calcium that arrives there
comes from cell body. As it needs to travel through cell body and branch to
leave the cell, the probabilistic effect is no longer evident.

Figure 19: Instantaneous calcium flux in the cell membrane where stimulus
occurs. Calcium flux in left end membrane of the branch where the stimulus
occurs is outflux, so it is negative. When there is ER in this branch, the
flux is larger and more oscillatory. This oscillation is due to probabilistic
nature of channels opening. When there is no ER, diffusion dissipates this
probabilistic effect, no longer being perceived.

The calcium outflux in the branch end where the stimulus does not occur, right
branch, can be observed in Fig. 20. This figure shows two types of behavior.
The first behavior corresponds to configurations containing endoplasmic
reticulum in the right branch. The second behavior corresponds to
configurations without reticulum in this branch. The ER close to output causes
output stream to start earlier. This is because IP3 quickly hits this ER,
releasing calcium. When there is no reticulum, calcium leaving the branch end
comes from the ER of cell body. As calcium diffusion is a tenth of IP3
diffusion, it takes longer to come out in this configuration.

Figure 20: Instantaneous calcium flux at right branch outlet. The fastest
dynamics correspond to settings where there is right branch endoplasmic
reticulum, L1R1 and L0R1. The differences between configurations with or
without left branch reticulum are due to IP3 degradation by this ER, reducing
IP3 amount in ER of cell body. The conclusion is that ER in branches away from
stimulus sites accelerates calcium advancement in one individual astrocyte, or
in a network of them.

When calculating the accumulated flux over time, the sequence of calcium
outflux value is L0R1$>$L0R0$>$ L1R1$>$L1R0. Although calcium outflux occurs
quicker when ER is near of the exit, the total calcium amount leaving the cell
is larger when there is no ER in the other branch (Fig. 21). This is because
the left branch reticulum, in addition to releasing calcium, promotes IP3
degradation. In this way, IP3 available to reach right branch reticulum is
smaller, releasing less calcium. From these results, it is concluded that the
configuration allowing greater amount of calcium advance is the one without ER
in the first branch. This ER incorporation reduces about 17% of calcium amount
that would advance either in astrocyte network or in cell itself. Now, with
respect to the ER near to the exit, it increases only 1% the amount of calcium
coming out.

Figure 21: Accumulated calcium flux at right branch outlet. The accumulated
calcium flux at the astrocyte outlet can be seen in two aspects, quantity and
velocity. In all configurations IP3 amount produced by glutamatergic stimulus
is the same. Existence of endoplasmic reticulum in the branch near the
stimulus, L1R1 and L1R0 configurations, causes IP3 degradation, reducing the
amount that reaches the ER in central boddy of the cell. This decreases
calcium total amount that leaves the astrocyte at the stimulus opposite end.
In cell and reticulum here used dimensions, this reduction is approximately
14.5%. By contrast, when observing velocity, reticulum existence at the exit
end, L1R1 and L0R1 configurations, increases the leaving calcium flux
velocity. The reason is clear: IP3’s higher rate of diffusion allows it to
reach the last ER and release calcium before the arrival of calcium generated
in central ER.

There are two ways to analyze the effect of endoplasmic reticulum positioned
in the branches of astrocytes: one with respect to the amount of calcium and
IP3 available to advance in the cell itself or in a network of astrocytes, and
the other with respect to the speed of progress of these compounds. The
positioning of ER in the branches decreases the amount of IP3 that advances in
the cell. As a result, there is less calcium release in the ERs that are
further ahead. By contrast, the ER positioned in the branches ahead makes the
calcium dynamics faster, because IP3, moving quickly, reaches them earlier,
releasing calcium at the site.

One hypothesis about organelles location and function such as ER in astrocytes
branches involve diffusion barriers formation (Oheim et al., 2018). It is
possible to evaluate the calcium barrier effect when analyzing flux in
connection between branches and cell body. Fig. 22 shows flux through
connection between left branch and cell body. To understand the graph is
important to consider that positive fluxes indicate left to right and negative
fluxes indicate right to left direction. Because this flux is negative, the
calcium that passes through the connection comes from the ER of cell body.
This happens because IP3 quickly gets there and the central ER is closer to
the connection than the branch ER. Comparing the L1R0 and L0R0 configurations,
red and blue lines, respectively, it is possible to visualize the calcium
barrier, because the calcium flux towards the branch decreases considerably.

Given the calcium barrier existence, the flux through the connection between
cell body and right branch would be expected to be higher in L1R0 situation
than flux in L0R0 configuration. However, the flux in L0R0 is higher (Fig.
23). This occurs because less IP3 is made available to cell body when there is
reticulum in left branch.

Figure 22: Flux through connection between left branch and cell body. The
lines are coincident, differentiated only by the existence or not of
endoplasmic reticulum in left branch. L1R1 and L1R0 configurations have a
lower calcium retrograde motion than L0R1 and L0R0 configurations. This is
possibly for two reasons, whose effects add up. The ER in the left branch
degrades IP3 causing smaller amounts to reach central ER and can at the same
time produce a calcium barrier that reduces calcium return by diffusion.
Figure 23: Flux through connection between right branch and cell body. Black
and green lines show calcium barrier effect. IP3 reaches the right branch ER
and releases calcium. Calcium moves by diffusion in both directions. Movement
in retrograde direction (right branch to cell) obstructs the front advancement
of calcium originating in central ER. This effect is observed at approximately
0.2s and does not exist in configurations represented by blue and red lines.
The effect of left branch ER is highlighted in black and red lines. Calcium
flux reaches a smaller total amount and ends earlier. The reason is that left
branch ER degrades IP3, reducing the amount that reaches central ER. Finally,
blue and green lines differ only by calcium barrier effect produced by right
branch ER.

To verify if calcium barrier really causes the flow decrease in the connection
between left branch and cell body, the L1R0 simulation was performed, but the
ER only degraded IP3 and did not release calcium. This simulation was called
L1R0*. As can be seen in Fig. 24, the difference between the blue (L0R0) and
pink (L1R0*) curves shows that IP3 degradation actually decreases this flow.
By contrast, the difference between pink (L1R0*) and red (L1R0) shows the
calcium barrier effect on this flow.

Figure 24: Flux through connection between left branch and cell body compared
to the situation where the left branch ER does not release calcium. The red
and pink lines show the flow of calcium in situations in which the left branch
ER respectively releases or not calcium. It is possible to observe the calcium
barrier effect. In situation L1R0 the flow is smaller. Red and blue lines
difference reveals the effect in calcium flow by left branch ER IP3
degradation. In left branch ER absence more IP3 reaches the cell body ER
causing greater calcium release.

In order to evaluate if calcium barrier promotes calcium wave advancement, one
can observe Fig. 25, which compares calcium fluxes in the connection between
right branch and cellular body. It is possible to observe no distinction
between L1R0 and L1R0* configurations. This shows the calcium barrier effects
are dissipated along the calcium pathway. In this way, the barrier contributes
to local, but not global, dynamics.

Figure 25: Flux through connection between right branch and cell body compared
to the situation where the left branch reticulum does not release calcium. The
graph shows that while there may be a calcium barrier formed by the left
branch ER calcium release, its effects are dissipated along the calcium
pathway. Therefore, in the connection between right branch and cell body,
there is no distinction between pink and red lines. At this point, the only
effect observed on calcium flow is IP3 degradation promoted by the left branch
ER.

Based on these results, it is possible to notice that several processes
influence calcium signaling but there is a clear distinction between local and
global effects (Di Castro et al., 2011), as well as between dynamic and
quantitative effects. ER on branches participates primarily in the local
effects and contributes mainly to the local process dynamics. The ER, when
disposed in vicinity of where the stimuli occur, contribute to local calcium
signaling. This may be interesting for the cell to speed up the release of
vesicles, which depends on increased intracellular calcium concentration
(Araque et al., 2014).

Figure 26: Instantaneous calcium flux at right branch outlet the solid lines
represent the L1R1 configurations and the dashed lines represent the L1R0
configurations. The times for the maximum points of each curve, with
dimensions of the right branch of L=6$\mu m$, 2L and 3L, are, respectively and
in seconds, 0.530, 0.560 and 0.700 for configurations L1R1 and 1.890, 2.790
and 3.840 for configurations L1R0.

ER allocated in the branches alter the calcium signaling dynamics, as they
allow the calcium concentration to increase more rapidly at the end of the
distant branches from where the stimulus occurred. This is because IP3 moves
very quickly, reaching ERs located in distant branches long before the calcium
wave produced by the ER located in the cell body reaches them. In order to
better visualize this effect, a new structure was created, similar to that
described in Fig. 18, a central cell and two branches of length L. Simulations
were performed with three sets of two configurations L1R0 and L1R1. The length
of the right branch was measured L, 2L and 3L. It was possible to notice an
increase in the time difference between the calcium wave at the right end, in
situations in which this branch was with or without ER. The results are shown
in Fig. 26. Here, the solid lines represent the L1R1 configurations and the
dashed lines represent the L1R0 configurations. The times for the maximum
points of each curve, with dimensions of the right branch of L, 2L and 3L,
are, respectively and in seconds, 0.530, 0.560 and 0.700 for configurations
L1R1 and 1.890, 2.790 and 3.840 for configurations L1R0. The differences in
times, between configurations L1R1 and L1R0, for dimensions L, 2L and 3L are,
respectively, 1.360, 2.230 and 3.140. These points have a linear behavior and
are represented by the line: $y=(5+10Lx)D_{C}/D_{I}$. On the other hand, Fig.
23 makes it clear that, in the configuration used in this work, the final
value of the accumulated calcium flow, at the exit of the right branch, does
not change with the existence of ER in this branch, just contributes to
accelerate calcium wave dynamics. These results suggest that we are not
dealing with calcium waves, but with IP3 waves with occasional geyser bursts
of calcium.

##  5 Conclusions

The present work presents three contributions to the scientific community: a
mathematical model and its computational implementation, a parametric
sensitivity analysis and aan analysis of the effects of ER on astrocytes
branches.

Regarding the first contribution, the model was computationally constructed
and implemented, to guarantee analytical flexibility of calcium signaling in
astrocytes. Astrocyte geometry can be quite extensive and branched, depending
only on available computational resources. In addition, the mechanisms can be
easily turned on or off, allowing an analysis of their effects on specific
studied phenomenon. The mathematical model also considers parameters that can
be experimentally determined, including channel density and mass transfer
coefficients. All these aspects are important when comparing model against
experimental results; this is because there is a huge diversity of
experimental results, as calcium signaling depends on location and timing of
calcium fluxes (Volterra et al., 2014). This type of information is easily fed
into the model, allowing a more accurate comparison.

It is important to point out that mathematical models should be made available
as the scientific community improves them. Manninen et al. (2017) reported
that several models are not reproducible and do not provide all the necessary
information. To avoid this, the authors of the present work made available the
program code so that all those who work in this area can use it, modify it,
and generate the most varied results, contributing to the discussion of the
role of astrocytes in brain information processing.

Regarding the parametric sensitivity analysis, more than providing a
reflection about the influence of parameters and mechanisms on calcium
signaling, it allowed two very pertinent conclusions. The first one is that
several sets of parameters lead to the same calcium dynamics in astrocytes.
This behavior is quite common in biological systems, because the system as a
whole must respond in a similar way in different types of situations to which
it is subjected (Prinz et al., 2004). The second is that parameters of
mathematical models cannot be regarded as absolute and require constant
revision, especially when they do not come from experimental results. Several
parameters used in the present work and in a sequence of literature available
works were initially adjusted by Keizer and De Young (1992). Several of them
were obtained from optimization of cell behavior as a whole, that is,
depending on the type of experiment that was performed. In addition, the
authors themselves emphasized that they chose only one of the possible sets
for the assessed situation. The wide use of these parameters has led the
scientific community to believe that there is no need to revise their values.
According to Manninen et al. (2018), the models for signaling calcium in
astrocytes are similar and use basically the same parameters. These values are
not discussed and they are used indiscriminately. The proposition here is,
given the experimental development that currently exists, specific experiments
should be made for each mechanism. These new parameters can be easily
aggregated into the our model, allowing for more accurate results.

Finally, the present work contributed to the analysis of ER effects on
astrocytes branches. The location and function of ER in astrocytes still
creates divergences of opinion in the scientific community (Golovina and
Blaustein, 2000; Golovina, 2005; Holtzclaw et al., 2002; Oheim et al., 2018).
Our results indicate that the ER has positive and negative effects when
allocated in the branches, depending of the branches position in the cell.
When placed next to the stimulus, they contribute to local calcium signaling;
however, they affect negatively the global signaling due to the IP3
degradation they promote. It is important to emphasize that the damage to
global calcium signaling caused by the ER located next to the stimulus is not
sufficient to prevent the signaling from advancing. By contrast, when
allocated far from the stimulus site, they allow calcium signaling to quickly
reach the most distant points. Thus, given brain processes speed, it is
possible that astrocyte allocate ER in all branches, to accelerate the advance
of calcium signaling. The consequence of the local signaling favored by the ER
at branches extremities where the stimulus occurs still has to be better
explored in future works.

It is important to emphasize that models for calcium signaling in astrocytes
require substantial sophistication in their contribution to the understanding
of astrocyte functions in cerebral information processing. Collaborative work
allows us to achieve greater goals in less time. Therefore, by making our
software available for the community, we welcome experimental and theoretical
scientists to contribute to the development of this model and to the
estimation of more appropriate parameters so all of us can make use of this
versatile and useful tool.

###### Acknowledgements.

The authors thank the Brazilian Federal Agencies, Coordination for Improvement
of Higher Education Personnel, CAPES, and National Council for Scientific and
Technological Development, CNPq, by scholarships for the first of them.

## References

  * Alberts et al. (2017) Alberts B, Johnson A, Lewis J, Morgan D, Raff M, Roberts K, Walter P (2017) Molecular Biology of the Cell. Garland Science
  * Allam et al. (2012) Allam SL, Ghaderi VS, Bouteiller JMC, Legendre A, Ambert N, Greget R, Bischoff S, Baudry M, Berger TW (2012) A computational model to investigate astrocytic glutamate uptake influence on synaptic transmission and neuronal spiking. Frontiers in Computational Neuroscience 6(70)
  * Amiri et al. (2011) Amiri M, Bahrami F, Janahmadi M (2011) Functional modeling of astrocytes in epilepsy: a feedback system perspective. Neural Computing and Applications 20(8):1131–1139
  * Araque et al. (2014) Araque A, Carmignoto G, Haydon PG, Oliet SH, Robitaille R, Volterra A (2014) Gliotransmitters travel in time and space. Neuron 81(4):728–739
  * Bazargani and Attwell (2016) Bazargani N, Attwell D (2016) Astrocyte calcium signaling: the third wave. Nature Neuroscience 19(2):182–189
  * Bear et al. (2007) Bear MF, Connors BW, Paradiso MA (2007) Neuroscience, vol 2. Lippincott Williams & Wilkins
  * Bernardinelli et al. (2014) Bernardinelli Y, Muller D, Nikonenko I (2014) Astrocyte-synapse structural plasticity. Neural Plasticity 2014:ID 232105
  * Berridge et al. (2003) Berridge MJ, Bootman MD, Roderick HL (2003) Calcium signalling: dynamics, homeostasis and remodelling. Nature Reviews Molecular Cell Biology 4(7):517–529
  * Bindocci et al. (2017) Bindocci E, Savtchouk I, Liaudet N, Becker D, Carriero G, Volterra A (2017) Three-dimensional Ca2+ imaging advances understanding of astrocyte biology. Science 356(6339):eaai8185
  * Bootman et al. (2002) Bootman MD, Berridge MJ, Roderick HL (2002) Calcium signalling: more messengers, more channels, more complexity. Current Biology 12(16):R563–R565
  * Budisantoso et al. (2013) Budisantoso T, Harada H, Kamasawa N, Fukazawa Y, Shigemoto R, Matsui K (2013) Evaluation of glutamate concentration transient in the synaptic cleft of the rat calyx of held. The Journal of Physiology 591(1):219–239
  * Clapham (1995) Clapham DE (1995) Calcium signaling. Cell 80(2):259–268
  * Di Castro et al. (2011) Di Castro MA, Chuquet J, Liaudet N, Bhaukaurally K, Santello M, Bouvier D, Tiret P, Volterra A (2011) Local Ca2+ detection and modulation of synaptic release by astrocytes. Nature Neuroscience 14(10):1276–84
  * Di Garbo et al. (2007) Di Garbo A, Barbi M, Chillemi S, Alloisio S, Nobile M (2007) Calcium signalling in astrocytes and modulation of neural activity. Biosystems 89(1):74–83
  * Durkee and Araque (2019) Durkee CA, Araque A (2019) Diversity and specificity of astrocyte-neuron communication. Neuroscience 396(1):73–78
  * Eisenberg (1999) Eisenberg R (1999) From structure to function in open ionic channels. Journal of Membrane Biology 171(1):1–24
  * Falcke (2003) Falcke M (2003) On the role of stochastic channel behavior in intracellular Ca2+ dynamics. Biophysical Journal 84(1):42–56
  * Fiacco and McCarthy (2018) Fiacco TA, McCarthy KD (2018) Multiple lines of evidence indicate that gliotransmission does not occur under physiological conditions, dual perspectives. The Journal of Neuroscience 38(1):3–13
  * Golovina (2005) Golovina VA (2005) Visualization of localized store-operated calcium entry in mouse astrocytes. close proximity to the endoplasmic reticulum. The Journal of Physiology 564(3):737–749
  * Golovina and Blaustein (2000) Golovina VA, Blaustein MP (2000) Unloading and refilling of two classes of spatially resolved endoplasmic reticulum Ca2+ stores in astrocytes. Glia 31(1):15–28
  * Golovina et al. (2003) Golovina VA, Song H, James PF, Lingrel JB, Blaustein MP (2003) Na+ pump $\alpha$ 2-subunit expression modulates Ca2+ signaling. American Journal of Physiology-Cell Physiology 284(2):C475–C486
  * Halassa et al. (2007a) Halassa MM, Fellin T, Haydon PG (2007a) The tripartite synapse: roles for gliotransmission in health and disease. Trends in Molecular Medicine 13(2):54–63
  * Halassa et al. (2007b) Halassa MM, Fellin T, Takano H, Dong JH, Haydon PG (2007b) Synaptic islands defined by the territory of a single astrocyte. Journal of Neuroscience 27(24):6473–6477
  * Höfer et al. (2002) Höfer T, Venance L, Giaume C (2002) Control and plasticity of intercellular calcium waves in astrocytes: a modeling approach. Journal of Neuroscience 22(12):4850–4859
  * Holtzclaw et al. (2002) Holtzclaw LA, Pandhit S, Bare DJ, Mignery GA, Russell JT (2002) Astrocytes in adult rat brain express type 2 inositol 1, 4, 5-trisphosphate receptors. Glia 39(1):69–84
  * Howarth (2014) Howarth C (2014) The contribution of astrocytes to the regulation of cerebral blood flow. Frontiers in Neuroscience 8(103)
  * Kanemaru et al. (2014) Kanemaru K, Sekiya H, Xu M, Satoh K, Kitajima N, Yoshida K, Okubo Y, Sasaki T, Moritoh S, Hasuwa H, et al. (2014) In vivo visualization of subtle, transient, and local activity of astrocytes using an ultrasensitive Ca2+ indicator. Cell Reports 8(1):311–318
  * Keizer and De Young (1992) Keizer J, De Young GW (1992) Two roles of Ca2+ in agonist stimulated Ca2+ oscillations. Biophysical Journal 61(3):649–660
  * Koh et al. (1995) Koh DS, Geiger J, Jonas P, Sakmann B (1995) Ca2+-permeable ampa and nmda receptor channels in basket cells of rat hippocampal dentate gyrus. The Journal of Physiology 485(2):383–402
  * Lallouette et al. (2018) Lallouette J, De Pittà M, Berry H (2018) Astrocyte networks and intercellular calcium propagation. bioRxiv pp 1–33
  * Li and Rinzel (1994) Li YX, Rinzel J (1994) Equations for InsP3 receptor-mediated [Ca2+]i oscillations derived from a detailed kinetic model: a hodgkin-huxley like formalism. Journal of Theoretical Biology 166(4):461–473
  * Losi et al. (2017) Losi G, Mariotti L, Sessolo M, Carmignoto G (2017) New tools to study astrocyte Ca2+ signal dynamics in brain networks in vivo. Frontiers in Cellular Neuroscience 11(134)
  * Luby-Phelps (1999) Luby-Phelps K (1999) Cytoarchitecture and physical properties of cytoplasm: volume, viscosity, diffusion, intracellular surface area. International Review of Cytology 192:189–221
  * MacVicar (1984) MacVicar BA (1984) Voltage-dependent calcium channels in glial cells. Science 226:1345–1348
  * Malli et al. (2007) Malli R, Frieden M, Hunkova M, Trenker M, Graier W (2007) Ca2+ refilling of the endoplasmic reticulum is largely preserved albeit reduced Ca2+ entry in endothelial cells. Cell Calcium 41(1):63–76
  * Manninen et al. (2017) Manninen T, Havela R, Linne ML (2017) Reproducibility and comparability of computational models for astrocyte calcium excitability. Frontiers in Neuroinformatics 11(11):1–18
  * Manninen et al. (2018) Manninen T, Havela R, Linne ML (2018) Computational models for calcium-mediated astrocyte functions. Frontiers in Computational Neuroscience 12(14):1–27
  * Matyash and Kettenmann (2010) Matyash V, Kettenmann H (2010) Heterogeneity in astrocyte morphology and physiology. Brain Research Reviews 63(1):2–10
  * McCarthy and Salm (1991) McCarthy KD, Salm AK (1991) Pharmacologically-distinct subsets of astroglia can be identified by their calcium response to neuroligands. Neuroscience 41(2):325–333
  * McDougal et al. (2017) McDougal RA, Morse TM, Carnevale T, Marenco L, Wang R, Migliore M, Miller PL, Shepherd GM, Hines ML (2017) Twenty years of modeldb and beyond: building essential modeling tools for the future of neuroscience. Journal of Computational Neuroscience 42(1):1–10
  * Means et al. (2006) Means S, Smith AJ, Shepherd J, Shadid J, Fowler J, Wojcikiewicz RJ, Mazel T, Smith GD, Wilson BS (2006) Reaction diffusion modeling of calcium dynamics with realistic er geometry. Biophysical Journal 91(2):537–557
  * Min et al. (2012) Min R, Santello M, Nevian T (2012) The computational power of astrocyte mediated synaptic plasticity. Frontiers in Computational Neuroscience 6(93)
  * Myung and Pitt (2002) Myung IJ, Pitt MA (2002) Mathematical modeling. Stevens’ Handbook of Experimental Psychology IV(Methodology):429–459
  * Nadkarni and Jung (2005) Nadkarni S, Jung P (2005) Synaptic inhibition and pathologic hyperexcitability through enhanced neuron-astrocyte interaction: a modeling study. Journal of Integrative Neuroscience 4(02):207–226
  * Navarrete et al. (2012) Navarrete M, Perea G, Maglio L, Pastor J, García de Sola R, Araque A (2012) Astrocyte calcium signal and gliotransmission in human brain tissue. Cerebral Cortex 23(5):1240–1246
  * Neary et al. (1988) Neary J, Van Breemen C, Forster E, Norenberg L, Norenberg M (1988) Atp stimulates calcium influx in primary astrocyte cultures. Biochemical and Biophysical Research Communications 157(3):1410–1416
  * O’Hayre et al. (2013) O’Hayre M, Vázquez-Prado J, Kufareva I, Stawiski EW, Handel TM, Seshagiri S, Gutkind JS (2013) The emerging mutational landscape of g-proteins and g-protein coupled receptors in cancer. Nature Reviews Cancer 13(6):412–424
  * Oheim et al. (2018) Oheim M, Schmidt E, Hirrlinger J (2018) Local energy on demand: Are ‘spontaneous’ astrocytic Ca2+-microdomains the regulatory unit for astrocyte-neuron metabolic cooperation? Brain Research Bulletin 136:54–64
  * Oschmann et al. (2018) Oschmann F, Berry H, Obermayer K, Lenk K (2018) From in silico astrocyte cell models to neuron-astrocyte network models: a review. Brain Research Bulletin 136:76–84
  * Patrushev et al. (2013) Patrushev I, Gavrilov N, Turlapov V, Semyanov A (2013) Subcellular location of astrocytic calcium stores favors extrasynaptic neuron–astrocyte communication. Cell Calcium 54(5):343–349
  * Pawelczyk and Matecki (1997) Pawelczyk T, Matecki A (1997) Structural requirements of phospholipase C$\delta$1 for regulation by spermine, sphingosine and sphingomyelin. The FEBS Journal 248(2):459–465
  * Perea et al. (2009) Perea G, Navarrete M, Araque A (2009) Tripartite synapses: astrocytes process and control synaptic information. Trends in Neuroscience 32(8):421–431
  * Perez-Alvarez et al. (2014) Perez-Alvarez A, Navarrete M, Covelo A, Martin ED, Araque A (2014) Structural and functional plasticity of astrocyte processes and dendritic spine interactions. Journal of Neuroscience 34(38):12738–12744
  * Pirttimaki et al. (2017) Pirttimaki TM, Sims RE, Saunders G, Antonio SA, Codadu NK, Parri HR (2017) Astrocyte-mediated neuronal synchronization properties revealed by false gliotransmitter release. Journal of Neuroscience 37(41):9859–9870
  * Pivneva et al. (2008) Pivneva T, Haas B, Reyes-Haro D, Laube G, Veh RW, Nolte C, Skibo G, Kettenmann H (2008) Store-operated Ca2+ entry in astrocytes: different spatial arrangement of endoplasmic reticulum explains functional diversity in vitro and in situ. Cell Calcium 43(6):591–601
  * Poskanzer and Molofsky (2018) Poskanzer KE, Molofsky AV (2018) Dynamism of an astrocyte in vivo: perspectives on identity and function. Annual Review of Physiology 80:143–157
  * Prinz et al. (2004) Prinz AA, Bucher D, Marder E (2004) Similar network activity from disparate circuit parameters. Nature Neuroscience 7(12):1345–1352
  * Rusakov (2015) Rusakov DA (2015) Disentangling calcium-driven astrocyte physiology. Nature Reviews Neuroscience 16(4):226–233
  * Sasaki et al. (2014) Sasaki T, Ishikawa T, Abe R, Nakayama R, Asada A, Matsuki N, Ikegaya Y (2014) Astrocyte calcium signalling orchestrates neuronal synchronization in organotypic hippocampal slices. The Journal of Physiology 592(13):2771–2783
  * Savtchouk and Volterra (2018) Savtchouk I, Volterra A (2018) Gliotransmission: Beyond black-and-white, dual perspectives. The Journal of Neuroscience 38(1):14–25
  * Shigetomi et al. (2016) Shigetomi E, Patel S, Khakh BS (2016) Probing the complexities of astrocyte calcium signaling. Trends in Cell Biology 26(4):300–312
  * Shtrahman et al. (2017) Shtrahman E, Maruyama D, Olariu E, Fink C, Zochowski M (2017) Understanding spatial and temporal patterning of astrocyte calcium transients via interactions between network transport and extracellular diffusion. Physical Biology 14(1):016001
  * Simpson et al. (1998) Simpson PB, Holtzclaw LA, Langley DB, Russell JT (1998) Characterization of ryanodine receptors in oligodendrocytes, type 2 astrocytes, and o-2a progenitors. Journal of Neuroscience Research 52(4):468–482
  * Skupin et al. (2008) Skupin A, Kettenmann H, Winkler U, Wartenberg M, Sauer H, Tovey SC, Taylor CW, Falcke M (2008) How does intracellular Ca2+ oscillate: by chance or by the clock? Biophysical Journal 94(6):2404–2411
  * Suffczynski et al. (2006) Suffczynski P, Wendling F, Bellanger JJ, Da Silva FL (2006) Some insights into computational models of (patho) physiological brain activity. Proceedings of the IEEE 94(4):784–804
  * Taheri et al. (2017) Taheri M, Handy G, Borisyuk A, White JA (2017) Diversity of evoked astrocyte Ca2+ dynamics quantified through experimental measurements and mathematical modeling. Frontiers in Systems Neuroscience 11(79):1–19
  * Thore et al. (2005) Thore S, Dyachok O, Gylfe E, Tengholm A (2005) Feedback activation of phospholipase c via intracellular mobilization and store-operated influx of Ca2+ in insulin-secreting $\beta$-cells. Journal of Cell Science 118(19):4463–4471
  * Traynelis et al. (2010) Traynelis SF, Wollmuth LP, McBain CJ, Menniti FS, Vance KM, Ogden KK, Hansen KB, Yuan H, Myers SJ, Dingledine R (2010) Glutamate receptor ion channels: structure, regulation, and function. Pharmacological Reviews 62(3):405–496
  * Tu et al. (2005a) Tu H, Wang Z, Bezprozvanny I (2005a) Modulation of mammalian inositol 1, 4, 5-trisphosphate receptor isoforms by calcium: a role of calcium sensor region. Biophysical Journal 88(2):1056–1069
  * Tu et al. (2005b) Tu H, Wang Z, Nosyreva E, De Smedt H, Bezprozvanny I (2005b) Functional characterization of mammalian inositol 1, 4, 5-trisphosphate receptor isoforms. Biophysical Journal 88(2):1046–1055
  * Ullah et al. (2006) Ullah G, Jung P, Cornell-Bell AH (2006) Anti-phase calcium oscillations in astrocytes via inositol (1, 4, 5)-trisphosphate regeneration. Cell Calcium 39(3):197–208
  * Venance et al. (1997) Venance L, Stella N, Glowinski J, Giaume C (1997) Mechanism involved in initiation and propagation of receptor-induced intercellular calcium signaling in cultured rat astrocytes. Journal of Neuroscience 17(6):1981–1992
  * Verkhratsky and Parpura (2014) Verkhratsky A, Parpura V (2014) Store-operated calcium entry in neuroglia. Neuroscience Bulletin 30(1):125–133
  * Verkhratsky et al. (2012) Verkhratsky A, Rodríguez JJ, Parpura V (2012) Calcium signalling in astroglia. Molecular and Cellular Endocrinology 353(1):45–56
  * Volterra et al. (2014) Volterra A, Liaudet N, Savtchouk I (2014) Astrocyte Ca2+ signalling: an unexpected complexity. Nature Reviews Neuroscience 15(5):327–335
  * Wu et al. (2014) Wu YW, Tang X, Arizono M, Bannai H, Shih PY, Dembitskaya Y, Kazantsev V, Tanaka M, Itohara S, Mikoshiba K, et al. (2014) Spatiotemporal calcium dynamics in single astrocytes and its modulation by neuronal activity. Cell Calcium 55(2):119–129