# A Theoretical Study on the Role of Astrocytic Activity in Neuronal Hyperexcitability by a Novel Neuron-Glia Mass Model

Aurélie Garnier, Alexandre Vidal, Habib Benali

► To cite this version:

Aurélie Garnier, Alexandre Vidal, Habib Benali. A Theoretical Study on the Role of Astrocytic Activity in Neuronal Hyperexcitability by a Novel Neuron-Glia Mass Model. Journal of Mathematical Neuroscience, 2016, 6 (1), pp.10. 10.1186/s13408-016-0042-0 . hal-01435048

HAL Id: hal-01435048

https://hal.sorbonne-universite.fr/hal-01435048

Submitted on 13 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# A Theoretical Study on the Role of Astrocytic Activity in Neuronal Hyperexcitability by a Novel Neuron-Glia Mass Model

Aurélie Garnier<fnref n="1" /> · Alexandre Vidal<fnref n="2" /> ·

Habib Benali<fnref n="1" />

Received: 28 January 2016 / Accepted: 26 November 2016 / Published online: 21 December 2016
© The Author(s) 2016. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

**Abstract** Recent experimental evidence on the clustering of glutamate and GABA transporters on astrocytic processes surrounding synaptic terminals pose the question of the functional relevance of the astrocytes in the regulation of neural activity. In this perspective, we introduce a new computational model that embeds recent findings on neuron-astrocyte coupling at the mesoscopic scale intra- and inter-layer local neural circuits. The model consists of a mass model for the neural compartment and an astrocyte compartment which controls dynamics of extracellular glutamate and GABA concentrations. By arguments based on bifurcation theory, we use the model to study the impact of deficiency of astrocytic glutamate and GABA uptakes on neural activity. While deficient astrocytic GABA uptake naturally results in increased neuronal inhibition, which in turn results in a decreased neuronal firing, deficient glutamate uptake by astrocytes may either decrease or increase neuronal firing either transiently or permanently. Given the relevance of neuronal hyperexcitability (or lack thereof) in the brain pathophysiology, we provide biophysical conditions for the onset identifying different physiologically relevant regimes of operation for astrocytic uptake transporters.

**Keywords** ODE model in neuroscience · Qualitative analysis of dynamical systems · Bifurcations · Neuron-astrocyte interactions · GABAergic and

<figure><img src="image_4.png" /></figure>

<footnote>A. Vidal
alexandre.vidal@univ-evry.fr</footnote>

<footnote>A. Garnier
aurelie.garnier@lib.upmc.fr</footnote>

<footnote>H. Benali
habib.benali@lib.upmc.fr</footnote>

<footnote n="1">¹ Laboratoire d'Imagerie Biomédicale (LIB), Sorbonne Universités, UPMC Univ Paris 06, CNRS, INSERM, Paris, 75013, France</footnote>

<footnote n="2">² Laboratoire de Mathématiques et Modélisation d'Évry (LaMME), CNRS UMR 8071, Université d'Évry-Val-d'Essonne, Évry, 91000, France</footnote>
glutamatergic neurotransmissions · Excitability modulation · Neuronal hyperexcitability

# 1 Introduction

Neural activity in the brain results from an interplay between neuronal excitation and inhibition and is subjected to the supply of nutrients by the cerebral blood flow. At the subcellular (microscopic) scale, neural activity propagates through axons and, mostly, chemical synapses which transmit electrical signals from one (presynaptic) neuron to another (postsynaptic) one, by release of neurotransmitters, like glutamate and GABA for example, in the extracellular space (or synaptic cleft) between pre- and postsynaptic terminals [1, 2]. These neurotransmitters diffuse in the synaptic cleft, binding and activating postsynaptic receptors, which can regulate the depolarization state of the postsynaptic neuron, triggering firing thereof. In parallel, spillout of these neurotransmitters from the synaptic cleft to the perisynaptic extracellular space may also activate receptors and transporters on perisynaptic astrocytic processes, with the potential to trigger glutamate release by the latter glutamate into the extra-synaptic space [3, 4].

The predominant expression of high affinity astrocytic glutamate and GABA transporters in the proximity of synaptic terminals [5–7] raises the question of the functional relevance of the astrocytic uptake in the regulation of synaptic transmission and thus in the tone of neural excitability in the brain pathophysiology [8–11]. In this perspective, computational modeling provides a key tool to aid diagnostic and possible treatment of brain pathophysiology, providing a rationale to interpret electrophysiological data when, at the current state of technology, dynamical recordings of intracranial electroencephalography (EEG) cannot provide any information on underlying concentrations of neurotransmitters. Over the past decade, modeling efforts made in this regard have exploited several directions. There exists indeed a wealth of published models for neuron-astrocyte interactions which consider neurotransmitter and ion dynamics at different scales: subcellular [12–16], cellular [4, 17–19] and supracellular, including a hemodynamic compartment [20–22], and network [23, 24]. Despite the different scales of analysis, consensus has grown from these studies as to the existence of a plausible causal link between astrocyte activation and neuronal (hyper)excitability [13, 15, 16, 23, 25]. With few exceptions [4, 16, 17], however, a main limitation of these models is that they are hardly mathematically tractable and the computational demand for scaling them up to meso/macroscopic scales requires dedicated computing platforms. A possible alternative which brings together computational portability with mathematical tractability is represented by neuron-astrocyte mass modeling of dynamical interactions between neurons and astrocytes populations. Indeed the mass approach is well suited to reproduce mesoscopic data obtained by current EEG technology [26]. The aggregated yet biologically significant parameters of a mass model can in fact efficiently be chosen so as to fit experimental data [27]. Moreover, such models are usually low dimensional and thus amenable to analytical treatment aimed at interpreting experimental observations as well as identifying parameter sets for different dynamical regimes of the model [28, 29]. In
particular, we benefit from the bifurcation analysis of the Neural Mass Model (NMM) performed in [29] to characterize in this study, the impact of astrocyte on neuronal activity. More exactly, our focus for the present study is on Noise Induced Spiking (NIS), which is characterized by quiescent phases of neuronal activity separated by isolated spikes resulting from synchronous neuronal activations. NIS typically mimics local field potentials (LFP) recordings obtained during pre-ictal phases of epilepsy as shown in [29] and may thus be used as a predictor for the insurgence of neuronal hyperexcitability.

A novel computational model based on the neural mass approach that focuses on astrocyte dynamics at the mesoscopic scale has recently been proposed [30]. This model links the LFP signal representative of neural activity measured by intracranial EEG to the cerebral blood flow dynamics measured by laser doppler recordings of astrocyte activity. Although this model incorporates astrocytic recycling of glutamate and GABA it does not take into account the modulation of neuronal excitability by extracellular neurotransmitter concentrations. In the following, we present an extension of this model that includes such modulation in the context of bilateral neuron-astrocyte interactions. In particular, we pursue an analytical characterization of the impact of the key model parameters of astrocytic activity on neuronal hyperexcitability.

Analytical investigation of our model allows us to make several predictions on the effect of deficiencies of GABA and glutamate astrocytic uptake on neural activity. In particular, we identify three possible neural regimes in the presence of deficiency of astrocytic glutamate uptake, respectively, consistent with reduced activity, transient or permanent hyperexcitability. Such a spectrum of responses is substantiated by the analysis of the dynamical structure of the model, which allows us to derive explicit conditions on the parameters involved in the astrocytic feedback corresponding to each type of regime and put emphasis on the delicate balance between neuronal excitation and inhibition and its sensitivity to extracellular concentrations of neurotransmitters.

## 2 Bilaterally Coupled Neuron-Astrocyte Mass Modeling

Physiologically, a pyramidal neuron (resp. an interneuron) releases glutamate (resp. GABA) in the synaptic cleft from where it binds to receptors on the postsynaptic neuron, which stands for the main mechanism of synaptic transmission. In the NMM, the role of neurotransmitter fluxes in this mechanism is implicit since the neuronal activity is reproduced at the scale of populations. Henceforth, the tonic impacts of neurotransmitters (GABA and glutamate) are embedded in conversion processes from average pulse density into postsynaptic potentials. Uptake processes of the neurotransmitters by the local astrocyte and presynaptic neuron regulates their concentration in the extracellular space (Fig. 1). In the presynaptic neuron, the uptake completes the stock whereas the uptake by the astrocyte triggers a cascade of reactions linked with the modulation of synaptic transmission (differentially according to the type of neurotransmitters) and the hemodynamics. This indirect mechanism involves concentrations represented by explicit variables of the glial compartment of
<figure><img src="image_5.png" /><figcaption>**Fig. 1** Scheme of neurotransmission mechanisms and neurotransmitter uptake. Red circles: glutamatergic vesicles. Blue circles: GABAergic vesicles. Red arrows: exchanges of glutamate. Blue arrows: exchanges of GABA. A: astrocytes</figcaption></figure>

our neuron-astrocyte mass model. Hence, feedbacks from the glial compartment onto the neural dynamics stand for the modulations of excitability (or more precisely the activation threshold of each population) by the extracellular concentration of each neurotransmitters.

In the following, we set to describe our mass modeling approach to account for the effect of extracellular glutamate and GABA concentrations on neural activity. With this aim, we start by recalling the NMM setup originally studied in [29, 30], its properties and time-series pattern that it produces. We then introduce our modeling arguments to develop bilateral coupling between the astrocytic and neuronal compartments of the model. Finally, we leverage on the novel dynamical features brought forth by inclusion of this bilateral coupling with respect to the original model in [29, 30].

## 2.1 Neural Mass Model

The NMM represents the dynamical interactions between two neural populations at a mesoscopic scale: a main population of pyramidal neurons (P) and a population of inhibitory interneurons (I). For the excitatory feedback of pyramidal neurons, two approaches have been considered in the literature: a direct link from their output to their input and an indirect track through synaptic coupling with distant neurons. The latter way of modeling amounts to considering an intermediary population of pyra-
<figure><img src="image_6.png" /><figcaption>Fig. 2 Two schematic representations of the NMM with double excitatory feedbacks. P: main population of pyramidal neurons. I: interneuron population. P': secondary population of pyramidal neurons. Red (resp. green) arrows in (a): excitatory (resp. inhibitory) interactions. Box hE (resp. hI): second-order process converting action potentials into excitatory (resp. inhibitory) postsynaptic potential. Box $\hat{S}$: process converting average membrane potential into average action potential density discharge by neurons of populations P, P', and I, respectively. Cᵢ for 1 ≤ i ≤ 4: coupling gain parameters depending on the maximum number C of synaptic connections between two populations. G: direct feedback coupling gain. p(t): external input. y₀, y₁, y₂: state variables. x₀, x₁, x₂: intermediary variables</figcaption></figure>

midal neurons P' interacting with P through synaptic connections. Both direct and indirect approaches model the excitatory synaptic interaction between neighbor principal cells, usually named collateral excitation. From the modeling perspective, we cannot privilege one type of feedback over the other, since both these couplings are physiologically relevant and can co-exist, a very local one and a more or less distant one. The NMM used in this article, proposed and studied in [29], includes both feedback circuits. Hence, it includes three feedback loops on population P activity: an inhibitory feedback through the interneuron population I, a direct excitatory feedback of P onto itself (referred to as “direct feedback”) and an indirect excitatory feedback (referred to as “indirect feedback”) involving the population P' (Fig. 2(a)).

The conversion process of average pulse density into excitatory and inhibitory postsynaptic potential, respectively, are based on the following α-functions introduced by Van Rotterdam et al. [31]:

$$h_E(t) = A a t e^{-a t},$$

$$h_I(t) = B b t e^{-b t}.$$

These are fundamental solutions (also known as Green's functions) for the differential operators $\mathcal{F}_E$ and $\mathcal{F}_I$, respectively, defined by

$$\mathcal{F}_E(h_E) = \frac{1}{A} \left( \frac{1}{a} \frac{d^2 h_E}{dt^2} + 2h'_E + ah_E \right), \qquad (1a)$$

$$\mathcal{F}_I(h_I) = \frac{1}{B} \left( \frac{1}{b} \frac{d^2 h_I}{dt^2} + 2h'_I + bh_I \right). \qquad (1b)$$
In this framework, the parameter A (resp. B) stands for the average excitatory (resp. inhibitory) synaptic gain and tunes the amplitude of excitatory (resp. inhibitory) post-synaptic potentials. Additionally, $\frac{1}{a}$ (resp. $\frac{1}{b}$) represents the time constant of excitatory (resp. inhibitory) postsynaptic potentials representative of the kinetics of synaptic connections and delays introduced by circuitry of the dendritic tree [31–33]. Following Freeman’s work [32], the functions converting the average membrane potential into an average pulse density can be approximated by sigmoids of the type

$$S(x, x_{\text{th}}, r_{\text{sl}}) = \frac{1}{1 + e^{r_{\text{sl}}(x_{\text{th}} - x)}}.$$

Yet, for the sake of compactness of the NMM presentation, we introduce an auxiliary parameterization after [32]:

$$\hat{S}(x, v) = 2e_0 S(x, v, r) = \frac{2e_0}{1 + e^{r(v-x)}}$$

where $2e_0$ represents the maximum discharge rate, $v$ the excitability threshold and $r$ the sigmoid slope at the inflection point. Finally, the NMM receives an excitatory input $p(t)$ standing for the action on population P of neural populations in other areas through long-range synaptic connections. This input can be either deterministic or stochastic, and in this case, being modeled for example by a Gaussian process.

Now we can write the dynamics for the intermediary variables $x_0$, $x_1$ and $x_2$, which represent the outputs of the population P, the sum of population $P'$ output and the external input $p(t)$, and the population I, respectively (Fig. 2(b)):

$$\frac{d^2 x_0}{dt^2} = A a \hat{S}(x_1 + G x_0 - x_2, v_P) - 2a \frac{dx_0}{dt} - a^2 x_0, \quad (2a)$$

$$\frac{d^2 x_1}{dt^2} = A a C_2 \hat{S}(C_1 x_0, v_{P'}) - 2a \frac{dx_1}{dt} - a^2 x_1 + A a p(t), \quad (2b)$$

$$\frac{d^2 x_2}{dt^2} = B b C_4 \hat{S}(C_3 x_0, v_I) - 2b \frac{dx_2}{dt} - b^2 x_2. \quad (2c)$$

Parameters $C_i$ ($i \in \{1, 2, 3, 4\}$) represent the average number of synapses between two populations. Following [34], each $C_i$ is proportional to the maximum number C of synapses between two populations. The excitation of P by its own output, resulting from the intra-population synaptic connections, is weighted by the coupling gain G.

For the sake of comparison, we use a variable change to obtain the same state variables as in the Jansen–Rit model [35]: the excitatory output ($y_0 = x_0$) and the excitatory ($y_1 = x_1 + Gx_0$) and inhibitory ($y_2 = x_2$) inputs of the main population P. The output $y_0$ acts on the secondary pyramidal neuron population $P'$ and on the interneuron population I. To analyze the model, we write the dynamics of the state variables $y_0$, $y_1$, and $y_2$ as a system of first order differential equations:

$$\frac{dy_0}{dt} = y_3, \quad (3a)$$

$$\frac{dy_1}{dt} = y_4, \quad (3b)$$
$$ \frac{dy_2}{dt} = y_5, \qquad (3c) $$

$$ \frac{dy_3}{dt} = Aa\hat{S}(y_1 - y_2, v_P) - 2ay_3 - a^2y_0, \qquad (3d) $$

$$ \frac{dy_4}{dt} = AaC_2\hat{S}(C_1y_0, v_{P'}) + AaG\hat{S}(y_1 - y_2, v_P) - 2ay_4 - a^2y_1 + Aap(t), \qquad (3e) $$

$$ \frac{dy_5}{dt} = BbC_4\hat{S}(C_3y_0, v_I) - 2by_5 - b^2y_2. \qquad (3f) $$

In this study, we consider the local field potential LFP(t) = y₁(t) - y₂(t) [33] as the main output of the model. It is important to note that, generally, studies of neural mass models, such as Jansen-Rit model, only considered the case with the same constant excitability thresholds for all populations, i.e.

$$ v_P = v_{P'} = v_I = v_0. $$

## 2.2 Astrocyte Model: Glutamate and GABA Concentration Dynamics

For reproducing the astrocyte activity, we use the model introduced in [30]. It focuses on the dynamics of glutamate and GABA concentrations, which are the main neurotransmitters of the central nervous system. In [30], the neuron-astrocyte coupling is feedforward: the astrocyte dynamics is driven by the neural activity, generated by the Jansen-Rit model, but it does not impact the neural compartment. The model considers the dynamics of glutamate and GABA concentrations, locally to the main population P of pyramidal neurons, at different stages of the recycling mechanism. The local nature of this interaction implies that the firing rate of the secondary population P' of pyramidal neurons does not impact the astrocyte dynamics associated with the neighboring astrocytes of the main population P. The mechanism is as follows (Fig. 3): excited pyramidal neurons (resp. interneurons) release glutamate (resp. GABA) in the extracellular space (synaptic cleft). Astrocytes and presynaptic neurons uptake the neurotransmitters. Astrocytes recycle or consume the neurotransmitters, while the presynaptic neurons capture them to complete their stock.

Following [30] the astrocyte compartment is built on the firing rate of the pyramidal neuron population and the firing rate of the interneuron population. The state variables are (i) $J_G$ and $J_γ$: the fluxes of glutamate and GABA from neurons to extracellular space, (ii) [Glu]e and [GABA]e: the neurotransmitter concentrations in the extracellular space, (iii) [Glu]a and [GABA]a: the quantity of neurotransmitters recycled and consumed by the astrocytes. The dynamics governing $J_G$ and $J_γ$ is driven by second-order differential operators similar to those for synaptic transfer dynamics introduced in (1a)-(1b) [31, 36], i.e.:

$$ \mathcal{F}_G(h_G) = \frac{1}{W} \left( \frac{1}{w_1} \frac{d^2h_G}{dt^2} + \frac{w_1 + w_2}{w_1} \frac{dh_G}{dt} + w_2h_G \right), $$

$$ \mathcal{F}_\gamma(h_\gamma) = \frac{1}{Z} \left( \frac{1}{z_1} \frac{d^2h_\gamma}{dt^2} + \frac{z_1 + z_2}{z_1} \frac{dh_\gamma}{dt} + z_2h_\gamma \right). $$
<figure><img src="image_7.png" /><figcaption>**Fig. 3** Neuron-astrocyte model with astrocyte feedback. *P* and *P*': main and secondary populations of pyramidal neurons. *I*: interneuron population. *p*(*t*): external input on population *P*. [Glu]*e* and [GABA]*e*: glutamate and GABA extracellular concentrations. [Glu]*a* and [GABA]*a*: glutamate and GABA astrocyte concentrations. Red arrows: *P* → *P*, *P* → *I* and *P* → *P'* couplings. Orange arrow: *P'* → *P* coupling. Green arrow: *I* → *P* coupling. Cyan arrows: glutamate and GABA release by populations *P* and *I* into extracellular space (fluxes *J*<sub>*γ*</sub><sup>*en*</sup> and *J*<sub>*γ*</sub>). Purple arrows: astrocyte and neural uptakes of neurotransmitters. Red dashed arrows: glutamate feedbacks on populations *P* and *I*. Brown dashed arrow: GABA feedback on population *P*</figcaption></figure>

Like for the synaptic transfer functions in (1a)-(1b), the parameter *W* (resp. *Z*) tunes the peak amplitude of glutamate (resp. GABA) concentrations and the parameters *w*₁ and *w*₂ (resp. *z*₁ and *z*₂) tune the rise and decay times of glutamate (resp. GABA) release transfer function. This dynamics is well suited to reproduce qualitative and quantitative properties of rise and decay of neurotransmitter concentrations as shown in [30].

Astrocytic uptake of glutamate and GABA from the extracellular space may be modeled by Michaelis–Menten kinetics [30], using the following function:

$$ \mathcal{H}(x, k) = \frac{x}{x + k}. \qquad (4) $$

The dynamics of the extracellular concentrations ([Glu]*e* and [GABA]*e*) result from neurotransmitter fluxes *J*<sub>*G*</sub> and *J*<sub>*γ*</sub> originating from synaptic release fluxes as just described. Intracellular concentrations of glutamate and GABA in the astrocyte instead mirror uptake dynamics and a nonspecific first order (linear) kinetics degradation [27].

To summarize, the feedforward model obtained by coupling the neural mass model defined by (3a)-(3f) with astrocytic dynamics introduced in [30] reads

$$ \frac{d^2 y_0}{dt^2} = A a \hat{S}(y_1 - y_2, v_P) - 2a \frac{dy_0}{dt} - a^2 y_0, \qquad (5a) $$
$$ \frac{d^2 y_1}{dt^2} = A a C_2 \hat{S}(C_1 y_0, v_{P'}) + A a G \hat{S}(y_1 - y_2, v_P) \\ - 2a \frac{dy_1}{dt} - a^2 y_1 + A a p(t), \qquad (5b) $$

$$ \frac{d^2 y_2}{dt^2} = B b C_4 \hat{S}(C_3 y_0, v_I) - 2b \frac{dy_2}{dt} - b^2 y_2, \qquad (5c) $$

$$ \frac{d^2 J_G}{dt^2} = W w_1 \hat{S}(y_1 - y_2, v_P) - (w_1 + w_2) \frac{d J_G}{dt} - w_1 w_2 J_G, \qquad (5d) $$

$$ \frac{d[\text{Glu}]_e}{dt} = J_G - (V_G^{\text{ae}} + V_G^{\text{ne}}) S([\text{Glu}]_e, s_g, r_g), \qquad (5e) $$

$$ \frac{d[\text{Glu}]_a}{dt} = V_G^{\text{ae}} S([\text{Glu}]_e, s_g, r_g) - V_G^{\text{ne}} [\text{Glu}]_a, \qquad (5f) $$

$$ \frac{d^2 J_\gamma}{dt^2} = Z z_1 \hat{S}(C_3 y_0, v_I) - (z_1 + z_2) \frac{d J_\gamma}{dt} - z_1 z_2 J_\gamma, \qquad (5g) $$

$$ \frac{d[\text{GABA}]_e}{dt} = J_\gamma - V_\gamma^{\text{ae}} \mathcal{H}([\text{GABA}]_e, K_\gamma^{\text{ae}}) - V_\gamma^{\text{ne}} \mathcal{H}([\text{GABA}]_e, K_\gamma^{\text{ne}}), \qquad (5h) $$

$$ \frac{d[\text{GABA}]_a}{dt} = V_\gamma^{\text{ae}} \mathcal{H}([\text{GABA}]_e, K_\gamma^{\text{ae}}) - V_\gamma^{\text{ne}} [\text{GABA}]_a. \qquad (5i) $$

In the sigmoidal functions for [Glu]e dynamics (5e), parameters $V_G^{\text{ae}}$ and $V_G^{\text{ne}}$ are the maximum rates of glutamate uptakes by the astrocytes and the neurons, respectively, $s_g$ represents the activation threshold and $r_g$ the sigmoidal slope at the inflection point. Parameters $V_\gamma^{\text{ae}}$ and $K_\gamma^{\text{ae}}$ (resp. $V_\gamma^{\text{ne}}$ and $K_\gamma^{\text{ne}}$) in Eq. (5h) are, respectively, the maximum rate and concentration for astrocyte (resp. neural) GABA transporter. Finally, $V_G^{\text{c}}$ and $V_\gamma^{\text{c}}$ are the glutamate and GABA degradation rates in astrocytes. We refer the reader to [30] for a detailed explanation of the dynamics.

System (5a)-(5i) is built as a feedforward coupling of the neural compartment onto the astrocyte one. Hence, in this model, the neural compartment is not impacted by the neurotransmitter concentrations in the extracellular space. As mentioned in the introduction, these concentrations have been proven to modulate the local neuron excitability and this feedback has been identified in recent studies [1] to be an essential mechanism of several pathologies triggered by astrocytic uptake deficiencies. Consequently, our aim is to include such feedback in the model in order to study the effects of different astrocyte dysfunctioning on the neuronal activity.

## 2.3 Astrocytic Feedback and Neuron-Astrocyte Mass Model

The concentrations of neurotransmitters in a synaptic cleft act on the excitability threshold of the postsynaptic neuron. In the neuron-astrocyte model (5a)-(5i) the alteration of this neural excitability threshold can be reproduced by dynamical changes in $v_P$, $v_{P'}$ and $v_I$. In the following, we describe how we model the modulation of the neuron excitability in each population by the neurotransmitter concentrations in the extracellular space basing ourselves on biological knowledge.
Extracellular concentrations of neurotransmitters have a thresholded impact on neural activity [1]. Precisely, on one hand, the impact of neurotransmitter concentrations on neural activity is implicitly taken into account in the neuronal compartment, thus the astrocyte feedback steps in only when the concentrations become larger than physiological ones. On the other hand, the postsynaptic neurons are saturated when these concentrations become too large and, consequently, the neural excitability remains bounded. As explained in the introduction, quantitative experimental data of the impact of neurotransmitter concentrations on neural excitability do not exist up to now. For fixing ideas, we consider sigmoidal functions to model the astrocyte feedback on neural excitability which is a natural choice for aggregating the qualitative experimental knowledge. Yet the upcoming mathematical analysis can easily be extended to any bounded increasing functions with a unique inflection point.

We introduce three sigmoidal functions to model the components of the astrocyte feedback:

(a) $m_G^P S([Glu]_e, v_G, r_G)$ for glutamate feedback on pyramidal neurons,

(b) $m_G^I S([Glu]_e, v_G, r_G)$ for glutamate feedback on interneurons,

(c) $m_\gamma S([GABA]_e, v_\gamma, r_\gamma)$ for GABA feedback on pyramidal neurons.

Note that glutamate binding mechanisms are the same for both pyramidal neurons and interneurons since both cell types express the same type of transporters [10, 37]. Thus, only parameters $m_G^P$ and $m_G^I$ representing the maximum coupling gains of the glutamate-dependent component of the astrocyte feedback discriminate between the coupling functions $m_G^P S([Glu]_e, v_G, r_G)$ and $m_G^I S([Glu]_e, v_G, r_G)$, insofar as the number of transporters is cell-specific and so is the ensuing uptake rate. Model parameters are summarized in Table 1 and their values used in the simulations of this study were chosen to qualitatively reproduce typical physiological data.

We previously mentioned that astrocyte feedback acts on the excitability thresholds of neurons. More specifically, if there is an excess of neurotransmitter in a synapse from a neuron $n_1$ of population $N_1$ to a neuron $n_2$ of population $N_2$, the extracellular concentration of neurotransmitter acts on the postsynaptic neuron $n_2$ by changing its excitability threshold. In system (5a)-(5i) the excitability threshold of neurons, that is, a parameter at the single cell scale, does not appear explicitly. However, when the excitability of the population $N_2$ neurons changes at the individual scale, the number of neurons activated in this population by a given input changes as well and so does the output of this population. Following the mass approach, we scale this feature to the population, and, consequently, we choose to change parameter $v_{N_2}$ in the equation corresponding to the output of population $N_2$ since this parameter represents a modulation of the threshold of the sigmoidal function $\hat{S}$.

Let us now describe how we build the feedbacks on the dynamics of the neural compartment using the sigmoidal functions of the neurotransmitter concentrations. The dashed arrows in Fig. 3 illustrate these feedbacks. We need to consider separately each type of synapse in the NMM, and the variables $x_0, x_1$, and $x_2$ for the feedbacks building. The NMM embeds five types of synaptic connections between populations: (i) $P' \to P$, (ii) $P \to P'$, (iii) $P \to I$, (iv) $I \to P$, (v) $P \to P$. In the following
<table><caption>**Table 1** Descriptions and values of the neuron-astrocyte mass model parameters</caption><thead><tr><th>Parameter</th><th>Interpretation</th><th>Value</th></tr></thead><tbody><tr><td colspan="3"><strong>Neuronal compartment</strong></td></tr><tr><td>A</td><td>Average excitatory synaptic gain</td><td>3.25 mV</td></tr><tr><td>B</td><td>Average inhibitory synaptic gain</td><td>22 mV</td></tr><tr><td>1/a</td><td>Time constant of excitatory postsynaptic potentials</td><td>1/100 s</td></tr><tr><td>1/b</td><td>Time constant of inhibitory postsynaptic potentials</td><td>1/50 s</td></tr><tr><td>e<sub>0</sub></td><td>Half of the maximum discharge rate of a neuronal population</td><td>2.5 s<sup>-1</sup></td></tr><tr><td>v<sub>0</sub></td><td>Basic excitability threshold for neurons</td><td>6 mV</td></tr><tr><td>r</td><td>Stiffness of neuronal excitability</td><td>0.56 mV<sup>-1</sup></td></tr><tr><td>C<sub>1</sub></td><td>Strength of the synaptic connections from P to P'</td><td>135</td></tr><tr><td>C<sub>2</sub></td><td>Strength of the synaptic connections from P' to P</td><td>108</td></tr><tr><td>C<sub>3</sub></td><td>Strength of the synaptic connections from P to I</td><td>33.75</td></tr><tr><td>C<sub>4</sub></td><td>Strength of the synaptic connections from I to P</td><td>33.75</td></tr><tr><td>G</td><td>Gain of the direct excitatory feedback from P to itself</td><td>40</td></tr><tr><td colspan="3"><strong>Glial compartment</strong></td></tr><tr><td>W</td><td>Tunes the peak amplitude of glutamate concentrations</td><td>53.6 µM · s<sup>-1</sup></td></tr><tr><td>Z</td><td>Tunes the peak amplitude of GABA concentrations</td><td>53.6 µM · s<sup>-1</sup></td></tr><tr><td>w<sub>1</sub></td><td>Tune the rise and decay times of glutamate release transfer function</td><td>90 s<sup>-1</sup></td></tr><tr><td>w<sub>2</sub></td><td>Tune the rise and decay times of GABA release transfer function</td><td>33 s<sup>-1</sup></td></tr><tr><td>z<sub>1</sub></td><td>Tune the rise and decay times of GABA release transfer function</td><td>90 s<sup>-1</sup></td></tr><tr><td>z<sub>2</sub></td><td>Tune the rise and decay times of GABA release transfer function</td><td>33 s<sup>-1</sup></td></tr><tr><td>V<sup>ne</sup><sub>G</sub></td><td>Maximal rate of glutamate uptake by neurons</td><td>0.5 µM · s<sup>-1</sup></td></tr><tr><td>V<sup>ae</sup><sub>G</sub></td><td>Maximal rate of glutamate uptake by astrocytes</td><td>4.5 µM · s<sup>-1</sup></td></tr><tr><td>s<sub>g</sub></td><td>Activation threshold of sigmoid glutamate uptakes</td><td>6 µM</td></tr><tr><td>r<sub>g</sub></td><td>Stiffness of sigmoid glutamate uptakes</td><td>0.9 µM<sup>-1</sup></td></tr><tr><td>V<sup>ae</sup><sub>γ</sub></td><td>Maximal rate of astrocyte GABA uptake</td><td>2 µM · s<sup>-1</sup></td></tr><tr><td>K<sup>ae</sup><sub>γ</sub></td><td>Maximal concentration for Hill dynamics of astrocyte GABA uptake</td><td>8 µM</td></tr><tr><td>V<sup>ne</sup><sub>γ</sub></td><td>Maximal rate of neuronal GABA uptake</td><td>5 µM · s<sup>-1</sup></td></tr><tr><td>K<sup>ne</sup><sub>γ</sub></td><td>Maximal concentration for Hill dynamics of neuronal GABA uptake</td><td>24 µM</td></tr><tr><td>V<sup>c</sup><sub>G</sub></td><td>Rate of glutamate degradation by astrocytes</td><td>9 µM · s<sup>-1</sup></td></tr><tr><td>V<sup>c</sup><sub>γ</sub></td><td>Rate of GABA degradation by astrocytes</td><td>9 µM · s<sup>-1</sup></td></tr><tr><td colspan="3"><strong>Neuron excitability modulations by neurotransmitter concentrations (feedbacks)</strong></td></tr><tr><td>v<sub>G</sub></td><td>Excitability threshold of glutamate feedback function</td><td>30 µM</td></tr><tr><td>r<sub>G</sub></td><td>Stiffness of sigmoid glutamate feedback function induced by glutamate</td><td>0.15 µM<sup>-1</sup></td></tr><tr><td>m<sup>P</sup><sub>G</sub></td><td>Maximal coupling gain of glutamate feedback on pyramidal neurons</td><td>2.5 mV</td></tr><tr><td>m<sup>I</sup><sub>G</sub></td><td>Maximal coupling gain of glutamate feedback on interneurons</td><td>1 mV</td></tr><tr><td>v<sub>γ</sub></td><td>Excitability threshold of GABA feedback function</td><td>25 µM<sup>-1</sup></td></tr><tr><td>r<sub>γ</sub></td><td>Stiffness of sigmoid GABA feedback function</td><td>0.12 µM<sup>-1</sup></td></tr><tr><td>m<sub>γ</sub></td><td>Maximal coupling gain of GABA feedback</td><td>1 mV</td></tr></tbody></table>
we detail the modulation of neural intermediary variables for each kind of synapse separately, then we gather these changes to specify the coupling terms reproducing the astrocyte feedback.

In the framework of the local neuron-astrocyte mass model, the astrocyte feedback does not impact the synaptic connections of type $P' \to P$ or $P \to P'$. As a matter of fact, the astrocyte compartment only takes into account neurotransmitters released locally by neurons of populations $P$ and $I$, whereas population $P'$ is nonlocal to population $P$. Hence, extracellular concentrations of neurotransmitters in the vicinity of $P'$ have no impact on the neuronal activity of $P$ and the concentrations in the neighborhood of $P$ and $I$ do not influence the postsynaptic neurons of population $P'$.

A synaptic connection of type $P \to I$ concerns the variable $x_2$. As shown in [1, 38], in the case of glutamate excess in the extracellular space, the postsynaptic neuron is more excitable. Consequently, more neurons are activated in the population $I$. We model this mechanism at the mesoscopic scale by introducing a dependency of population $I$ excitability threshold $v_I$ on the extracellular glutamate concentration and set in Eq. (2c):

$$v_I = v_0 - m_G^I S([Glu]_e, v_G, r_G).$$

While experiments have evidenced impacts of extracellular glutamate concentrations on the excitability of GABAergic neurons [38], we lack experimental data and conclusions on the way GABA acts on the excitability of each subtypes of GABAergic neurons. In particular, observations raise the question of whether all subtypes are subject to mechanisms of GABA-mediated excitability modulation in the same way [39, 40]. Since our model is not subtype-specific for the sake of compactness, we chose not to embed any GABA feedback on interneuron population dynamics.

On the one hand, a synaptic connection of type $I \to P$ is concerned by extracellular concentrations of GABA since it involves GABAergic interneurons. In case of a GABA excess in the extracellular space, the inhibition of the postsynaptic neuron is strengthened (see e.g. [1]), i.e. less neurons are activated in population $P$, which is translated in the NMM by an increase of the threshold of the sigmoidal term in the $x_0$ dynamics. On the other hand, a synaptic connection of type $P \to P$ is impacted by the extracellular concentration of glutamate implying a modulation of variable $x_0$ dynamics as well. In the case of an excess of glutamate in this kind of synapse, the postsynaptic neuron is more excitable [1]. Hence, more neurons are activated in population $P$ which can be reproduced by a decrease in the threshold parameter appearing in (3d). Gathering both modulations impacting the excitability of population $P$, we set in Eq. (3d)

$$v_P = v_0 + m_\gamma S([GABA]_e, v_\gamma, r_\gamma) - m_G^P S([Glu]_e, v_G, r_G).$$

The new neuron-astrocyte mass model embedding the astrocyte feedback is obtained from model (5a)-(5i) by considering the dynamical entries $v_I$ and $v_P$ mentioned above. Accordingly, the sigmoidal functions appearing in Eqs. (5a), (5b), (5c), (5d), and (5g) become

$$\hat{S}(y_1 - y_2, v_P) = \hat{S}(y_1 - y_2, v_0 + m_\gamma S([GABA]_e, v_\gamma, r_\gamma) - m_G^P S([Glu]_e, v_G, r_G)),$$
$$ \hat{S}(C_1 y_0, v_{P'}) = \hat{S}(C_1 y_0, v_0), $$

$$ \hat{S}(C_3 y_0, v_I) = \hat{S}(C_3 y_0, v_0 - m_G^I S([\text{Glu}]_e, v_G, r_G)). $$

## 2.4 Bifurcation-Based Characterization of the Neural Activity: The Case of Noise-Induced Spiking

The behavior of NMMs can be deduced from the bifurcation diagram according to the value $p(t) = p$ considered as a parameter, as it has been performed in [28] on the Jansen-Rit model. In [29], we have classified the types of time-series patterns generated by model (3a)-(3f) and the associated bifurcation structures according to the strengths of the different excitatory inputs to population P. Let us recall the bifurcation diagram underlying the predominant type of generated time series, which we will consider in this article.

Model (3a)-(3f) has the following useful features [29]. First, for a fixed value of parameter $p$, the $y_0$ value of a singular point suffices to have explicit expressions of all the other components. Second, for a given $y_0$ value, there exists a unique value of $p$ such that $y_0$ corresponds to a singular point. In other terms, the set of singular points obtained for the different values of $p$ is a graph over $y_0$. Hence, we can visualize the shape of the singular point locus in the plane $(p, y_0)$: in the case presented here (see top panel of Fig. 4), this curve of singular points is S-shaped. In the following description, for a given bifurcation “X” according to parameter $p$, we denote by $p_X$ the bifurcation value and, if the bifurcation involves a singular point, we denote by $y_X$ the corresponding $y_0$ value.

Two saddle-node bifurcations $SN_1$ and $SN_2$ split the curve of singular points into three branches. We name “lower branch”, “middle branch”, and “upper branch” the sets of singular points, respectively, satisfying $y_0 < y_{SN_1}$, $y_{SN_1} < y_0 < y_{SN_2}$, and $y_0 > y_{SN_2}$. Singular points on the lower branch are stable (blue) and those on the middle branch are unstable (cyan). Singular points on the upper branch are unstable (green) for $p < p_{H1}$ and stable (blue) otherwise. At $p = p_{H1}$ the system undergoes a supercritical Hopf bifurcation $H_1$ giving birth to a stable limit cycle for $p < p_{H1}$ that persists until $p = p_{SN1} = p_{SNIC}$ where it disappears by a saddle-node on invariant circle (SNIC) bifurcation (dashed orange line). The existence of the SNIC bifurcation is essential because it implies the appearance of a large amplitude stable limit cycle with large period. Hence, according to the value of $p$, the system alternates between oscillatory phases (for $p > p_{SNIC}$) and quiescent phases (for $p < p_{SNIC}$). In other terms, $p_{SNIC}$ plays the role of an activation threshold for the neural compartment.

This means that, given a constant input $p(t) = \bar{p}$ such that $\bar{p} < p_{SNIC}$, we expect the NMM to reproduce quiescent neuronal activity. On the other hand, if we consider Gaussian noise as input that is normally distributed with $\bar{p}$ as mean value and large enough variance $\sigma$, the NMM output results instead in an alternation of quiescent phases and isolated LFP spikes (see Fig. 4), consistent with emergence of Noise-Induced Spiking (NIS) [29]. Importantly, NIS frequency depends on the input noise statistics (i.e. mean and variance in the case of Gaussian noise) as reflected by the confidence interval [$\bar{p} - \sigma, \bar{p} + \sigma$] reported by the gray bar above the bifurcation diagram on the left panel of Fig. 4.
<figure><img src="image_8.png" /><figcaption>**Fig. 4** Bifurcation diagram according to *p* (left) and associated LFP and *y*₀ time series (right). Blue curves: stable singular points. Cyan (resp. green) curves: singular points with one (resp. two) eigenvalues with positive real parts. Black curves: *y*₀ extrema along stable limit cycles. Black points (SN₁ and SN₂): saddle-node bifurcations. Red point (H₁): supercritical Hopf bifurcation. Dashed orange line: Saddle-Node on Invariant Circle (SNIC) bifurcation. Horizontal gray bar: confidence interval of the Gaussian variable *p*(*t*) used to generate the time series i.e. [*p*̄ − *σ*, *p*̄ + *σ*] where *p*̄ and *σ* are the mean and the variance of the associated normal distribution, respectively</figcaption></figure>

Such LFP activity, i.e. sparse large amplitude spikes, corresponds to episodic syn-
chronization of the neuron activities among the populations. This pattern of activity
is symptomatic of a strong excitability of the neuronal system that can turn into hy-
perexcitability during pathological crisis. In [29] we have fitted the NMM parameter
values for reproducing interictal and pre-ictal spiking activities by tuning the coupling
gain of the indirect excitatory feedback. We have compared these outputs to exper-
imental data recorded from epileptic mice (Mesial Temporal Lobe Epilepsy mouse
model). It appeared that the NIS behavior is well suited for reproducing the high am-
plitude sharp waves characterizing the activity pattern recorded between pathological
crises. For these fixed parameter values of the NMM, the activity is stable in the sense
that the oscillation frequency does not change much with time. The values of the pa-
rameters associated with the NMM used in this study and given in Table 1 have been
chosen so as to reproduce NIS behavior using the analysis in [29]. In the following,
<figure><img src="image_9.png" /><figcaption>**Fig. 5** Comparison of the effects of a GABA bolus on the feedforward and bilaterally coupled neuron-astrocyte mass models. Bifurcation diagrams according to *p* (top panels) computed for the model without feedback for all *t* (a1) and for the model with feedback at *t* = 0 s (b1) and *t* = 30 s (b2). Time series (bottom panels) corresponding to LFP, glutamate and GABA extracellular concentrations and time variations of *p*<sub>SNIC</sub> generated by the models without feedback (left) and with feedback (right). *The purple lines* on the bifurcation diagrams materialize the fixed value of input *p*. *The dark red arrows* above each time series materialize the time of GABA bolus injection</figcaption></figure>

we study the variations of the activity when the neural dynamics is altered by the
surrounding activity, i.e. the astrocyte feedback.

## 2.5 NMM Dynamics: Example from GABA Injection

To illustrate the impact of the astrocytic feedback on neural excitability and intro-
duce the dynamical tools that we are using in the remainder of the article, we set to
study LFP dynamics generated by our model both with and without astrocytic feed-
back, leveraging on the role of the SNIC bifurcation occurring at *p* = *p*<sub>SNIC</sub> on the
emergence of different regimes of neuronal firing.

With this aim, we simulate injection of GABA (20 AU) into the extracellular space
at t = 0. Assuming p(t) = p ≳ p<sub>SNIC</sub> (Fig. 5, panels (a1) and (b1)). Either in the
presence of astrocytic feedback or lack thereof, the model generates low-frequency
NIS right after injection time. Nonetheless it may be noted that, in the presence of astrocytic feedback, even if extracellular GABA concentration is low—and the corresponding sigmoidal feedback is consequently very weak—this is sufficient to affect neural activity—as reflected by the difference in spike frequency between the two LFP time series.

In the time series generated by the model without feedback (left panels of Fig. 5), the neural activity and the glutamate concentration dynamics remain unchanged (Fig. 5(a)) after the instantaneous increase of [GABA]$_{e}$ at $t = 0$, which mimics GABA injection. In contrast, in the model with astrocyte feedback (right panels of Fig. 5), the strong increase of [GABA]$_{e}$ halts neuronal activity, which, in turn reduces synaptic release of glutamate. Once [GABA]$_{e}$ has become sufficiently low, neural activity re-emerges and both glutamate and GABA concentrations come back to their respective basal values, oscillating with neural activity.

These phenomena can be well accounted by arguments based on bifurcation analysis (top panels of Fig. 5). In the model without astrocytic feedback, neural dynamics is decoupled from astrocytic dynamics. Hence, the bifurcation diagram of the NMM remains unchanged during all simulated time (Fig. 5(a1)). In this fashion, the excitability threshold $p_{\text{SNIC}}$ of the neural compartment is constant in time (Fig. 5, bottom panel). On the other hand, in the system with feedback, $v_{P}$, $v_{P'}$ and $v_{I}$ depend on GABA and glutamate concentrations. Consequently, the bifurcation diagram for the neural compartment is time dependent, and so is the value of $p_{\text{SNIC}}$, which changes along with the astrocytic variables. Thus, right after GABA injection, $p_{\text{SNIC}}$ is larger than the input value $p$ (Fig. 5(b2)), and because $p_{\text{SNIC}}$ plays the role of activation threshold for the neural population, then as long as GABA concentration remains high, the neural population remains quiescent. One may easily see by a simple algebraic calculation (not shown) that the astrocytic compartment admits a unique singular stable point for any neural state. Hence, starting from a high extracellular GABA concentration, this concentration decreases toward the attractive state slowly decreasing $p_{\text{SNIC}}$, as shown in the bottom panel of Fig. 5, case (b). Once, GABA is low enough, $p_{\text{SNIC}}$ becomes smaller than $p$, and the system oscillates again.

This analysis shows how the model with feedback can take into account changes in glutamate and GABA dynamics to modify all the dynamics of the system and illustrates the interest of embedding the astrocyte feedback in such neuron-astrocyte model. Our model allows us to study the effects of variations in glutamate or GABA dynamics on neural activity. In the following, we study the effects of deficiencies in the uptake of neurotransmitters by the astrocytes both on extracellular concentrations and neural activity. For both types of deficiency (GABA and glutamate uptake), we first describe the biologic context and mechanisms and their outcomes, then we provide a mathematical analysis of the underlying dynamical mechanisms to explain the effects that can be expected in a biological system.

## 3 Deficiency of Astrocytic GABA Uptake

We now set to investigate some predictions of our model, starting from the consideration of the scenario of deficiency of astrocytic GABA transporters whose experimental correlate is a low capacity by astrocytes to uptake extracellular GABA [7].
Because a single astrocyte may ensheath hundreds of thousands of synapses [41, 42], the increase of extracellular GABA resulting from slow, defective astrocytic uptake is expected to promote neuronal inhibition and thus reduces synaptic release in line with what is elucidated in the previous Sect. 2.5. Thus, once we consider many defective astrocytes, all neurons in the proximity of these astrocytes are likely to be inhibited [43].

To simulate deficiency of astrocytic GABA uptake in our model, we decrease the maximum (saturated) rate of astrocytic GABA uptake from the extracellular space ($V_Y^{ae}$) as this parameter directly depends on the uptake efficiency of astrocytic GABA transporters. At the neuronal population level we are interested instead in deriving the dependence of $p_{SNIC}$ on astrocytic feedback. With this aim, because glutamate binding by neuronal receptors is independent of cell types, we may set:

$$ 
\begin{aligned}
m_G^I S([Glu]_e, v_G, r_G) &\equiv v_1, \\
m_G^P S([Glu]_e, v_G, r_G) &\equiv \frac{m_G^P}{m_G^I} v_1, \\
m_\gamma S([GABA]_e, v_\gamma, r_\gamma) &\equiv v_2,
\end{aligned}
$$

where $v_1$ and $v_2$ are defined within the boundaries of $m_G^I S([Glu]_e, v_G, r_G)$ and $m_\gamma S([GABA]_e, v_\gamma, r_\gamma)$, that is, $0 \le v_1 \le m_G^I$ and $0 \le v_2 \le m_\gamma$.

With these new notations, the dynamical excitability thresholds $v_P$, $v_{P'}$, and $v_I$ of populations P, P', and I become

$$ 
\begin{aligned}
v_P &= v_0 + v_2 - \frac{m_G^P}{m_G^I} v_1, \\
v_{P'} &= v_0, \\
v_I &= v_0 - v_1.
\end{aligned}
$$

By these new parameters, an increase (decrease) of extracellular GABA or glutamate are represented by an increase (decrease) of $v_2$ or $v_1$, respectively. Accordingly, because a deficiency of astrocytic GABA uptake increases extracellular GABA concentration, then we aim at characterizing the dependency of $p_{SNIC}$ on the value of $v_2$. In particular we formalize our first model prediction in Proposition 3.1 and exemplify it by numerical continuation of the SNIC bifurcation for varying $v_2$ in Fig. 6 (left panel). In carrying out our following analysis, we assume $v_1$ constant, leaving justification of such assumption a posteriori at the end of the section.

**Proposition 3.1** The threshold for neuronal activation, $p_{SNIC}$, linearly increases with the GABA-induced modulation of neural excitability $v_2$.

*Proof* The set of the system singular points obtained for the different values of parameter $p$ can be explicitly expressed according to $y_0$, $v_1$, and $v_2$ all other parameters being fixed. Using system (3a)-(3f), a singular point satisfies $y_3 = y_4 = y_5 = 0$. Under these conditions, annihilation of the r.h.s. of Eq. (3f) allows expressing $y_2$ with
<figure><img src="image_10.png" /><figcaption>**Fig. 6** The dependency of the SNIC value on $v_2$ explains the impact of an alteration of the GABA astrocytic uptake on the neural activity. Variation of $p_{SNIC}$ value according to $v_2$ (left). Time series corresponding to LFP, [GABA]$_e$, [Glu]$_e$ and $v_2 = m_γS([GABA]_e, v_γ, r_γ)$ (right from top to bottom) for $p(t)$ a Gaussian variable. At $t = 40$ s, the GABA astrocytic uptake is artificially altered by setting $V_γ^{ae} = 0$ for any subsequent time. The gray time window highlights the transient toward the new behavior</figcaption></figure>

respect to $y_0$. Then annihilation of the r.h.s. of (3d) provides $y_1$ as function of $y_0$. Finally, by annihilation of Eq. (3e), we find that the $y_0$ components of the singular points for given values of $p, v_1$ and $v_2$ are characterized as solutions of

$$p = f(y_0, v_1, v_2), \qquad (6)$$

where

$$f(y_0, v_1, v_2) = \frac{a}{A} \left( v_0 - \frac{m_G^P}{m_G^I} v_1 + v_2 \right) - \frac{a}{Ar} \ln \left( \frac{2Ae_0 - ay_0}{ay_0} \right) - \frac{aG}{A} y_0 - C_2 \hat{S}(C_1 y_0, v_0) + \frac{aB}{bA} C_4 \hat{S}(C_3 y_0, v_0 - v_1). \qquad (7)$$

All the other components of a given singular point can be explicitly expressed as functions of the $y_0$ component. The calculation is straightforward from the annihilation of the vector field (5a)-(5i) as explained above. Since these components are not needed for the following analysis and are only minor calculus details, we do not specify these heavy expressions. We rewrite Eq. (7) as follows:

$$f(y_0, v_1, v_2) = \frac{a}{A} v_2 + q(y_0, v_1). \qquad (8)$$

Obviously the two saddle-node bifurcation values $p_{SN_1}$ and $p_{SN_2}$ are relative extrema of $f(y_0, v_1, v_2)$. In particular $p_{SN_1} = p_{SNIC}$ coincide with the local maximum of $f(y_0, v_1, v_2)$ and, as such, is of the solution of

$$p = f(y_0, v_1, v_2), \qquad (9a)$$

$$\frac{\partial f}{\partial y_0}(y_0, v_1, v_2) = 0, \qquad (9b)$$

$$\frac{\partial^2 f}{\partial y_0^2}(y_0, v_1, v_2) \le 0. \qquad (9c)$$
Since $\frac{\partial f}{\partial y_0}(y_0, v_1, v_2)$ is independent on $v_2$, so is $y_{SNIC}$ and it can be considered as a parameter in Eq. (9a). From (8) and (9a)-(9c) we obtain the following expression for $p_{SNIC}$:

$$p_{SNIC} = \frac{a}{A}v_2 + q(y_{SNIC}, v_1). \quad \square$$

Let us consider the model generating an oscillatory output with a fixed value of $p$ ($p_{SNIC} < p$). If the extracellular concentration of GABA increases (e.g. by an injection of GABA as in Fig. 5), the value of $v_2$ increases and Proposition 3.1 asserts that the value of $p_{SNIC}$ also increases. As already explained, the closest $p_{SNIC}$ is to $p$, with $p_{SNIC} < p$, the largest is the limit cycle period, thus the oscillation frequency of the outputs decreases. If $p_{SNIC}$ increases enough such that $p_{SNIC} > p$, the stable limit cycle of the system disappears, and the neural compartment becomes quiescent.

In the case of deficiency of GABA uptake, the extracellular concentration of GABA increases, and we can use Proposition 3.1 to explain the following effects. For that purpose, we use the following *in silico* protocol: we initialize the neuron-astrocyte model in an oscillatory phase with a low oscillation frequency and consider $p(t)$ a Gaussian input. At $t = 40$ s, we turn off the GABA astrocyte uptake by setting $V_y^{ae} = 0$ (Fig. 6). The result is an increase in GABA extracellular concentration implying an increase in $p_{SNIC}$. As $p_{SNIC}$ increases, the probability for $p(t)$ to overcome $p_{SNIC}$ along the associated Brownian motion decreases, and also does the oscillation frequency (Fig. 6). Consequently, we observe a decrease in the oscillation frequency after $t > 40$ s. In the time series, the oscillation frequency decreases gradually during a transient ($40 s < t < 60$ s) until reaching its minimum. This can be explained by the slow increase of GABA extracellular concentration that reaches its new baseline at $t = 60$ s.

To summarize, our hitherto analysis reveals that a deficiency of astrocytic GABA uptake causes a decrease of neural firing, leaving extracellular glutamate concentration close to its baseline. This emphasizes that, under such conditions, changes in glutamate-induced modulation of neuronal excitability ($v_1$) are overall negligible insofar as neural activity is sufficiently inhibited by high extracellular GABA.

## 4 Deficiency of Astrocytic Glutamate Uptake

We now focus on the impact of deficient astrocytic glutamate uptake on neuronal activation. This scenario results in an increase of extracellular glutamate which in principle is expected to make neurons more excitable and thus more prone to fire. However, such increase in excitability may be (at least partly) counteracted by the resulting boost in interneurons firing which triggers synaptic release of GABA, overall increasing extracellular levels of this species and promoting, in turn, neural inhibition. This suggests that the possible balance between glutamatergic neuronal excitation and GABAergic inhibition may lead to nontrivial dynamics regimes for neural activity. Considering the rationale outlined in the previous section, we now set out to characterize the changes in the SNIC bifurcation value $p_{SNIC}$ (i.e. the threshold for neural activation) with $v_1$ value. The latter variable, introduced in the beginning of Sect. 3, quantifies the instantaneous impacts of glutamate-induced modulations on
both interneuron and pyramidal neuron excitabilities. Hence, it depends directly on the astrocytic glutamate uptake. We therefore analyze the various consequence on the neural activity of a deficiency of astrocytic glutamate uptake simulated by choosing small values for the maximum uptake rate $V_G^{ae}$.

We note in this regard, however, that, for certain parameter sets of the whole model, the SNIC bifurcation may be lost as $v_1$ changes. Due to the high dimensionality of our parameter space, it is difficult to estimate the region of existence of the SNIC for any value of $v_1$, yet by the arguments exposed in [29] we may safely state that this region is large enough to be always found for our choice of parameters (see Table 1). Accordingly, in the following, we assume that $p_{SNIC}$ always exists and the associated saddle-node bifurcation is not degenerate for any $v_1$ such as $0 \le v_1 \le m_G^I$, that is the maximum interval of values taken by $v_1 = m_G^I S([Glu]_e, v_G, r_G)$.

We recall that $p_{SNIC}$ can be written as follows:

$$p_{SNIC} = f(y_{SNIC}, v_1, v_2), \quad (10)$$

where $f$ is given by (7). Since we consider $v_2$ fixed we introduce the function

$$g(y_0, v_1) \equiv f(y_0, v_1, v_2)|_{v_2=\text{const}}.$$

As explained above, for each $v_1$, there exists a unique bifurcation value $p_{SNIC}$ occurring at the non-hyperbolic (saddle-node) singular point characterized by $y_{SNIC}$, which is defined by

$$\frac{\partial g}{\partial y_0}(y_{SNIC}, v_1) = 0,$$

$$\frac{\partial^2 g}{\partial y_0^2}(y_{SNIC}, v_1) < 0.$$

This value satisfies $p_{SNIC} = g(y_{SNIC}, v_1)$. We cannot find the explicit expressions of $y_{SNIC}(v_1)$ and $p_{SNIC}(v_1)$. Thus, for characterizing the variations of $p_{SNIC}$ with $v_1$, we take advantage of the implicit definitions above and focus on localizing the extrema of $p_{SNIC}(v_1)$.

**Proposition 4.1** Assume that, for all $v_1$ such that $0 \le v_1 \le m_G^I$, $p_{SNIC}$ exists and the associated saddle-node bifurcation is not degenerate. Note $\chi = \frac{Be_0rC_4}{2b}$. Then:

1. if $\frac{m_G^P}{m_G^I} \ge \chi$, $p_{SNIC}(v_1)$ has no local extremum,

2. if $0 < \frac{m_G^P}{m_G^I} < \chi$, $p_{SNIC}(v_1)$ may admit two local extrema: a minimum at $v_1^*$ and a maximum at $v_1^{**}$. If both exist, then $v_1^* < v_1^{**}$.

*Proof* Let us search for local extrema of the function $p_{SNIC}(v_1)$, which is implicitly defined by (9a)-(9c). Hence, we are interested in the following constrained optimization problems:

$$\max \left\{ g(y_0, v_1) \mid \frac{\partial g}{\partial y_0}(y_0, v_1) = 0 \right\} \quad \left( \min \left\{ g(y_0, v_1) \mid \frac{\partial g}{\partial y_0}(y_0, v_1) = 0 \right\} \right). \quad (11)$$
We introduce the associated Lagrangian function

$$L(y_0, v_1, \lambda) = g(y_0, v_1) - \lambda \frac{\partial g}{\partial y_0}(y_0, v_1).$$

The necessary condition for the existence of an extremum of g under the constraint $\frac{\partial g}{\partial y_0} = 0$ is

$$\vec{\nabla} L(y_0, v_1, \lambda) = 0,$$

that is,

$$\frac{\partial g}{\partial y_0}(y_0, v_1) - \lambda \frac{\partial^2 g}{\partial y_0^2}(y_0, v_1) = 0, \qquad (12a)$$

$$\frac{\partial g}{\partial v_1}(y_0, v_1) - \lambda \frac{\partial^2 g}{\partial v_1 \partial y_0}(y_0, v_1) = 0, \qquad (12b)$$

$$\frac{\partial g}{\partial y_0}(y_0, v_1) = 0. \qquad (12c)$$

By assumption, the saddle-node bifurcation associated with the SNIC bifurcation is not degenerate, i.e. every solution of (11) such as $0 \le v_1 \le m_G^l$ satisfies $\frac{\partial^2 g}{\partial y_0^2}(y_0, v_1) \ne 0$. Thus, system (12a)-(12c) reads

$$\lambda = 0, \qquad (13a)$$

$$\frac{\partial g}{\partial v_1}(y_0, v_1) = 0, \qquad (13b)$$

$$\frac{\partial g}{\partial y_0}(y_0, v_1) = 0. \qquad (13c)$$

Therefore, if the problem under constraint admits an extremum, at this extremum it is $\frac{\partial g}{\partial v_1} = 0$. Following the assumption that a SNIC bifurcation occurs for any value of $v_1$ between 0 and $m_G^l$, then Eq. (13c) admits a solution for any $v_1$ as well. Hence, if our constrained problem admits an extremum, this corresponds to the SNIC bifurcation occurring at $(y_0, v_1)$ such that

$$\frac{\partial g}{\partial y_0}(y_0, v_1) = 0.$$

From (7), we obtain

$$\frac{\partial g}{\partial v_1}(y_0, v_1) = -\frac{a}{A} \left( \frac{m_G^P}{m_G^l} + \frac{B}{b} C_4 \frac{\partial \hat{S}}{\partial v} (C_3 y_0, v_0 - v_1) \right).$$

Using the facts that, for any fixed value of $y_0$, function $v_1 \to \frac{\partial g}{\partial v_1}(y_0, v_1)$ is bell-shaped and its maximum value does not depend on $y_0$ (see Fig. 7), one finds that the
<figure><img src="image_11.png" /><figcaption>**Fig. 7** Interval of $m_G^P/m_G^I$ values corresponding to a realized local minimum of neural excitability threshold. Graphic representation of the function $h$ defined by (19) and interval $[I_1, I_2]$ for which $v_1^* \in [0, m_G^I]$</figcaption></figure>

function $\frac{\partial g}{\partial v_1}(y_0, v_1)$ vanishes in $v_1$ if

$$0 < \frac{m_G^P}{m_G^I} < \chi. \qquad (14)$$

If $\frac{m_G^P}{m_G^I} \ge \chi$, the function $\frac{\partial g}{\partial v_1}(y_0, v_1)$ admits no zero, which proves the first item of Proposition 4.1.

Now, we assume that condition (14) is fulfilled and we search the values of $v_1$ satisfying $\frac{\partial g}{\partial v_1}(y_0, v_1) = 0$, i.e.

$$\frac{m_G^P}{m_G^I} + \frac{B}{b} C_4 \frac{\partial \hat{S}}{\partial v}(C_3 y_0, v_0 - v_1) = 0.$$

Noting $X = e^{r(v_0-v_1-C_3y_0)} = (S(C_3y_0, v_0 - v_1, r))^{-1} - 1$, this equation reads

$$\frac{m_G^P}{m_G^I} X^2 + 2 \left( \frac{m_G^P}{m_G^I} - 2\chi \right) X + \frac{m_G^P}{m_G^I} = 0. \qquad (15)$$

Setting

$$V_{\pm} = 2\chi \frac{m_G^I}{m_G^P} - 1 \pm 2\sqrt{\chi \frac{m_G^I}{m_G^P} \left( \chi \frac{m_G^I}{m_G^P} - 1 \right)}. \qquad (16)$$

we obtain the two solutions $v_1^* < v_1^{**}$ of $\frac{\partial g}{\partial v_1}(y_0, v_1) = 0$:

$$v_1^* = v_0 - C_3 y_0 - \frac{1}{r} \ln(V_+), \qquad (17)$$
$$v_1^{**} = v_0 - C_3 y_0 - \frac{1}{r} \ln(V_-). \qquad (18)$$

Note that $v_1^*$ (resp. $v_1^{**}$) corresponds to the extremum when the saddle-node $SN_1$ (resp. $SN_2$) crosses the fold of the surface $g(y_0, v_1) = p$. We consider $v_1 = v_1^*$ and we note $y_0^*$ the value of $y_0$ corresponding to the SNIC connection for this value of $v_1$, i.e. the solution of

$$\frac{\partial g}{\partial y_0}(y_0, v_1^*) = 0.$$

$$\frac{\partial^2 g}{\partial y_0^2}(y_0, v_1^*) < 0.$$

To prove that $p_{SNIC}$ reaches a local minimum at $v_1 = v_1^*$, we introduce the bordered Hessian matrix $\bar{H}$ associated with the Lagrangian function at its singular point $(y_0, v_1, \lambda) = (y_0^*, v_1^*, 0)$ (solution of system (13a)-(13c)):

$$\bar{H}(y_0^*, v_1^*, 0) = \begin{pmatrix} 0 & \frac{\partial^2 g}{\partial y_0^2} & \frac{\partial^2 g}{\partial v_1 \partial y_0} \\ \frac{\partial^2 g}{\partial y_0^2} & \frac{\partial^2 \bar{L}}{\partial y_0^2} & \frac{\partial^2 \bar{L}}{\partial v_1 \partial y_0} \\ \frac{\partial^2 g}{\partial v_1 \partial y_0} & \frac{\partial^2 \bar{L}}{\partial v_1 \partial y_0} & \frac{\partial^2 \bar{L}}{\partial v_1^2} \end{pmatrix}_{|(y_0^*, v_1^*, 0)}$$

$$= \begin{pmatrix} 0 & \frac{\partial^2 g}{\partial y_0^2} & \frac{\partial^2 g}{\partial v_1 \partial y_0} \\ \frac{\partial^2 g}{\partial y_0^2} & \frac{\partial^2 g}{\partial y_0^2} & \frac{\partial^2 g}{\partial v_1 \partial y_0} \\ \frac{\partial^2 g}{\partial v_1 \partial y_0} & \frac{\partial^2 g}{\partial v_1 \partial y_0} & \frac{\partial^2 g}{\partial v_1^2} \end{pmatrix}_{|(y_0^*, v_1^*, 0)}.$$

The determinant of $\bar{H}(y_0^*, v_1^*, 0)$ is given by

$$\det \bar{H}(y_0^*, v_1^*, 0) = -\frac{\partial^2 g}{\partial y_0^2}(y_0^*, v_1^*) \left[ \frac{\partial^2 g}{\partial y_0^2}(y_0^*, v_1^*) \frac{\partial^2 g}{\partial v_1^2}(y_0^*, v_1^*) - \left( \frac{\partial^2 g}{\partial v_1 \partial y_0}(y_0^*, v_1^*) \right)^2 \right].$$

On one hand, the saddle-node associated with the SNIC bifurcation is not degenerate and is a local maximum of $g(y_0, v_1)$, thus $\frac{\partial^2 g}{\partial y_0^2}(y_0^*, v_1^*) < 0$. On the other hand, for any $y_0, v_1 \to \frac{\partial g}{\partial v_1}(y_0, v_1)$ is increasing at $(y_0, v_1^*)$ (see Fig. 7), thus $\frac{\partial^2 g}{\partial v_1^2}(y_0^*, v_1^*) > 0$. Finally

$$\det \bar{H}(y_0^*, v_1^*, 0) < 0$$

and $(y_0^*, v_1^*)$ corresponds to a local minimum of $p_{SNIC}$. A similar argument proves that $(y_0^{**}, v_1^{**})$ corresponds to a local maximum of $p_{SNIC}$ (where $y_0^{**}$ is the $y_0$ value corresponding to $SN_2$ bifurcation for $v_1 = v_1^{**}$). $\square$
The above proposition can be interpreted as a necessary condition for having a change in the sense of variations of $p_{SNIC}$ when $v_1$ varies in $[0, m_G^I]$. We now derive a necessary and sufficient condition so that $v_1^*$ actually lies in $[0, m_G^I]$. Since $v_1^*$ satisfies $\frac{\partial g}{\partial v_1}(y_0^*, v_1^*) = 0$, one obtains, from Eq. (15),

$$ \frac{m_G^P}{m_G^I} = 4\chi(1 - S(C_3y_0^*, v_0 - v_1^*, r))S(C_3y_0^*, v_0 - v_1^*, r) = 4\chi h(v_1^*). \quad (19) $$

For any $y_0$, function $h: v_1 \mapsto (1 - S(C_3y_0, v_0 - v_1, r))S(C_3y_0, v_0 - v_1, r)$ is strictly increasing over $[0, m_G^I]$, and $0 \le v_1^* \le m_G^I$ if and only if $I_1 \le \frac{m_G^P}{m_G^I} \le I_2$ where

$$ I_1 = 4\chi(1 - S(C_3y_0^*, v_0, r))S(C_3y_0^*, v_0, r), \quad (20a) $$

$$ I_2 = 4\chi(1 - S(C_3y_0^*, v_0 - m_G^I, r))S(C_3y_0^*, v_0 - m_G^I, r). \quad (20b) $$

Figure 7 gives an illustration of the function $h$ as well as how $I_1$ and $I_2$ are built.

In conclusion, for a fixed value of $v_2$, $p_{SNIC}$ reaches a local minimum at a value $v_1^*$ such that $0 \le v_1^* \le m_G^I$ if and only if

$$ \frac{m_G^P}{m_G^I} \in (0, \chi) \cap [I_1, I_2]. $$

Moreover, in Sect. 3 about the astrocyte GABA uptake deficiency, we proved that, for a fixed value of $v_1$, $p_{SNIC}$ is linear and increasing with $v_2$. Both results allow us to predict that there exist three shapes of $p_{SNIC}(v_1, v_2)$ according to the value of $\frac{m_G^P}{m_G^I}$.

(a) If $\frac{m_G^P}{m_G^I} < I_1$ then $v_1^* < 0$ and $p_{SNIC}$ strictly increases with $v_1$ and $v_2$.

(b) If $\frac{m_G^P}{m_G^I} > I_2$, then $v_1^* > m_G^I$ and $p_{SNIC}$ strictly decreases when $v_1$ increases (for $v_2$ fixed) and strictly increases with $v_2$ (for $v_1$ fixed).

(c) If $I_1 \le \frac{m_G^P}{m_G^I} \le I_2$, then $0 \le v_1^* \le m_G^I$ and $p_{SNIC}$ decreases when $v_1$ increases in $[0, v_1^*]$ and increases with $v_1 > v_1^*$ (for $v_2$ fixed).

In the following we illustrate the three qualitative types of neural activity resulting from an astrocyte deficiency to capture glutamate using the following values:

$$ (a) \quad \frac{m_G^P}{m_G^I} = 1.7 < I_1, $$

$$ (b) \quad \frac{m_G^P}{m_G^I} = 3.2 > I_2, $$

$$ (c) \quad \frac{m_G^P}{m_G^I} = 2.43 \in [I_1, I_2]. $$
For each case, we provide simulations representing the value of $p_{\text{SNIC}}$ in ($v_1, v_2$) space and time series generated by the model when the glutamate astrocytic uptake is altered (Fig. 8).

For $\frac{m_G^P}{m_I^P} = 1.7 < I_1$ (panel (a) of Fig. 8), $v_1^*$ is negative and thus $p_{\text{SNIC}}$ increases with $v_1$. Hence in this scenario, when $v_1$ increases as a consequence of the astrocytic glutamate uptake deficiency, we expect to observe a decrease in the oscillation frequency in the neural activity. This is possible because, although a reduction of astrocytic glutamate uptake is paired with a strong decrease of neural activity, it limits synaptic glutamate release allowing deficient astrocytic uptake to still stabilize extracellular glutamate levels and thus neural activity.

Panel (b) of Fig. 8 is obtained with $\frac{m_G^P}{m_I^P} = 3.2 > I_2$. Since $v_1^* > 1$, $p_{\text{SNIC}}$ decreases as $v_1$ increases. Thereby, an astrocyte glutamate uptake deficiency triggers an increase of the $v_1$ value, and we observe an increase in the oscillation frequencies in the neural time series. Since the astrocytic uptake is reduced, glutamate accumulates in the extracellular space and the corresponding concentration baseline increases drastically.

Finally, the intermediate scenario where $\frac{m_G^P}{m_I^P} = 2.43 \in [I_1, I_2]$ is shown in panel (c) of Fig. 8. In this case $p_{\text{SNIC}}$ decreases for $v_1 < v_1^*$ and increases otherwise. Therefore, in the presence of deficient astrocytic glutamate uptake, we observe only a transient increase of the frequency of LFP spikes. This is the scenario presciently mentioned at the beginning of the section where the effect of increased extracellular glutamate levels on neuronal excitability are counterbalanced by increased extracellular GABA levels resulting from increased interneurons firing and accounts for the delay observed in our simulations in reaching steady frequency of LFP oscillations.

This scenario is physiologically relevant insofar as it is conceivable that an excess of extracellular glutamate concentration is regulated after a delay, triggering a decrease of neural activity after the initial increase. Moreover, the frequency after the regulation delay can be greater or lower than the initial one, depending on the value of the ratio $\frac{m_G^P}{m_I^P}$. Note that this value can be tuned to obtain $v_1^*$ small enough and $p_{\text{SNIC}}$ large enough so that the frequency after regulation is equal or lower than the one before uptake deficiency. This property offers the possibility of fitting the model outputs to experimental data and allows us to propose hypotheses about physiological and pathological mechanisms.

# 5 Discussion and Conclusions

We have introduced a novel neuron-astrocyte mass model, which combines two previously studied models [29, 30] by bidirectional neuron-astrocyte coupling based on astrocytic uptake-mediated modulations of extracellular glutamate and GABA concentrations and their effect on neural activity in agreement with experimental observations [10, 37, 44, 45].

Using analytical arguments of bifurcation theory and contained optimization, we have characterized different types of change in the neural activity behavior under the
<figure><img src="image_12.png" /><figcaption>**Fig. 8** Alteration of astrocytic glutamate uptake: (a) lessening the excitability, (b) resulting in sustained hyperexcitability, (c) resulting in transient hyperexcitability. In each case, *the colormap on the left* displays the value of *p*<sub>SNIC</sub> in (v<sub>1</sub>, v<sub>2</sub>) plane, and the time series on *the right* correspond to LFP, [GABA]<sub>e</sub>, [Glu]<sub>e</sub> and v<sub>1</sub> = *m*<sub>G</sub><sup>I</sup>S([Glu]<sub>e</sub>, v<sub>G</sub>, r<sub>G</sub>). *The black curves on the colormap* are the trace of (*m*<sub>G</sub><sup>I</sup>S([Glu]<sub>e</sub>, v<sub>G</sub>, r<sub>G</sub>), *m*<sub>γ</sub>S([GABA]<sub>e</sub>, v<sub>γ</sub>, r<sub>γ</sub>)) along the associated orbits of the model. At *t* = 20 s, we alter the glutamate astrocytic uptake by setting V<sub>G</sub><sup>ac</sup> = 0. The three cases are obtained with: (a) *m*<sub>G</sub><sup>P</sup>/*m*<sub>G</sub><sup>I</sup> = 1.7, (b) *m*<sub>G</sub><sup>P</sup>/*m*<sub>G</sub><sup>I</sup> = 3.2, (c) *m*<sub>G</sub><sup>P</sup>/*m*<sub>G</sub><sup>I</sup> = 2.43. All other parameters are the same in the three cases and given in Table 1</figcaption></figure>
impact of the neurotransmitters-mediated modulation of neuron excitability. Based on the interpretation of the aggregated—yet biophysically significant parameters of our model—we have reproduced *in silico* both glutamate and GABA uptake deficiency in astrocytes and illustrated by numerical simulations the different types of change in the neural activity resulting from such deficiency.

Our study leads to several predictions. First is the observation that deficiency of GABA uptake by astrocytes increases the threshold for neuronal activation in a linear fashion. In other words, neurons tend to fire more the slower astrocyte GABA uptake is, or equivalently, the larger extracellular GABA concentration is, which is consistent with experimental findings [46, 47]. The second prediction comes instead from the analysis of the neuronal response in the presence of deficient astrocytic glutamate uptake. In this case, neural activity may either be reduced or enhanced or, alternatively, may display a transient of high activity before stabilizing around a new regime whose firing frequency is close to the one measured in the absence of astrocyte deficiency. A prominent feature of our model is its mathematical tractability, despite the fact that, differently from the sole neural mass model in [29], an explicit formula for the SNIC bifurcation underpinning the rise of neural activity in the presence of bilateral neuron-astrocyte coupling cannot be obtained. It is nonetheless possible to recast the characterization of variations of neuronal excitability upon astrocytic uptake as an optimization problem under an equality constraint resulting from the implicit characterization of the SNIC bifurcation. In particular, we identify by the quantity $m_G^P/m_G^I$, which denotes the ratio between maximally affecting concentrations of extracellular glutamate impacting on excitability of pyramidal neurons and interneurons, respectively, the biophysical conditions for the occurrence of different neural activity regimes caused by deficient astrocytic glutamate uptake. In this regard, our analysis suggests that a neural population may counteract the effect of excess extracellular glutamate resulting from deficient astrocytic uptake only if the value of $\frac{m_G^P}{m_G^I}$ lies in a well-defined interval $[I_1, I_2]$ which depends on inherent biophysical and cellular properties of the population under consideration. Outside this interval, for values of $\frac{m_G^P}{m_G^I} < I_1$, the frequency of neural activity decreases and may enter quiescence while for values of $\frac{m_G^P}{m_G^I} > I_2$, neural firing drastically and persistently increases leading to epileptiform activity in agreement with experiment [48]. On the other hand, for $I_1 < \frac{m_G^P}{m_G^I} < I_2$ neural activity may recover to close-to-baseline levels after a high frequency transient, regardless of large extracellular concentrations of glutamate and GABA. We note that while both $I_1$ and $I_2$ depend on the strength of inhibition of pyramidal neurons by interneurons only $I_2$ further depends on $m_G^I$. Hence, one might interpret the latter dependence as an estimator of the balance between glutamatergic astrocytic feedbacks on pyramidal neurons and correlated interneuron hyperexcitability that the system should seek to recover from astrocytic uptake deficits and avoid hyperexcitability.

A further interesting prediction of our model is that, when $\frac{m_G^P}{m_G^I} > I_1$, an increase of extracellular glutamate due to deficits in astrocytic uptake develops in parallel with an increase of extracellular GABA which may even be larger than the one expected
in an alternative scenario of faulty GABAergic uptake by astrocytes. Dynamically, this GABA increase is necessary for the system to reach a steady low-frequency activity but it may not be sufficient to properly counterbalance the positive feedback on neuronal activation exerted by the faulty astrocytic glutamate uptake. This suggests GABA increases induced by excess extracellular glutamate and the related neuronal activity as a putative biophysical correlate to distinguish between different mechanisms for rise of neural hyperexcitability. In light of the current technological limits whereby monitoring of extracellular neurotransmitters dynamics in humans is not available, these results could have important translational implications, helping to design either new experimental protocols to characterize and identify epileptogenic foci as well as to define more accurate treatments of neuropathogenic states associated with neuronal hyperexcitability or lack thereof.

## Competing Interests

The authors declare that they have no competing interests.

## Authors' Contributions

All authors designed the study and built the model; AV designed the codes for the simulations; AG and AV fitted the parameters and performed the simulations; AG and AV produced the figures; AV stated and proved the results of the propositions; all authors participated in writing the manuscript. All authors read and approved the final manuscript.

**Acknowledgements** This work was performed within the Labex SMART (ANR-11-LABX-65) supported by French state funds managed by the ANR within the Investissements d'Avenir programme under reference ANR-11-IDEX-0004-02. We thank H. Berry for helpful comments on this work. We also thank S. Blanchard and F. Wendling for allowing us to use their astrocyte activity model.

## References

1. Araque A, Parpura V, Sanzgiri RP, Haydon PG. Glutamate-dependent astrocyte modulation of synaptic transmission between cultured hippocampal neurons. Eur J Neurosci. 1998;10(6):2129-42.

2. Parpura V, Haydon PG. Physiological astrocytic calcium levels stimulate glutamate release to modulate adjacent neurons. Proc Natl Acad Sci USA. 2000;97(15):8629-34.

3. Bazargani N, Attwell D. Astrocyte calcium signaling: the third wave. Nat Neurosci. 2016;19(2):182-9.

4. De Pittà M, Brunel N, Volterra A. Astrocytes: orchestrating synaptic plasticity? Neuroscience. 2016;323:43-61.

5. Arizono M, Bannai H, Nakamura K, Niwa F, Enomoto M, Matsu-ura T, Miyamoto A, Sherwood MW, Nakamura T, Mikoshiba K. Receptor-selective diffusion barrier enhances sensitivity of astrocytic processes to metabotropic glutamate receptor stimulation. Sci Signal. 2012;5(218):ra27.

6. Sarup A, Bonde C, Zimmer J, Danbolt N, Schousboe A, Gegelashvili G. Expression and dynamic regulation of glutamate transporters GLT1 and GLAST in primary astroglial and hippocampal slice cultures: effects of oxygen/glucose deprivation (OGD) and GDNF. In: GLIA. New York: Wiley; 2002. p. 65-6.

7. Schousboe A, Madsen KK, White HS. GABA transport inhibitors and seizure protection: the past and future. Future Med Chem. 2011;3(2):183-7.

8. Wang D, Bordey A. The astrocyte odyssey. Prog Neurobiol. 2008;86(4):342-67.
9. Murphy-Royal C, Dupuis JP, Varela JA, Panatier A, Pinson B, Baufreton J, Groc L, Oliet SH. Surface diffusion of astrocytic glutamate transporters shapes synaptic transmission. Nat Neurosci. 2015;18(2):219–26.

10. Huang YH, Sinha SR, Tanaka K, Rothstein JD, Bergles DE. Astrocyte glutamate transporters regulate metabotropic glutamate receptor-mediated excitation of hippocampal interneurons. J Neurosci. 2004;24(19):4551–9.

11. Tzingounis AV, Wadiche JI. Glutamate transporters: confining runaway excitation by shaping synaptic transmission. Nat Rev Neurosci. 2007;8(12):935–47.

12. Nadkarni S, Jung P. Modeling synaptic transmission of the tripartite synapse. Phys Biol. 2007;4(1):1–9.

13. Postnov DE, Koreshkov RN, Brazhe NA, Brazhe AR, Sosnovtseva OV. Dynamical patterns of calcium signaling in a functional model of neuron-astrocyte networks. J Biol Phys. 2009;35(4):425–45.

14. Postnov DE, Ryazanova LS, Sosnovtseva OV. Functional modeling of neural-glia interaction. Biosystems. 2007;89(1–3):84–91.

15. Volman V, Bazhenov M, Sejnowski TJ. Computational models of neuron-astrocyte interaction in epilepsy. Front Comput Neurosci. 2012;6:58.

16. Volman V, Ben-Jacob E, Levine H. The astrocyte as a gatekeeper of synaptic information transfer. Neural Comput. 2007;19(2):303–26.

17. De Pittà M, Volman V, Berry H, Ben-Jacob E. A tale of two stories: astrocyte regulation of synaptic depression and facilitation. PLoS Comput Biol. 2011;7(12):1002293.

18. Gruetter R, Seaquist ER, Ugurbil K. A mathematical model of compartmentalized neurotransmitter metabolism in the human brain. Am J Physiol: Endocrinol Metab. 2001;281(1):100–12.

19. Silchenko AN, Tass PA. Computational modeling of paroxysmal depolarization shifts in neurons induced by the glutamate release from astrocytes. Biol Cybern. 2008;98(1):61–74.

20. Aubert A, Costalat R. A model of the coupling between brain electrical activity, metabolism, and hemodynamics: application to the interpretation of functional neuroimaging. NeuroImage. 2002;17(3):1162–81.

21. Aubert A, Costalat R, Magistretti PJ, Pellerin L. Brain lactate kinetics: modeling evidence for neuronal lactate uptake upon activation. Proc Natl Acad Sci USA. 2005;102(45):16448–53.

22. Aubert A, Costalat R. Interaction between astrocytes and neurons studied using a mathematical model of compartmentalized energy metabolism. J Cereb Blood Flow Metab. 2005;25(11):1476–90.

23. Savin C, Triesch J, Meyer-Hermann M. Epileptogenesis due to glia-mediated synaptic scaling. J R Soc Interface. 2009;6(37):655–68.

24. Volman V, Bazhenov M, Sejnowski TJ. Divide and conquer: functional segregation of synaptic inputs by astrocytic microdomains could alleviate paroxysmal activity following brain trauma. PLoS Comput Biol. 2013;9(1):1002856.

25. Nadkarni S, Jung P. Synaptic inhibition and pathologic hyperexcitability through enhanced neuron-astrocyte interaction: a modeling study. J Integr Neurosci. 2005;4(2):207–26.

26. Nunez PL, Silberstein RB. On the relationship of synaptic activity to macroscopic measurements: does co-registration of EEG with fMRI make sense? Brain Topogr. 2000;13(2):79–96.

27. Keener JP, Sneyd J. Mathematical physiology. vol. 1. Berlin: Springer; 1998.

28. Touboul J, Wendling F, Chauvel P, Faugeras O. Neural mass activity, bifurcations, and epilepsy. Neural Comput. 2011;23(12):3232–86.

29. Garnier A, Vidal A, Huneau C, Benali H. A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics. Neural Comput. 2015;27:329–64.

30. Blanchard S, Saillet S, Ivanov A, Benquet P, Benar C, Pelegrini-Issac M, Benali H, Wendling F. A new computational model for neuro-glio-vascular coupling: astrocyte activation can explain cerebral blood flow nonlinear response to interictal events. PLoS ONE. 2016. doi:10.1371/journal.pone.0147292.

31. Van Rotterdam A, Lopes da Silva FH, van den Ende J, Viergever M, Hermans A. A model of the spatial-temporal characteristics of the alpha rhythm. Bull Math Biol. 1982;44(2):283–305.

32. Freeman W. Mass action in the nervous system. New York: Academic Press; 1975.

33. Jansen BH, Zouridakis G, Brandt ME. A neurophysiologically-based mathematical model of flash visual evoked potentials. Biol Cybern. 1993;68(3):275–83.

34. Braitenberg V, Schüz A. Cortex: statistics and geometry of neuronal connectivity. Berlin: Springer; 1998.

35. Jansen BH, Rit VG. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol Cybern. 1995;73(4):357–66.
36. Molaee-Ardekani B, Márquez-Ruiz J, Merlet I, Leal-Campanario R, Gruart A, Sánchez-Campusano R, Birot G, Ruffini G, Delgado-García J-M, Wendling F. Effects of transcranial Direct Current Stimulation (tDCS) on cortical activity: a computational modeling study. Brain Stimul. 2013;6(1):25-39.

37. Pittenger C, Bloch MH, Williams K. Glutamate abnormalities in obsessive compulsive disorder: neuropsychiology, pathophysiology, and treatment. Pharmacol Ther. 2011;132(3):314-32.

38. Hertz L, Rodrigues TB. Astrocytic-neuronal-astrocytic pathway selection for formation and degradation of glutamate/GABA. Lausanne: Frontiers E-books; 2014.

39. Patel AB, de Graaf RA, Mason GF, Rothman DL, Shulman RG, Behar KL. The contribution of GABA to glutamate/glutamine cycling and energy metabolism in the rat cortex in vivo. Proc Natl Acad Sci. 2005;102(15):5588-93.

40. Liang SL, Carlson GC, Coulter DA. Dynamic regulation of synaptic GABA release by the glutamate-glutamine cycle in hippocampal area CA1. J Neurosci. 2006;26(33):8537-48.

41. Oberheim NA, Takano T, Han X, He W, Lin JH, Wang F, Xu Q, Wyatt JD, Pilcher W, Ojemann JG, Ransom BR, Goldman SA, Nedergaard M. Uniquely hominid features of adult human astrocytes. J Neurosci. 2009;29(10):3276-87.

42. Bushong EA, Martone ME, Jones YZ, Ellisman MH. Protoplasmic astrocytes in CA1 stratum radiale occupy separate anatomical domains. J Neurosci. 2002;22(1):183-92.

43. Henneberger C, Papouin T, Oliet SH, Rusakov DA. Long-term potentiation depends on release of D-serine from astrocytes. Nature. 2010;463(7278):232-6.

44. Bellone C, Lüscher C, Mameli M. Mechanisms of synaptic depression triggered by metabotropic glutamate receptors. Cell Mol Life Sci. 2008;65(18):2913-23.

45. Niswender CM, Conn PJ. Metabotropic glutamate receptors: physiology, pharmacology, and disease. Annu Rev Pharmacol Toxicol. 2010;50:295-322.

46. Brickley SG, Mody I. Extrasynaptic GABA(A) receptors: their function in the CNS and implications for disease. Neuron. 2012;73(1):23-34.

47. Lüscher B, Fuchs T, Kilpatrick CL. GABAA receptor trafficking-mediated plasticity of inhibitory synapses. Neuron. 2011;70(3):385-409.

48. Campbell SL, Hablitz JJ. Decreased glutamate transport enhances excitability in a rat model of cortical dysplasia. Neurobiol Dis. 2008;32(2):254-61.