Emergence of Relaxation Oscillations in Neurons Interacting With Non-stationary Ambient GABA.

Dynamics of a homogeneous neural population interacting with active extracellular medium were considered. The corresponding mathematical model was tuned specifically to describe the behavior of interneurons with tonic GABA conductance under the action of non-stationary ambient GABA. The feedback provided by the GABA mediated transmembrane current enriched the repertoire of population activity by enabling the oscillatory behavior. This behavior appeared in the form of relaxation oscillations which can be considered as a specific type of brainwaves.

1. Introduction

Historically, the focus of experimental and theoretical studies of brain signaling was almost exceptionally on neurons and their networks. Being the only electrically excitable cells in the nervous system, neurons are able to communicate by receiving, processing and generating electrical signals in the form of spike trains (Nicholls et al., 2001). All other structures constituting the nervous tissue such as glial cells and extracellular matrix (ECM) until very recent decades were not taken into account in the mechanisms of information processing.

Glial cells and various exctracellular structures were primarily thought to perform a number of auxiliary functions such as trophic, supportive and immune (Allen et al., 2009). The comprehension of inalienability of glia and ECM to the neuronal signaling came with the discovery of chemical synaptic transmission machinery (Krnjević, 1974) and secretory function of astrocytes (Martin, 1992). Glia turned out to be a gigantic chemical factory of the nervous system, governing neurons and using the extracellular space as an intermediary (Barres, 1991).

At present, there has been a great number of theoretical and experimental studies devoted to neuron-glia interaction (Bezzi and Volterra, 2001). One of the most prominent concepts in the field was that of the tripartite synapse (Araque et al., 1999). Glial cells, particularly astrocytes, can effectively influence and modulate the synaptic transmission. Many aspects of such modulations were discussed in a number of computational studies (Postnov et al., 2007; Gordleeva et al., 2012; Kazantsev et al., 2012; Volman et al., 2012; Lazarevich et al., 2017).

Besides several glial cell types, the extracellular space itself can be an important player in neuronal signaling. It serves as an interstitial transport system mediating cell-to-cell communications by means of numerous active chemicals (Sykovaá and Nicholson, 2008). This type of communications is called “volume” transmission and is characterized by signal diffusion in a three-dimensional fashion within the brain extracellular fluid (Agnati et al., 1995). The “volume” transmission depends crucially upon the actual geometry of the ECS (Syková, 2004) which has great relevance for pharmacokinetics and actions of neuropsychoactive drugs (Zoli et al., 1999).

One of the major neurotransmitters in the CNS is γ-Aminobutyric acid (GABA) (Webster, 2001). It mediates intercellular communications by participating in both “wiring” and “volume” transmission (Semyanov et al., 2004). The “wiring” action of GABA is through mediating the synaptic transmission by activating the postsynaptic (phasic) GABA-receptors. The “volume” transmission is, in its turn, carried out by “overspilled” ambient GABA which regulates neuronal excitability by creating the extra transmembrane current through extrasynaptic (tonic) GABA-receptors.

GABA was reported to maintain the fast neuronal oscillations (gamma, 20–80 Hz) in inhibitory interneuron networks (Whittington et al., 1995; Bartos et al., 2007). In the computational study (Wang and Buzsáki, 1996), GABA synaptic transmission was shown to provide a suitable mechanism for synchronized gamma oscillations in a sparsely connected network of fast-spiking interneurons. Incidentally, GABA was reported to enhance collective behavior in neuronal axons (Traub et al., 2003). Specifically, gamma-frequency oscillations were demonstrated to coexist with phasic high-frequency oscillations (>90 Hz) in principal cell axon populations.

Primarily, GABA was considered to be the main inhibitory neurotransmitter in the brain until it was shown experimentally to be able to perform the bi-directional regulation of neuronal spiking activity (Song et al., 2011). Based on this experimental finding, a number of mathematical models describing the action of ambient GABA on the excitability properties of interneurons were suggested (Adamchik et al., 2015; Adamchik and Kazantsev, 2017). In Adamchik et al. (2015), the behavior of a single interneuron embedded in the extracellular space with constant ambient GABA concentration was studied. It was shown that depending on the parameters of tonic current, such as tonic conductance density and GABA reversal potential, the interneuron demonstrated different behavioral modes including self-oscillations. The impact of stationary GABA at the population level was studied subsequently in Adamchik and Kazantsev (2017). Specifically, it led to bistability between asynchronous firing and zero-activity state.

In this paper, we study the effects of non-stationary, activity dependent GABA upon population dynamics of interneurons. To this end, we propose a mathematical model accounting for the feedback between interneurons and ambient GABA (section 2). The origin of the feedback has the following explanation. Extracellular GABA creates the additional transmembrane current through activation of extrasynaptic (tonic) GABA-receptors. This current further changes the firing properties of interneurons (Adamchik et al., 2015), which immediately affects the synaptic release of GABA (Destexhe et al., 1994). Since extracellular GABA concentration depends, among others, on spillover, i.e., the diffusion of the neurotransmitter out of the synaptic cleft (Semyanov et al., 2004), it changes, which futher affects tonic conductance and provides the respective feedback.

Based on these considerations, we proposed a mathematical model using the following assumptions. First, we considered a particular case of a homogeneous population of interneurons which allowed us to describe their collective behavior using the simple rate-based formalism. Second, we neglected any spatial gradient of neurotransmitter, considering its concentration to be uniformly distributed over the entire extracellular space. This assumption allowed us to build the minimal model of the feedback avoiding dealing with an explicit model of spatiotemporal GABA dynamics. The model consisted of two coupled equations, one of which described the time-course of population activity (section 2.1) while the other–the concentration of ambient GABA (section 2.2). The model was explored both numerically and analytically (section 3). The results including the appearance of relaxation oscillations were discussed in (section 4).

2. Materials and methods

2.1. Population dynamics

Within the framework of rate-based formalism, a homogenous population of neurons is described by a single variable, e.g., the population activity, A. The rate of change of A is determined by the so called gain fuction gλ(I), which is unique for each cell type. The respective equation reads:

where τ is the membrane time constant and I is the total input current an arbitrary neuron receives from the entire network. The latter is linearly dependent on population activity, I = JA, where the proportionality factor, J, is called coupling strength (Gerstner et al., 2014).

The gain gλ(I) is primarily a function of input current, I, but can also depend on a number of factors, collectively denoted here by λ. In our case, these are tonic conductance density, G, and GABA reversal potential, E, i.e., λ = (E, G).

The exact form of the gain function can be derived analytically only for a few simple neuron models, such as, for example, the quadratic integrate-and-fire (QIF), whose dimensionless normal form reads: . Using separation of variables and integration over infinite potential bounds, one can get:

where the dimensionless parameter κ depends both on input current I and other factors. The absolute refractory period, τ, is added to the period of oscillations to prevent the firing frequency from taking an arbitrarily large value. Note, that formula (Equation 2) is valid only for positive κ; when κ < 0 no oscillations occur and, as a result, gλ(I) = 0.

The gain function (Equation 2) describes qualitatively the responce of Class I excitability neurons. In these neurons, the transition from resting to spiking occurs via saddle-node on invariant circle bifurcation (SNIC), that allows them to fire with arbitrarily small frequency (Izhikevich, 2007). The Wang-Buzsáki interneuron (Wang and Buzsáki, 1996) belongs exactly to this type of neurons. In Adamchik and Kazantsev (2017), the original conductance-based model, modified in a way to account for the additional transmembrane tonic current, was reduced to the QIF neuron. The dimensionless parameter κ took the following form:

where E, G and k are constants listed in Table 1. The details of the reduction can be found in Adamchik and Kazantsev (2017).

Table 1

Model constants.

Constant	Value	Description
τm	8.925 ms	Membrane time constant
τr	0.627 ms	Absolute refractory period
Gm	0.112 mS · cm−2	Conductance density at threshold
Em	−60.414 mV	Halfway between resting and
		threshold membrane potentials
k	0.0155 μA · cm−2 mV−2	Proportionality factor
α	5 mmol−1 ms−1	Forward rate
β	0.18 ms−1	Backward rate

Equation (1) along with the relations (Equations 2, 3) describes the time-course of population activity. It contains parameters such as coupling strength, J, and GABA reversal potential, E, which can take arbitrary values but remain unchanged. Tonic conductance density, G, is on the contrary a variable, which depends on local GABA concentration, C. The form of this dependence can be determined using a common kinetic formalism (Destexhe et al., 1998). According to a simplified kinetic scheme of the GABA-receptor, which is assumed to exist in two conformations, open (O) and closed (C), one can get:

where α and β are forward (activation) and backward (deactivation) rates, respectively, and Ḡ is maximum conductance density. Rate values were taken from Koch and Segev (1998) and are given in Table 1 for reference.

Note, that equation (4) does not determine a momentary but rather a stable-state value of G. Nevertheless, we may use it because conductance relaxation time τ = 1/(αC + β) ≤ 1/β ≈ 5 ms, which is far less than the operating time of ambient GABA concentration which amounts to hundreads of milliseconds (Semyanov et al., 2004).

2.2. Ambient GABA dynamics

According to Semyanov et al. (2004), extracellular GABA concentration is regulated by uptake, non-synaptic release and spillover. Uptake is carried out by GABA transporters which decrease the concentration by binding and removing GABA molecules from the extracellular space. Ambient GABA can originate from various sources. It can escape from synaptic cleft (spillover) and can be released via non-vesicular mechanism by neurons and glia. Both spillover and non-synaptic release increase ambient GABA concentration but do it differently. Unlike non-synaptic release, spillover depends crucially on synaptic dynamics and, as a result, on population activity. These general considerations allowed us to write a governing equation for ambient GABA concentration:

The first term is supposed to describe the mutual action of uptake and non-vesicular release, which counterbalance each other by maintaining an optimal background GABA level, denoted here by C0. Spillover is described, in its turn, by the second term, S(A), representing the production function of GABA and depending explicitly on population activity. The exact form of the production function was chosen in a way to describe qualitatively correctly the properties of synaptic neurotransmitter release. In the most common case, the production function reads as follows:

where τ is GABA production time constant and Q is the maximum production rate.

The exact form of the production function can be derived based on the following consideration. Let δC be the amount of GABA released in response to a single spike. Then, due to exhaustion of synaptic vesicle pools, the next spike, coming τ time units after the first one, will evoke the release of a lesser amount of neurotransmitter, precisely δC = δC[1 − exp(−τ/τ)]. For a Poisson spike train, the interspike interval distribution (ISI) with the mean firing rate equal to A reads: P(τ) = A·exp(Aτ). Then, the average amount of GABA released in response to a spike from the spike train will be . The product A·〈δC〉 gives the required production rate (cf. 6), where Q = δC/τ.

2.3. Dynamical system

Two coupled ordinary differential equations (ODE) (1, 5) with relations (2–4, 6) form a 2D dynamical system. Its state variables are population activity, A, and ambient GABA concentration, C. Besides some constants (see Table 1), the equations contain a number of free parameters, which can roughly be split into two distinct groups. The first group, (E, Ḡ, J), consists of the parameters controlling population activity, while the second one, (C0, Q), determines ambient GABA concentration. Our task is to reveal how the dynamics of equations (1, 5) depend on all these parameters. Some preliminary considerations concerning the matter are the following.

In absense of the second equation (Equation 5), the system reduces to a simple 1D phase line, corresponding to the case of stationary external medium. This particular case was a subject of our previous study (Adamchik and Kazantsev, 2017). It was shown then that introduction of tonic current did not lead to any new dynamical effects compared to the reference case, (G = 0), characterized by a simple stable-state dynamics. It resulted, however, in appearance of a monostable regime of asynchronous firing once tonic current parameters, E and G, were properly tuned, specifically, in a way that the point (E, G) was located above a certain curve on the parameter plane. This curve was shown previously to be a border of the self-oscillatory mode in the model of a tonically driven single neuron (Adamchik et al., 2015). The impact of coupling strength, J, on population activity consisted, in its turn, in inducing bistability, i.e., coexistence of resting and asynchronous firing, and expanding the bistability region at the cost of trivial (zero-activity) monostable solutions.

In this study our particular interest is on the feedback-induced dynamical effects. For this purpose we specifically focus on the parameters governing ambient GABA concentration, i.e., C0 and Q. Note that such parameters as baseline concentation, C0, and maximum production rate, Q, can be controlled pharmacologically in experiments.

2.4. Semi-explicit model

To verify the validity of our population model prediction as well as to visualize neural oscillations we built a semi-explicit computational model of the respective spiking network. To do this we replaced the equation for the population activity (Equation 1) with an explicit spiking neural network model but preserved (Equation 5), which describes the dynamics of ambient GABA.

We considered specifically a network of N = 100 interneurons randomly coupled with probability p = 0.1. Each neuron was described by the original conductance-based model (Wang and Buzsáki, 1996) with the additional tonic current term: I = G(u − E), where u is the membrane potential. Tonic conductance density, G, depended here on extracellular GABA concentration just in accordance with equation (4). Synaptic transmission was mediated by phasic GABA-receptors. The total synaptic current to an arbitrary neuron was determined as the normalized sum over the contributions of all its presynaptic neighbors:

where M is the mean number of presynaptic inputs: M = Np. The fraction of the receptors in the open state, r, obeyed the kinetic equation

while the synaptic concentration of GABA, T, strictly followed the potential at the presynapse:

Besides transmembrane and synaptic components, the total current to each neuron included also a constant one, I0, which served to regulate the level of depolarization. The parameters of synaptic transmission are listed in Table 2.

Table 2

Model synapse.

Parameter	Value	Description
Gsyn	0.1 mS · cm−2	Maximum phasic (synaptic) conductance density
Tmax	1 mmol	Peak synaptic cleft concentration of GABA
Θ	0 mV	Threshold of the activation
Δ	2 mV	Width of the transition area
I0	−1.4 μA cm−2	Depolarization current

To couple the explicit spiking network model with the equation describing ambient GABA dynamics we calculated at each time step the instantaneous population activity, using for averaging the α-function (Dayan and Abbott, 2001):

Here, S(t − s; t − s + τ) is the total number of spikes, the entire network generates within the respective time window; τ = 20 ms and τ = 1 ms are averaging and sliding windows, respectively.

3. Results

3.1. Numerical and phase plane analysis

First, we performed numerical analysis of the dynamical system (Equations 1, 5). We carried out numerous simulation trials for a wide range of biologically relevant parameter values. In each trial, the system proceeded with the same initial conditions (A| = 0, C| = C0) corresponding to zero population activity and baseline GABA concentration, respectively. We found out that, depending on parameters, the system demonstrated either oscillatory or stationary behavior (see Figure 1). The parameters of the simulations are listed in Table 3.

Figure 1

The time-course of population activity (solid) and ambient GABA concentration (dashed) along with the respective phase plane trajectories (thick solid) for oscillatory (A,B) and stationary (C,D) modes. On the phase portraits (B,D), thin gray lines correspond to A- (solid) and C- (dashed) nullclines, respectively. Separatrices of the saddle are represented by a pair of thin black lines with arrows indicating the flow direction.

Table 3

Model parameters.

Parameter	Value	Description
J	50 ms ·μA cm−2	Coupling strength
E	−50 mV	GABA reversal potential
Ḡ	1 mS · cm−2	Maximum tonic conductance density
τC	100 ms	GABA relaxation time constant
τP	100 ms	GABA production time constant
C0	0.05 mmol	Baseline ambient GABA concentration
Q	0.02 mmol · ms−1 (oscillations)	Maximum GABA production rate
	0.01 mmol · ms−1 (stationary)

Different kinds of behavior can be accounted for using the phase plane. In case of oscillations, the trajectory first makes a big loop before converging to the limit cycle (see Figure 1B). The latter is intersected by the A-nullcline, coinciding for the small A's with the border between zero and non-zero gain: gλ(0) = 0 (see Figure 2B). It means that the system in the oscillatory mode sequentially visits the region of excitatory GABA action. The oscillations would evidently not occur if the baseline GABA concentration, C0, exceeded the borderline value between inhibition and excitation, C+. In other words, if the system was placed into the region of inhibitory GABA, it would never leave it. In the stationary mode, the trajectory moves up slower than the A-nullcline does (see Figure 1D), so it converges to the fixed point corresponding to stationary asynchronous firing. This scenario realizes if the maximun production rate, Q, is lower than a certain threshold value. These considerations helped us subsequently to determine conditions for oscillations (section 3.2).

Figure 2

Stationary relaxation oscillations in the model of the feedback between a population of interneurons and GABA-containing extracellular medium: (A) the time-course of the dynamical variables during one period of oscillations–population activity (solid) and ambient GABA concentration (dashed); (B) the corresponding phase plane trajectory in the form of a limit cycle. The unstable fixed point (denoted by the empty circle) lies close to the limit cycle.

Let us now describe the biophysical mechanism underlying periodical oscillatory solutions (see Figure 1B). If the baseline ambient GABA concentration, C0, is high enough to make neurons fire but not too high to inhibit them by shunting, the initially silent neurons start firing. Non-zero population activity makes activity dependent ambient GABA concentration steadily grow up through synaptic release and spillover (OM). If parameters, governing ambient GABA dynamics, are properly tuned, then, at some point, tonic GABA switches from excitation to inhibition. On the phase plane, it corresponds to the intersection of the trajectory with the A-nullcline (M). As soon as the intersection occurs, the gain becomes zero and population activity starts decreasing to its steady-state (zero) value with a time constant of the membrane, τ. While population activity goes down, the concentration keeps growing but its grow rate gradually slows down (MN). Eventually, the rate of change of C becomes negative and the trajectory moves down with a time constant of concentration, τ (NO′). At some point (O′), the trajectory re-enters the region of excitatory GABA and the entire process starts from the scratch. Note, that concentration does not reach its baseline level, C0, so the magnitude of the limit cycle is less than that of the initial loop.

To avoid the trajectory from making a loop before converging to the limit cycle, we took the initial conditions exactly at the upper border between inhibitory and excitatory GABA: A| = 0, C| = C+. The condition for the border follows directly from equation (1) as: g(0) > 0, g(0) = 0. Based on the explicit analytical expression for the gain function (Equations 2, 3) as well as on the relation between tonic conduction density, G, and GABA concentration, C, (Equation 4), one can get:

where

The minus-subscripted concentration, C−, corresponds to the transition from inhibition to excitation as we move upwards the C-axis and is given here just for reference. Its value for the actual choice of parameters (see Table 3) is negligible and its existence does not play any substantial role for oscillations. Oscillations occur essentially at the border between excitatory and inhibitory GABA and not vice versa.

Note, that the solution exists only if GABA reversal potential, E, lies above a certain threshold, whose value is determined by zero determinant condition (see equation (11)): mV, which exceeds the resting membrane potential (≈−64 mV) by 9.2 mV.

The time-course of the dynamical variables as well as the shape of the limit cycle (see Figure 2) are typical for relaxation oscillations. In the excitatory region, the population activity relaxes to the value determined by the gain function but as soon as it leaves it, A starts the exponential decay to zero. The concentration follows the population activity with a delay caused by the difference between the time constants of membrane, τ, and of concentration, τ.

3.2. Conditions for oscillations and their characterisctics

Our next step was to determine the region in the parameter plane (C0, Q) where the system had periodical solutions. For this purpose we implemented the following calculations. For each point (C0, Q) from the rectangular (0, C) × (0, Q) we traced the time evolution of the dynamical system Equations (1, 6) under the same initial conditions: A| = 0, C| = C+. Depending on the parameters, the system either remained at the starting point or relaxed to the upper stable state, or oscillated periodically. The periodicity was ascertained based on multiple crossings of the trajectory with the A-nullcline. The border between the regions of oscillatory and transient solutions is depicted in Figure 3 for two different values of GABA reversal potential. Specifically, Figure 3A corresponds to highly depolarizing GABA, while Figure 3B corresponds to GABA, whose reversal potential is slighly above (precisely by 1.8 mV) the oscillatory threshold.

Figure 3

Oscillatory region on the plane (gray) (C0, Q) for (A) E = −50 mV, (B) E = −55 mV.

The oscillatory region for highly depolarizing GABA has roughly the shape of a semi-infinite strip, (0, C+) × (Q−, ∞), where the border values C+ and Q− are both dependent on the reversal potential, E. As far as we get closer to the oscillatory threshold, the region shrinks until collapsing at E = −56.8 mV. The shape of the oscillatory region for high E's implies that oscillations occur if both the baseline GABA concentration and the GABA production rate are located below and above their respective threshold values: C < C+, Q > Q−. Note, that the upper bound C+ corresponds exactly to the concentration value at the starting point of simulation. The shape of the oscillatory region has the straightforward phase-plane interpretation (see section 3.1). Although, we managed to find the explicit analytical expression for C+ (see Equations 10, 11), there was no way to obtain such for Q−(E), the more so it depends not only on E but on C+ as well, which is illustrated in Figure 3B.

Having found the oscillatory region we looked for the magnitude and period of oscillations as functions of C0 and Q. To this end we started time-course simulations from the same point at the limit cycle as we did before in oder to pass the transition phase. We defined the period of oscillations as the time before two subsequent intersections with the A-nullcline with the same sign of the slope. The magnitude of oscillations was determined as the maximum value of population activity and concentration, respectively. The results are depicted in Figure 4. First, we fixed the maximum production rate, Q, and found numerically the dependence of both period and magnitude on the baseline concentration of GABA, C0 (see Figures 4A,B). Next, we fixed C0 and found the respective dependencies on the parameter Q; the corresponding graphs are depicted in Figures 4C,D. Note, that the period grows infinitely as soon as we approach the boundaries of the oscillatory region. The dependence of the magnitude on the production rate is more drastic than that on the background GABA level. The explanation is quite simple – the higher Q is, the faster the periodic trajectory intersects the A-nullcline, the lower is the magnitude.

Figure 4

(A,C) Period and (B,D) magnitude of relaxation oscillations (A–solid, C–dash) vs. baseline GABA concentration and maximal GABA production rate. The parameters are: E = −50 mV, (A,B) Q [mmol · ms−1]: 0.02 (black), 0.06 (dim gray), 0.1 (silver); (C,D) C0 [mmol]: 0 (black), 0.03 (dim gray), 0.06 (silver).

3.3. Spiking network simulation

The results of our computer simulation are depicted in Figure 5. Note that they quite correctly reproduce those obtained using the original rate model. In case of periodicity, after an initial burst of population activity, the network demonstrates stationary oscillatory behavior (Figure 5C) just in accordance with the prediction (Figure 1A). The stationary-like behavior in the spiking network model was another option that we can verify at the network level (see Figure 5D).

Figure 5

Spike raster plots and time-course of population activity for (A,C) oscillatory and (B,D) stationary modes. The parameters are: C0 = 0.05 mmol, Q [mM/ms]: (A,C) 0.01, (B,D) 0.003.

4. Discussion

We proposed a self-consistent model of interneurons interacting with extracellular, activity dependent GABA. The model represented two coupled nonlinear ODEs describing the dynamics of population activity and ambient GABA concentration, respectively. To write the first equation we used the well-known Wilson-Cowan formalism describing the low-pass behavior of a neural ensemble (Gerstner et al., 2014). The gain function was chosen in a way to properly mimic the behavior of the interneuron with tonic GABA conductance (Adamchik and Kazantsev, 2017). The dynamics of ambient neurotransmitter were quantitatively accounted for on the basis of empirical evidence about the sources and sinks of extracellular GABA (Semyanov et al., 2004). Mathematically, the model was a continuous-time dynamical system on a plane. It was shown to admit both stationary and periodic solutions depending on the parameters governing neurotransmitter concentration. Unlike the stationary-like behavior, periodicity was a feedback-induced feature with clear biophysical explanation. In oscillatory mode, the system evolved between the regions of excitatory and inhibitory GABA. In each of these regions the dynamical variables relaxed to their stable-state values, so that the type of the oscillations was essentially relaxational. Such a pattern of synchronized population activity can be regarded as a specific type of brainwave.

We determined the conditions for oscillations and their characteristics such as period and magnitude as a function of GABA parameters. In particular, we found out that oscillations were possible only for strongly depolarizing GABA. For interneurons, GABA reversal potential had to exceed the RMP by at least 9.2 mV for oscillations to occur, which is above the reported values (Michelson and Wong, 1991; Verheugen et al., 1999; Chavas and Marty, 2003; Banke and McBain, 2006). The oscillatory region on the plane of baseline GABA concentration and maximum GABA production rate had roughly the shape of a semi-infinite stripe, i.e., there existed an upper background level of GABA and a lower intensity of its production, beyond which no oscillations occured. As a consequence, neural oscillations could be induced or suppressed pharmacologically, by changing GABA control parameters. This might have some biomedical implications since extrasynaptic GABA is believed to contribute to epileptic or schizophrenic brain activity (Brickley and Mody, 2012).

In conclusion, we need to mention the limitations of our present consideration. The suggested model of the feedback between neurons and extracellular GABA is minimal in the sense that it does not account for many key features of real neural networks and their environment. For example, when discussing ambient GABA dynamics (section 2.2) we assumed GABA sources, sinks, and receptors to be co-local. This allowed us to describe the time evolution of ambient GABA concentration with a simple ODE instead of building a detailed model accounting for the actual geometry of the extracellular space. Futher, we assumed uptake and non-vesicular release independent on population activity, although there is experimental evidence that this traditional view was too simplistic (Richerson and Wu, 2003). In addition, we considered the special case of a homogeneous neural network which is a rough representation of real neuronal ensembles. Taking into view all these considerations, we must admit that our conclusions can offer only primary insights into the feedback-induced dynamics of GABA-driven interneurons. At the same time, they can be regarded as reference point for future studies applying more sophisticated methods.

Author contributions

DA designed the model and the computational framework, carried out the implementation and wrote the manuscript. VM contributed to the interpretation of the results. VK conceived the study and was in charge of overall direction and planning. All authors provided critical feedback and helped shape the research, analysis and manuscript.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer CB and handling Editor declared their shared affiliation.