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 τ_m 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: v˙=v2+κ. 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, τ_r, 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_m, G_m and k are constants listed in Table 1. The details of the reduction can be found in Adamchik and Kazantsev (2017).

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_A-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 τ_G = 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).

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 C₀. 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 τ_P 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_m 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_m[1 − exp(−τ/τ_P)]. 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 〈δC〉=∫0∞δC·P(τ)dτ=δCm/(1+AτP). The product A·〈δC〉 gives the required production rate (cf. 6), where Q = δC_m/τ_P.

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, (C₀, 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., C₀ and Q. Note that such parameters as baseline concentation, C₀, and maximum production rate, Q, can be controlled pharmacologically in experiments.

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_GABA = 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_A-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, I₀, which served to regulate the level of depolarization. The parameters of synaptic transmission are listed in Table 2.

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 + τ_s) is the total number of spikes, the entire network generates within the respective time window; τ_w = 20 ms and τ_s = 1 ms are averaging and sliding windows, respectively.

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|_{t = 0} = 0, C|_{t = 0} = C₀) 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.

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, C₀, 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).

Let us now describe the biophysical mechanism underlying periodical oscillatory solutions (see Figure 1B). If the baseline ambient GABA concentration, C₀, 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, τ_m. 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, τ_C (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, C₀, 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|_{t = 0} = 0, C|_{t = 0} = C₊. The condition for the border follows directly from equation (1) as: g_C+−0(0) > 0, g_C++0(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)): E*=Em+Gm2k≈-56.8 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, τ_m, and of concentration, τ_C.

Our next step was to determine the region in the parameter plane (C₀, Q) where the system had periodical solutions. For this purpose we implemented the following calculations. For each point (C₀, Q) from the rectangular (0, C_max) × (0, Q_max) we traced the time evolution of the dynamical system Equations (1, 6) under the same initial conditions: A|_{t = 0} = 0, C|_{t = 0} = 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.

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 C₀ 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, C₀ (see Figures 4A,B). Next, we fixed C₀ 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.

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

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.