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

HAL
open science

# Biologically realistic mean field model of spiking neural networks with fast and slow inhibitory synapses

Claudio Di Geronimo, Alain Destexhe, Matteo Di Volo

► To cite this version:

Claudio Di Geronimo, Alain Destexhe, Matteo Di Volo. Biologically realistic mean field model of spiking neural networks with fast and slow inhibitory synapses. Journal of Computational Neuroscience, 2025, 53 (1), pp.289-303. 10.1007/s10827-025-00904-7. hal-05084112

HAL Id: hal-05084112

https://hal.science/hal-05084112v1

Submitted on 26 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Biologically realistic mean field model of spiking
neural networks with fast and slow inhibitory
synapses

Claudio Di Geronimo<sup>1,2</sup>, Alain Destexhe<sup>3</sup>, Matteo di Volo<sup>1*</sup>

<sup>1</sup>Université Claude Bernard Lyon 1, Institut National de la Santé et de la Recherche Médicale, Stem Cell and Brain Research Institute U1208, Bron, France.

<sup>2</sup>Dipartimento di Fisica, Universita di Firenze, Via G. Sansone 1,
I-50019 Sesto Fiorentino (FI), Italy.

<sup>3</sup>Paris-Saclay University, CNRS, Institute of Neuroscience (NeuroPSI),
Saclay, France.

*Corresponding author(s). E-mail(s): matteo.di-volo@univ-lyon1.fr;

Abstract

We present a mean field model for a spiking neural network of excitatory and inhibitory neurons with fast GABAA and nonlinear slow GABAB inhibitory conductance-based synapses. This mean field model can predict the spontaneous and evoked response of the network to external stimulation in asynchronous irregular regimes. The model displays theta oscillations for sufficiently strong GABAB conductance. Optogenetic activation of interneurons and an increase of GABAB conductance caused opposite effects on the emergence of gamma oscillations in the model. In agreement with direct numerical simulations of neural networks and experimental data, the mean field model predicts that an increase of GABAB conductance reduces gamma oscillations. Furthermore, the slow dynamics of GABAB synapses regulates the appearance and duration of transient gamma oscillations, namely gamma bursts, in the mean field model. Finally, we show that nonlinear GABAB synapses play a major role to stabilize the network from the emergence of epileptic seizures.
# 1 Introduction

Neural networks in the cerebral cortex exhibit a wide range of dynamical regimes, from asynchronous dynamics to synchronous irregular regimes characterized by oscillations at the scale of a population of neurons [1, 2]. Oscillations at the population scale, detected through mesoscopic signals such as Local Field Potentials (LFPs), are typically characterized by irregular and low firing activity of neurons at the microscopic scale [3] and have been shown to play an important role in brain computation. In particular, gamma oscillations (40-120 Hz) have been found to be implicated in cognitive functions, more specifically in attention [4, 5]. At the same time, slower oscillations such as theta oscillations (4-12 Hz) [6] have been shown to coordinate attention selection and memory [7, 8]. Moreover, the cross-frequency coupling between these different rhythms controls the functional connectivity and the large scale interactions in the brain [9]. Nevertheless, it is still unclear today what are the circuit mechanisms behind the different oscillations observed in cortical networks' dynamics.

Low dimensional mean field models are powerful tools to propose potential circuit mechanisms regulating neural networks’ dynamics and stability. Characterized by few differential equations for population variables, mean field models define a formal link between the microscopic structure and the mesoscopic emergent dynamics of neural networks. Moreover, being low dimensional, mean field models also serve as the building block to model large scale dynamics at the whole brain scale [10]. Since the pioneer work by Wilson and Cowan [11] there has been intense research in the area of mean field models applied to the cortical dynamics [12–16]. While the first mean field models were purely phenomenological, many efforts have been made to develop mean field models that actually represent the population dynamics of a spiking neural network. This propriety is important because it allows us to validate mean field model’s predictions. Nevertheless, given the complexity and nonlinearities of networks of spiking neurons, such dimensional reduction procedure is mathematically challenging and mean field models need to rely on assumptions on the microscopic structure of the network. These assumptions often limit the inclusion of biologically realistic ingredients such as network’s structure, neurons diversity or conductance-based interactions. Recently, an exact mean field model has been introduced for Quadratic integrate and fire neurons [17]. This approach has been shown to be able to model gamma oscillations and theta-gamma coupling [18], including the bursting features of gamma oscillations [19]. Nevertheless, this approach is limited to QIF neurons, all-to-all networks, and Cauchy distributions of heterogeneities. Recent works have shown that new dynamical regimes emerge instead in sparse networks [20, 21] or in the presence of non-Cauchy distributions [22].

In this paper we consider instead networks with biologically realistic ingredients, such as sparse network connections, conductance-based interactions and different spike shape between regular spiking excitatory and fast spiking inhibitory neurons. We refer to our model as biologically realistic, because it integrates these ingredients. In order to derive a mean field model of this network we employ a different approach, based on master equation methodology [14] that has been shown to be extendable to different neurons types (from Hodgkin Huxley to Morris-Lecar) [23], to account for any distribution of heterogeneity [24] and to predict the population dynamics of a sparse
network of spiking neurons [25]. This approach has been successfully employed for modeling asynchronous dynamics, slow oscillations [24] and recently gamma oscillations with the addition of synaptic delay [26]. Nevertheless, these previous mean field models consider only two types of synapses, fast AMPA and $GABA_A$ receptors.

In this work we are interested in the role of slow and nonlinear inhibitory synapses such as $GABA_B$ receptors. We thus developed a new mean field model including the effect of nonlinear slow $GABA_B$ receptors. We have employed an adiabatic approximation, already tested effective in this model for spike frequency adaptation [25], and we have performed mean field approximation to account for the experimentally observed nonlinearity of $GABA_B$ receptors. This method allows to account for nonlinear slow inhibitory synapses. We show that this mean field model can reproduce the network response to external stimulation and that it predicts the emergence of theta oscillations at 4-5 Hz appearing for sufficiently strong $GABA_B$ conductance. Interestingly, this new mean-field model predicts distinct effects between the direct stimulation of inhibitory neurons and the enhancement of $GABA_B$ receptors in the emergence of gamma oscillations. Indeed, while the stimulation of inhibitory neurons does increases the power of gamma oscillations, an increase of $GABA_B$ conductance eliminates gamma oscillations. This results shows that increasing inhibition can have different effects, depending on how inhibition is increased in a network of neurons. These predictions are consistent with direct network simulations of spiking neural networks and with experimental observations. Indeed, $GABA_B$ agonist baclofen has been shown to suppress gamma oscillations in hippocampal slices [27]. By employing electrical stimulation, it is in fact found that synaptically released GABA was capable of modulating gamma frequency oscillations via activation of $GABA_B$ receptors [27]. On the contrary, optogenetic stimulation of parvalbumin-positive (PV) interneurons has been shown to increase the power of gamma oscillations [28, 29]. We also employ our mean field model to model gamma bursts, i.e. periods of few gamma cycles separated by asynchronous activity [30]. We show, through our mean field model, that $GABA_B$ dynamics is strongly linked with gamma oscillatory dynamics. In fact, the appearance and the duration of bursts of gamma oscillations are driven by the temporal evolution and the decay of $GABA_B$ receptors. Finally, we show that the nonlinearity of $GABA_B$ synapses (activating only when a certain amount of presynaptic spikes arrive in synchrony [31]) allows the network to stabilize and avoid epileptic seizures.

In summary, this mean field model is a good candidate to model different dynamical regimes in mesoscopic signals, from asynchronous dynamics to slow and fast oscillations. Importantly, it can be used to test pharmacological action on different inhibitory receptors at the mesoscopic scale and, when integrated in large scale models of the whole brain, it can produce valuable predictions on the effect of such manipulation for large scale activity in the brain and epileptic spreading.

# 2 Materials and Methods

## 2.1 Spiking neural network model

We consider a network of $N=10000$ Adaptive Exponential Integrate and Fire neurons. $N_e = 8000$ are excitatory Regular Spiking (RS) neurons and the remaining $N_i = 2000$
are inhibitory Fast Spiking (FS). Neurons are randomly connected over a sparse network with a probability of connection between two neurons of p = 5%. The temporal evolution of neurons’ membrane potential is described by the following differential equations:

$$C_m \frac{dv_k(t)}{dt} = g_L(E_L - v_k(t)) + \Delta e^{\frac{v_k(t)-v_{thr}}{\Delta}} - w_k(t) + I_{syn}(t) + I_{ext} \quad (1)$$

$$\frac{dw_k(t)}{dt} = -\frac{w_k(t)}{\tau_w} + b \sum_{t_{sp}(k)} \delta(t - t_{sp}(k)) + a(v_k(t) - E_L) \quad (2)$$

At spike time ($v_k > v_{thr}$): $v_k \to v_{rest}, w_k \to w_k + b$,

where $C_m = 200pF$ is the membrane capacity, $v_k$ is the voltage of neuron k, and whenever $v_k > v_{thr} = -30mV$, $v_k$ is reset to the resetting voltage $v_{rest} = -65mV$ and fixed to that value for a refractory time of $T_{refr} = 5ms$. The variable $w_k$ is a self generated current that simulates the adaptation of a neuron after a spike. Due to physiological insight only excitatory neurons have adaptation with values of $a = 0, b = 20pA, \tau_w = 500ms$ ($a = b = 0$ for all inhibitory neurons). The leak conductance used is $g_l = 10nS$ and $E_l = -65mV$ is the leakage reversal potential. The exponential term is slightly different for RS ($\Delta = 2mV$) and FS neurons ($\Delta = 0.5mV$) while the effective threshold is the same for all neurons: $v_{thr} = -50mV$. $I_{ext}$ describes an external current that is set to zero if not stated otherwise. We will employ this in Sec. 3.5 to manually activate inhibitory neurons as in optogenetic experiments. $I_{syn}$ describes all types of presynaptic inputs. In particular, in conductance-based models, such a current is generated by changes in conductances when presynaptic neurons fire:

$$I_k^{syn} = (E_{ampa} - v_k)G_{ampa,k} + (E_{gabaa} - v_k)G_{gabaa,k} + (E_{gabab} - v_k)G_{gabab,k}. \quad (3)$$

In our model we will consider 3 types of receptors that can change the conductances of the membrane with the opening of associated ionic channels. $G_{ampa}$ is increased by incoming excitatory spikes, while $G_{gabaa}$ and $G_{gabab}$ are increased by spikes of presynaptic inhibitory neurons. Each channel has its own reversal potential: $E_{ampa} = 0mV$, $E_{gabaa} = -80mV$ and $E_{gabab} = -90mV$. The time evolution of individual conductances is modeled as a decaying exponential function increased of a fixed amount $Q_x$ at each presynaptic spike:

$$G_{x,k}(t) = g_x \sum_{spike,pre} \Theta(t - t_{sp} - \tau_D) e^{-\frac{t-t_{sp}-\tau_D}{\tau_x}}, \quad (4)$$

with $x = ampa, gabaa$ or $gabab$, $\Theta$ is the Heaviside function and $\tau_D$ the synaptic delay. In this work we consider either instantaneous synapses (Section 3.1-3.4), i.e. with a delay of synaptic transmission $\tau_D$ set to 0ms, either homogeneous delay $\tau_D = 3ms$ (Section 3.5-3.6). It is known that $GABA_B$ receptor have non linear activation. Indeed,
receptors are fully open only when a sufficient number of presynaptic spikes occur within short time intervals, typically 4-10 spikes in 100ms [31]. In order to include such nonlinearity, we consider that the conductance $g_{gabab}$ depends on the amount of inhibitory inputs received in the last 100ms, that we approximate as $0.1K_i\nu_i$, where $K_i = pN_I$ is the average number of inhibitory presynaptic connections and $\nu_i$ ($\nu_e$) the inhibitory (excitatory) population firing rate. More specifically, following [31] we consider the fraction of open receptors $\mathcal{G}$ as

$$ \mathcal{G}(\nu_i) = \frac{1}{1 + e^{-\frac{0.1K_i\nu_i - 7}{1.4}}}, \qquad (5) $$

and the conductance of $GABA_B$ synapses is thus $g_{gabab} \cdot \mathcal{G}(\nu_i)$, instead of simply $g_{gabab}$. We determine the population activity of one cell type as

$$ \nu_{e/i}(t) = \frac{1}{\Delta t} \frac{n_{act}^{e/i}(t, t + \Delta t)}{N_{e/i}}, \qquad (6) $$

where $N_{e/i}$ represents the number of neurons within the population. The term $n_{act}^{e/i}(t, t + \Delta t)$ signifies the cumulative spike count across all neurons in the group within the time window $(t, t + \Delta t)$. In our simulations we use $\Delta t = 10$ms. Synaptic parameters here employed are the same as in [25] (apart from those of $GABA_B$ synapses). The $GABA_B$ conductance $g_{gabab}$ varies in our simulations. The constants of decay for the conductances are $\tau_{ampa} = \tau_{gabaa} = 5$ms and $\tau_{gabab} = 180$ms if not stated otherwise. The quantal conductances are $g_{ampa} = 1$nS for RS to FS coupling, $g_{ampa} = 1.05$nS for RS to RS coupling and $g_{gabaa} = 5$nS for both FS to RS and FS to FS coupling. Each neuron receives an external drive, modeled by considering $K_e = 400$ RS neurons generating action potentials with a Poissonian statistic at a fixed frequency rate $\nu_{ext} = 2.5$Hz. In Section 3.3 we have instead employed a time varying $\nu_{ext}$. First we used a rising and decaying input of the form:

$$ \nu_{ext}(t) = \nu_{baseline} + A \left( \Theta(t_{end} - t) e^{-\frac{(t-t_{end})^2}{2T_1^2}} + \Theta(t - t_{rise}) e^{-\frac{(t-t_{rise})^2}{2T_2^2}} \right), \quad (7) $$

with $\nu_{baseline} = 2.5$Hz, $T_1 = T_2 = 20$ms and $A = 6$Hz. Then we have also employed a simple sinusoidal wave form for $\nu_{ext}$ (see Sec. 3.3).

## 2.2 Mean Field Model

In order to derive the mean field model for this spiking neural network we adopt the master equation formalism first introduced in [32], based on a markovian assumption such that the network dynamics is memoryless over a time scale $T$. Nevertheless, as shown in [25] this assumption can be employed also in spiking neural networks where there are slower mean field processes such as spike frequency adaptation, by employing an adiabatic approximation. In our network we have two relatively slow processes: spike frequency adaptation and synaptic $GABA_B$ conductances. We thus employ an adiabatic approximation for both processes. The equation for the population average $GABA_B$ conductance, namely $G_b$, is obtained by taking the derivative of Eq. (4) and
then by applying the population average. All together, the first order equations for this mean field model are the following:

$$T \partial_t \nu_e = F_e(\nu_e(t - \tau_D), \nu_i(t - \tau_D), W, G_b(t - \tau_D)) - \nu_e \quad (8a)$$

$$T \partial_t \nu_i = F_i(\nu_e(t - \tau_D), \nu_i(t - \tau_D), W, G_b(t - \tau_D)) - \nu_i \quad (8b)$$

$$\partial_t W = -\frac{W}{\tau_w} + b\nu_e + a(\mu_v^e - E_L) \quad (8c)$$

$$\partial_t G_b = -\frac{G_b}{\tau_{gabab}} + K_i \nu_i \mathcal{G}(\nu_i) g_{gabab}, \quad (8d)$$

where $K_i$ is the average number of presynaptic FS neurons, $F_{e,i}$ are respectively the stationary transfer functions of RS and FS neurons. They express the output firing rate as a function of the rates of incoming excitatory and inhibitory neurons, the level of adaptation $W$ and the level of $GABA_B$ conductance $G_b$. The impact of synaptic delay has been introduced following [26]. In the next chapter we will describe how to evaluate these transfer functions $F_{e,i}$, their dependence on $GABA_B$ conductance $G_b$ and the mean population voltage $\mu_v^e$ of excitatory/inhibitory neurons using a semianalytical derivation. In the following analysis we will use $T = 10ms$ if not stated otherwise. We will often refer to stationary version of this mean field (i.e. stationary mean field), which is the version where we do not consider the temporal evolution of $GABA_B$ synapses and $G_b(t) = \tau_{gabab} K_i g_{gabab} \nu_i(t)$, just like we do for fast $GABA_A$ synapses.

## 2.3 Neurons transfer function (TF)

The theoretical analysis of the TF of RS and FS neurons is based on the shot-noise theory in [33], adapted to RS and FS neurons in [34, 35]. The assumption behind this approach is that it is possible to write the output firing rate of a neuron as a function of the subthreshold voltage statistic. We will estimate the average subthreshold voltage $\mu_v$, the standard deviation $\sigma_v$ and the time correlation decay $\tau_v$ of neurons’ membrane potential. First we will report how to evaluate these quantities, and then how to relate them to the firing rate.

### 2.3.1 From Input Rates to Subthreshold Voltage Moments

The mean membrane potential is obtained following the Methods described in [33], by estimating the stationary conductance given by an average synaptic Poissonian input rate ($\nu_e, \nu_i$):

$$\mu_{ampa}(\nu_e) = \nu_e K_e \tau_{ampa} g_{ampa} \\ \mu_{gabaa}(\nu_i) = \nu_i K_i \tau_{gabaa} g_{gabaa}, \quad (9)$$

where $K_\mu = pN_\mu$ represents the number of incoming synapses from excitatory ($\mu = e$) or inhibitory ($\mu = i$) neurons. These stationary conductances, together with the
instantaneous and time-varying GABA<sub>B</sub> conductance G<sub>b</sub>, will define the effective time constant τ<sub>m</sub><sup>eff</sup> of the membrane:

$$ \mu_G(\nu_e, \nu_i) = \mu_{ampa} + \mu_{gabaa} + G_b + g_L \quad (10) $$

$$ \tau_m^{\text{eff}}(\nu_e, \nu_i) = \frac{C_m}{\mu_G}. \quad (11) $$

Given a specific value of W and G<sub>b</sub>, we can estimate the mean membrane potential as:

$$ \mu_V(\nu_e, \nu_i, W, G_b) = \frac{\mu_{ampa} E_{ampa} + \mu_{gabaa} E_{gabab} + G_b E_{gabab} + g_L E_L - W}{\mu_G}. \quad (12) $$

The calculation of σ<sub>v</sub> and of τ<sub>v</sub> is based on the assumption of Poissonian spiking statistics [35]. Each conductance gives an additional term to σ<sub>v</sub> and τ<sub>v</sub>, so the effect of GABA<sub>B</sub> synapses can be taken into account:

$$ \sigma_V(\nu_e, \nu_i) = \sqrt{\sum_s K_s \nu_s \frac{(U_s \cdot \tau_s)^2}{2(\tau_m^{\text{eff}} + \tau_s)}}, \quad (13) $$

$$ \tau_V(\nu_e, \nu_i) = \left( \frac{\sum_s (K_s \nu_s (U_s \cdot \tau_s)^2)}{\sum_s (K_s \nu_s (U_s \cdot \tau_s)^2 / (\tau_m^{\text{eff}} + \tau_s))} \right), \quad (14) $$

where s = ampa, gabaa, gabab and U<sub>s</sub> = g<sub>s</sub>/μ<sub>G</sub>(E<sub>s</sub> - μ<sub>V</sub>).

### 2.3.2 From Subthreshold Voltage Moments to the Output Firing Rate

We will now estimate the output firing rate of a neuron as a function of (μ<sub>V</sub>, σ<sub>V</sub>, τ<sub>V</sub>), using the following formula:

$$ \nu_{out} = \frac{1}{2\tau_V} \cdot \text{Erfc} \left( \frac{V_{thre}^{eff} - \mu_V}{\sqrt{2}\sigma_V} \right). \quad (15) $$

This formula is based on the hypothesis that the membrane potential follows a Gaussian distribution evolving over a timescale of τ<sub>V</sub>. We thus define the firing rate as the area of such distribution beyond the threshold V<sub>thre</sub><sup>eff</sup>. In this way we obtain an analytical framework to predict the firing rates as a function of (μ<sub>V</sub>, σ<sub>V</sub>, τ<sub>V</sub>). We then fit the function V<sub>thre</sub><sup>eff</sup>(μ<sub>V</sub>, σ<sub>V</sub>, τ<sub>V</sub>) to single neurons' simulations. This procedure overcomes the problem of finding a theoretical value for the threshold, which is not possible due to the nonlinear terms in the membrane potential equation, such as the exponential term in the membrane evolution and the conductance based input. The phenomenological threshold was taken as a second-order polynomial in the following way:

$$ V_{thre}^{eff}(\mu_V, \sigma_V, \tau_V^N) = P_0 + \sum_{x \in \{\mu_V, \sigma_V, \tau_V^N\}} P_x \cdot \left( \frac{x - x^0}{\delta x^0} \right) $$
<table><caption>**Table 1** Fit Parameters (Expressed in mV)</caption><thead><tr><th>Cell Type</th><th>P<sub>0</sub></th><th>P<sub>μV</sub></th><th>P<sub>σV</sub></th><th>P<sub>τ<sub>V</sub><sup>N</sup></sub></th><th>P<sub>μ<sub>V</sub><sup>2</sup></sub></th><th>P<sub>σ<sub>V</sub><sup>2</sup></sub></th><th>P<sub>(τ<sub>V</sub><sup>N</sup>)<sup>2</sup></sub></th><th>P<sub>μVσV</sub></th><th>P<sub>μVτ<sub>V</sub><sup>N</sup></sub></th><th>P<sub>σVτ<sub>V</sub><sup>N</sup></sub></th></tr></thead><tbody><tr><td>FS</td><td>-49.47</td><td>0.01</td><td>7.81</td><td>0.66</td><td>8.90</td><td>-0.53</td><td>1.43</td><td>0.85</td><td>13.68</td><td>0.25</td></tr><tr><td>RS</td><td>-47.04</td><td>1.66</td><td>-4.81</td><td>1.08</td><td>7.59</td><td>0.10</td><td>-0.63</td><td>1.88</td><td>24.66</td><td>0.40</td></tr></tbody></table>

$$ + \sum_{x,y \in \{\mu_V, \sigma_V, \tau_V^N\}^2} P_{xy} \cdot \left( \frac{x - x^0}{\delta x^0} \right) \left( \frac{y - y^0}{\delta y^0} \right), \quad (16) $$

where we used the dimensionless timescale $\tau_V^N = \tau_V G_l / C_m$. The fit of $V_{thre}^{eff}$ is evaluated using a simulation of single neuron with a Poissonian external input fixed and correcting the predicted output until it match the single neuron simulation. We used the following values to expand $V_{thre}^{eff}$: $\mu_V^0 = -50\text{mV}$, $\sigma_V^0 = 0.003\text{mV}$, $(\tau_V^N)^0 = 1$, $\delta\mu_V^0 = 0.01\text{mV}$, $\delta\sigma^0 = 0.006\text{mV}$ and $\delta(\tau^N)^0 = 0.2$. In order to evaluate the values of P, we performed fitting using single neuron simulations with different values of $g_{gabab}$. Results of the fit are reported in Table 1.

# 3 Results

## 3.1 Transfer function estimation

The building block of the mean field model in Eqs. (8) is the transfer function F. The transfer function represents the output firing rate of one RS or FS neuron when receiving excitatory and inhibitory Poissonian spikes at a rate $\nu_e$ and $\nu_i$. The transfer function is thus a propriety of one neuron type, not of a recurrent neural network (that we will treat in next section). In the formalism developed here the transfer function F depends on the rate of incoming excitatory and inhibitory Poissonian spikes $\nu_e$ and $\nu_i$, on the level of adaptation W and on the level of $GABA_B$ conductance $G_b$. In stationary conditions, those where the transfer function is estimated, adaptation and $GABA_B$ conductance reach the value $W^* = b\tau_w\nu_e$ and $G_b^* = \tau_{gabab}\nu_i K_i g_{gabab}$, which depends only on $\nu_e$ and $\nu_i$. In Fig. 1a we show the transfer function for RS and FS neurons in absence of $GABA_B$ synapses ($g_{gabab} = 0\text{nS}$). Dots stand for results of simulation of single neurons and the continuous line is the prediction based on the semi-analytical approach we employed to estimate F. We observe a very good agreement, as it was already proven in previous works [25, 34]. Notice that this semi-analytical approach allows us to predict the output firing rate of a neuron for different parameters, without the need to perform again the fitting. The fitting is performed only once for one parameter setup, and thanks to the semi-analytic approach, we can theoretically predict isolated neurons firing rate for different parameters’ values. Indeed, even if the transfer function has been fitted for a specific parameter setup, we can test its capacity to predict neurons’ output activity in other parameters setup (i.e. for different values of $GABA_B$ quantal conductance) synapses. This basically corresponds to testing the efficacy of our semi-analytic approach (see Section 2.3.1 of the Methods) in presence of time varying conductance based synapses. In Fig. 1b we
report the transfer function for an increased level of $GABA_B$ conductance, showing that the same analytical transfer function works also in presence of $GABA_B$ synapses. We then extended this analysis by increasing systematically the amount of $GABA_B$ conductance $g_{gabab}$. In Fig. 1c we report the firing rate of RS and FS neurons for fixed incoming spike rates $v_e$ and $v_i$ but for different values of $GABA_B$ quantal conductance $g_{gabab}$.

<figure><img src="image_3.png" /><figcaption>**Fig. 1** (a) Theoretical predictions (continuous lines) based on Eq. (15) compared with numerical simulations of the firing rates of RS (FS) neurons, represented as blue (red) dots, for different excitatory input rates $v_e$ and a fixed inhibitory input of $v_i = 8$Hz. Quantal conductance of $GABA_B$ synapses is $g_{gabab} = 0$pS. (b) Same as (a), but for $g_{gabab} = 25$pS. (c) RS (FS) neurons firing rates for fixed excitatory inputs $v_e = 4$Hz and inhibitory input $v_i = 8$Hz, shown as blue (red) dots for simulations and continuous line for theoretical predictions, as a function of $g_{gabab}$. (d) Average membrane potential of RS neurons as function of $g_{gabab}$ (same simulations as panel c).</figcaption></figure>

As expected, the firing rate decreases increasing $g_{gabab}$. Importantly, we observe a very good agreement between numerical simulations and theoretical predictions. Moreover, in Fig. 1d we observe that our theory (in particular Eq. 12, see Methods) allows us to predict the mean membrane potential of neurons for different level of $GABA_B$ quantal conductance $g_{gabab}$. We can thus now employ the transfer function to predict the emergent dynamics of recurrent neural networks of RS and FS neurons through the mean field Eq.s (8). Notice that even if we obtain a good theoretical transfer function for isolated neurons, this does not necessarily imply that we can predict the population dynamics of a sparsely connected network of interacting neurons. We show in the next chapters that the complete mean field Eq.s (8) here presented, are able to correctly predict sparse network dynamics.

## 3.2 Spontaneous activity in Asynchronous irregular regimes

Here, we study the emergent dynamics of a recurrently connected network composed of RS and FS neurons (see Methods section for details) in the absence of synaptic delay
($\tau_D$ = 0ms). In Fig. 2a we show the raster plot for RS neurons. Notice that we considered the presence of $GABA_B$ synapses, with a quantal conductance $g_{gabab} = 25pS$. We observe an asynchronous irregular regime characterized by neurons with Poissonian firing statistics (we have verified that the coefficient of variation $CV$ of neurons firing statistics is around $CV = 1$ for both RS and FS neurons, as expected for Poisson processes). Since the network dynamics is asynchronous, it should be represented by a fixed point in the mean-field model of Eq.s (8).

<figure><img src="image_4.png" /><figcaption>**Fig. 2** (a) Raster plot of excitatory neurons spiking activity for $g_{gabab} = 25pS$. The network simulation consists of 10,000 neurons, but we report only a subsample of 1,000 excitatory neurons. (b) Function $\mathcal{F}(\nu_e)$ derived from the mean field equations (8) (see main text) is reported as a continuous blue line. The bisector is the orange line, while the dots represent the corresponding excitatory neurons firing rate in numerical simulations of the spiking neural network. (c) Comparison between the fixed point predicted by the mean field equations (continuous lines) and the values obtained from numerical simulations (dots), for different values of $g_{gabab}$.</figcaption></figure>

We thus report in Fig. 2b the consistency map for the theoretical excitatory neurons firing rate $\nu_e$ (see the continuous blue line). This map is obtained by estimating, for each value of excitatory input $\nu_e$, the corresponding stationary value of inhibitory neurons activity $\nu_i^*$, the adaptation variable $W^*$, the $GABA_B$ conductance $G_b^*$ and corresponding excitatory firing rate $\mathcal{F}(\nu_e) = F_e(\nu_e, \nu_i^*, W^*, G_b^*)$. A consistent solution (fixed point) is obtained when $\mathcal{F}(\nu_e) = \nu_e$, i.e. when the blue line crosses the bisector (see the orange dashed line in Fig. 2a). As shown in Fig. 2b the mean field model predicts a stable fixed point, with a value of $\sim 2$ Hz. When reporting the values of $\nu_e$ at different time steps from numerical simulation of the neural network (see blue dots in Fig. 2b) we observe values with a mean of $\sim 2.3$ Hz characterized by natural fluctuations over time. These fluctuations are normally distributed [25] and we have verified that they vanish in the limit of infinite number of neurons. Moreover, the value of the mean rate in network simulations does not change for sufficiently large networks (i.e. more than 1,000 neurons). From these results we conclude that the mean field predictions are consistent with numerical simulations in terms of the mean activity of population rates. Notice that the mean formalism employed here can be extended to include equations predicting the dynamics of second order quantities, e.g. the amount
of fluctuations in finite size networks [14, 25]. We limit here to discuss the prediction
of the mean field in terms of the mean activity of population rates. We then asked if
the prediction of the mean field model for stationary firing activity of excitatory and
inhibitory neurons are consistent with network simulations increasing the strength of
$GABA_B$ conductances $g_{gabab}$. In Fig. 2c we observe a very good agreement for both
excitatory and inhibitory activity over a wide range of values of $g_{gabab}$. In summary
the mean field model correctly predicts stationary population activity for different
values of $GABA_B$ synaptic conductances.

## 3.3 Evoked activity in Asynchronous irregular regimes

We then investigated whether the mean field model correctly predicts network response
to transient external stimulation. We thus introduced a time varying excitatory input
frequency $v_{ext}$ and observed the response of the network.

In Fig. 3a and b we report the raster plot of excitatory neurons and the population
firing rate $v_e(t)$ in response to a rising stimuli at time $t_{rise} = 3$s, that goes back to
the baseline value at time $t_{end} = 5$s, as described by Eq. (7) (see Methods section
and the dashed line in the inset of Fig. 3b). Notice that we employ a network with
excitatory and inhibitory neurons as described in the method section. We observe a
very good agreement between the complete mean field model of Eq.s (8) (that includes
the dynamics of $GABA_B$ synapses) and the network simulation. Notice that before
stimulation at $t_{rise} = 3$s the mean field predicts a rate ~2Hz and network simulation
show a rate of ~2.3Hz with intrinsic fluctuations of ±0.5Hz. The fluctuations observed
in network simulations are due to the finite effects in our sparse network. We also
report the results from a stationary version of this mean field model (dashed line).
This stationary version is obtained by considering instantaneous $GABA_B$ synapses,
i.e. without considering their slow decay in time but only their tonic stationary value,
$G_b^* = \tau_{gabab} K_i v_i g_{gabab}$. This is the same approach we employ for fast $GABA_A$ and
$AMPA$ receptors. In the complete mean field instead, $G_b$ follows its own differential
equation, see Eq. (8d). Notice that when the population dynamics is stationary (i.e.
non-evolving in time), the stationary mean field predicts the same fixed values for
population rates to the complete mean field. In fact, before the stimulus onset at $t=3$s,
the stationary mean field (dashed line) and the complete mean field (continuous line)
are superimposed, and they predict very well the network activity. We observe instead
that the complete mean field captures much better the transient response at stimulus
onset, while the stationary mean field underestimates the response to the stimulus
arrival at $t = 3$s and overestimates network activity at stimulus offset at time $t = 5$s.
This is a direct consequence of the time scale of $GABA_B$ synapses, that is taken into
account in our mean field model. We have also tested a sinusoidal input at a frequency
$f = 2$Hz. The results are shown in Fig.3c, where we observe that the complete mean
field is capable to correctly reproduce the amplitude of the response entrained by the
oscillatory input, at variance with the previous stationary mean field model. As we
have shown in a previous work, the adiabatic approximation here proposed induced
discrepancies between the mean field model predictions and the network simulation
for too fast time scales [25]. This implies that whenever $\tau_{gabab}$ is too fast (typically
faster that 100ms) and when the input frequency is too high (typically faster that 10
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Comparison between mean field predictions, network simulations, and a stationary mean field where $G_b$ is considered always equal to its stationary value. (a) Raster plot of excitatory neurons spiking activity for $g_{gabab}$ = 25pS in response to an external input. The input is reported as a dashed line in panel b, and follows the Eq. (7), see Methods. The network simulation consists of 8,000 excitatory and 2,000 inhibitory neurons, we show here only a subsample of 1,000 excitatory neurons. (b) Corresponding RS neurons population rate $v_e(t)$ (blue noisy line). The shape of the stimuli is represented by the black dotted line below the graph. The thick blue line is the mean field model prediction from Eq.s (8) while the dashed line is the stationary mean field model prediction (i.e. without $GABA_B$ dynamics). (c) As (b) but for a sinusoidal input of frequency $f = 2$Hz and amplitude $A = 1$Hz (see Methods section).</figcaption></figure>

Hz), the mean field and the network simulation show a quantitative discrepancy due to the adiabatic approximation.

## 3.4 Theta oscillations with $GABA_B$ receptors

We have then employed the mean field model to predict if new dynamical regimes can be induced by the presence of $GABA_B$ synapses. In Fig. 4a we report the standard deviation (over time) of the excitatory population firing rate, namely $\Sigma(v_e)$. Red dots are obtained for the value of the external drive $v_{ext} = 2.5$Hz used in Fig. 2 and
a

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

b

<figure><img src="image_7.png" /><figcaption>**Fig. 4** (a) The standard deviation (over time) $\Sigma(\nu_e)$ of the excitatory neurons population rate $\nu_e(t)$ estimated from the mean field model in Eq.s (8) as a function of $g_{gabab}$ for $\nu_{ext} = 2.5$ Hz (red dots) and $\nu_{ext} = 1.5$ Hz (grey dots). (b) Excitatory $\nu_e(t)$ (top panel) and inhibitory $\nu_i(t)$ (middle panel) population rate predicted by the mean field model when $g_{gabab}$ increases (see bottom panel, from zero to $g_{gabab} = 0.08$ nS ) for $\nu_{ext} = 1.5$ Hz. The power spectrum of $\nu_e(t)$ is reported in the inset of the top panel. Dashed line correspond to the stationary mean field model (i.e. with stationary $GABA_B$). The light solid and noisy lines correspond to network simulations with 10,000 neurons.</figcaption></figure>

Fig. 3. We observe that $\Sigma(\nu_e)$ is always zero, meaning that the excitatory firing rate $\nu_e(t)$ does not change over time, i.e. the mean field model displays a stable fixed point corresponding to asynchronous network dynamics. On the contrary, for lower external drive $\nu_{ext} = 1.5$Hz we observe a transition to a non-zero $\Sigma(\nu_e)$ for sufficiently large values of $GABA_B$ conductances, indicating that the excitatory firing rate $\nu_e(t)$ changes over time. In Fig. 4b we report the time trace of neural population rates for values of $GABA_B$ conductances in the regime of large $\Sigma(\nu_e)$ (i.e. $\nu_{ext} = 1.5$Hz and $g_{gabab} = 0.08$ nS). We observe periodic oscillations, at a frequency of $f_{osc} = 4$ Hz in the theta range. In the same figure we report the same simulation for the stationary mean field model, that does not account for the dynamics of $GABA_B$ synapses (dashed lines). Importantly, we observe that the stationary mean field (with stationary $GABA_B$ synapses) does not reproduces the emergence of theta oscillations, at variance with the new mean field model. We have verified that the frequency of oscillations depends on the value of the synaptic time scale $\tau_{gabab}$: the longer $\tau_{gabab}$
the slower network oscillations. Notice that the mean excitatory neurons firing rate (over time) is around 1.3 Hz, which is lower than the frequency of oscillations (around 5Hz). This means that neurons do not spike at every oscillation cycle. We have also estimated the number of spikes per neuron during one oscillation by integrating the excitatory neurons firing rate, and we have found that one neuron fires around 0.27 spikes during one cycle of oscillations (i.e. one spike every 4 cycles on average). These results confirm that these oscillations appear in a sparsely firing regime. We then tested the predictions of the new mean field models against direct numerical simulations. We observe that theta oscillations emerge for sufficiently large values of $GABA_B$ also in direct neural network simulations. The amplitude of these oscillations is lower in network simulations with respect to mean field predictions (see light noisy line in Fig. 4b) but the frequency of oscillations is almost the same as the one predicted by the mean field model (~ 5 Hz). In Fig. 5 we report the raster plot (Fig. 5a-b) and the power spectrum of population activity $v_e$ for network simulation with and without $GABA_B$ (Fig. 5c), confirming mean field predictions.

<figure><img src="image_8.png" /><figcaption>**Fig. 5** Numerical simulations to confirm mean field predictions of Fig. 4. Panel (a) shows the raster plot obtained from numerical simulation of the neural network with $g_{gabab} = 0$ pS and panel (b) for $g_{gabab} = 0.08$ nS and $v_{ext} = 1.5$ Hz. Panel (c) shows the power spectrum of the excitatory population firing rate $v_e(t)$ for the two conditions.</figcaption></figure>

## 3.5 Distinct roles of increasing inhibition for gamma oscillations

Experimental data suggest that increasing inhibition in cortical circuits can have different impact on population dynamics, and in particular on gamma oscillations. Optogenetically activating interneurons has been found to increase the power of gamma oscillations [28, 29]. On the other side, increasing inhibition by increasing the conductance of $GABA_B$ synapses has seems to have opposite effect, i.e. $GABA_B$ conductance was found to eliminate gamma oscillations in rat hippocampal slices [27]. In this section we show that our mean field model reproduces and predicts this differential role of gamma oscillations. We have thus equipped our mean field model with synaptic delay $\tau_D = 3$ ms, while all previous simulations were obtained for instantaneous synapses ($\tau_D = 0$ ms). In this version of the model, following the recent work published in [26], we consider an additive noise modeled as an Ornstein-Ulhenbeck (OU)
process to Eq.s (8). Accordingly, the external drive $v_{ext}$ becomes a time-dependent variable $v_{ext}(t) = v_{ext} + \xi(t)$, where $\xi(t)$ is an (OU) process evolving according to $d\xi = -\xi\theta dt + \sigma dW_t$ with $\theta = 200$Hz, $\sigma = 5$ and where $W_t$ is the standard Brownian motion. In Fig. 6 we simulate such mean field model with $\tau_D = 3$ms (see the circle in panel a, and the corresponding time trace in the top of panel b). We observe the emergence of gamma oscillations at around $f_{osc} = 80$Hz, see top panel in Fig. 6b. Gamma oscillations are not observed in the mean field model for $\tau_D = 0$ms (dashed line, no delay). From this control parameter setup, we have then independently increased the external drive to inhibitory neurons $I_{ext}$ and the conductance of slow $GABA_B$ synapses in the mean field model. In order to evaluate the amplitude of gamma oscillations we report in Fig. 6a the standard deviation over time of inhibitory firing rate $\Sigma(\nu_i)$. Note that this quantity measures the degree of synchrony among interneurons, which contributes to LFP power. The LFP power can indeed be estimated as the sum of net inhibitory and excitatory inputs to pyramidal cells. [36]. In the control condition (circle) we observe a moderate amplitude of gamma oscillations, see the time trace in the top panel of Fig. 6b. By manually activating inhibitory neurons through increasing of $I_{ext}$ the power of oscillations increases (see the star and the middle panel of Fig. 6b). On the contrary, increasing $GABA_B$ conductance kills gamma oscillations (see the square and the bottom panel of Fig. 6b.

We have then verified these predictions in numerical simulation of the neural network. In Fig. 6c we report the raster plot for the parameters in panel b (circle, star and square) of the same figure. We observe that indeed $I_{ext}$ and $GABA_B$ have a differential role in the generation of gamma oscillations: increase in inhibitory network synchronization for $I_{ext}$ and asynchronous dynamics by increasing $GABA_B$. These results are in agreement with our mean field model predictions and with experimental data.

## 3.6 Gamma bursts are regulated by GABAB synapses

Oscillations in brain recordings such as LFP are not stationary limit cycles but typically appear as bursts, i.e. short interval of time characterized by few cycles of oscillatory behavior, separated by asynchronous activity. This has been observed for beta [37] but also for gamma oscillations [30, 38]. Oscillations have thus a non-trivial dynamics, characterized by the appearance and the duration of gamma bursts. Our model has shown that $GABA_B$ synapses affect the amplitude of gamma oscillations. In this section we investigate the impact of $GABA_B$ synapses on the appearance of gamma bursts. Also here we consider an additive noise modeled as an Ornstein–Uhlenbeck (OU) process to Eq.s (8), see previous section for details and parameters. In Fig. 7 we report the time-frequency analysis of excitatory neurons’ firing rate predicted by the mean field model. We then superimpose the mean field $GABA_B$ conductance $G_b$. Different panels are obtained for different values of $GABA_B$ synaptic time scale $\tau_{gabab}$. We observe that for sufficiently fast inhibition, gamma oscillations are very powerful and sustained over time, resembling the limit cycles observed in Fig. 6b for low $GABA_B$ conductance, i.e. a situation where inhibition is dominated by fast $GABA_A$ synapses. Increasing the time scale of $\tau_{gabab}$ we observe that gamma bursts become more and more rare and less powerful (lower amplitude). By looking at the trace of
<figure><img src="image_9.png" /><figcaption>**Fig. 6** Panel (a) shows the standard deviation over time of the inhibitory firing rate $\Sigma(\nu_i)$ for different values of $g_{gabab}$ and $I_{drive}$ estimated from simulations of the mean field model with $\tau_D = 3$ms. We used $T = 2$ms for the mean field simulations. Circle, star and square indicate a set of parameter values employed in panel b (from top to bottom), where we report $\nu_i$ over time for the mean field with $\tau_D = 3$ms (continuous line) and $\tau_D = 0$ms (dashed line). In panel c we report results from numerical simulations of the neural network (we report the raster plots) for the parameter values of circle, star and square, from left to right.</figcaption></figure>

$G_b$ we observe that gamma bursts are appearing by the disinhibition from $GABA_B$ release and they terminate after a sufficient level of $GABA_B$ is accumulated by the activation of inhibitory neurons in a gamma burst. We thus conclude that the dynamics of $GABA_B$ receptors and its time scales regulates the appearance and the duration of gamma bursts.

## 3.7 Nonlinear $GABA_B$ synapses block epileptic seizures

We now employ our excitatory and inhibitory network of neurons and its mean field model to investigate the stability of network dynamics. In particular, we focus on the nonlinearity of $GABA_B$ synapses, i.e. the fact that synapses are activated only if there is a sufficient number of synchronous incoming spikes. The function $G(\nu_i)$ follows a sigmoid that increases when $K_i\nu_i > 1$, that corresponds to more than 10 presynaptic spikes in 100ms. In our previous simulations we considered $K_i = 100$ and we worked in setups such that the inhibitory neurons firing rate $\nu_i$ was usually much higher than 0.01 Hz, thus $G$ was always equal to 1. We have indeed verified
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Power spectrograms of the excitatory population rate $v_e(t)$ from the mean field model 8 for different values of $GABA_B$ conductance time scale $\tau_{gabab}$ (panels a-b-c). We use the short-time Fourier transform (STFT) subroutine from the signal package of the SciPy library [39] to obtain the Fourier transform of $v_e(t)$ within a running time window of length $\Delta_t$ at time t. We perform STFT using a 90% overlap, with a window length of $\Delta_t = 0.05$ s. The spectrogram was visualized using a colormap, where the color code represents the normalized power spectral density. The white trace represents the $GABA_B$ conductance $G_b$. For this simulation $v_{ext} = 7$ Hz and $g_{gabab} = 0.01$ nS.</figcaption></figure>

that the results of previous simulations are unchanged if we consider linear synapses with $G = 1$. Nevertheless, the situation changes if we consider very sparse networks with low activity. In the simulations of Fig. 8 we considered a sparser network with $p = 0.025\%$, i.e. mean inhibitory (excitatory) inputs per postsynaptic neuron $K_i = 50$ ($K_e = 200$) and a lower external drive $v_{baseline} = 1.5$ Hz. This allows to model a sparse network with low activity, where the amount of open $GABA_B$ receptors is low during spontaneous activity. In Fig. 8a-d we have performed a simulation on this network, by including a perturbation at time $t_{ep}=3$ s, where excitatory AMPA conductance are suddenly increased (see panel d). Before $t_{ep} = 3$ s we observe a sparse asynchronous firing of neurons (see panel a and b) and a very low fraction $G \sim 1\%$ of activated $GABA_B$ synapses (see panel c). In panel a we considered linear synapses, i.e. $G$ does not depend on $v_i$ and is fixed to $G \sim 1\%$. We observe that when $g_{ampa}$ increases at time $t_{ep}$ there is a transition to high frequency activity (around 150 Hz) characterized by hypersynchronous bursts (see the inset of panel a). This corresponds to a pathological hyperactivity typical of epileptic seizures. In panel b we considered instead the nonlinear $GABA_B$ synapses. We observe that the network does not transit to hyperactivity but to relatively slow oscillations driven by the dynamics of the fraction of activated $GABA_B$ receptors $G$ (see blue continuous line in panel c). We have then performed the same simulations in the mean field model observing very similar results. In this case, we considered a lower value of $g_{ampa}$ for the transition to epilepsy. Data are reported in panels e-g. We observe that linear $GABA_B$ yields
<figure><img src="image_11.png" /><figcaption>**Fig. 8** Panels a-d refer to neural network simulations. Panel a shows the raster plot of excitatory neurons for linear $GABA_B$ synapses, i.e. $G \sim 1\%$, while panel b shows the raster plot for nonlinear $GABA_B$ synapses, see Methods. In panel c we report the values of $G$ for the simulation in panel a (green dashed line) and panel b (blue continuous line). In panel d we show that the values of $g_{ampa}$, at $t = t_{ep} = 3s$ increases, leading to an epileptic seizure for the linear model (panel a). In the inset of the panel a we report a zoom for a duration of 20ms where we observe hypersynchronous bursts. Panels e-g refer to mean field model simulations. Top panel shows the firing rate of excitatory neurons $v_e$ for the linear (green dashed line) and nonlinear (blue continuous line) $GABA_B$ conductance model. Panel f shows the values of $G$ and panel g the values of $g_{ampa}$. In these simulations $K_i = 50$, $v_{baseline} = 1.5Hz$ and $b = 2pA$.</figcaption></figure>

hyperactivity at 150Hz (green dashed line in panel e), while nonlinear $GABA_B$ drive the system to healthy slow oscillations driven by the dynamics of $G$. Notice that our
mean field approach cannot reproduce oscillations as fast as (or even faster) than the time scale of the Markovian process $T \sim 10$ms. As a result, it cannot reproduce the hyper-synchronous oscillations that we observe in the network during epileptic state (see the inset of Fig.8a). In the mean field model, this state is represented by an abnormally high population firing rate $\nu_e \sim 150$Hz. These results show that our mean field model can capture important features of the effect of nonlinear $GABA_B$ synapses, such as their potential role in stabilizing network activity to avoid epileptic seizures.

# 4 Discussion

In this paper we have developed a new mean field model of sparse conductance-based spiking neural networks including the slow dynamics of $GABA_B$ receptors. The predictions of the mean field models have been directly compared to numerical simulations of a spiking neural network of excitatory and inhibitory neurons. We have found a very good agreement between the mean field model and the network's simulations, both in spontaneous activity and in the response to external stimulation.

A first important outcome of this model is the emergence of theta oscillations at 4-5 Hz. The mean field model predicts that these oscillations can be locally produced by the slow time scale of $GABA_B$ synapses. We have verified that a decrease in $GABA_B$ time scale increased the frequency of these oscillations, covering the theta range observed in rats (4-12)Hz and in Humans (4-8)Hz [40]. The important point of this analysis is that the stationary mean field model (that does not account for the dynamics of slow $GABA_B$ synapses) cannot predict the emergence of these oscillations.

A second important outcome from this mean field model is its capacity to predict different role of increasing inhibition in neural networks. Inhibition can be increased either by directly activating inhibitory neurons (such as in experiments of optogenetics) or by increasing $GABA_B$ conductances (such as by pharmacological agonists). These two manipulations yield different results. Experimental data suggests that $GABA_B$ reduces gamma oscillations, at variance with optogenetic stimulations of interneurons. Our mean field model reproduces gamma oscillations and the different role of these manipulations, thus suggesting that these are nontrivial emergent properties of neuronal networks captured by our "simple" and low dimensional mean field model.

In our approach we have considered $GABA_B$ receptors having a slow decaying activity, at variance with $GABA_A$ synapses. This is a main difference between these two receptor types, but not the only one. In fact, as mentioned in the method section, $GABA_B$ receptors have a nonlinear activation at variance with $GABA_A$ receptors. They get activated only when there are 3-10 repetitive presynaptic inputs in a short time interval of around 100ms [31]. This property could play a crucial role to stabilize the network and avoid epileptic seizures, as our model predicts. Nevertheless, this happens only if we consider very sparse networks with very low firing rates. In our mean field model (and in direct numerical simulations of neural networks) nonlinear $GABA_B$ synapses seem to prevent epileptic seizures by strongly moving to a fully open $GABA_B$ receptor state, characterized by slow oscillations of few Hz driven by the nonlinear opening and closing of $GABA_B$ receptors. As these nonlinear receptors
are strongly activated by transient events of strongly active interneurons, it will be
interesting in future studies to investigate the impact of bursting interneurons on the
emergent dynamics of neural networks through the nonlinearity of GABA<sub>B</sub> synapses.
In this view, bursting interneurons could serve as openers of GABA<sub>B</sub> synapses, putting
the network in a safe mode to avoid epilepsy. This could generalize to any type of high
firing activity, which would likely activate GABA<sub>B</sub> receptors [31].

It has been largely reported that dysfunction of the GABAergic system has a
crucial impact on the genesis and the propagation of epileptic seizures [41]. Many
studies have shown an impact of GABA<sub>A</sub> receptors for the genesis of epileptic seizure
in mice models [42]. Our model proposes instead a major role of GABA<sub>B</sub> receptors,
due to their nonlinearity and slower dynamics. Interestingly, some experimental work
points to GABA<sub>B</sub> receptors as principal actors for the absence of seizures, presumably
in the thalamic nuclei [43, 44]. It is important to notice that the data in this regard
are quite complex as they show that GABA<sub>B</sub> receptor involvment can be pro or anti
epilepsy depending on the nature of the pathological neuronal networks involved [45].
Interestingly, our mechnistic model suggests a major role of GABA<sub>B</sub> receptors, but
only for sufficiently sparsely connected networks with low firing activity in healthy
conditions. Finally, let us point out that on top of GABAergic involvment in epilepsy,
many experimental works suggest a major role of glial cells and plasticity [46]. Future
computational works should study the distinct or synergetic role of these mechanisms
for the stability of network dynamics.

In our model we considered only one type of interneurons, fast spiking inhibitory
neurons resembling parvalbumin (PV) cells. Nevertheless, it is known that different
interneuons types are present in the cortex [47] and in the hippocampus [48]. Also,
recent works seem to suggest that different interneurons types play synergetic role in
network dynamics and stability [49]. Future works should aim to develop mean field
models able to include the diversity of these interneurons. Our approach, capable to
model different neurons proprieties (such as duration of action potential and spike
frequency adaptation) seem suitable for this goal. Another important aspect is the
ratio between excitatory and inhibitory neurons. We have here considered 80-20% for
the sake of generality, but this ratio is variable across brain areas. Future works could
employ our model to study the effect of GABA<sub>B</sub>, for example on seizure susceptibility,
on different regions depending on their specific E/I ratio.

Furthermore, the adiabatic approach employed here to model slow inhibitory
synapses can be employed also for slow excitatory NMDA receptors. A similar
approach has been employed in the context of the exact reduction of QIF neurons
[50], by taking also into account the additional nonlinearity of NMDA receptor. This
approach could be employed, combined with the one presented here, to include NMDA
receptors in the context of sparse conductance based networks such as those employed
in this work. This model could be also employed to investigate the effects of NMDA
receptors to regulate gamma oscillations [51].

Finally, this mean field approach can be employed to study large scale interac-
tion between brain areas through cross frequency coupling, and the role of receptor
specific agonist on interneurons synapses. This approach can show differential role of
pharmacological action (namely $GABA_A$ or $GABA_B$ agonists) in the functional con-
nectivity at large scale by using different species of connectomes, following the large
scale modeling developed in [10] for the mean field model developed in [25].

**Acknowledgements.** We thank B. Delord, E. Procyk and F. Tahvili for useful com-
ments and discussions on this work. This work was supported by the French Ministry
of Higher Education (Ministère de l’Enseignement Supérieur) and the project LABEX
CORTEX (ANR-11-LABX-0042) of Université Claude Bernard Lyon 1 operated by
the ANR. A. Destexhe received support from the European Union (Virtual Brain Twin
project 101137289) and ANR (FLAG-ERA BrainAct project).

References

[1] Dehghani, N., Peyrache, A., Telenczuk, B., Le Van Quyen, M., Halgren, E., Cash, S.S., Hatsopoulos, N.G., Destexhe, A.: Dynamic balance of excitation and inhibition in human and monkey neocortex. Scientific reports **6**, 23176 (2016)

[2] Buzsáki, G.: Rhythms of the brain. Oxford university press (2006)

[3] Buzsáki, G., Wang, X.-J.: Mechanisms of gamma oscillations. Annual review of neuroscience **35**, 203–225 (2012)

[4] Wang, X.-J.: Neurophysiological and computational principles of cortical rhythms in cognition. Physiological reviews **90**(3), 1195–1268 (2010)

[5] Jensen, O., Colgin, L.L.: Cross-frequency coupling between neuronal oscillations. *Trends in cognitive sciences* **11**(7), 267–269 (2007)

[6] Zhang, H., Jacobs, J.: Traveling theta waves in the human hippocampus. Journal of Neuroscience **35**(36), 12477–12487 (2015)

[7] Belluscio, M.A., Mizuseki, K., Schmidt, R., Kempter, R., Buzsáki, G.: Cross-frequency phase–phase coupling between theta and gamma oscillations in the hippocampus. Journal of Neuroscience **32**(2), 423–435 (2012)

[8] Colgin, L.L.: Theta-gamma coupling in the entorhinal-hippocampal system. *Current opinion in neurobiology* **31**, 45–50 (2015)

[9] Fries, P.: Rhythms for cognition: communication through coherence. Neuron **88**(1), 220–235 (2015)

[10] Goldman, J.S., Kusch, L., Aquilue, D., Yalçınkaya, B.H., Depannemaecker, D., Ancourt, K., Nghiem, T.-A.E., Jirsa, V., Destexhe, A.: A comprehensive neural simulation of slow-wave sleep and highly responsive wakefulness dynamics. *Frontiers in Computational Neuroscience* **16**, 1058957 (2023)

[11] Wilson, H.R., Cowan, J.D.: Excitatory and inhibitory interactions in localized populations of model neurons. *Biophysical journal* **12**(1), 1–24 (1972)
[12] Ohira, T., Cowan, J.D.: Master-equation approach to stochastic neurodynamics. Phys. Rev. E **48**, 2259–2266 (1993) https://doi.org/10.1103/PhysRevE.48.2259

[13] Ginzburg, Sompolinsky: Theory of correlations in stochastic neural networks. Physical review. E, Statistical physics, plasmas, fluids, and related interdisciplinary topics **50** 4, 3171–3191 (1994)

[14] El Boustani, S., Destexhe, A.: A master equation formalism for macroscopic modeling of asynchronous irregular activity states. Neural computation **21**(1), 46–100 (2009)

[15] Buice, M.A., Cowan, J.D., Chow, C.C.: Systematic fluctuation expansion for neural network activity equations. Neural Computation **22**(2), 377–426 (2010) https://doi.org/10.1162/neco.2009.02-09-960
https://doi.org/10.1162/neco.2009.02-09-960. PMID: 19852585

[16] Schwalger, T., Deger, M., Gerstner, W.: Towards a theory of cortical columns: From spiking neurons to interacting neural populations of finite size. PLoS computational biology **13**(4), 1005507 (2017)

[17] Montbrió, E., Pazó, D., Roxin, A.: Macroscopic description for networks of spiking neurons. Physical Review X **5**(2), 021028 (2015)

[18] Segneri, M., Bi, H., Olmi, S., Torcini, A.: Theta-nested gamma oscillations in next generation neural mass models. Frontiers in computational neuroscience **14**, 47 (2020)

[19] Nandi, M.K., Valla, M., Di Volo, M.: Bursting gamma oscillations in neural mass models. Frontiers in Computational Neuroscience **18**, 1422159 (2024)

[20] Di Volo, M., Torcini, A.: Transition from asynchronous to oscillatory dynamics in balanced spiking networks with instantaneous synapses. Physical review letters **121**(12), 128301 (2018)

[21] Goldobin, D.S., Di Volo, M., Torcini, A.: Reduction methodology for fluctuation driven population dynamics. Physical Review Letters **127**(3), 038301 (2021)

[22] Pyragas, V., Pyragas, K.: Mean-field equations for neural populations with q-gaussian heterogeneities. Physical Review E **105**(4), 044402 (2022)

[23] Carlu, M., Chehab, O., Dalla Porta, L., Depannemaecker, D., Héricé, C., Jedynak, M., Köksal Ersöz, E., Muratore, P., Souihel, S., Capone, C., et al.: A mean-field approach to the dynamics of networks of complex neurons, from nonlinear integrate-and-fire to hodgkin–huxley models. Journal of neurophysiology **123**(3), 1042–1051 (2020)

[24] Di Volo, M., Destexhe, A.: Optimal responsiveness and information flow in
networks of heterogeneous neurons. Scientific reports 11(1), 17611 (2021)

[25] Di Volo, M., Romagnoni, A., Capone, C., Destexhe, A.: Biologically realistic mean-field models of conductance-based networks of spiking neurons with adaptation. Neural computation 31(4), 653-680 (2019)

[26] Tahvili, F., Destexhe, A.: A mean-field model of gamma-frequency oscillations in networks of excitatory and inhibitory neurons. Journal of Computational Neuroscience 52(2), 165-181 (2024)

[27] Brown, J.T., Davies, C.H., Randall, A.D.: Synaptic activation of gabab receptors regulates neuronal network activity and entrainment. European Journal of Neuroscience 25(10), 2982-2990 (2007)

[28] Cardin, J.A., Carlén, M., Meletis, K., Knoblich, U., Zhang, F., Deisseroth, K., Tsai, L.-H., Moore, C.I.: Driving fast-spiking cells induces gamma rhythm and controls sensory responses. Nature 459(7247), 663-667 (2009)

[29] Sohal, V.S., Zhang, F., Yizhar, O., Deisseroth, K.: Parvalbumin neurons and gamma rhythms enhance cortical circuit performance. Nature 459(7247), 698-702 (2009)

[30] Douchamps, V., Di Volo, M., Torcini, A., Battaglia, D., Goutagny, R.: Gamma oscillatory complexity conveys behavioral information in hippocampal networks. Nature Communications 15(1), 1849 (2024)

[31] Thomson, A.M., Destexhe, A.: Dual intracellular recordings and computational models of slow inhibitory postsynaptic potentials in rat neocortical and hippocampal slices. Neuroscience 92(4), 1193-1215 (1999)

[32] Destexhe, A.: Self-sustained asynchronous irregular states and up-down states in thalamic, cortical and thalamocortical networks of nonlinear integrate-and-fire neurons. Journal of computational neuroscience 27(3), 493 (2009)

[33] Kuhn, A., Aertsen, A., Rotter, S.: Neuronal integration of synaptic input in the fluctuation-driven regime. Journal of Neuroscience 24(10), 2345-2356 (2004)

[34] Zerlaut, Y., Chemla, S., Chavane, F., Destexhe, A.: Modeling mesoscopic cortical dynamics using a mean-field model of conductance-based networks of adaptive exponential integrate-and-fire neurons. Journal of computational neuroscience 44(1), 45-61 (2018)

[35] Zerlaut, Y., Teleńczuk, B., Deleuze, C., Bal, T., Ouanounou, G., Destexhe, A.: Heterogeneous firing rate response of mouse layer v pyramidal neurons in the fluctuation-driven regime. The Journal of physiology 594(13), 3791-3808 (2016)

[36] Mazzoni, A., Lindén, H., Cuntz, H., Lansner, A., Panzeri, S., Einevoll, G.T.:
Computing the local field potential (lfp) from integrate-and-fire network models. PLoS computational biology **11**(12), 1004584 (2015)

[37] Lundqvist, M., Miller, E.K., Nordmark, J., Liljefors, J., Herman, P.: Beta: bursts of cognition. Trends in Cognitive Sciences (2024)

[38] Lundqvist, M., Rose, J., Herman, P., Brincat, S.L., Buschman, T.J., Miller, E.K.: Gamma and beta bursts underlie working memory. Neuron **90**(1), 152–164 (2016)

[39] Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cour- napeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., et al.: Scipy 1.0: fundamental algorithms for scientific computing in python. Nature methods **17**(3), 261–272 (2020)

[40] Jacobs, J.: Hippocampal theta oscillations are slower in humans than in rodents: implications for models of spatial navigation and memory. Philosophical Trans- actions of the Royal Society B: Biological Sciences **369**(1635), 20130304 (2014)

[41] Sperk, G., Furtinger, S., Schwarzer, C., Pirker, S.: Gaba and its receptors in epilepsy. Recent advances in epilepsy research, 92–103 (2004)

[42] DeLorey, T., Handforth, A., Anagnostaras, S., Homanics, G., Minassian, B., Asa- tourian, A., Fanselow, M., Delgado-Escueta, A., Ellison, G., Olsen, R.: Mice lacking the β3 subunit of the gabaa receptor have the epilepsy phenotype and many of the behavioral characteristics of angelman syndrome. Journal of Neuroscience **18**(20), 8505–8514 (1998)

[43] Hosford, D.A., Clark, S., Cao, Z., Wilson Jr, W.A., Lin, F.-h., Morrisett, R.A., Huin, A.: The role of gabab receptor activation in absence seizures of lethargic (lh/lh) mice. Science **257**(5068), 398–401 (1992)

[44] Liu, Z., Vergnes, M., Depaulis, A., Marescaux, C.: Involvement of intrathalamic gabab neurotransmission in the control of absence seizures in the rat. Neuroscience **48**(1), 87–93 (1992)

[45] Joshi, K., Cortez, M.A., Snead, O.C.: Targeting the gaba b receptor for the treatment of epilepsy. GABAB receptor, 175–195 (2016)

[46] Heinemann, U., Gabriel, S., Jauch, R., Schulze, K., Kivi, A., Eilers, A., Kovacs, R., Lehmann, T.-N.: Alterations of glial cell function in temporal lobe epilepsy. Epilepsia **41**, 185–189 (2000)

[47] Tremblay, R., Lee, S., Rudy, B.: Gabaergic interneurons in the neocortex: from cellular properties to circuits. Neuron **91**(2), 260–292 (2016)

[48] Klausberger, T., Somogyi, P.: Neuronal diversity and temporal dynamics: the unity of hippocampal circuit operations. Science **321**(5885), 53–57 (2008)
[49] Tahvili, F., Vinck, M., Di Volo, M.: A cortical microcircuit model reveals distinct inhibitory mechanisms of network oscillations and stability. bioRxiv, 2025-02 (2025)

[50] Sheheitli, H., Jirsa, V.: Incorporating slow nmda-type receptors with nonlinear voltage-dependent magnesium block in a next generation neural mass model: derivation and dynamics. Journal of Computational Neuroscience 52(3), 207-222 (2024)

[51] Susin, E., Destexhe, A.: A network model of the modulation of γ oscillations by nmda receptors in cerebral cortex. Eneuro 10(11) (2023)