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

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.

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.

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

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

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].

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 pyramidal neurons 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 (Fig. 2(a)).

Fig. 2

Two schematic representations of the NMM with double excitatory feedbacks. P: main population of pyramidal neurons. I: interneuron population. : secondary population of pyramidal neurons. Red (resp. green) arrows in (a): excitatory (resp. inhibitory) interactions. Box (resp. ): second-order process converting action potentials into excitatory (resp. inhibitory) postsynaptic potential. Box : process converting average membrane potential into average action potential density discharge by neurons of populations P, , and I, respectively. for : coupling gain parameters depending on the maximum number C of synaptic connections between two populations. G: direct feedback coupling gain. : external input. , , : state variables. , , : intermediary variables

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]: These are fundamental solutions (also known as Green’s functions) for the differential operators and , respectively, defined by 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) postsynaptic potentials. Additionally, (resp. ) 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 Yet, for the sake of compactness of the NMM presentation, we introduce an auxiliary parameterization after [32]: where 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 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 , and , which represent the outputs of the population P, the sum of population output and the external input , and the population I, respectively (Fig. 2(b)): Parameters () represent the average number of synapses between two populations. Following [34], each 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 () and the excitatory () and inhibitory () inputs of the main population P. The output acts on the secondary pyramidal neuron population and on the interneuron population I. To analyze the model, we write the dynamics of the state variables , , and as a system of first order differential equations: In this study, we consider the local field potential [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.

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

Fig. 3

Neuron-astrocyte model with astrocyte feedback. P and : main and secondary populations of pyramidal neurons. I: interneuron population. : external input on population P. and : glutamate and GABA extracellular concentrations. and : glutamate and GABA astrocyte concentrations. Red arrows: , and couplings. Orange arrow: coupling. Green arrow: coupling. Cyan arrows: glutamate and GABA release by populations P and I into extracellular space (fluxes and ). 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

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)  and : the fluxes of glutamate and GABA from neurons to extracellular space, (ii)  and : the neurotransmitter concentrations in the extracellular space, (iii)  and : the quantity of neurotransmitters recycled and consumed by the astrocytes. The dynamics governing and is driven by second-order differential operators similar to those for synaptic transfer dynamics introduced in (1a)-(1b) [31, 36], i.e.:

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 and (resp. and ) 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: The dynamics of the extracellular concentrations ( and ) result from neurotransmitter fluxes and 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 In the sigmoidal functions for dynamics (5e), parameters and are the maximum rates of glutamate uptakes by the astrocytes and the neurons, respectively, represents the activation threshold and the sigmoidal slope at the inflection point. Parameters and (resp. and ) in Eq. (5h) are, respectively, the maximum rate and concentration for astrocyte (resp. neural) GABA transporter. Finally, and 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.

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 , and . 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: 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 and representing the maximum coupling gains of the glutamate-dependent component of the astrocyte feedback discriminate between the coupling functions and , 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.

Table 1

Descriptions and values of the neuron–astrocyte mass model parameters

Parameter	Interpretation	Value
Neuronal compartment
A	Average excitatory synaptic gain	3.25mV
B	Average inhibitory synaptic gain	22mV
1a	Time constant of excitatory postsynaptic potentials	1100s
1b	Time constant of inhibitory postsynaptic potentials	150s
e_0	Half of the maximum discharge rate of a neuronal population	2.5s^−1
v_0	Basic excitability threshold for neurons	6mV
r	Stiffness of neuronal excitability	0.56mV^−1
C_1	Strength of the synaptic connections from P to P^′	135
C_2	Strength of the synaptic connections from P^′ to P	108
C_3	Strength of the synaptic connections from P to I	33.75
C_4	Strength of the synaptic connections from I to P	33.75
G	Gain of the direct excitatory feedback from P to itself	40
Glial compartment
W	Tunes the peak amplitude of glutamate concentrations	53.6μM⋅s^−1
Z	Tunes the peak amplitude of GABA concentrations	53.6μM⋅s^−1
w_1	Tune the rise and decay times of glutamate release transfer function	90s^−1
w_2		33s^−1
z_1	Tune the rise and decay times of GABA release transfer function	90s^−1
z_2		33s^−1
VGne	Maximal rate of glutamate uptake by neurons	0.5μM⋅s^−1
VGae	Maximal rate of glutamate uptake by astrocytes	4.5μM⋅s^−1
s_g	Activation threshold of sigmoid glutamate uptakes	6μM
r_g	Stiffness of sigmoid glutamate uptakes	0.9μM^−1
Vγae	Maximal rate of astrocyte GABA uptake	2μM⋅s^−1
Kγae	Maximal concentration for Hill dynamics of astrocyte GABA uptake	8μM
Vγne	Maximal rate of neuronal GABA uptake	5μM⋅s^−1
Kγne	Maximal concentration for Hill dynamics of neuronal GABA uptake	24μM
VGc	Rate of glutamate degradation by astrocytes	9μM⋅s^−1
Vγc	Rate of GABA degradation by astrocytes	9μM⋅s^−1
Neuron excitability modulations by neurotransmitter concentrations (feedbacks)
v_G	Excitability threshold of glutamate feedback function	30μM
r_G	Stiffness of sigmoid glutamate feedback function induced by glutamate	0.15μM^−1
mGP	Maximal coupling gain of glutamate feedback on pyramidal neurons	2.5mV
mGI	Maximal coupling gain of glutamate feedback on interneurons	1mV
v_γ	Excitability threshold of GABA feedback function	25μM^−1
r_γ	Stiffness of sigmoid GABA feedback function	0.12μM^−1
m_γ	Maximal coupling gain of GABA feedback	1mV

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 of population to a neuron of population , the extracellular concentration of neurotransmitter acts on the postsynaptic neuron 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 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 in the equation corresponding to the output of population 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 dashed arrows in Fig. 3 illustrate these feedbacks. We need to consider separately each type of synapse in the NMM, and the variables , , and for the feedbacks building. The NMM embeds five types of synaptic connections between populations: (i) , (ii) , (iii) , (iv) , (v) . 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 or . As a matter of fact, the astrocyte compartment only takes into account neurotransmitters released locally by neurons of populations P and I, whereas population is nonlocal to population P. Hence, extracellular concentrations of neurotransmitters in the vicinity of 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 .

A synaptic connection of type concerns the variable . 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 on the extracellular glutamate concentration and set in Eq. (2c):

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 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 dynamics. On the other hand, a synaptic connection of type is impacted by the extracellular concentration of glutamate implying a modulation of variable 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)

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

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 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 value of a singular point suffices to have explicit expressions of all the other components. Second, for a given value, there exists a unique value of p such that corresponds to a singular point. In other terms, the set of singular points obtained for the different values of p is a graph over . Hence, we can visualize the shape of the singular point locus in the plane : 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 the bifurcation value and, if the bifurcation involves a singular point, we denote by the corresponding value.

Fig. 4

Bifurcation diagram according to p (left) and associated LFP and 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: extrema along stable limit cycles. Black points ( and ): saddle-node bifurcations. Red point (): supercritical Hopf bifurcation. Dashed orange line: Saddle-Node on Invariant Circle (SNIC) bifurcation. Horizontal gray bar: confidence interval of the Gaussian variable used to generate the time series i.e. where p̄ and σ are the mean and the variance of the associated normal distribution, respectively

Two saddle-node bifurcations and 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 , , and . 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 and stable (blue) otherwise. At the system undergoes a supercritical Hopf bifurcation giving birth to a stable limit cycle for that persists until 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 ) and quiescent phases (for ). In other terms, plays the role of an activation threshold for the neural compartment.

This means that, given a constant input such that , 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 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 reported by the gray bar above the bifurcation diagram on the left panel of Fig. 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 [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 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 with time. The values of the parameters 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, we study the variations of the activity when the neural dynamics is 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 on the emergence of different regimes of neuronal firing.

With this aim, we simulate injection of GABA () into the extracellular space at . Assuming (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.

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 (b1) and (b2). Time series (bottom panels) corresponding to LFP, glutamate and GABA extracellular concentrations and time variations of 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

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 at , which mimics GABA injection. In contrast, in the model with astrocyte feedback (right panels of Fig. 5), the strong increase of halts neuronal activity, which, in turn reduces synaptic release of glutamate. Once 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 of the neural compartment is constant in time (Fig. 5, bottom panel). On the other hand, in the system with feedback, , and depend on GABA and glutamate concentrations. Consequently, the bifurcation diagram for the neural compartment is time dependent, and so is the value of , which changes along with the astrocytic variables. Thus, right after GABA injection, is larger than the input value p (Fig. 5(b2)), and because 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 , as shown in the bottom panel of Fig. 5, case (b). Once, GABA is low enough, 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.

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 () 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 on astrocytic feedback. With this aim, because glutamate binding by neuronal receptors is independent of cell types, we may set: where and are defined within the boundaries of and , that is, and .

With these new notations, the dynamical excitability thresholds , , and of populations P, , and I become

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

Fig. 6

The dependency of the SNIC value on explains the impact of an alteration of the GABA astrocytic uptake on the neural activity. Variation of value according to (left). Time series corresponding to LFP, , and (right from top to bottom) for a Gaussian variable. At , the GABA astrocytic uptake is artificially altered by setting for any subsequent time. The gray time window highlights the transient toward the new behavior

Proposition 3.1

The threshold for neuronal activation, , linearly increases with the GABA-induced modulation of neural excitability .

Proof

The set of the system singular points obtained for the different values of parameter p can be explicitly expressed according to , , and all other parameters being fixed. Using system (3a)-(3f), a singular point satisfies . Under these conditions, annihilation of the r.h.s. of Eq. (3f) allows expressing with respect to . Then annihilation of the r.h.s. of (3d) provides as function of . Finally, by annihilation of Eq. (3e), we find that the components of the singular points for given values of p, and are characterized as solutions of where All the other components of a given singular point can be explicitly expressed as functions of the 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: Obviously the two saddle-node bifurcation values and are relative extrema of . In particular coincide with the local maximum of and, as such, is of the solution of Since is independent on , so is and it can be considered as a parameter in Eq. (9a). From (8) and (9a)-(9c) we obtain the following expression for :  □

Let us consider the model generating an oscillatory output with a fixed value of p (). If the extracellular concentration of GABA increases (e.g. by an injection of GABA as in Fig. 5), the value of increases and Proposition 3.1 asserts that the value of also increases. As already explained, the closest is to p, with , the largest is the limit cycle period, thus the oscillation frequency of the outputs decreases. If increases enough such that , 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 a Gaussian input. At , we turn off the GABA astrocyte uptake by setting (Fig. 6). The result is an increase in GABA extracellular concentration implying an increase in . As increases, the probability for to overcome along the associated Brownian motion decreases, and also does the oscillation frequency (Fig. 6). Consequently, we observe a decrease in the oscillation frequency after . In the time series, the oscillation frequency decreases gradually during a transient () until reaching its minimum. This can be explained by the slow increase of GABA extracellular concentration that reaches its new baseline at .

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 () 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 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 (i.e. the threshold for neural activation) with 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 .

We note in this regard, however, that, for certain parameter sets of the whole model, the SNIC bifurcation may be lost as 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 , 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 always exists and the associated saddle-node bifurcation is not degenerate for any such as , that is the maximum interval of values taken by .

We recall that can be written as follows: where f is given by (7). Since we consider fixed we introduce the function As explained above, for each , there exists a unique bifurcation value occurring at the non-hyperbolic (saddle-node) singular point characterized by , which is defined by This value satisfies . We cannot find the explicit expressions of and . Thus, for characterizing the variations of with , we take advantage of the implicit definitions above and focus on localizing the extrema of .

Proposition 4.1

Assume that, for all such that , exists and the associated saddle-node bifurcation is not degenerate. Note . Then:

if , has no local extremum,

if , may admit two local extrema: a minimum at and a maximum at . If both exist, then .

Let us search for local extrema of the function , which is implicitly defined by (9a)-(9c). Hence, we are interested in the following constrained optimization problems: We introduce the associated Lagrangian function The necessary condition for the existence of an extremum of g under the constraint is that is, By assumption, the saddle-node bifurcation associated with the SNIC bifurcation is not degenerate, i.e. every solution of (11) such as satisfies . Thus, system (12a)-(12c) reads Therefore, if the problem under constraint admits an extremum, at this extremum it is . Following the assumption that a SNIC bifurcation occurs for any value of between 0 and , then Eq. (13c) admits a solution for any as well. Hence, if our constrained problem admits an extremum, this corresponds to the SNIC bifurcation occurring at such that From (7), we obtain Using the facts that, for any fixed value of , function is bell-shaped and its maximum value does not depend on (see Fig. 7), one finds that the function vanishes in if If , the function admits no zero, which proves the first item of Proposition 4.1.

Fig. 7

Interval of values corresponding to a realized local minimum of neural excitability threshold. Graphic representation of the function h defined by (19) and interval for which

Now, we assume that condition (14) is fulfilled and we search the values of satisfying , i.e. Noting , this equation reads Setting we obtain the two solutions of : Note that (resp. ) corresponds to the extremum when the saddle-node (resp. ) crosses the fold of the surface . We consider and we note the value of corresponding to the SNIC connection for this value of , i.e. the solution of To prove that reaches a local minimum at , we introduce the bordered Hessian matrix H̅ associated with the Lagrangian function at its singular point (solution of system (13a)-(13c)): The determinant of is given by On one hand, the saddle-node associated with the SNIC bifurcation is not degenerate and is a local maximum of , thus . On the other hand, for any , is increasing at (see Fig. 7), thus . Finally and corresponds to a local minimum of . A similar argument proves that corresponds to a local maximum of (where is the value corresponding to bifurcation for ). □

The above proposition can be interpreted as a necessary condition for having a change in the sense of variations of when varies in . We now derive a necessary and sufficient condition so that actually lies in . Since satisfies , one obtains, from Eq. (15), For any , function is strictly increasing over , and if and only if where Figure 7 gives an illustration of the function h as well as how and are built.

In conclusion, for a fixed value of , reaches a local minimum at a value such that if and only if Moreover, in Sect. 3 about the astrocyte GABA uptake deficiency, we proved that, for a fixed value of , is linear and increasing with . Both results allow us to predict that there exist three shapes of according to the value of . In the following we illustrate the three qualitative types of neural activity resulting from an astrocyte deficiency to capture glutamate using the following values: For each case, we provide simulations representing the value of in space and time series generated by the model when the glutamate astrocytic uptake is altered (Fig. 8).

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 in plane, and the time series on the right correspond to LFP, , and . The black curves on the colormap are the trace of (, ) along the associated orbits of the model. At , we alter the glutamate astrocytic uptake by setting . The three cases are obtained with: (a) , (b) , (c) . All other parameters are the same in the three cases and given in Table 1

If then and strictly increases with and .

If , then and strictly decreases when increases (for fixed) and strictly increases with (for fixed).

If , then and decreases when increases in and increases with (for fixed).

For (panel (a) of Fig. 8), is negative and thus increases with . Hence in this scenario, when 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 . Since , decreases as increases. Thereby, an astrocyte glutamate uptake deficiency triggers an increase of the 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 is shown in panel (c) of Fig. 8. In this case decreases for 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 . Note that this value can be tuned to obtain small enough and 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.

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 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 [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 , 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 lies in a well-defined interval 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 , neural firing drastically and persistently increases leading to epileptiform activity in agreement with experiment [48]. On the other hand, for 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 and depend on the strength of inhibition of pyramidal neurons by interneurons only further depends on . 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 , 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.