Aurélie Garnier

Alexandre Vidal
email: alexandre.vidal@univ-evry.fr

Habib Benali

Benali A Habib

A theoretical

Keywords: ODE Model in Neuroscience, Qualitative analysis of dynamical systems, Bifurcations, Neuronastrocyte interactions, GABAergic and glutamatergic neurotransmissions, Excitability modulation, Neuronal hyperexcitability AMS classification: 34C15, 34C23, 34C28, 34C46, 34C60, 34H20, 92C20, 92B25

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 post-synaptic terminals [START_REF] Araque | Glutamate-dependent astrocyte modulation of synaptic transmission between cultured hippocampal neurons[END_REF][START_REF] Parpura | Physiological astrocytic calcium levels stimulate glutamate release to modulate adjacent neurons[END_REF]. 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 these latter glutamate into the extra-synaptic space [START_REF] Bazargani | Astrocyte calcium signaling: the third wave[END_REF][START_REF] De Pittà | Astrocytes: orchestrating synaptic plasticity?[END_REF].

The predominant expression of high affinity astrocytic glutamate and GABA transporters in proximity of synaptic terminals [START_REF] Arizono | Receptor-selective diffusion barrier enhances sensitivity of astrocytic processes to metabotropic glutamate receptor stimulation[END_REF][START_REF] Sarup | 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[END_REF][START_REF] Schousboe | GABA transport inhibitors and seizure protection: the past and future[END_REF] 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 [START_REF] Wang | The astrocyte odyssey[END_REF][START_REF] Murphy-Royal | Surface diffusion of astrocytic glutamate transporters shapes synaptic transmission[END_REF][START_REF] Huang | Astrocyte glutamate transporters regulate metabotropic glutamate receptor-mediated excitation of hippocampal interneurons[END_REF][START_REF] Tzingounis | Glutamate transporters: confining runaway excitation by shaping synaptic transmission[END_REF]. 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 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 [START_REF] Nadkarni | Modeling synaptic transmission of the tripartite synapse[END_REF][START_REF] Postnov | Dynamical patterns of calcium signaling in a functional model of neuron-astrocyte networks[END_REF][START_REF] Postnov | Functional modeling of neural-glial interaction[END_REF][START_REF] Volman | Computational models of neuron-astrocyte interaction in epilepsy[END_REF][START_REF] Volman | The astrocyte as a gatekeeper of synaptic information transfer[END_REF], cellular [START_REF] De Pittà | A tale of two stories: astrocyte regulation of synaptic depression and facilitation[END_REF][START_REF] Gruetter | A mathematical model of compartmentalized neurotransmitter metabolism in the human brain[END_REF][START_REF] Silchenko | Computational modeling of paroxysmal depolarization shifts in neurons induced by the glutamate release from astrocytes[END_REF][START_REF] De Pittà | Astrocytes: orchestrating synaptic plasticity?[END_REF] and supracellular, including a hemodynamic compartment [START_REF] Aubert | A model of the coupling between brain electrical activity, metabolism, and hemodynamics: application to the interpretation of functional neuroimaging[END_REF][START_REF] Aubert | Brain lactate kinetics: Modeling evidence for neuronal lactate uptake upon activation[END_REF][START_REF] Aubert | Interaction between astrocytes and neurons studied using a mathematical model of compartmentalized energy metabolism[END_REF], and network [START_REF] Savin | Epileptogenesis due to glia-mediated synaptic scaling[END_REF][START_REF] Volman | Divide and conquer: functional segregation of synaptic inputs by astrocytic microdomains could alleviate paroxysmal activity following brain trauma[END_REF]. Despite the different scales of analysis, consensus has grown from these studies for the existence of plausible causal link between astrocyte activation and neuronal (hyper)excitability [START_REF] Volman | The astrocyte as a gatekeeper of synaptic information transfer[END_REF][START_REF] Postnov | Dynamical patterns of calcium signaling in a functional model of neuron-astrocyte networks[END_REF][START_REF] Volman | Computational models of neuron-astrocyte interaction in epilepsy[END_REF][START_REF] Nadkarni | Synaptic inhibition and pathologic hyperexcitability through enhanced neuronastrocyte interaction: a modeling study[END_REF][START_REF] Savin | Epileptogenesis due to glia-mediated synaptic scaling[END_REF]. With few exceptions however [START_REF] De Pittà | Astrocytes: orchestrating synaptic plasticity?[END_REF][START_REF] Volman | The astrocyte as a gatekeeper of synaptic information transfer[END_REF][START_REF] De Pittà | A tale of two stories: astrocyte regulation of synaptic depression and facilitation[END_REF], 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 [START_REF] Nunez | On the relationship of synaptic activity to macroscopic measurements: does co-registration of EEG with fMRI make sense?[END_REF]. The aggregated yet biologically significant parameters of a mass model can in fact efficiently be chosen so as to fit experimental data [36]. 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 [START_REF] Touboul | Neural mass activity, bifurcations, and epilepsy[END_REF][START_REF] Garnier | A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics[END_REF]. In particular, we benefit from the bifurcation analysis of the Neural Mass Model (NMM) performed in [START_REF] Garnier | A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics[END_REF] 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 [START_REF] Garnier | A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics[END_REF] 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 [START_REF] Blanchard | A new computational model for neuro-glio-vascular coupling: Astrocyte activation can explain cerebral blood flow nonlinear response to interictal events[END_REF]. 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 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 analytical characterization of the impact of 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 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.

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 (Figure 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 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 what follows, 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 [START_REF] Blanchard | A new computational model for neuro-glio-vascular coupling: Astrocyte activation can explain cerebral blood flow nonlinear response to interictal events[END_REF][START_REF] Garnier | A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics[END_REF], 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 [START_REF] Blanchard | A new computational model for neuro-glio-vascular coupling: Astrocyte activation can explain cerebral blood flow nonlinear response to interictal events[END_REF][START_REF] Garnier | A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics[END_REF].

GABA

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. This latter way of modeling amounts to considering an intermediary population of pyramidal 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 [START_REF] Garnier | A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics[END_REF], 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 (Figure 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. [START_REF] Van Rotterdam | A model of the spatial-temporal characteristics of the alpha rhythm[END_REF]:

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 F E and F I respectively defined as:

F E (h E ) = 1 A 1 a d 2 h E dt 2 + 2 h E + a h E , (1a)
F I (h I ) = 1 B 1 b d 2 h I dt 2 + 2 h I + b h I . (1b)
In this framework, parameter A (resp. B) stands for the average excitatory (resp. inhibitory) synaptic gain and tunes the amplitude of excitatory (resp. inhibitory) postsynaptic potentials. Additionally, 1 a (resp. 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 [START_REF] Freeman | Mass Action in the Nervous System[END_REF][START_REF] Van Rotterdam | A model of the spatial-temporal characteristics of the alpha rhythm[END_REF][START_REF] Jansen | A neurophysiologically-based mathematical model of flash visual evoked potentials[END_REF]. Following Freeman's work [START_REF] Freeman | Mass Action in the Nervous System[END_REF], the functions converting the average membrane potential into an average pulse density can be approximated by sigmoids of the type: S(x, x th , r sl ) = 1 1 + e r sl (x th -x) . Yet, for sake of compactness of the NMM presentation, we introduce an auxiliary parameterization after [START_REF] Freeman | Mass Action in the Nervous System[END_REF]:

+ + + _ (a) (b) p(t) C 1 C 2 C 3 C 4 G y 0 y 1 y 2 + x 0 x 1 x 2 h e h E
Ŝ(x, v) = 2 e 0 S(x, v, r) = 2 e 0 1 + e r (v-x)
where 2 e 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 (Figure 2(b)):

d 2 x 0 dt 2 = A a Ŝ(x 1 + G x 0 -x 2 , v P ) -2 a dx 0 dt -a 2 x 0 , (2a)
d 2 x 1 dt 2 = A a C 2 Ŝ(C 1 x 0 , v P ) -2 a dx 1 dt -a 2 x 1 + A a p(t), (2b)
d 2 x 2 dt 2 = B b C 4 Ŝ(C 3 x 0 , v I ) -2 b dx 2 dt -b 2 x 2 . ( 2c
)
Parameters C i (i ∈ {1, 2, 3, 4}) represent the average number of synapses between two populations. Following [START_REF] Braitenberg | Cortex: Statistics and Geometry of Neuronal Connectivity[END_REF], 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 [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF]: the excitatory output (y 0 = x 0 ) and the excitatory (y 1 = x 1 + G x 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:

dy 0 dt = y 3 , (3a)
dy 1 dt = y 4 , (3b)
dy 2 dt = y 5 , (3c)
dy 3 dt = A a Ŝ(y 1 -y 2 , v P ) -2 a y 3 -a 2 y 0 , (3d)
dy 4 dt = A a C 2 Ŝ(C 1 y 0 , v P ) + A a G Ŝ(y 1 -y 2 , v P ),
-2 a y 4a 2 y 1 + A a p(t), (3e)

dy 5 dt = B b C 4 Ŝ(C 3 y 0 , v I ) -2 b y 5 -b 2 y 2 . ( 3f
)
In this study, we consider the local field potential LFP(t) = y 1 (t)y 2 (t) [START_REF] Jansen | A neurophysiologically-based mathematical model of flash visual evoked potentials[END_REF] 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 .

Astrocyte model: glutamate and GABA concentration dynamics

For reproducing the astrocyte activity, we use the model introduced in [START_REF] Blanchard | A new computational model for neuro-glio-vascular coupling: Astrocyte activation can explain cerebral blood flow nonlinear response to interictal events[END_REF]. It focuses on the dynamics of glutamate and GABA concentrations, which are the main neurotransmitters of the central nervous system. In [START_REF] Blanchard | A new computational model for neuro-glio-vascular coupling: Astrocyte activation can explain cerebral blood flow nonlinear response to interictal events[END_REF], the neuronastrocyte coupling is feedforward: the astrocyte dynamics is driven by the neural activity, generated by the Jansen-Rit model, but 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 (Figure 3): excited pyramidal neurons (resp. interneurons) release glutamate (resp. GABA) in the extracellular space (synaptic cleft). Astrocytes and pre-synaptic neurons uptake the neurotransmitters. Astrocytes recycle or consume the neurotransmitters while the presynaptic neurons capture them to complete their stock. Following [START_REF] Blanchard | A new computational model for neuro-glio-vascular coupling: Astrocyte activation can explain cerebral blood flow nonlinear response to interictal events[END_REF] 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. Dynamics governing J G and J γ are driven by second-order differential operators similar to those for synaptic transfer dynamics introduced in (1) [START_REF] Molaee-Ardekani | Effects of transcranial Direct Current Stimulation (tDCS) on cortical activity : A computational modeling study[END_REF][START_REF] Van Rotterdam | A model of the spatial-temporal characteristics of the alpha rhythm[END_REF], i.e.:

F G (h G ) = 1 W 1 w 1 d 2 h G dt 2 + w 1 + w 2 w 1 dh G dt + w 2 h G , F γ (h γ ) = 1 Z 1 z 1 d 2 h γ dt 2 + z 1 + z 2 z 1 dh γ dt + z 2 h γ .
Likewise synaptic transfer functions in (1), parameter W (resp. Z) tunes the peak amplitude of glutamate (resp. GABA) concentrations and parameters w 1 and w 2 (resp. z 1 and z 2 ) tune the rise and decay times of glutamate (resp. GABA) release transfer function. These dynamics are well-suited to reproduce qualitative and quantitative properties of rise and decay of neurotransmitter concentrations as shown in [START_REF] Blanchard | A new computational model for neuro-glio-vascular coupling: Astrocyte activation can explain cerebral blood flow nonlinear response to interictal events[END_REF].

Astrocytic uptake of glutamate and GABA from the extracellular space may be modeled by Michaelis-Menten kinetics [START_REF] Blanchard | A new computational model for neuro-glio-vascular coupling: Astrocyte activation can explain cerebral blood flow nonlinear response to interictal events[END_REF], using the following function: The dynamics of the extracellular concentrations ([Glu] e and [GABA] e ) result from neurotransmitter fluxes J G and J γ 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 [36].

H(x, k) = x x + k (4) p(t) A h I h E h E [GABA] e [Glu] e [Glu] a [GABA] a P P ′ I Ŝ Ŝ Ŝ J G J γ
To summarize, the feedforward model obtained by coupling the neural mass model defined by (3) with astrocytic dynamics introduced in [START_REF] Blanchard | A new computational model for neuro-glio-vascular coupling: Astrocyte activation can explain cerebral blood flow nonlinear response to interictal events[END_REF] reads:

d 2 y 0 dt 2 = A a Ŝ(y 1 -y 2 , v P ) -2 a dy 0 dt -a 2 y 0 , (5a)
d 2 y 1 dt 2 = A a C 2 Ŝ(C 1 y 0 , v P ) + A a G Ŝ(y 1 -y 2 , v P ) -2 a dy 1 dt -a 2 y 1 + A a p(t), (5b)
d 2 y 2 dt 2 = B b C 4 Ŝ(C 3 y 0 , v I ) -2 b dy 2 dt -b 2 y 2 , (5c)
d 2 J G dt 2 = W w 1 Ŝ(y 1 -y 2 , v P ) -(w 1 + w 2 ) dJ G dt -w 1 w 2 J G , (5d)
d[Glu] e dt = J G -(V ae G + V ne G ) S([Glu] e , s g , r g ), (5e)
d[Glu] a dt = V ae G S([Glu] e , s g , r g ) -V c G [Glu] a , (5f)
d 2 J γ dt 2 = Z z 1 Ŝ(C 3 y 0 , v I ) -(z 1 + z 2 ) dJ γ dt -z 1 z 2 J γ , (5g)
d[GABA] e dt = J γ -V ae γ H([GABA] e , K ae γ ) -V ne γ H([GABA] e , K ne γ ), (5h)
d[GABA] a dt = V ae γ H([GABA] e , K ae γ ) -V c γ [GABA] a . (5i)
In the sigmoidal functions for [Glu] e dynamics (5e), parameters V ae G and V ne G 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 ae γ and K ae γ (resp. V ne γ and K ne γ ) in equation (5h) are respectively the maximum rate and concentration for astrocyte (resp. neural) GABA transporter. Finally, V c G and V c γ are the glutamate and GABA degradation rates in astrocytes. We refer the reader to [START_REF] Blanchard | A new computational model for neuro-glio-vascular coupling: Astrocyte activation can explain cerebral blood flow nonlinear response to interictal events[END_REF] for a detailed explanation of the dynamics.

System ( 5) 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 [START_REF] Araque | Glutamate-dependent astrocyte modulation of synaptic transmission between cultured hippocampal neurons[END_REF] 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.

Astrocytic feedback and neuron-astrocyte mass model

The concentrations of neurotransmitters in a synaptic cleft act on the excitability threshold of the post-synaptic neuron. In the neuron-astrocyte model ( 5) 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 [START_REF] Araque | Glutamate-dependent astrocyte modulation of synaptic transmission between cultured hippocampal neurons[END_REF]. 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 be easily 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 P G S([Glu] e , v G , r G ) for glutamate feedback on pyramidal neurons, (b) m I G S([Glu] e , v G , r G )
for glutamate feedback on interneurons, (c) m γ S([GABA] e , v γ , r γ ) 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 [START_REF] Pittenger | Glutamate abnormalities in obsessive compulsive disorder: neurobiology, pathophysiology, and treatment[END_REF][START_REF] Huang | Astrocyte glutamate transporters regulate metabotropic glutamate receptor-mediated excitation of hippocampal interneurons[END_REF]. Thus, only parameters m P G and m I G representing the maximum coupling gains of the glutamate-dependent component of the astrocyte feedback discriminate between the coupling functions m P G S([Glu] e , v G , r G ) and m I G 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 (5) 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 N2 in the equation corresponding to the output of population N 2 since this parameter represents a modulation of the threshold of the sigmoidal function Ŝ.

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 dotted arrows in Figure 3 illustrate these feedbacks. We need to consider separately each type of synapse in the NMM, and the variables x 0 , x 1 an x 2 for the feedbacks building. The NMM embeds five types of synaptic connections between populations: (i) P → P , (ii) P → P , (ii) P → I, (iv) I → P , (v) P → P . In the following 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 → P or P → 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 non local 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 postsynaptic neurons of population P .

A synaptic connection of type P → I concerns the variable x 2 . As shown in [START_REF] Araque | Glutamate-dependent astrocyte modulation of synaptic transmission between cultured hippocampal neurons[END_REF][START_REF] Hertz | Astrocytic-neuronal-astrocytic Pathway Selection for Formation and Degradation of glutamate/GABA[END_REF], 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 equation (2c):

v I = v 0 -m I G S([Glu] e , v G , r G ).
While experiments have evidenced impacts of extracellular glutamate concentrations on the excitability of GABAergic neurons [START_REF] Hertz | Astrocytic-neuronal-astrocytic Pathway Selection for Formation and Degradation of glutamate/GABA[END_REF], 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 [START_REF] Patel | The contribution of GABA to glutamate/glutamine cycling and energy metabolism in the rat cortex in vivo[END_REF][START_REF] Liang | Dynamic regulation of synaptic GABA release by the glutamateglutamine cycle in hippocampal area CA1[END_REF]. 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 → 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. [START_REF] Araque | Glutamate-dependent astrocyte modulation of synaptic transmission between cultured hippocampal neurons[END_REF]), 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 → 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 [START_REF] Araque | Glutamate-dependent astrocyte modulation of synaptic transmission between cultured hippocampal neurons[END_REF]. 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 equation (3d)

v P = v 0 + m γ S([GABA] e , v γ , r γ ) -m P G S([Glu] e , v G , r G ).
The new neuron-astrocyte mass model embedding the astrocyte feedback is obtained from model (5) by considering the dynamical entries v I and v P mentioned above. Accordingly, the sigmoidal functions appearing in equations (5a), (5b), (5c), (5d) and (5g) become

Ŝ(y 1 -y 2 , v P ) = Ŝ(y 1 -y 2 , v 0 + m γ S([GABA] e , v γ , r γ ) -m P G S([Glu] e , v G , r G )), Ŝ(C 1 y 0 , v P ) = Ŝ(C 1 y 0 , v 0 ), Ŝ(C 3 y 0 , v I ) = Ŝ(C 3 y 0 , v 0 -m I G S([Glu] e , v G , r G )).

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 [START_REF] Touboul | Neural mass activity, bifurcations, and epilepsy[END_REF] on the Jansen-Rit model. In [START_REF] Garnier | A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics[END_REF], we have classified the types of time series patterns generated by model ( 3) 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 (3) has the following useful features [START_REF] Garnier | A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics[END_REF]. 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 Figure 4), this curve of singular points is S-shaped. In the following description, for a given bifurcation "X" according to parameter p, we denote p X the bifurcation value and, if the bifurcation involves a singular point, we denote 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 SN1 , y SN1 < y 0 < y SN2 and y 0 > y SN2 . 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) = p such that 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 p as mean value and large enough variance σ, the NMM output results instead in an alternation of quiescent phases and isolated LFP spikes (see Figure 4), consistent with emergence of Noise-Induced Spiking (NIS) [START_REF] Garnier | A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics[END_REF]. 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 [pσ, p + σ] reported by the grey bar above the bifurcation diagram on the left panel of Figure 4.

Such LFP activity, i.e. sparse large amplitude spikes, corresponds to episodic synchronization 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 hyperexcitability during pathological crisis. In [START_REF] Garnier | A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics[END_REF] 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 experimental data recorded from epileptic mice (Mesial Temporal Lobe Epilepsy mouse model). It appeared that the NIS behavior is well-suited for reproducing the high amplitude 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 along time. The values of the parameters associated with the NMM used in this study and given in Table 1 have been chosen to reproduce NIS behavior using the analysis in [START_REF] Garnier | A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics[END_REF]. In the following, we study the variations of the activity when the neural dynamics are altered by the surrounding activity, i.e. the astrocyte feedback.

NMM dynamics: example from GABA injection

To illustrate the impact of the astrocytic feedback on neural excitability and introduce 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 feedback, leveraging on the role of the SNIC bifurcation occurring at p = p SNIC 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 SNIC (Figures 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 -10 -20 -30

(µM) (µM) (mV) (mV) (s -1 ) p (s -1 ) p (s -1 ) (s -1 ) t (s) t (s)
Figure 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 = 0s (b1) and t = 30s (b2). Time series (bottom panels) corresponding to LFP, glutamate and GABA extracellular concentrations and time variations of p SNIC 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 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 Figure 5), the neural activity and the glutamate concentration dynamics remain unchanged (Figure 5(a)) after the instantaneous increase of [GABA] e at t = 0 that mimics GABA injection. In contrast, in the model with astrocyte feedback (right panels of Figure 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 phenomenons can be well accounted by arguments based on bifurcation analysis (top panels of Figure 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 (Figure 5(a1)). In this fashion, the excitability threshold p SNIC of the neural compartment is constant in time (Figure 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 SNIC which changes along with astrocytic variables. Thus, right after GABA injection, p SNIC is larger than the input value p (Figure 5(b2)), and because p 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 towards the attractive state slowly decreasing p SNIC , as shown in the bottom panel of Figure 5 case (b). Once, GABA is low enough, p 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 the biological system.

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 [START_REF] Schousboe | GABA transport inhibitors and seizure protection: the past and future[END_REF]. Because a single astrocyte may ensheathe hundreds of thousands of synapses [START_REF] Oberheim | Uniquely hominid features of adult human astrocytes[END_REF][START_REF] Bushong | Protoplasmic astrocytes in CA1 stratum radiatum occupy separate anatomical domains[END_REF], 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 Section 2.5. Thus, once we consider many defective astrocytes, all neurons in the proximity of these astrocytes are likely to be inhibited [START_REF] Henneberger | Long-term potentiation depends on release of d-serine from astrocytes[END_REF].

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 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:

m I G S([Glu] e , v G , r G ) ≡ v 1 , m P G S([Glu] e , v G , r G ) ≡ m P G m I G v 1 , m γ S([GABA] e , v γ , r γ ) ≡ v 2 .
where v 1 and v 2 are defined within the boundaries of

m I G S([Glu] e , v G , r G ) and m γ S([GABA] e , v γ , r γ ), that is 0 ≤ v 1 ≤ m I G and 0 ≤ v 2 ≤ m γ .
With these new notations, the dynamical excitability thresholds v P , v P and v I of populations P , P and I become:

v P = v 0 + v 2 - m P G m I G v 1 , v P = v 0 , v I = v 0 -v 1 .
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 the following Proposition 3.1 and exemplify it by numerical continuation of the SNIC bifurcation for varying v 2 in Figure 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. (AU)

(µM) (µM) (mV) v 2 (mV) (s -1 )
t (s) 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 (3), a singular point satisfies y 3 = y 4 = y 5 = 0. Under these conditions, annihilation of the r.h.s. of equation (3f) allows expressing y 2 with 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 equation (3e), we obtain that 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 ) (6)
where

f (y 0 , v 1 , v 2 ) = a A (v 0 - m P G m I G v 1 + v 2 ) - a A r ln 2 A e 0 -a y 0 a y 0 - a G A y 0 -C 2 Ŝ(C 1 y 0 , v 0 ) + a B b A C 4 Ŝ(C 3 y 0 , v 0 -v 1 ). ( 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 (5) as explained above. Since these components are not needed for the following analysis and are only minor calculus details, we don't specify these heavy expressions. We rewrite equation ( 7) as follows

f (y 0 , v 1 , v 2 ) = a A v 2 + q(y 0 , v 1 ). ( 8
)
Obviously the two saddle-node bifurcation values p SN1 and p SN2 are relative extrema of f (y 0 , v 1 , v 2 ). In particular p SN1 = 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 ), ( 9a
) ∂f ∂y0 (y 0 , v 1 , v 2 ) = 0, (9b)
∂ 2 f ∂y 2 0 (y 0 , v 1 , v 2 ) 0. ( 9c
)
Since ∂f ∂y0 (y 0 , v 1 , v 2 ) is independent on v 2 , so is y SNIC and it can be considered as a parameter in equation (9a). From ( 8) and ( 9) we obtain the following expression for p SNIC

p SNIC = a A v 2 + q(y SNIC , v 1 ).
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 Figure 5), the value of v 2 increases and Proposition 3.1 asserts that 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, 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 = 40s, we turn off the GABA astrocyte uptake by setting V ae γ = 0 (Figure 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 (Figure 6). Consequently, we observe a decrease in the oscillation frequency after t > 40s. In the time series, the oscillation frequency decreases gradually during a transient (40s < t < 60s) until reaching its minimum. This can be explained by the slow increase of GABA extracellular concentration that reaches its new baseline at t = 60s.

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, in 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.

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 non trivial dynamics regimes for neural activity. Along the rationale outlined in the previous section, we now set to characterize the changes in the SNIC bifurcation value p SNIC (i.e. the threshold for neural activation) with v 1 value. This latter variable, introduced in the beginning of Section 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 ae G . We note on 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 [START_REF] Garnier | A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics[END_REF] 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 ≤ v 1 ≤ m I G , that is the maximum interval of values taken by v 1 = m I G 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 ) ( 10
)
where f is given by ( 7). Since we consider v 2 fixed we introduce the function

g(y 0 , v 1 ) ≡ f (y 0 , v 1 , v 2 )| v2=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

∂g ∂y 0 (y SNIC , v 1 ) = 0, ∂ 2 g ∂y 2 0 (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 ≤ v 1 ≤ m I G , p SNIC exists and the associated saddle-node bifurcation is not degenerate. Note χ = B e0 r C4

2 b . Then 1. if m P G m I G χ, p SNIC (v 1 ) has no local extremum, 2. if 0 < m P G m I G < χ, 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 function p SNIC (v 1 ) which is implicitly defined by [START_REF] Murphy-Royal | Surface diffusion of astrocytic glutamate transporters shapes synaptic transmission[END_REF]. Hence, we are interested in the following constrained optimization problems:

max g(y 0 , v 1 ) | ∂g ∂y 0 (y 0 , v 1 ) = 0 min g(y 0 , v 1 ) | ∂g ∂y 0 (y 0 , v 1 ) = 0 . (11a)
We introduce the associated Lagrangian function

L(y 0 , v 1 , λ) = g(y 0 , v 1 ) -λ ∂g ∂y 0 (y 0 , v 1 ).
The necessary condition for the existence of an extremum of g under the constraint ∂g ∂y0 = 0 is

- → ∇L(y 0 , v 1 , λ) = 0, that is ∂g ∂y 0 (y 0 , v 1 ) -λ ∂ 2 g ∂y 2 0 (y 0 , v 1 ) = 0, (12a)
∂g ∂v 1 (y 0 , v 1 ) -λ ∂ 2 g ∂v 1 ∂y 0 (y 0 , v 1 ) = 0, (12b)
∂g ∂y 0 (y 0 , v 1 ) = 0. ( 12c
)
By assumption, the saddle-node bifurcation associated with the SNIC bifurcation is not degenerate, i.e. every solution of (11a) such as

0 ≤ v 1 ≤ m I G satisfies ∂ 2 g ∂y 2 0
(y 0 , v 1 ) = 0. Thus, system [START_REF] Nadkarni | Modeling synaptic transmission of the tripartite synapse[END_REF] reads

λ = 0, (13a) ∂g ∂v 1 (y 0 , v 1 ) = 0, (13b)
∂g ∂y 0 (y 0 , v 1 ) = 0. (13c)
Therefore, if the problem under constraint admits an extremum, at this extremum it is ∂g ∂v1 = 0. Following the assumption that a SNIC bifurcation occurs for any value of v 1 between 0 and m I G , then equation (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 ∂g ∂y 0 (y 0 , v 1 ) = 0.

From [START_REF] Schousboe | GABA transport inhibitors and seizure protection: the past and future[END_REF], we obtain

∂g ∂v 1 (y 0 , v 1 ) = - a A m P G m I G + B b C 4 ∂ Ŝ ∂v (C 3 y 0 , v 0 -v 1 )
.

Using the facts that, for any fixed value of y 0 , function v 1 → ∂g ∂v1 (y 0 , v 1 ) is bell-shaped and its maximum value does not depend on y 0 (see Figure 7), one obtains that function ∂g ∂v1 (y 0 , v 1 ) vanishes in v 1 if

0 < m P G m I G < χ. (14)
If

m P G m I G
χ, function ∂g ∂v1 (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 ∂g ∂v1 (y 0 , v 1 ) = 0, i.e.

h(v 1 ) m P G m I G v * * 1 v * 1 m I G v 0 C 3 y * 0 - 0 I 1 I 2 v 1 (mV) χ
m P G m I G + B b C 4 ∂ Ŝ ∂v (C 3 y 0 , v 0 -v 1 ) = 0.
Noting X = e r (v0-v1-C3 y0) = (S(C 3 y 0 , v 0v 1 , r)) -1 -1, this equation reads:

m P G m I G X 2 + 2 m P G m I G -2χ X + m P G m I G = 0. (15)
Setting

V ± = 2 χ m I G m P G -1 ± 2 χ m I G m P G χ m I G m P G -1 (16)
we obtain the two solutions v * 1 < v * * 1 of ∂g ∂v1 (y 0 , v 1 ) = 0:

v * 1 = v 0 -C 3 y 0 - 1 r ln (V + ) , (17)
v * * 1 = v 0 -C 3 y 0 - 1 r ln (V -) . (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 ∂g ∂y 0 (y 0 , v * 1 ) = 0,

∂ 2 g ∂y 2 0 (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 H associated with the lagrangian function at its singular point (y 0 , v 1 , λ) = (y * 0 , v * 1 , 0) (solution of system ( 13)):

H(y * 0 , v * 1 , 0) =     0 ∂ 2 g ∂y 2 0 ∂ 2 g ∂v1∂y0 ∂ 2 g ∂y 2 0 ∂ 2 L ∂y 2 0 ∂ 2 L ∂v1∂y0 ∂ 2 g ∂v1∂y0 ∂ 2 L ∂v1∂y0 ∂ 2 L ∂v 2 1     | (y * 0 ,v * 1 ,0) =     0 ∂ 2 g ∂y 2 0 ∂ 2 g ∂v1∂y0 ∂ 2 g ∂y 2 0 ∂ 2 g ∂y 2 0 ∂ 2 g ∂v1∂y0 ∂ 2 g ∂v1∂y0 ∂ 2 g ∂v1∂y0 ∂ 2 g ∂v 2 1     | (y * 0 ,v * 1 ,0)
The determinant of

H(y * 0 , v * 1 , 0) is given by det H(y * 0 , v * 1 , 0) = - ∂ 2 g ∂y 2 0 (y * 0 , v * 1 ) ∂ 2 g ∂y 2 0 (y * 0 , v * 1 ) ∂ 2 g ∂v 2 1 (y * 0 , v * 1 ) - ∂ 2 g ∂v 1 ∂y 0 (y * 0 , v * 1 ) 2 .
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 ∂ 2 g ∂y 2 0 (y * 0 , v * 1 ) < 0. On the other hand, for any y 0 , v 1 → ∂g ∂v1 (y 0 , v 1 ) is increasing at (y 0 , v * 1 ) (see Figure 7), thus ∂ 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 I G ]. We now derive a necessary and sufficient condition so that v * 1 actually lies in [0, m I G ]. Since v * 1 satisfies ∂g ∂v1 (y * 0 , v * 1 ) = 0, one obtains, from equation [START_REF] Volman | Computational models of neuron-astrocyte interaction in epilepsy[END_REF],

m P G m I G = 4χ (1 -S(C 3 y * 0 , v 0 -v * 1 , r)) S(C 3 y * 0 , v 0 -v * 1 , r) = 4χ h(v * 1 ). (19)
For any y 0 , function h :

v 1 → (1 -S(C 3 y 0 , v 0 -v 1 , r)) S(C 3 y 0 , v 0 -v 1 , r) is strictly increasing over [0, m I G ], and 0 ≤ v * 1 ≤ m I G if and only if I 1 ≤ m P G m I G ≤ I 2
where

I 1 = 4χ (1 -S(C 3 y * 0 , v 0 , r)) S(C 3 y * 0 , v 0 , r), (20a)
I 2 = 4χ 1 -S(C 3 y * 0 , v 0 -m I G , r) S(C 3 y * 0 , v 0 -m I G , r). (20b)
Figure 7 gives an illustration of 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

≤ v * 1 ≤ m I G if and only if m P G m I G ∈ (0, χ) ∩ [I 1 , I 2 ] .
Moreover, in section 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

c) If I 1 ≤ m P G m I G ≤ I 2 , then 0 ≤ v * 1 ≤ m I G 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) m P G m I G = 1.7 < I 1 , (b) m
P G m I G = 3.2 > I 2 , (c) m P G m I G = 2.43 ∈ [I 1 , I 2 ] .
For each case, we provide simulations representing the value of p SNIC in (v 1 , v 2 ) space and time series generated by the model when the glutamate astrocytic uptake is altered (Figure 8). For 8), v * 1 is negative and thus p 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 and thus neural activity.

m P G m I G = 1.7 < I 1 (Panel (a) of Figure
Panel (b) of Figure 8 is obtained with

m P G m I G = 3.2 > I 2 . Since v * 1 > 1,
p SNIC decreases as v 1 increases. Thereby, an astrocyte glutamate uptake deficiency triggers an increase of 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

m P G m I G = 2.43 ∈ [I 1 , I 2 ]
is shown in panel (c) of Figure 8. In this case p 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 prescienced 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

m P G m I G
. Note that this value can be tuned to obtain v * 1 small enough and p 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. All other parameters are the same in the three cases and given in Table 1.

p SNIC v 1 v 2 p SNIC v 1 v 2 p SNIC v 1 v 2 (a) (b) (c) LFP [Glu] e [GABA] e (AU) v 1 LFP [Glu] e [GABA] e (AU) v 1 LFP [Glu] e [GABA] e (AU) (µM) (µM) (mV) (µM) (µM) (mV) (µM) (µM) (mV) (mV) (mV) (mV) (mV) (mV) (mV) t (s) t (s) t (s)

Discussion and Conclusions

We have introduced a novel neuron-astrocyte mass model that combines two previously studied models [START_REF] Blanchard | A new computational model for neuro-glio-vascular coupling: Astrocyte activation can explain cerebral blood flow nonlinear response to interictal events[END_REF][START_REF] Garnier | A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics[END_REF] 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 [START_REF] Bellone | Mechanisms of synaptic depression triggered by metabotropic glutamate receptors[END_REF][START_REF] Huang | Astrocyte glutamate transporters regulate metabotropic glutamate receptor-mediated excitation of hippocampal interneurons[END_REF][START_REF] Niswender | Metabotropic glutamate receptors: physiology, pharmacology, and disease[END_REF][START_REF] Pittenger | Glutamate abnormalities in obsessive compulsive disorder: neurobiology, pathophysiology, and treatment[END_REF]. Using analytical arguments of bifurcation theory and contained optimization, we have characterized different types of change in the neural activity behavior under the 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 tends to fire more the slower astrocyte GABA uptake is, or equivalently, the larger extracellular GABA concentration is, which is consistent with experimental findings [START_REF] Brickley | Extrasynaptic GABA(A) receptors: their function in the CNS and implications for disease[END_REF][START_REF] Luscher | GABAA receptor trafficking-mediated plasticity of inhibitory synapses[END_REF]. 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 [START_REF] Garnier | A neural mass model with direct and indirect excitatory feedback loops: identification of bifurcations and temporal dynamics[END_REF], 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 P G /m I G , 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. With 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 G . Hence, one might interpret this 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 hyper-excitability.

A further interesting prediction of our model is that when

m P G m I G > I 1
, an increase of extracellular glutamate due to deficits in astrocytic uptake develops in parallel with an increase of extracellular GABA which may be even 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 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 possibly define more accurate treatments of neuropathogenic states associated with neuronal hyperexcitability or lack thereof.

Figure 1 :

 1 Figure 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.

Figure 2 :

 2 Figure 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 h E (resp. h I ): second order process converting action potentials into excitatory (resp. inhibitory) post-synaptic potential. Box Ŝ: process converting average membrane potential into average action potential density discharge by neurons of populations P , P and I respectively. C i 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 0 , y 1 , y 2 : state variables. x 0 , x 1 , x 2 : intermediary variables.

Figure 3 :

 3 Figure 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 en γ and J γ ). 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 .

Figure 4 :

 4 Figure 4: Bifurcation diagram according to p (left) and associated LFP and y 0 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 0 extrema along stable limit cycles. Black points (SN 1 and SN 2 ): saddle-node bifurcations. Red point (H 1 ): supercritical Hopf bifurcation. Dashed orange line: Saddle-Node on Invariant Circle (SNIC) bifurcation. Horizontal grey 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.

Figure 6 :

 6 Figure 6: The dependency of the SNIC value on v 2 explains the impact of an alteration of 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 grey time window highlights the transient towards the new behavior.

Figure 7 :

 7 Figure 7: Interval of m P G /m I values corresponding to a realized local minimum of neural excitability threshold.Graphic representation of function h defined by (19) and interval [I 1 , I 2 ] for which v * 1 ∈ [0, m I G ].

1 )

 1 > 0. Finally det 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 ).

< I 1 then v * 1 <> I 2 ,

 12 0 and p SNIC strictly increases with v 1 and v 2 . then v * 1 > m I G and p SNIC strictly decreases when v 1 increases (for v 2 fixed) and strictly increases with v 2 (for v 1 fixed).

(

Figure 8 :

 8 Figure 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 SNIC in (v 1 , v 2 ) plane, and the time series on the right correspond to LFP, [GABA] e , [Glu] e and v 1 = m I G S([Glu] e , v G , r G ). The white curves on the colormap are the traces of (m I G S([Glu] e , v G , r G ), m γ S([GABA] e , v γ , r γ )) along the associated orbits of the model. At t = 20s, we alter the glutamate astrocytic uptake by setting V ae G = 0. The three cases are obtained with : (a) m P G m I G

1 , 2

 12 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 the frequency of neural activity decreases and may enter quiescence while for values of and persistently increases leading to epileptiform activity in agreements with experiments[START_REF] Campbell | Decreased glutamate transport enhances excitability in a rat model of cortical dysplasia[END_REF]. On the other hand, for I 1 < 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 I

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.

Tables