Astroglial cells express all major ion channel types, including potassium (K⁺), sodium (Na⁺), and Ca²⁺ channels, and also various types of anion and chloride (Cl⁻) channels, water channels (aquaporins), transient receptor potential (TRP) channels, and non-selective channels. The ion distribution is also somewhat different from neurons: intracellular concentrations of K⁺ and Ca²⁺ are similar to neurons, but the concentrations of Na⁺ and especially Cl⁻ are higher compared to neurons. Astrocytes have a rather negative resting membrane potential (around −80 to −90 mV) because of the predominance of K⁺ conductance. Electrical depolarization of astroglia does not produce regenerative action potentials as in neurons.

Ca²⁺-mediated signals have been proposed to be the main mediator of communication between astrocytes and other cellular elements in the brain (Nimmerjahn, 2009; Volterra et al., 2014; Bazargani and Attwell, 2016). Transient Ca²⁺ increases restricted to single cells are called Ca²⁺ oscillations. In isolated astrocytes, intracellular Ca²⁺ oscillations have been shown to depend mainly on the Ca²⁺-induced Ca²⁺ release (CICR) via IP₃Rs from the ER to the cytosol (see e.g., Agulhon et al., 2008). CICR has been shown to depend on Ca²⁺ with or without the influence of IP₃. Ca²⁺ influx via voltage-gated Ca²⁺ channels (VGCCs) and other types of Ca²⁺ influxes from the extracellular space have been linked with Ca²⁺ oscillations as well (Aguado et al., 2002).

The membrane transporters are particularly important for astroglia because they control movements of various substances, including ions, neurotransmitters, and metabolic substrates. Astroglial transporters include adenosine and adenosine triphosphate (ATP)-dependent transporters such as the Na⁺/K⁺-ATPase (NKA, also called Na⁺/K⁺ pump) and Ca²⁺-ATPase [also called Ca²⁺ pump or plasma membrane Ca²⁺-ATPase (PMCA)] on the plasma membrane, in addition to sarco/ER Ca²⁺-ATPase (SERCA) located on the ER membrane. They also contain so-called secondary transporters, such as glutamate transporters [excitatory amino acid transporters (EAATs)], gamma-aminobutyric acid (GABA) transporters, glycine transporters, Na⁺/Ca²⁺ exchangers (NCXs), Na⁺/hydrogen (H⁺) exchangers, Na⁺/bicarbonate (HCO3-) cotransporters, Na⁺/K⁺/Cl⁻ cotransporters (NKCC1), and some others. Although, for example, glutamate transporters are expressed by all cell types in the brain, astrocytes are the main cell type responsible for glutamate uptake. Astrocytes have enzymes that convert both glutamate and GABA into glutamine. Glutamine is then released into the extracellular space and taken up by the presynaptic terminal, and can be converted to glutamate or GABA.

The NKCC1 cotransporter specifically contributes to the regulation of extracellular K⁺ homeostasis in the central nervous system. During excessive neuronal firing, the local extracellular K⁺ concentration can increase markedly and lead to neural dysfunction. Uptaking extracellular K⁺, transporting it intracellularly and releasing to areas with lower concentration are some of the most important functions of astrocytes. This process of intracellularly transporting K⁺ from extracellular areas of high concentration to areas of low concentration, via active uptake and release, is called spatial K⁺ buffering. K⁺ can also diffuse in extracellular space which also helps in lowering the local high concentration. Furthermore, astrocytes also buffer other ions than K⁺.

Astrocytes are able to sense with the help of their receptors local and more distant environments, including the neural activity and the synaptically released neurotransmitters, both in in vitro cell cultures and brain slices as well as in in vivo (Glaum et al., 1990; Dani et al., 1992; Porter and McCarthy, 1996; Hirase et al., 2004). Astroglial cells can express the same variety of receptors as neurons, both ionotropic and metabotropic. These include glutamate, GABA, glycine, acetylcholine, adrenergic, serotonin, histamine, cannabinoid, and neuropeptide receptors, and purinoceptors for adenosine and ATP. For example, extracellular ATP can bind to astroglial purinergic G-protein receptors, resulting in IP₃-mediated CICR, and, eventually, release of ATP to the extracellular space.

As described above, astroglia are capable of receiving signals from neurons via membrane receptors and converting the received information into Ca²⁺ excitability. The fact that astrocytes could release extracellular signaling molecules, regulated by this Ca²⁺ excitability, not only to the vascular system but also to the neuronal one implies a very active role for astroglia in the brain. The concept of regulated transmitter/modulator release from astrocytes to neurons is generally known as gliotransmission (Parpura et al., 1994; Bezzi and Volterra, 2001), and the released substances are referred to as gliotransmitters. Most common gliotransmitters are glutamate, D-serine, and ATP (Parpura et al., 1994; Parri et al., 2001). Astrocytes are, in general, thought to release transmitters, and probably also modulators, by several different mechanisms, which can be broadly divided into exocytotic release, diffusional release through membrane pores, and transporter mediated release. It is, however, not known which of these mechanisms are used in certain astroglial functions.

These terms “gliotransmission” and “gliotransmitter” may be somewhat misleading as the very same transmitters are also released by neurons. Moreover, it is assumed that a mechanism similar to vesicular release from neuronal synaptic terminals exists also in astrocytes. Some studies have indeed detected vesicle-type structures in astrocytes in vitro, and thus it has been proposed that gliotransmitter release might be purely vesicular. It is, however, important to bear in mind that the existence of molecular machinery for packing gliotransmitters into vesicles and for arranging the vesicles into a readily releasable pool has not so far been supported by experimental evidence in vivo (see e.g., Fujita et al., 2014; Sloan and Barres, 2014). More evidence on the release mechanism, using improved experimental model systems and techniques that allow studies at deeper resolution in physiological conditions, is required (Li et al., 2013; Bazargani and Attwell, 2016; Fiacco and McCarthy, 2018; Savtchouk and Volterra, 2018).

In our evaluation of models, we use the term “gliotransmission” for all biophysical and phenomenological mechanisms that were modeled to take into account the release of substances from astrocytes and targeting neurons. The reason for this is that the term “gliotransmission” is often used in the original modeling publications.

In addition, glutamate released from astrocytes can activate extrasynaptic N-methyl-D-aspartate receptor (NMDAR)-dependent currents, often called NMDAR-dependent slow inward current (SIC). In modeling studies, SIC is many times modeled similarly to, for example, the modulating current (I_astro) presented by Nadkarni and Jung (2003).

It is not just neurons that form networks but also astrocytes. A fundamental difference between neuronal and astroglial networks is that astrocytes are connected, through gap junctions composed mainly of connexin 43 hemichannels, to form a functional cellular syncytium in the central nervous system. In their open state, these channels are permeable to large hydrophilic solutes with molecular mass of several hundred Daltons, and are permeable to small solutes in their closed state (see e.g., Bao et al., 2007). The gap junction connectivity is instrumental for astrocytes' functions, including generation of Ca²⁺ waves, water transport, K⁺ buffering, and control of vascular system, and is one of the most important mechanisms to be modeled in astrocyte networks.

Three different pathways have been discovered so far to induce Ca²⁺ waves in astroglial networks. The first route depends on the transfer of IP₃ via gap junctions (Giaume and Venance, 1998). Transported IP₃ via gap junctions triggers CICR in the coupled astrocytes and induces Ca²⁺ wave propagation in astroglial syncytium. The second route to induce Ca²⁺ waves depends on the extracellular diffusion of ATP (see e.g., Newman and Zahs, 1997; Guthrie et al., 1999 and section 2.1.4). The third route has been shown to depend on extracellularly applied potassium chloride, causing, among others, a pathophysiological phenomenon called cortical spreading depression (Peters et al., 2003). The regulation of gap junction communication within the astroglial syncytium is a complex process and is intensively studied.

Most of the above described biophysical and biochemical mechanisms have been modeled in some detail in astrocytes. Below we address altogether 106 models developed until the end of 2017 and describe their capacity to represent the dynamics of astrocyte biophysics and biochemistry.

Half of the single astrocyte models were so-called generic, meaning that they did not describe astrocytes in any specific anatomical brain location. Others, however, were specified to model astrocytes in the cerebrum (Farr and David, 2011; Witthoft and Karniadakis, 2012), cerebral cortex (Diekman et al., 2013; Witthoft et al., 2013; Mesiti et al., 2015b; Kenny et al., 2018), cortex (De Pittà et al., 2009b; Toivari et al., 2011), hippocampus (Riera et al., 2011a,b; Chander and Chakravarthy, 2012), as well as the visual cortex (Gibson et al., 2007; Bennett et al., 2008b) and somatosensory cortex (Bennett et al., 2008b; Taheri et al., 2017). One third of the single astrocyte models took into account neurotransmitters in a simplistic way just as a stimulus, having either the neurotransmitter as a constant, step function, or something similar (see e.g., Larter and Craig, 2005; Gibson et al., 2007; Bennett et al., 2008b; De Pittà et al., 2009a; Dupont et al., 2011; Toivari et al., 2011; Witthoft and Karniadakis, 2012; Hadfield et al., 2013; Witthoft et al., 2013; Kenny et al., 2018). Only two models (Chander and Chakravarthy, 2012; Oschmann et al., 2017) actually modeled the amount of neurotransmitter with a differential equation. The stimulus to the astrocyte model by Oschmann et al. (2017) was taken from the model by Tsodyks and Markram (1997). In addition, Mesiti et al. (2015b) modeled the presynaptic neuron. We included these three models (Chander and Chakravarthy, 2012; Mesiti et al., 2015b; Oschmann et al., 2017) under single astrocyte models, because these models did not have bidirectional communication between astrocytes and neurons. The characteristics of single astrocyte models can be found in Table 2.

Most of the single astrocyte models studied Ca²⁺ oscillations, of which a few models specifically focused on modeling only spontaneous Ca²⁺ oscillations (see Table 2). All the other models had components for CICR and SERCA pump except the model by Montaseri and Yazdanpanah (2014). Furthermore, all the other models except the models by López-Caamal et al. (2014) and Montaseri and Yazdanpanah (2014) modeled leak from the ER into the cytosol. Half of the models had influx of Ca²⁺ from outside of the astrocyte or efflux of Ca²⁺ to outside of the astrocyte. About one third of the models took into account Ca²⁺ buffers and astrocytic release of signaling molecules. None of the models had gap junctions, because these were single astrocyte models. Thus, these models had similar core structure with small variations. As an example, six modeled capacitive Ca²⁺ entry (CCE) which is mediated via store-operated Ca²⁺ (SOC) channels (Riera et al., 2011a,b; Toivari et al., 2011; Handy et al., 2017; Taheri et al., 2017; Ding et al., 2018) and two modeled VGCCs (Zeng et al., 2009; Ding et al., 2018), as well as one had Ca²⁺ flux via ryanodine receptors (RyRs) (López-Caamal et al., 2014) and three had Ca²⁺ influx via TRP vanilloid-related channel 4 (TRPV4) (Witthoft and Karniadakis, 2012; Witthoft et al., 2013; Kenny et al., 2018). Only Diekman et al. (2013) and Komin et al. (2015) modeled mitochondrial Ca²⁺ unitransporters (MCUs). Two modeled mitochondrial NCX (Diekman et al., 2013; Komin et al., 2015) and one modeled NCX via plasma membrane (Oschmann et al., 2017). Skupin et al. (2010), Komin et al. (2015), and Ding et al. (2018) used detailed IP₃R model by De Young and Keizer (1992). Six models had diffusion either in the cytosol, in the ER, or in the extracellular space (Roth et al., 1995; Gibson et al., 2007; Bennett et al., 2008b; Skupin et al., 2010; López-Caamal et al., 2014; Komin et al., 2015).

The first models developed in this category were the models by Roth et al. (1995) and Larter and Craig (2005). Roth et al. (1995) extended the model by Li and Rinzel (1994) by adding Ca²⁺ diffusion in the cytosol and ER. They also divided their cell to alternating active and passive compartments. In passive compartments, Ca²⁺ signal propagated by diffusion, whereas in active compartments there were also additional Ca²⁺ fluxes propagating the signal. Roth et al. (1995) showed with their model that different parts of the astrocyte, such as the ER membrane, produced different frequencies of Ca²⁺ oscillations if the diffusion constant had a value close to the physiological range. If the diffusion constant had unphysiologically high value, then the different parts of the astrocyte all oscillated with the same frequency.

Larter and Craig (2005) presented the only model in this category that included a differential equation for extracellular glutamate ([Glu]_ext). Their equation for extracellular glutamate included the rate of glutamate vesicle release from the astrocyte, the input glutamate rate arriving from nearly synapses, and two terms describing the rates of glutamate uptake by neurons and astrocytes. Larter and Craig (2005) tested their model by setting the extracellular glutamate to zero during the simulations (meaning both glutamate release from the astrocyte and glutamate arriving from the nearby synapses were zero) and showed that the astrocytic Ca²⁺ was not oscillating in this case. They further tested their model when the astrocyte was not releasing glutamate, but they still had the input glutamate rate arriving from nearly synapses. In this case, the astrocytic Ca²⁺ was oscillating. They also tested how the number of glutamate in the astrocytic vesicles affected the model behavior. They showed that when increasing the number of glutamate in the astrocytic vesicles, the Ca²⁺ oscillations transformed into Ca²⁺ bursting oscillations. Thus, their message was that the presence of glutamate, but not the positive astrocytic glutamate feedback, was required for astrocytic Ca²⁺ oscillations. In addition, the positive glutamate feedback process made it possible for the astrocyte to have bursting oscillations when the input glutamate rate arriving from nearby synapses had low values.

Only six of the single astrocyte models had K⁺ concentration modeled in the perisynaptic space, perivascular space [in Table 2 called as extracellular space (ext)], or intracellular space of the astrocyte. These are marked in Table 2 as [K⁺]_N, [K⁺]_ext, and [K⁺], respectively. The model by Diekman et al. (2013) had the most simple K⁺ model, they modeled only astrocytic K⁺ that depended on NKA. Oschmann et al. (2017) modeled both the astrocytic and extracellular K⁺ concentrations. The astrocytic K⁺ concentration depended on three transmembrane transporters; glutamate transporter, NKA, and NCX. The extracellular K⁺ concentration was a function of astrocytic K⁺ concentration. Farr and David (2011) and Witthoft and Karniadakis (2012) modeled the K⁺ concentration both in the perisynaptic space [named as synaptic cleft and synaptic space by Farr and David (2011) and Witthoft and Karniadakis (2012)] and in the perivascular space. They both used the same simple model for the K⁺ concentration in the perisynaptic space. K⁺ was released from the neurons into the perisynaptic space and removed from the perisynaptic space via the astrocytic NKA (Farr and David, 2011). However, Witthoft and Karniadakis (2012) named the removal so that it included both astrocytic NKA and astrocytic inwardly rectifying K⁺- and voltage-gated K⁺ (KIR) channel. Farr and David (2011) and Witthoft and Karniadakis (2012) also used the same simple model for the K⁺ concentration in the perivascular space. The K⁺ concentration in the perivascular space depended on astrocytic voltage- and Ca²⁺-operated K⁺ (BK) channels and smooth muscle cell's KIR channels. They also had a K⁺-dependent decay term in the perivascular space. Witthoft et al. (2013) modeled in more detail the K⁺ concentration in the perisynaptic space (named as synaptic space by Witthoft et al., 2013), intracellular space of the astrocyte, and perivascular space. In the model by Witthoft et al. (2013), the K⁺ concentration in the perisynaptic space depended on K⁺ released from the neurons and removed via the astrocytic NKCC1 and NKA, in addition to (bidirectional) astrocytic KIR channels and a decay term. The astrocytic K⁺ concentration depended on K⁺ entering from the perisynaptic space via NKA and NKCC1, in addition to astrocytic KIR channels on the perisynaptic side as well as astrocytic BK and KIR channels on the perivascular side of the astrocyte. The equation also had a decay term. The K⁺ concentration in the perivascular space depended on astrocytic BK and KIR channels as well as arteriolar KIR channels and a decay term. Kenny et al. (2018) modeled the K⁺ concentration in the perisynaptic space (named as synaptic cleft by Kenny et al., 2018), intracellular space of the astrocyte, perivascular space, intracellular space of the smooth muscle cell, and extracellular space. In the model by Kenny et al. (2018), the K⁺ concentration in the perisynaptic space depended on K⁺ released from the neuron and removed via the astrocytic K⁺/Cl⁻ cotransporter (KCC1), NKCC1, and NKA, in addition to K⁺ diffusion between extracellular space and perisynaptic space as well as astrocytic K⁺ channels. The astrocytic K⁺ concentration depended on K⁺ entering from the perisynaptic space via KCC1, NKCC1, and NKA, in addition to K⁺ channels on the perisynaptic side and BK channels on the perivascular side of the astrocyte. The K⁺ concentration in the perivascular space depended on astrocytic BK channels and smooth muscle cell's KIR channels. In conclusion, only the model by Witthoft et al. (2013) took into account spatial K⁺ buffering.

Some of the most recent models developed in this category were the models by Komin et al. (2015), Handy et al. (2017), and Taheri et al. (2017). Komin et al. (2015) presented two models, a reaction-diffusion model and a reaction model. With both models they tested if the temperature-dependent SERCA activity was the reason for the differences in Ca²⁺ activity. They showed that their reaction-diffusion model behaved similarly to the experimental data, thus increased SERCA activity (higher temperature) led to decreased Ca²⁺ activity. On the other hand, their reaction model showed the opposite. Thus, they claimed that spatiality was needed to be taken into account to get biologically correct results. However, since the core models were different in the reaction-diffusion and reaction models, it would be interesting to see how the results would look like if the same core model was tested with and without diffusion.

Handy et al. (2017) and Taheri et al. (2017) used the same model but explored somewhat different parameter spaces. They studied the role of SOC channels as well as the PMCA and SERCA pumps in Ca²⁺ activity. They specifically tested which form the Ca²⁺ response had with different parameter values of the channel and pumps (single peak, multiple peaks, plateau, or long-lasting response). They found out that SOC channels were necessary for plateau and long-lasting responses as well as for stable oscillations with multiple peaks. Stable oscillations disappeared when the SERCA pump was partially blocked, but plateau and long-lasting responses were still present. The likelihood of having multiple peaks increased when the PMCA pump was blocked. Taheri et al. (2017) also did Ca²⁺ imaging on cortical astrocytes in mice. They applied ATP on acute brain slices and recorded the Ca²⁺ responses from different subcompartments of the astrocytes, from soma as well as from large and short processes, and categorized the results into four different types of responses named above. Their conclusion was that the variability mainly stemmed from differences in IP₃ dynamics and Ca²⁺ fluxes via SOC channels. To take into account the experimental variability between the different subcompartments, Taheri et al. (2017) ran simulations with different parameter values of the SOC channel and the PMCA and SERCA pumps together with the input IP₃ kinetics. Next, they chose the parameter subspaces that matched their experimental data the best. Their simulations suggested that the soma had higher PMCA and lower SERCA flux rates as well as shorter rise duration for the IP₃ transient than the small and large processes.

Half of the astrocyte network models were so-called generic. Others, however, were specified to model astrocytes in the cerebrum (Iacobas et al., 2006; Ghosh et al., 2010), cortex (Goldberg et al., 2010; Wallach et al., 2014), neocortex (Li et al., 2012), visual cortex and somatosensory cortex (Bennett et al., 2008a), hippocampus (Goto et al., 2004; Ullah et al., 2006), retina (Edwards and Gibson, 2010), spinal cord (Bennett et al., 2006; Gibson et al., 2008), as well as the striatum (Höfer et al., 2002). One fourth of the astrocyte network models took into account neurotransmitters in a simplistic way just as a stimulus, having either the glutamate as a constant, step function, or something similar (see e.g., Goto et al., 2004; Ullah et al., 2006; Bennett et al., 2008a; Kang and Othmer, 2009; MacDonald and Silva, 2013). Only Wallach et al. (2014) actually modeled the amount of neurotransmitter glutamate with a differential equation. The stimulus to the astrocyte model by Wallach et al. (2014) was taken from the model by Tsodyks and Markram (1997). We included this model under astrocyte models because this model was not bidirectional between astrocytes and neurons. The characteristics of astrocyte network models can be found in Table 3.

All the astrocyte network models studied Ca²⁺ waves and few models specifically addressed spontaneous Ca²⁺ waves and vascular events (see Table 3). All the models except the model by Iacobas et al. (2006) had the components for all three; CICR, leak from the ER into the cytosol, and the SERCA pump. About one fourth of the models took into account Ca²⁺ buffering. About one third of the models had either influx of Ca²⁺ from outside of the astrocyte or efflux of Ca²⁺ to outside of the astrocyte, or both. About half of the models took into account astrocytic release of signaling molecules. Thus, the models had equations mostly for extracellular ATP, but one considered equations for extracellular glutamate (Bellinger, 2005). However, none of the models presented a detailed mechanistic description of how the release occurs. More than half of the models took into account diffusion, and, especially, almost half of the models studied the ATP diffusion in the extracellular space. Three quarters of the astrocyte network models had gap junctions for IP₃ but some models had them also for Ca²⁺. Thus, these models had similar core structure with small variations. As an example, Li et al. (2012) were the only ones that modeled K⁺ concentration, both in astrocytic and extracellular spaces, and VGCCs. Goto et al. (2004) were the only ones to use the detailed IP₃R model by De Young and Keizer (1992). Höfer et al. (2002), Bellinger (2005), Ullah et al. (2006), Kazantsev (2009), Ghosh et al. (2010), and Matrosov and Kazantsev (2011) modeled CCE.

The first model developed in this category was the model by Höfer et al. (2002). Höfer et al. (2002) showed with their two-dimensional (19 × 19) astrocyte network model that IP₃ permeability in gap junctions was a more important factor in intercellular Ca²⁺ waves than Ca²⁺ permeability. When blocking the IP₃ permeability, intercellular Ca²⁺ wave propagation was prevented. However, intercellular Ca²⁺ wave propagation was not prevented when the Ca²⁺ permeability was blocked. In the model by Höfer et al. (2002), the IP₃ concentration depended on two distinct production terms via phospholipase C (PLC), one corresponding to PLC isotype β (PLCβ) and the other to PLC isotype δ (PLCδ). Höfer et al. (2002) showed that PLCδ was needed to be modeled to get the downstream cells to respond to the stimulus with a Ca²⁺ increase.

Two of the newest models developed in this category were the models by Lallouette et al. (2014) and Wallach et al. (2014). Lallouette et al. (2014) simplified the astrocyte network model by Goldberg et al. (2010) to be able to simulate the function of a three-dimensional (11 × 11 × 11) astrocyte network. With this network of more than a thousand astrocytes, Lallouette et al. (2014) studied how the variability in the topology of gap-junction coupled astrocytes affected the intercellular Ca²⁺ wave propagation. They tested five different coupling rules and found out that these different coupling rules can be used to reproduce the variation in the experimental data. They showed that dense connectivity or having long-distance gap-junction coupled astrocytes reduced the intercellular Ca²⁺ wave propagation. Wallach et al. (2014) continued the previous study by stimulating astrocyte network by Tsodyks and Markram (1997) model. In the present study, the model by Wallach et al. (2014) is listed in the category of astrocyte network models because the astrocytes did not have an effect on the neuron. Wallach et al. (2014) demonstrated through experimental and simulation studies that there was a threshold stimulation frequency when astrocytes started to respond with Ca²⁺ oscillations. However, this threshold frequency was different for different astrocytes and it increased with the number of astrocytes coupled.

Neuron-astrocyte synapse models include models which have only one astrocyte and one to multiple synapses. Of the models covered in our study, about half of the neuron-astrocyte synapse models were found to be so-called generic models, in other words they were developed with no specific brain area or cell in mind. Others, however, were specified to model neuron-astrocyte synapses in the cortex or neocortex (Nadkarni and Jung, 2003, 2004; Di Garbo et al., 2007; Volman et al., 2007; DiNuzzo et al., 2011; Valenza et al., 2011; Nazari et al., 2015b, 2017; Amiri et al., 2016; Li et al., 2017), hippocampus (Nadkarni and Jung, 2005, 2007; Nadkarni et al., 2008; Tewari and Majumdar, 2012a,b; Tang et al., 2013, 2016; Tewari and Parpura, 2013; Li et al., 2016b), thalamocortical networks (Amiri et al., 2011a), or brainstem (Oku et al., 2016).

The modeling approaches for neurons varied depending on the author. Half of the studied publications utilized relatively complex biophysical neuron models, namely Hodgkin and Huxley (1952) model (Nadkarni and Jung, 2003; Sotero and Martínez-Cancino, 2010; Tewari and Majumdar, 2012a,b; Tang et al., 2013, 2016; Guo et al., 2017; Li et al., 2017), Traub et al. (1991) model or its derivative Pinsky and Rinzel (1994) model (Nadkarni and Jung, 2004, 2005, 2007; Silchenko and Tass, 2008; Tewari and Parpura, 2013; Li et al., 2016a,b), or Olufsen et al. (2003) model (Di Garbo et al., 2007; Di Garbo, 2009). Simpler phenomenological models used in the studied publications were the FitzHugh-Nagumo (FitzHugh, 1961) model (Postnov et al., 2007), leaky integrate-and-fire (LIF, Gerstner and Kistler, 2002) model (Wade et al., 2011, 2012; Amiri et al., 2016; De Pittà and Brunel, 2016; Nazari et al., 2017; Liu et al., 2018), Izhikevich (2007) model (Valenza et al., 2011; Nazari et al., 2015a,b,c), Morris and Lecar (1981) model or its derivatives (Volman et al., 2007; Amiri et al., 2011b, 2013b), and Suffczynski et al. (2004) neuronal population model (Amiri et al., 2011a). The released neurotransmitter was modeled explicitly by Nadkarni and Jung (2005, 2007), Volman et al. (2007), Nadkarni et al. (2008), Silchenko and Tass (2008), Amiri et al. (2011b, 2013b), Wade et al. (2011), Tewari and Majumdar (2012a,b), Tang et al. (2013, 2016), Tewari and Parpura (2013), Nazari et al. (2015a,b,c), De Pittà and Brunel (2016), and Li et al. (2017), however not necessarily with a differential equation. Other models utilized, for example, phenomenological transfer functions between neuronal events and astrocytic IP₃ concentration. Several models had increased astrocytic IP₃ concentration after the neuronal membrane potential increased over some threshold (marked in Table 4 as V_m,N → [IP₃]). The details of the neuron-astrocyte synapse models can be found in Table 4.

The neuron-astrocyte synapse models were developed to explain many different biological events as can be seen in Table 4. Examples of the studied phenomena included Ca²⁺ dynamics, plasticity, hyperexcitability, information transfer, synchronization, and vascular events. All the other models except the models by Oku et al. (2016) and Guo et al. (2017) had all astrocytic components for CICR, leak from the ER into the cytosol, and the SERCA pump. About one third of the models had influx of Ca²⁺ from outside of the astrocyte and efflux of Ca²⁺ to outside of the astrocyte. The models had neither gap junction signaling because these models had only one astrocyte nor buffers. Thus, these models had similar core structure with small variations. As an example, two modeled CCE (Di Garbo et al., 2007; Di Garbo, 2009) and one modeled efflux via the PMCA pump and NCX (DiNuzzo et al., 2011) as well as diffusion (Silchenko and Tass, 2008).

The first model developed in this category was the model by Nadkarni and Jung (2003). Nadkarni and Jung (2003) studied Ca²⁺ oscillations with a model consisting of a single astrocyte and single Hodgkin-Huxley type neuron. The astrocyte model was based on the model by Li and Rinzel (1994) but they added a differential equation for IP₃. When the membrane potential of the neuron was higher than certain threshold, IP₃ concentration was increased in the astrocyte. The model also included a degradation of IP₃. The membrane potential of the neuron depended on additional inward current (I_astro) that depended on astrocytic Ca²⁺. Nadkarni and Jung (2003) showed that with higher rate of IP₃ production, the astrocytic Ca²⁺ concentration was higher. In addition, with large enough rate of IP₃ production, the neuron oscillated spontaneously without any external stimulus. Their hypothesis was that a higher expression of metabotropic glutamate receptors (mGluRs) and the resulting spontaneous oscillations caused epilepsy.

As the exact mechanisms of signaling from astrocytes to neurons are still largely unknown or controversial, most of the models had phenomenological descriptions of gliotransmitter release. About one third of the models took into account gliotransmitter release with equations for extracellular glutamate (see Table 4). Other models used phenomenological transfer functions to relay the effect of gliotransmission to the target synaptic terminal (I_astro, part of I_ast, f, and G_m). Only a few of the studied models had detailed astrocytic vesicle release model, namely the models by Silchenko and Tass (2008), Tewari and Majumdar (2012a,b), and Tewari and Parpura (2013). Silchenko and Tass (2008) modeled five vesicle states which represented different stages in the transmitter vesicle cycle. The rate of prefusion complex formation depended on astrocytic Ca²⁺ concentration which was modeled by two steps. In the first step, they simplified the equation where Ca²⁺ activated Ca²⁺-binding soluble N-ethylmaleimide-sensitive factor attachment protein receptor (SNARE) proteins by assuming that the concentration of activated SNARE-proteins was considered stationary. In the second step, they simplified the equation for the fusion of vesicles leading to an irreversible exocytosis of glutamate. However, Silchenko and Tass (2008) did not provide all the details of the model which makes the reuse of the model difficult. The models by Tewari and Majumdar (2012a,b) and Tewari and Parpura (2013) assumed, based on experimental data on cultured hippocampal astrocytes, that the binding of three Ca²⁺ ions was required for gliotransmitter release. The fusion and recycling process of the synaptic-like micro-vesicle was modeled using two differential equations that both depended on the probability that the synaptic-like micro-vesicle was ready to be released. In addition to these more detailed vesicle release models, De Pittà and Brunel (2016) modeled astrocytic glutamate exocytosis in a phenomenological way with just a few equations. They assumed that a fraction of gliotransmitter resources was available for release at any time. Then, every time astrocytic Ca²⁺ increased beyond a certain threshold, the fraction of readily releasable astrocytic glutamate resources was released into the periastrocytic space.

Two of the newest models were provided by Li et al. (2016a, 2017). However, these studies contained, to the best of our understanding, fundamental errors in the biological terminology. Basically, the model by Li et al. (2016a) was the same as presented by Nadkarni and Jung (2004), but the neuronal membrane potential depended on astrocytic glutamate, as presented by Postnov et al. (2009), instead of astrocytic Ca²⁺, as presented by Nadkarni and Jung (2004). Li et al. (2017) developed a GABA-activated astrocyte model (which they, misleadingly, termed GABAergic). The model by Li et al. (2017) is similar to the model by Li et al. (2016a), but Li et al. (2017) added a more complex differential equation for IP₃ by taking into account both the GABA released by the interneuron and glutamate released by the astrocyte, somewhat similarly to Ullah et al. (2006), Nadkarni and Jung (2005), Volman et al. (2007), and others. The differential equations for the extracellular glutamate released by the astrocyte had similar forms as the IP₃ equations and were somewhat similar to the equation by Wade et al. (2012). Li et al. (2016a) showed how a higher equilibrium concentration of extracellular glutamate or glutamate degradation time constant predicted a higher neuronal firing frequency and existence of epileptic seizures. Li et al. (2017), on the other hand, presented using their GABA-activated astrocyte model (misleadingly termed GABAergic) that after a 0.5 s long GABA stimulation, astrocytic Ca²⁺ oscillations were long-lasting. After combining the GABA-activated astrocyte model (misleadingly termed GABAergic) and a neuronal seizure model, they concluded that in this model, the astrocyte, through stimulating pyramidal neurons and thus increasing excitatory activity, prevented the transition from seizure activity into a normal firing activity state, which GABA alone was capable of inducing by inhibiting pyramidal neuron activity.

Neuron-astrocyte network models include models that have several astrocytes in addition to neurons. Half of the neuron-astrocyte network models were so-called generic models. Others, however, were specified to model neuron-astrocyte interactions in the cortex (Allegrini et al., 2009; Liu and Li, 2013a; Chan et al., 2017; Tang et al., 2017; Yao et al., 2018), hippocampus (Amiri et al., 2012a, 2013a; Mesiti et al., 2015a; Li et al., 2016c), spinal cord (Yang and Yeo, 2015), or thalamocortical networks (Amiri et al., 2012b,c).

The modeling strategies for neurons varied depending on the author. Three of the studied publications utilized Hodgkin and Huxley (1952) model (Liu and Li, 2013b; Li et al., 2016c; Yao et al., 2018) and one utilized Traub et al. (1991) model's derivative Pinsky and Rinzel (1994) model (Mesiti et al., 2015a). Simpler phenomenological models used in the studied publications were the FitzHugh-Nagumo (FitzHugh, 1961) model (Postnov et al., 2009; Hayati et al., 2016), LIF (Gerstner and Kistler, 2002) model (Liu and Li, 2013a; Naeem et al., 2015), Izhikevich (2007) model (Allegrini et al., 2009; Haghiri et al., 2016, 2017; Tang et al., 2017), Morris and Lecar (1981) model or its derivatives (Amiri et al., 2012a, 2013a; Chan et al., 2017), and Suffczynski et al. (2004) neuronal population model (Amiri et al., 2012b,c). The released neurotransmitter was modeled explicitly by Amiri et al. (2012a, 2013a), Liu and Li (2013a), Yang and Yeo (2015), Li et al. (2016c), and Yao et al. (2018). Other models utilized phenomenological transfer functions between the neurotransmitter and astrocytic IP₃ concentration. The details of the neuron-astrocyte network models can be found in Table 5.

The neuron-astrocyte network models were developed to explain many different biological events as can be seen in Table 5. Examples included Ca²⁺ dynamics, synchronization, information transfer, plasticity, and hyperexcitability. All the other models except the model by Allegrini et al. (2009) had components for all three; CICR, leak from the ER into the cytosol, and the SERCA pump. More than half of the models had influx of Ca²⁺ from outside of the astrocyte and efflux of Ca²⁺ to outside of the astrocyte. About one third of the models took into account gliotransmitter release by modeling extracellular glutamate, and few were also modeling extracellular ATP. Other models used phenomenological transfer functions to relay the effect of gliotransmission to the target synaptic terminal (I_astro, I_syn, part of I_ast, and G_m). None of the studied models had a detailed astrocytic vesicle release model. Most of the models had gap junction signaling for IP₃, and some also for Ca²⁺. Thus, these models had a similar core structure with small variations. As an example, only Yao et al. (2018) modeled buffering as well as astrocytic and extracellular K⁺. Diffusion was taken into account in the models by Allegrini et al. (2009), Postnov et al. (2009), Mesiti et al. (2015a), Yang and Yeo (2015), Li et al. (2016c), and Yao et al. (2018). Yao et al. (2018) presented one of the available models for cortical spreading depression.

One of the first models developed in this category was the two-dimensional model by Postnov et al. (2009). They studied how different lengths of stimulus affected astrocytic Ca²⁺ and showed how short stimulus of less than 100 s did not induce Ca²⁺ wave propagation. However, a longer stimulus of 320 s showed Ca²⁺ wave propagation for a short distance and a stimulus of about 2,000 s showed Ca²⁺ wave propagation along the astrocyte network. They also tested how Ca²⁺ wave propagation was affected by different noise levels added to the model. They found out that the stronger the noise, the more accelerated was the Ca²⁺ wave propagation. With the largest noise level they tested, they found out that the spatially synchronized behavior was destroyed, and the model started to behave irregularly.

A few publications presented simplification of model complexity. Simplification is, in general, used to reduce the model order to allow cost-effective computation yet preserving the major, key dynamical behavior of the model. Soleimani et al. (2015), Haghiri et al. (2016, 2017), and Hayati et al. (2016) presented the original and simplified versions of the earlier published models by Postnov et al. (2007, 2009). However, most of the reduced astrocyte models were not detailed enough based on our criteria in section 2.2. In the future, it is important to put more emphasis on the model order reduction of the complex neuron-astrocyte interaction models to be able to simulate the behavior of large networks biologically more accurately (see e.g., Lehtimäki et al., 2017).

One of the newest models was the model by Chan et al. (2017). Their neuron-astrocyte network model was able to generate ultra-slow neuronal oscillatory activity patterns with frequencies less than 0.01 Hz, recorded by high-density microelectrode arrays and examined using methods reflecting conventional EEG analysis rather than conventional microelectrode array analysis. Their model showed that the frequency of these neuronal oscillations depended on the astrocytic Ca²⁺ oscillations. These astrocytic Ca²⁺ oscillations depended on the properties of IP₃Rs and lasted from a few seconds to minutes. The results of Chan et al. (2017) also suggested that astrocytes preferred asynchronous neuronal firing to synchronous neuronal firing.