Asynchronous and Coherent Dynamics in Balanced Excitatory-Inhibitory Spiking Networks

Dynamic excitatory-inhibitory (E-I) balance is a paradigmatic mechanism invoked to explain the irregular low firing activity observed in the cortex. However, we will show that the E-I balance can be at the origin of other regimes observable in the brain. The analysis is performed by combining extensive simulations of sparse E-I networks composed of N spiking neurons with analytical investigations of low dimensional neural mass models. The bifurcation diagrams, derived for the neural mass model, allow us to classify the possible asynchronous and coherent behaviors emerging in balanced E-I networks with structural heterogeneity for any finite in-degree K. Analytic mean-field (MF) results show that both supra and sub-threshold balanced asynchronous regimes are observable in our system in the limit N >> K >> 1. Due to the heterogeneity, the asynchronous states are characterized at the microscopic level by the splitting of the neurons in to three groups: silent, fluctuation, and mean driven. These features are consistent with experimental observations reported for heterogeneous neural circuits. The coherent rhythms observed in our system can range from periodic and quasi-periodic collective oscillations (COs) to coherent chaos. These rhythms are characterized by regular or irregular temporal fluctuations joined to spatial coherence somehow similar to coherent fluctuations observed in the cortex over multiple spatial scales. The COs can emerge due to two different mechanisms. A first mechanism analogous to the pyramidal-interneuron gamma (PING), usually invoked for the emergence of γ-oscillations. The second mechanism is intimately related to the presence of current fluctuations, which sustain COs characterized by an essentially simultaneous bursting of the two populations. We observe period-doubling cascades involving the PING-like COs finally leading to the appearance of coherent chaos. Fluctuation driven COs are usually observable in our system as quasi-periodic collective motions characterized by two incommensurate frequencies. However, for sufficiently strong current fluctuations these collective rhythms can lock. This represents a novel mechanism of frequency locking in neural populations promoted by intrinsic fluctuations. COs are observable for any finite in-degree K, however, their existence in the limit N >> K >> 1 appears as uncertain.

1. Introduction

Cortical neurons are subject to a continuous bombardment from thousands of presynaptic neurons, mostly pyramidal ones, evoking postsynaptic potentials of sub-millivolt or millivolt amplitudes (Destexhe and Paré, 1999; Bruno and Sakmann, 2006; Lefort et al., 2009). This stimulation would induce an almost constant depolarization of the neurons leading to a regular firing, However, cortical neurons fire quite irregularly and with low firing rates (Softky and Koch, 1993). This apparent paradox can be solved by introducing the concept of a balanced network, where excitatory and inhibitory synaptic currents are approximately balanced and the neurons are kept near their firing threshold crossing it at random times (Shadlen and Newsome, 1994, 1998). However, the balance should naturally emerge in the network without fine-tuning of the parameters and the highly irregular firing observed in vivo should be maintained also for a large number of connections (in-degree) K >> 1. This is possible by considering a sparse excitatory-inhibitory (E-I) neural network composed of N neurons and characterized by an average in-degree K << N and by synaptic couplings scaling as $1/\sqrt{K}$ (van Vreeswijk and Sompolinsky, 1996). This scaling as well as many other key predictions of the theory developed in (van Vreeswijk and Sompolinsky, 1996) have been recently confirmed by experiments on a neural culture (Barral and Reyes, 2016). Furthermore, Barral and Reyes (2016) have shown that the major predictions of the seminal theory (van Vreeswijk and Sompolinsky, 1996) also hold under conditions far from the asymptotic limits where K and N are large.

The dynamics usually observable in balanced neural networks is asynchronous and characterized by irregular neural firing joined to stationary firing rates (van Vreeswijk and Sompolinsky, 1996; Monteforte and Wolf, 2010; Renart et al., 2010; Litwin-Kumar and Doiron, 2012; Ullner et al., 2020). However, other asynchronous regimes characterized by sub-Poissonian and super-Poissonian statistics have been reported in balanced homogenous and heterogeneous networks (Lerchner et al., 2006; Ullner et al., 2020). Furthermore, regular and irregular collective oscillations (COs) have been shown to emerge in balanced networks composed of rate models (van Vreeswijk and Sompolinsky, 1996) and of spiking neurons (Brunel, 2000; Ostojic, 2014; di Volo and Torcini, 2018; Ullner et al., 2018; Bi et al., 2020). The balanced asynchronous irregular state has been experimentally observed both in vivo and in vitro (Shu et al., 2003; Haider et al., 2006) and dynamic balance of excitation and inhibition is observable in the neocortex across all states of the wake-sleep cycle, in both human and monkey (Dehghani et al., 2016). However, this is not the unique balanced state observable in neural systems. In particular, balancing of excitation and inhibition appears to be crucial for the emergence of cortical oscillations (Okun and Lampl, 2008; Isaacson and Scanziani, 2011; Le Van Quyen et al., 2016) as well as for the instantaneous modulation of the frequency of gamma oscillations in the hippocampus (Atallah and Scanziani, 2009). Moreover, balancing of excitation and inhibition is essential for the generation of respiratory rhythms in the brainstem (Ramirez and Baertsch, 2018) and the rhythmic activity of irregular firing motoneurons in the spinal cord of the turtle (Berg et al., 2007, 2019).

In this work, we characterize in detail the asynchronous regimes and the emergence of COs (population rhythms) in E-I balanced networks with structural heterogeneity. In particular, we consider sparse random networks of quadratic integrate-and-fire (QIF) neurons (Ermentrout and Kopell, 1986) pulse coupled via instantaneous post synaptic potentials. We compare numerical findings with analytical results obtained in the mean-field (MF) limit by employing an effective low-dimensional neural mass model recently developed for sparse QIF networks (Montbrió et al., 2015; di Volo and Torcini, 2018; Bi et al., 2020).

In the asynchronous regime, our analytical MF predictions are able to reproduce the mean membrane potentials and the population firing rates of the structurally heterogeneous network for any finite K value. Furthermore, in the limit N >> K >> 1, we analytically derive the asymptotic MF values of the population firing rates and the effective input currents. This analysis shows that the system always achieves balanced dynamics, whose supra or sub-threshold nature is determined by the model parameters. Detailed numerical investigations of the microscopic dynamics allow identifying three different groups of neurons, whose activity is essentially controlled by their in-degrees and by the effective input currents.

In the balanced network, we have identified three types of COs depending on the corresponding solution displayed by the neural mass model. The first type, termed O_P emerges in the MF via a Hopf bifurcation (HB) from a stable focus solution. These COs gives rise to collective chaos via a period-doubling sequence of bifurcations. Another type of CO, already reported for pure inhibitory networks (di Volo and Torcini, 2018), denoted as O_F corresponds in the MF to a stable focus characterized by relaxation oscillations toward the fixed point which in the sparse network become noise sustained oscillations due to fluctuations in the input currents. The last type of COs identified in the finite network are named O_S and characterized by abnormally synchronized dynamics among the neurons, and the high level of synchronization prevents their representation in the MF formulation (Montbrió et al., 2015).

O_P and O_S emerge as sustained oscillations in the network via a mechanism resembling that invoked for pyramidal-interneuron gamma (PING) rhythms (Whittington et al., 2011) despite the frequency of these oscillations is not restricted to the γ band. Excitatory neurons start to fire followed by the inhibitory ones and the peak of activity of the excitatory population precedes that of the inhibitory one of a time delay Δt. Furthermore, Δt tends to vanish when the amplitude of the current fluctuations in the network increases. Indeed, for O_F oscillations, which cannot emerge in absence of current fluctuations, no delay has been observed between the activation of excitatory and inhibitory populations. The last important question that we tried to address in our study was whether the COs, observable for any finite K, are still present in the limit N >> K >> 1.

The study is organized as follows. Section 2 is devoted to the introduction of the network model and of the corresponding effective neural mass model, as well as of the microscopic and macroscopic indicators employed to characterize the neural dynamics. In the same section, the stationary solutions for the balanced neural mass model are analytically obtained as finite in-degree expansion and their range of stability is determined. The macroscopic dynamical regimes emerging in our network are analyzed in section 3. In particular, we report bifurcation phase diagrams obtained from the neural mass model displaying the possible dynamical states and network simulations. The focus of this section is on the analysis of the asynchronous balanced state for structurally heterogeneous networks and the emergence of the different types of COs observable at finite in-degrees. A discussion of the obtained results and conclusions are reported in section 4.

2. Models and Dynamical Indicators

2.1. Network Model

We consider two sparsely coupled excitatory and inhibitory populations composed of N^(e) and N⁽ⁱ⁾ QIF neurons, respectively (Ermentrout and Kopell, 1986). The evolution equation for the membrane potentials $v_j^{(e)}$ and $v_j^{(i)}$ of the excitatory and inhibitory neurons can be written as:

$\begin{array}{l}{\tau }_m{\dot{v}}_j^{(e)}={\left(v_j^{(e)}\right)}^2+I^{(e)}+2{\tau }_m{\displaystyle [g^{(ee)}{\displaystyle \sum _{l|t_l^{(n)}<t}}{ϵ}_{jl}^{(ee)}\delta (t-t_l^{(n)})}\\ {\displaystyle -g^{(ei)}{\displaystyle \sum _{k|t_k^{(m)}<t}}{ϵ}_{jk}^{(ei)}\delta (t-t_k^{(m)})]}& (\mathrm{1a})\end{array}$

$\begin{array}{l}{\tau }_m{\dot{v}}_j^{(i)}={\left(v_j^{(i)}\right)}^2+I^{(i)}+2{\tau }_m{\displaystyle [g^{(ie)}{\displaystyle \sum _{l|t_l^{(n)}<t}}{ϵ}_{jl}^{(ie)}\delta (t-t_l^{(n)})}\\ {\displaystyle -g^{(ii)}{\displaystyle \sum _{k|t_k^{(m)}<t}}{ϵ}_{jk}^{(ii)}\delta (t-t_k^{(m)})]}& (\mathrm{1b})\end{array}$

where τ_m = 20 ms is the membrane time constant that we set identical for excitatory and inhibitory neurons, I^(e) (I⁽ⁱ⁾) is the external direct current (DC) acting on excitatory (inhibitory) population, g^(αβ) represents the synaptic coupling strengths between post synaptic neurons in the population α and pre synaptic ones in population β, with α, β ∈ {e, i}. The elements of the adjacency matrices ${ϵ}_{jk}^{(\alpha \beta )}$ are equal to 1 (0) if a connection from a pre synaptic neuron k of population β toward a post synaptic neuron j of population α, exists (or not). Furthermore, $k_j^{(\alpha \beta )}=\sum _k{ϵ}_{jk}^{(\alpha \beta )}$ is the number of pre synaptic neurons in the population β connected to neuron j in population α, or in other terms, its in-degree restricted to population β. The emission of the n-th spike emitted by neuron l of the population α occurs at time $t_l^{(n)}$ whenever the membrane potential $v_l^{(\alpha )}({t_l^{(n)}}^{-})\to \infty$, while the reset mechanism is modeled by setting $v_l^{(\alpha )}({t_l^{(n)}}^{+})\to -\infty$ immediately after the spike emission. The postsynaptic potentials are assumed to be δ-pulses and the synaptic transmissions to be instantaneous. The Equations (1) can be formally rewritten as

$\begin{array}{ll}{\tau }_m{\dot{v}}_j^{(e)}={\left(v_j^{(e)}\right)}^2+i_{eff,j}^{(e)},{\tau }_m{\dot{v}}_j^{(i)}={\left(v_j^{(i)}\right)}^2+i_{eff,j}^{(i)};& (2)\end{array}$

where $i_{eff,j}^{(e)}$ ($i_{eff,j}^{(i)}$) represents the instantaneous excitatory (inhibitory) effective currents, which include the external DC current and the synaptic currents due to the recurrent connections.

We consider the neurons within the excitatory and inhibitory populations as randomly connected, with in-degrees k^(αα) distributed according to a Lorentzian distribution

$\begin{array}{ll}P(k^{(\alpha \alpha )})=\frac{{\Delta }_k^{(\alpha \alpha )}}{{(k^{(\alpha \alpha )}-K^{(\alpha \alpha )})}^2+{{\Delta }_k^{(\alpha \alpha )}}^2}·\frac{1}{\pi }& (3)\end{array}$

peaked at K^(αα) and with a half-width half-maximum (HWHM) ${\Delta }_k^{(\alpha \alpha )}$, this latter parameter measures the level of structural heterogeneity in each population. For simplicity, we set K^(ee) = K⁽ⁱⁱ⁾ ≡ K. Furthermore, we assume that also neurons from a population α are randomly connected to neurons of a different population β ≠ α. However, in this case, we consider no structural heterogeneity with in-degrees fixed to a constant value K^(ei) = K^(ie) = K. We have verified that by considering Erdös-Renyi distributed in-degrees K^(ei) and K^(ie) with average K does not modify the observed dynamical behavior.

The DC current and the synaptic coupling are rescaled with the median in degree as $I^{(\alpha )}=\sqrt{K}{I}_0^{(\alpha )}$ and $g^{(\alpha \beta )}=g_0^{(\alpha \beta )}/\sqrt{K}$, as done in previous studies to obtain a self-sustained balanced dynamics for N >> K >> 1 (van Vreeswijk and Sompolinsky, 1996; Renart et al., 2010; Litwin-Kumar and Doiron, 2012; Kadmon and Sompolinsky, 2015). The structural heterogeneity parameters are rescaled as ${\Delta }_k^{(\alpha \alpha )}={\Delta }_0^{(\alpha \alpha )}\sqrt{K}$ in analogy to Erdös-Renyi networks. The choice of the Lorentzian distribution for the k^(αα) is needed in order to obtain an effective MF description for the microscopic dynamics (di Volo and Torcini, 2018; Bi et al., 2020) as detailed in the next section.

The microscopic activity can be analyzed by considering the inter-spike interval (ISI) distribution as characterized by the coefficient of variation cv_i for each neuron i, which is the ratio between the SD and the mean of the ISIs associated with the train of spikes emitted by the considered neuron. To characterize the macroscopic dynamics of each population α, we measure the average coefficient of variation $C{V}^{(\alpha )}=\sum _{i=1}^{N^{(\alpha )}}c{v}_i/N^{(\alpha )}$, the mean membrane potential $V^{(\alpha )}(t)=\sum _{i=1}^{N^{(\alpha )}}{v}_i^{(\alpha )}(t)/N^{(\alpha )}$, and the population firing rate R^(α)(t), corresponding to the number of spikes emitted within the population α per unit of time and per neuron.

Furthermore, the level of coherence in the neural activity of the population α can be quantified in terms of the following indicator (Golomb, 2007),

$\begin{array}{ll}{\rho }^{(\alpha )}={\left(\frac{{\sigma }_{V^{(\alpha )}}^2}{{\displaystyle {\sum }_{i=1}^{N^{(\alpha )}}}{\sigma }_i^2/N^{(\alpha )}}\right)}^{1/2}& (4)\end{array}$

where ${\sigma }_{V^{(\alpha )}}$ is the SD of the mean membrane potential, ${\sigma }_i^2=\langle {(v_i^{(\alpha )})}^2\rangle -\langle v_i^{(\alpha )}{\rangle }^2$ and 〈·〉 denotes a time average. A perfect synchrony corresponds to ρ^(α) = 1, while an asynchronous dynamics to a vanishing small ${\rho }^{(\alpha )}\approx O(1/\sqrt{N^{(\alpha )}})$.

The frequencies associated with collective motions can be identified by measuring the power spectra S(ν) of the mean membrane potentials V(t) of the whole network. In the case of a periodic motion, the position of the main peak ν_CO represents the frequency of the COs, while for quasi-periodic motions, the spectrum is characterized by many peaks that can be obtained as a linear combination of two fundamental frequencies (ν₁, ν₂). The spectra obtained in the present case, always exhibit a continuous background due to the intrinsic fluctuations present in the balanced network. The power spectra have been obtained by calculating the temporal Fourier transform of V(t) sampled at time intervals of 10 ms. Time traces composed of 10,000 consecutive intervals have been considered to estimate the spectra, which are obtained at a frequency resolution of Δν = 0.01 Hz. Finally, the power spectra have been averaged over five independent realizations of the random network.

The network dynamics are integrated by employing an Euler scheme with time step dt = 0.0001 ms, while time averages and fluctuations are usually estimated on time intervals T_s ≃ 100 s, after discarding transients T_t ≃ 10 s. Usually, we consider networks composed of N^(e) = 10, 000 excitatory and N⁽ⁱ⁾ = 2, 500 inhibitory neurons.

2.2. Effective Neural Mass Model

In this sub-section, we derive a low dimensional effective neural mass formulation for the spiking network (Equation 1) by following Montbrió et al. (2015). In such an article, the authors obtained an exact MF model for a globally coupled heterogeneous population of QIF neurons by generalizing to neural systems using a reduction methodology previously developed for phase-coupled oscillators by Ott and Antonsen (2008). In particular, the neural mass model can be obtained by performing a rigorous mathematical derivation from the original spiking network in the limit N → ∞ by assuming that the heterogeneity present in the network, which can be either neuronal excitabilities or synaptic couplings, are distributed as Lorentzians. This MF reduction methodology gives rise to a neural mass model written in terms of only two collective variables: the mean membrane potential V and the instantaneous population rate R. For sufficiently large network size, the agreement between the simulation results and the neural mass model is impressive as shown in Montbrió et al. (2015) and in several successive publications.

The detailed derivation of the neural mass models from the corresponding spiking networks can be found in Montbrió et al. (2015), in this study, we limit to report its expression for a fully coupled homogeneous network of QIF neurons with synaptic couplings randomly distributed according to a Lorentzian:

$\begin{array}{ll}{\tau }_mṘ=2RV+\frac{\Gamma }{\pi }R& (\mathrm{5a})\end{array}$

$\begin{array}{ll}{\tau }_m\dot{V}=V^2+I+ḡ{\tau }_{m}R-{(\pi {\tau }_{m}R)}^2& (\mathrm{5b})\end{array}$

where ḡ is the median and Γ the HWHM of the Lorentzian distribution of the synaptic couplings.

Such formulation can be applied to the random sparse network studied in this article, in this paper. Indeed, as shown for a single sparse inhibitory population (di Volo and Torcini, 2018; Bi et al., 2020), the quenched disorder associated to the in-degree distribution can be rephrased in terms of random synaptic couplings. Namely, each neuron i in population α is subject to currents of amplitude $g_0^{(\alpha \beta )}{k}_i^{(\alpha \beta )}{R}^{(\beta )}/(\sqrt{K})$ proportional to their in-degrees $k_i^{(\alpha \beta )}$, with β ∈ {e, i}. Therefore, we can consider the neurons as fully coupled but with random values of the couplings distributed as Lorentzian of median $g_0^{(\alpha \beta )}\sqrt{K}$ and HWHM $g_0^{(\alpha \beta )}{\Delta }_0^{(\alpha \beta )}$.

The neural mass model corresponding to the spiking network (Equation 1) can be written as follows:

$\begin{array}{ll}{\tau }_m{Ṙ}^{(e)}=R^{(e)}\left[2V^{(e)}+g_0^{(ee)}\frac{{\Delta }_0^{(ee)}}{\pi }\right]& (\mathrm{6a})\end{array}$

$\begin{array}{l}{\tau }_m{\dot{V}}^{(e)}={\left[V^{(e)}\right]}^2-{\left[\pi R^{(e)}{\tau }_m\right]}^2\\ +\sqrt{K}\left[I_0^{(e)}+(g_0^{(ee)}{R}^{(e)}-g_0^{(ei)}{R}^{(i)}){\tau }_m\right]& (\mathrm{6b})\end{array}$

$\begin{array}{ll}{\tau }_m{Ṙ}^{(i)}=R^{(i)}\left[2V^{(i)}+g_0^{(ii)}\frac{{\Delta }_0^{(ii)}}{\pi }\right]& (\mathrm{6c})\end{array}$

$\begin{array}{l}{\tau }_m{\dot{V}}^{(i)}={\left[V^{(i)}\right]}^2-{\left[\pi R^{(i)}{\tau }_m\right]}^2\\ +\sqrt{K}\left[I_0^{(i)}+(g_0^{(ie)}{R}^{(e)}-g_0^{(ii)}{R}^{(i)}){\tau }_m\right];& (\mathrm{6d})\end{array}$

where we have set ${\Delta }_0^{(ei)}={\Delta }_0^{(ie)}=0$, since we have assumed that the connections among neurons of different populations are random but with a fixed in-degree K^(ei) = K^(ie) = K.

2.2.1. Stationary Solutions

The stationary solutions $\left\{{\overline{V}}^{(e)},{\overline{V}}^{(i)},{\overline{R}}^{(e)},{\overline{R}}^{(i)}\right\}$ of Equation (6) can be explicitly obtained for the mean membrane potentials as

$\begin{array}{ll}{\overline{V}}^{(e)}=-\frac{g_0^{(ee)}{\Delta }_0^{(ee)}}{2\pi },{\overline{V}}^{(i)}=-\frac{g_0^{(ii)}{\Delta }_0^{(ii)}}{2\pi };& (7)\end{array}$

while the instantaneous population rates are the solutions of the following quadratic system

$\begin{array}{ll}{g}_0^{(ee)}{\overline{R}}^{(e)}{\tau }_m-g_0^{(ei)}{\overline{R}}^{(i)}{\tau }_m=-I_0^{(e)}+\epsilon \left\{{\left[\pi {\overline{R}}^{(e)}{\tau }_m\right]}^2-{\left[{\overline{V}}^{(e)}\right]}^2\right\}& (\mathrm{8a})\end{array}$

$\begin{array}{ll}{g}_0^{(ie)}{\overline{R}}^{(e)}{\tau }_m-g_0^{(ii)}{\overline{R}}^{(i)}{\tau }_m=-I_0^{(i)}+\epsilon \left\{{\left[\pi {\overline{R}}^{(i)}{\tau }_m\right]}^2-{\left[{\overline{V}}^{(i)}\right]}^2\right\}& (\mathrm{8b})\end{array}$

where $\epsilon =1/\sqrt{K}$ is a smallness parameter taking in to account finite in-degree corrections. It is interesting to notice that the parameters controlling the structural heterogeneity ${\Delta }_0^{(ii)}$ and ${\Delta }_0^{(ee)}$ fix the stationary values of the mean membrane potentials reported in Equation (7). The solutions of Equation (8) can be exactly obtained and the associated bifurcations analyzed by employing the software XPP AUTO developed for orbit continuation (Ermentrout, 2007).

For sufficiently large K, one can obtain analytic approximations of the solution of Equation (8) by expanding the population rates as follows:

$\begin{array}{ll}{\overline{R}}^{(\alpha )}={\overline{R}}_0^{(\alpha )}+\epsilon {\overline{R}}_1^{(\alpha )}+{\epsilon }^2{\overline{R}}_2^{(\alpha )}+{\epsilon }^3{\overline{R}}_3^{(\alpha )}+\dots \alpha \in \left\{e,i\right\},& (9)\end{array}$

by inserting these expressions in Equation (8), and finally by solving order by order in ε.

The solutions at any order can be written as follows:

$\begin{array}{l}{\overline{R}}_k^{(e)}{\tau }_m=\frac{N_k^{(e)}{g}_0^{(ii)}-N_k^{(i)}{g}_0^{(ei)}}{g_0^{(ei)}{g}_0^{(ie)}-g_0^{(ee)}{g}_0^{(ii)}},\\ {\overline{R}}_k^{(i)}{\tau }_m=\frac{N_k^{(e)}{g}_0^{(ie)}-N_k^{(i)}{g}_0^{(ee)}}{g_0^{(ei)}{g}_0^{(ie)}-g_0^{(ee)}{g}_0^{(ii)}};& (10)\end{array}$

where,

$\begin{array}{ll}{N}_0^{(\alpha )}=I_0^{(\alpha )},N_1^{(\alpha )}={\left[{\overline{V}}^{(\alpha )}\right]}^2-{\left[\pi {\overline{R}}_0^{(\alpha )}{\tau }_m\right]}^2& (\mathrm{11a})\end{array}$

$\begin{array}{ll}{N}_{2j}^{(\alpha )}=-2{\left[\pi {\tau }_m\right]}^2{\displaystyle \sum _{k=1}^j}\left[{\overline{R}}_{k-1}^{(\alpha )}{\overline{R}}_{2j-k}^{(\alpha )}\right]& (\mathrm{11b})\end{array}$

$\begin{array}{l}{N}_{2j+1}^{(\alpha )}=-2{\left[\pi {\tau }_m\right]}^2\left\{\left[{\displaystyle \sum _{k=1}^j}{\overline{R}}_{k-1}^{(\alpha )}{\overline{R}}_{2j+1-k}^{(\alpha )}\right]+\frac{1}{2}{\left[{\overline{R}}_j^{(\alpha )}\right]}^2\right\}\\ \text{for }j\ge 1& (\mathrm{11c})\end{array}$

The systems (Equation 10) with parameters given by Equation (11) can be resolved recursively for any order and the final solution obtained from the expression (Equation 9). The zeroth-order approximation, valid in the limit K → ∞, corresponds to the usual solution found for rate models in the balanced state (van Vreeswijk and Sompolinsky, 1996; Rosenbaum and Doiron, 2014), such solution is physical whenever one of the following inequalities is satisfied

$\begin{array}{ll}\frac{I_0^{(e)}}{I_0^{(i)}}>\frac{g_0^{(ei)}}{g_0^{(ii)}}>\frac{g_0^{(ee)}}{g_0^{(ie)}},\frac{I_0^{(e)}}{I_0^{(i)}}<\frac{g_0^{(ei)}}{g_0^{(ii)}}<\frac{g_0^{(ee)}}{g_0^{(ie)}};& (12)\end{array}$

which ensure the positive sign of ${\overline{R}}_0^{(e)}$ and ${\overline{R}}_0^{(i)}$. The zeroth-order solution does not depend on the structural heterogeneity, since the ratio Δ^(αα)/K vanishes in the limit K → ∞. It should be stressed that this ratio does not correspond to the coefficient of variation introduced in Landau et al. (2016) to characterize the in-degree distribution. This is because we are considering a Lorentzian distribution, where the average and the SD are not even defined. Moreover, already the first-order corrections depend on ${\Delta }_0^{(\alpha \alpha )}$.

To characterize the level of balance in the system, one usually estimates the values of the effective input currents $i_{eff,j}^{(e)}$ and $i_{eff,j}^{(i)}$ driving the neuron dynamics. These at a population level can be rewritten as

$\begin{array}{l}{I}_{eff}^{(e)}=\sqrt{K}\left[I_0^{(e)}+{\tau }_m(g_0^{(ee)}{R}^{(e)}-g_0^{(ei)}{R}^{(i)})\right],\\ I_{eff}^{(i)}=\sqrt{K}\left[I_0^{(i)}+{\tau }_m(g_0^{(ie)}{R}^{(e)}-g_0^{(ii)}{R}^{(i)})\right].& (13)\end{array}$

In a balanced state, these quantities should not diverge with K, instead, they should approach some constant value. For an asynchronous state we can estimate analytically, within our MF formulation, the values of the effective currents in the limit K → ∞. These read as

$\begin{array}{l}{I}_a^{(e)}={\tau }_m\left[g_0^{(ee)}{\overline{R}}_1^{(e)}-g_0^{(ei)}{\overline{R}}_1^{(i)}\right],\\ I_a^{(i)}={\tau }_m\left[g_0^{(ie)}{\overline{R}}_1^{(e)}-g_0^{(ii)}{\overline{R}}_1^{(i)}\right].& (14)\end{array}$

It should be noticed that these asymptotic values depend on the first-order corrections to the balanced solution (Equation 10). Therefore, they depend not only on the synaptic couplings $g_0^{(\alpha \beta )}$ and on the external DC currents but also on the parameters ${\Delta }_0^{(\alpha \alpha )}$ controlling the structural heterogeneities.

Depending on the parameter values, the currents $I_a^{(\alpha )}$ can be positive or negative, thus, indicating a balanced dynamics where most parts of the neurons are supra or below the threshold, respectively. Usually, in order to obtain a stationary state characterized by a low rate and a Poissonian statistic, as observed in the cortex, one assumes that the excitation and inhibition nearly cancel. So that the mean membrane potential remains slightly below the threshold, and the neurons can fire occasionally due to the input current fluctuations (van Vreeswijk and Sompolinsky, 1996; Brunel, 2000). However, as pointed out in Lerchner et al. (2006), this is not the only possible scenario for a balanced state. In particular, the authors have developed a self-consistent MF theory for balanced Erdös-Renyi networks made of heterogeneous Leaky Integrate-and-Fire (LIF) neurons. In this context, they have shown that Poisson-like dynamics are visible only at intermediate synaptic couplings. While mean driven dynamics are expected for low couplings, and at large couplings bursting behaviors appear in the balanced network. Recently, analogous dynamical behaviors have been reported also for a purely inhibitory heterogeneous LIF network (Angulo-Garcia et al., 2017). These findings are consistent with the results in Lerchner et al. (2006), where the inhibition is indeed predominant in the balanced regime.

2.2.2. Lyapunov Analysis

To analyze the linear stability of generic solutions of Equation (6), we have estimated the corresponding Lyapunov spectrum (LS) {λ_k} (Pikovsky and Politi, 2016). This can be done by considering the time evolution of the tangent vector δ = {δR^(e), δV^(e), δR⁽ⁱ⁾, δV⁽ⁱ⁾}, that is ruled by the linearization of the Equation (6), namely

$\begin{array}{ll}{\tau }_m\delta {Ṙ}^{(e)}=\left[2V^{(e)}+g_0^{(ee)}\frac{{\Delta }_0^{(ee)}}{\pi }\right]\delta R^{(e)}+2R^{(e)}\delta V^{(e)}& (\mathrm{15a})\end{array}$

$\begin{array}{l}{\tau }_m\delta {\dot{V}}^{(e)}=2V^{(e)}\delta V^{(e)}-2{(\pi {\tau }_m)}^2{R}^{(e)}\delta R^{(e)}\\ +\sqrt{K}{\tau }_m\left[g_0^{(ee)}\delta R^{(e)}-g_0^{(ei)}\delta R^{(i)}\right]& (\mathrm{15b})\end{array}$

$\begin{array}{ll}{\tau }_m\delta {Ṙ}^{(i)}=\left[2V^{(i)}+g_0^{(ii)}\frac{{\Delta }_0^{(ii)}}{\pi }\right]\delta R^{(i)}+2R^{(i)}\delta V^{(i)}& (\mathrm{15c})\end{array}$

$\begin{array}{l}{\tau }_m\delta {\dot{V}}^{(i)}=2V^{(i)}\delta V^{(i)}-2{(\pi {\tau }_m)}^2{R}^{(i)}\delta R^{(i)}\\ +\sqrt{K}{\tau }_m\left[g_0^{(ie)}\delta R^{(e)}-g_0^{(ii)}\delta R^{(i)}\right].& (\mathrm{15d})\end{array}$

In this case, the LS is composed of four Lyapunov exponents (LEs) {λ_k} with k = 1, …, 4, which quantify the average growth rates of infinitesimal perturbations along the orthogonal manifolds. The LEs can be estimated as follows:

$\begin{array}{ll}{\lambda }_k={\displaystyle \underset{t\to \infty }{lim}}\frac{1}{t}log\frac{|{\delta }_k(t)|}{|{\delta }_k(0)|},& (16)\end{array}$

where the tangent vectors δ_k are maintained ortho-normal during the time evolution by employing a standard technique introduced in Benettin et al. (1980). The autonomous system will be chaotic for λ₁ > 0, while a periodic (two frequency quasi-periodic) dynamics will be characterized by λ₁ = 0 (λ₁ = λ₂ = 0) and a fixed point by λ₁ < 0.

In order to estimate the LS for the neural mass model, we have integrated the direct and tangent space evolution with a Runge-Kutta 4th order integration scheme with dt = 0.01 ms, for a duration of 200 s, after discarding a transient of 10 s.

2.2.3. Linear Stability of Stationary Solutions

The linear stability of the stationary solutions $\left\{{\overline{V}}^{(e)},{\overline{V}}^{(i)},{\overline{R}}^{(e)},{\overline{R}}^{(i)}\right\}$ can be analyzed by solving the eigenvalue problem for the linear Equations (15) estimated for stationary values of the mean membrane potentials and of the population firing rates. This approach gives rise to a fourth-order characteristic polynomial of the complex eigenvalues ${\Lambda }^{(k)}={\Lambda }_R^{(k)}+i{\Lambda }_I^{(k)}$ with k = 1, …, 4. The stability of the fixed point is controlled by the maximal ${\Lambda }_R^{(k)}$, whenever it is positive (negative), the stationary solution is unstable (stable). The nature of the fixed point is determined by ${\Lambda }_I^{(k)}$, if the imaginary parts of the eigenvalues are all zero, we have a node, otherwise a focus. Due to the fact that the coefficients of the characteristic polynomial are real, the eigenvalues are real or if complex they appear in complex conjugates couples ${\Lambda }_R^{(j)}\pm i{\Lambda }_I^{(k)}$. Therefore, the relaxation toward the fixed point is characterized by one or two frequencies ${\nu }_k={\Lambda }_I^{(k)}/(2\pi )$. These latter quantities, as discussed in detail in the following, can give good predictions for the frequencies ν_CO of fluctuation driven COs observable for the same parameters in the network dynamics.

In the limit K >> 1, we can approximate the linear stability (Equations 15) as follows:

$\begin{array}{ll}{\tau }_m\delta {Ṙ}^{(e)}=2{\overline{R}}_0^{(e)}\delta V^{(e)}& (\mathrm{17a})\end{array}$

$\begin{array}{ll}{\tau }_m\delta {\dot{V}}^{(e)}=\sqrt{K}{\tau }_m\left[g_0^{(ee)}\delta R^{(e)}-g_0^{(ei)}\delta R^{(i)}\right]& (\mathrm{17b})\end{array}$

$\begin{array}{ll}{\tau }_m\delta {Ṙ}^{(i)}=2{\overline{R}}_0^{(i)}\delta V^{(i)}& (\mathrm{17c})\end{array}$

$\begin{array}{ll}{\tau }_m\delta {\dot{V}}^{(i)}=\sqrt{K}{\tau }_m\left[g_0^{(ie)}\delta R^{(e)}-g_0^{(ii)}\delta R^{(i)}\right];& (\mathrm{17d})\end{array}$

where we have considered the zeroth-order approximation for the population rates ${\overline{R}}_0^{(e)}$ and ${\overline{R}}_0^{(i)}$.

In this case, the complex eigenvalues Λ^(k) are given by the following expression:

$\begin{array}{l}{\left[{\Lambda }^{(k)}\right]}^2={\displaystyle \frac{\sqrt{K}}{{\tau }_m}[\left(g_0^{(ee)}{\overline{R}}_0^{(e)}-g_0^{(ii)}{\overline{R}}_0^{(i)}\right)}\\ {\displaystyle \pm \sqrt{{\left(g_0^{(ee)}{\overline{R}}_0^{(e)}+g_0^{(ii)}{\overline{R}}_0^{(i)}\right)}^2-4g_0^{(ei)}{g}_0^{(ie)}{\overline{R}}_0^{(e)}{\overline{R}}_0^{(i)}}]}.& (18)\end{array}$

From Equation (18), it is evident that Λ^(k) ∝ (K)^1/4, and by assuming $I_0^{(i)}\propto I_0^{(e)}$, as we will do in this study, we also have that ${\Lambda }^{(k)}\propto {(I_0^{(e)})}^{1/2}$. Therefore, for a focus solution, we will have the following scaling relation for the relaxation frequencies for sufficiently large K:

$\begin{array}{ll}{\nu }_k^R=\frac{{\Lambda }_I^{(k)}}{2\pi }\propto \sqrt{I_0^{(e)}{K}^{1/2}};& (19)\end{array}$

This scaling is analogous to that found for purely inhibitory QIF networks in di Volo and Torcini (2018). In van Vreeswijk and Sompolinsky (1996), it has been found that the eigenvalues, characterizing the stability of the asynchronous state, scale proportionally to $\sqrt{K}$, therefore, the convergence (divergence) from the stationary stable (unstable) solution is somehow slower with K in our model. This is due to the presence in our MF of an extra macroscopic variable, the mean membrane potential, with respect to the usual rate models.

3. Results

3.1. Phase Diagrams

In this sub-section, we will investigate the possible dynamical regimes emerging in our model by employing its neural mass formulation. In particular, the dynamics of the neural mass model (Equation 6) take place in a four-dimensional space {R^(e), V^(e), R⁽ⁱ⁾, V⁽ⁱ⁾} and it depends on nine parameters, namely on the four synaptic coupling strengths $\left\{g_0^{(ee)},g_0^{(ei)},g_0^{(ii)},g_0^{(ie)}\right\}$, the two external stimulation currents $\left\{I_0^{(e)},I_0^{(i)}\right\}$, the median in-degree K, and the HWHM of the two distributions of the in-degrees $\left\{{\Delta }_0^{(ee)},{\Delta }_0^{(ii)}\right\}$.

However, in order to reduce the space of parameters to investigate and at the same time to satisfy the inequalities (Equation 12) required for the existence of a balanced state in the large K limit, we fix the inhibitory DC current as $I_0^{(i)}=I_0^{(e)}/1.02$ and the synaptic couplings as $g_0^{(ee)}=0.27$, $g_0^{(ii)}=0.953939$, $g_0^{(ie)}=0.3$, and $g_0^{(ei)}=0.96286$ analogously to what was done in Monteforte and Wolf (2010). Therefore, we are left with four control parameters, namely ${\Delta }_0^{(ee)}$, ${\Delta }_0^{(ii)}$, $I_0^{(e)}$, and K, that we will vary to investigate the possible dynamical states.

Three bidimensional bifurcation diagrams for the neural mass model (Equation 6) are reported in Figure 1 for the couple of parameters $(I_0^{(e)},{\Delta }_0^{(ee)})$, $(K,{\Delta }_0^{(ee)})$, and $({\Delta }_0^{(ii)},{\Delta }_0^{(ee)})$. From the bifurcation analysis, we have identified five different dynamical states for the excitatory population: namely, (I) an unstable focus; (II) a stable focus coexisting with an unstable limit cycle (LC); (III) a stable node; (IV) a stable limit cycle coexisting with an unstable focus; and (V) a chaotic regime. For the analysis reported in the following, it is important to remark that the stable foci are usually associated with four complex eigenvalues arranged in complex conjugate couples, therefore, the relaxation toward a stable focus is characterized by two frequencies (ν₁, ν₂) corresponding to the complex parts of the eigenvalues. In region (III), the macroscopic fixed point is characterized by two real eigenvalues and a couple of complex conjugated ones. Thus, the relaxation toward the macroscopic node is, in this case, guided by a single relaxation frequency. The inhibitory population, reveals the same bifurcation structure as the excitatory one, apart from an important difference: the inhibitory population never displays stable nodes. Therefore, the region (III) for the inhibitory population is also a region of type (II).

FIGURE 1

Figure 1. Bifurcation diagrams of the neural mass model. The bifurcation diagrams concern the dynamical state exhibited by the excitatory population in the bidimensional parameter spaces $(I_0^{(e)},{\Delta }_0^{(ee)})$ (A), $(K,{\Delta }_0^{(ee)})$ (B), and $({\Delta }_0^{(ii)},{\Delta }_0^{(ee)})$ (C). The regions marked by Roman numbers correspond to the following collective solutions: (I) an unstable focus; (II) a stable focus coexisting with an unstable limit cycle (LC); (III) a stable node; (IV) an unstable focus coexisting with a stable LC; (V) a chaotic dynamics. The green solid line separates the regions with a stable node (III) and a stable focus (II). The blue solid (dashed) curve is a line of super-critical (sub-critical) Hopf bifurcations (HBs), and the red one of saddle-node (SN) bifurcations of LCs. The yellow curve denotes the period doubling (PD) bifurcation lines. In (C), we also report the coherence indicator ρ^(e) (Equation 4) estimated from the network dynamics with N^(e) = 10,000 and N⁽ⁱ⁾ = 2,500. The dashed lines in (A) indicate the parameter cuts we will consider in Figures 4, 5 (black) and Figure 7 (purple), while the open circles in (A,B) denote the set of parameters employed in Figure 11. In the three panels, the inhibitory DC current and the synaptic couplings are fixed to $I_0^{(i)}=I_0^{(e)}/1.02$, $g_0^{(ee)}=0.27$, $g_0^{(ii)}=$0.953939, $g_0^{(ie)}=0.3$, $g_0^{(ei)}=$0.96286; other parameters: (A) K=1,000, ${\Delta }_0^{(ii)}=0.3$, (B) $I_0^{(e)}=$0.001, ${\Delta }_0^{(ii)}=0.3$, (C) K= 1,000, and $I_0^{(e)}=0.1$.

As shown in Figures 1A,B, for fixed ${\Delta }_0^{(ii)}$ and for low values of the structural heterogeneity ${\Delta }_0^{(ee)}$ and of the excitatory DC current $I_0^{(e)}$, one observes a stable node (III) that becomes a stable focus (II) by increasing ${\Delta }_0^{(ee)}$, these transitions are signaled as green solid lines in Figure 1. By further increasing the degree of heterogeneity ${\Delta }_0^{(ee)}$, the stable focus gives rise to COs (IV) via a super-critical HB (blue solid lines). Depending on the values of K and $I_0^{(e)}$, one can have the emergence of chaotic behaviors (V) via a period doubling (PD) cascade (yellow solid lines). For sufficiently large ${\Delta }_0^{(ee)}$, the COs disappear via a saddle-node (SN) bifurcation of limit cycles (LC) (red solid lines) and above the SN line, the only remaining solution is an unstable focus (I).

As shown in Figure 1A, for fixed structural heterogeneities, the increase of $I_0^{(e)}$ leads to the disappearance of the stable focus (II) via a sub-critical HB (dashed blue line). The dependence of the observed MF solutions on the in-degree K is reported in Figure 1B for a current $I_0^{(e)}=0.001$, and it is not particularly dramatic apart from for the emergence of a chaotic region (V) from a CO regime (IV).

In order to observe the emergence of COs (IV) from the destabilization of a node solution (III), we should vary the structural inhibitory heterogeneity ${\Delta }_0^{(ii)}$, as shown in Figure 1C. Indeed, for sufficiently low ${\Delta }_0^{(ii)}$ and ${\Delta }_0^{(ee)}$, we can observe a super-critical bifurcation line from a node to a stable LC. From this analysis, it emerges that the excitatory heterogeneity has an opposite effect with respect to the inhibitory one, indeed by increasing ${\Delta }_0^{(ee)}$, the value of ρ^(e) increases indicating the presence of more synchronized COs. This effect is due to the fact that the increase of ${\Delta }_0^{(ee)}$ leads to more and more neurons with large $k_j^{(ee)}>>K$, therefore, receiving higher and higher levels of recurrent excitation. These neurons are definitely supra-threshold and drive the activity of the network toward coherent behaviors.

In order to understand the limits of our MF formulation, it is of particular interest to compare the network simulations with the MF phase diagram. To this aim, we report in Figure 1C, the coherence indicator ρ^(e) (Equation 4) estimated from the network dynamics. The indicator ρ^(e) reveals that no COs are present in the region (III), where the MF displays a stable node, however, COs emerge in all the other MF regimes for sufficiently low ${\Delta }_0^{(ii)}<1$. The presence of COs is expected from the MF analysis only in the regions (IV) and (V), but neither in (II) where the MF forecasts the existence of a stable focus nor in (I) where no stable solutions are envisaged. The origin of the discrepancies among the MF and the network simulations in the region (II) is due to the fact that the considered neural mass neglects the dynamical fluctuations in the input currents present in the original networks, which can give rise to noise induced COs (Goldobin et al., 2021). However, as shown in di Volo and Torcini (2018) and Bi et al. (2020) for purely inhibitory populations, the analysis of the neural mass model can still give relevant information on the network dynamics. In particular, the frequencies of the fluctuation induced COs observable in the network simulations can be well estimated from the frequencies (ν₁, ν₂) of the relaxation oscillations toward the stable MF focus. The lack of agreement between MF and network simulations in the region (I) is due to finite size effects, indeed in this case, the system tends to fully synchronize. Therefore, in the network, one observes highly synchronized COs characterized by population firing rates that diverge for increasing K and N and the MF is unable to reproduce these unrealistic solutions (Montbrió et al., 2015).

On the basis of these observations, we can classify the COs observable in the network in three different types accordingly to the corresponding MF solutions: O_P, when in the MF we observe periodic, quasi-periodic, or chaotic collective solutions in regions (IV) and (V); O_F, when the MF displays relaxation oscillations toward the stable focus in regions (II) and (III), that in the sparse network become noise sustained oscillations due to fluctuations in the input currents; O_S, when the MF fully synchronizes as in region (I).

In the following sub-sections, we will analyze the macroscopic dynamics of the E-I network of QIF neurons in order to test the predictions of the effective neural mass mode for asynchronous and coherent dynamics. In this latter case, we will focus on the three types of identified COs: namely, O_P, O_F, and O_S. These can manifest as periodic, quasi-periodic, and chaotic solutions as we will see by examining two main scenarios indicated as dashed horizontal lines in Figure 1A corresponding to the transition to chaos (black dashed line) and the emergence of abnormal synchronization from a stable focus (purple dashed line).

3.2. Asynchronous Regimes

We will first consider a situation where the network dynamics remains asynchronous for any value of the median in-degree K, this occurs for sufficiently high structural inhibitory heterogeneities ${\Delta }_0^{(ii)}$ and external DC currents as shown in Figures 1B,C for E-I networks and as reported in di Volo and Torcini (2018) for purely inhibitory populations. If the population dynamics are asynchronous, we expect that at an MF level, the system will converge toward a stationary state corresponding to a stable equilibrium. Therefore, we have compared the results of the network simulations with the stationary rates $({\overline{R}}^{(e)},{\overline{R}}^{(i)})$ solutions of Equation (6). As shown in Figures 2A,B, the macroscopic activity of the excitatory and inhibitory populations is well reproduced by the fixed point solutions (Equation 8) in a wide range of values of the in-degrees 10 ≤ K ≤ 10⁴. This is particularly true for the inhibitory population, while at low K < 100, the excitatory firing rate is slightly underestimated by the macroscopic solution ${\overline{R}}^{(e)}$. Due to our choice of parameters, the average inhibitory firing rate is larger than the excitatory one for K > 100. This is consistent with experimental data reported for the barrel cortex of behaving mice (Gentet et al., 2010) and other cortical areas (Mongillo et al., 2018). Moreover, the rates have a non-monotonic behavior with K with a maximum at K ≃ 450 (K ≃ 2,500) for excitatory (inhibitory) neurons. As expected, the balanced state solutions ${\overline{R}}_0^{(e)}=3.18$ Hz and ${\overline{R}}_0^{(i)}\simeq 11.28$ Hz (dashed horizontal lines) are approached only for sufficiently large K >> 1. In Figures 1A,B are reported also the first (second) order approximation ${\overline{R}}_0^{(e)}+\epsilon {\overline{R}}_1^{(e)}$ (${\overline{R}}_0^{(e)}+\epsilon {\overline{R}}_1^{(e)}+{\epsilon }^2{\overline{R}}_2^{(e)}$) given by Equation (10). These approximations reproduce quite well the complete solutions already at K ≥ 10⁴.

FIGURE 2

Figure 2. Asynchronous dynamics: Instantaneous population rate R^(e) (R⁽ⁱ⁾) of excitatory (inhibitory) neurons in function of the median in-degree K are shown in (A,B). The effective input currents $I_{eff}^{(e)}$ ($I_{eff}^{(i)}$) given by Equations (13) are reported in (C) and the fluctuations of the input currents $\Delta I_{eff}^{(e)}$ ($\Delta I_{eff}^{(i)}$), as obtained from Equations (20), in (D). Red (blue) color refers to excitatory inhibitory population. The solid continuous lines represent the value obtained by employing the exact MF solutions ${\overline{R}}^{(x)}$ of Equation (8), the dotted (dash-dotted) lines correspond to the first (second) order approximation ${\overline{R}}_0^{(x)}+\epsilon {\overline{R}}_1^{(x)}$ (${\overline{R}}_0^{(x)}+\epsilon {\overline{R}}_1^{(x)}+{\epsilon }^2{\overline{R}}_2^{(x)}$) and the dashed horizontal lines to the zeroth-order one ${\overline{R}}_0^{(x)}$ in (A,B,D), and to $I_a^{(x)}$ in (C) with x = e, i. The circles correspond to data obtained from numerical simulations of N^(e) = N⁽ⁱ⁾ = 10,000 neurons for K < 4,096, N^(e) = N⁽ⁱ⁾ = 20,000 for K = 4,096,8,192 and N^(e) = N⁽ⁱ⁾ = 30,000 for K > 8, 192, averaging the population rates over a window of T = 40 s, after discarding a transient of T = 60 s. The error bars in (A,B) are obtained as the SD (over the time window T) of the population rates, while the average CV of neurons is around 0.15 for all the reported simulations. Synaptic couplings and the ratio between the currents are fixed as stated in sub-section 3.1, other parameters are ${\Delta }_0^{(ii)}=1$, ${\Delta }_0^{(ee)}=2.5$, and $I_0^{(e)}=0.2$. The values of the asymptotic solutions (dashed lines) are : in (A,B) ${\overline{R}}_0^{(e)}=3.18$ Hz and ${\overline{R}}_0^{(i)}=11.28$ Hz, respectively; in (C) $I_a^{(e)}=$ 0.0284 and $I_a^{(i)}\simeq 0.4791$; in (D) $\Delta I_{eff}^{(e)}=$ 0.4623 and $\Delta I_{eff}^{(i)}=$ 0.4593.

Let us now consider the effective input currents (Equation 13), these are reported in Figure 2C vs. the median in-degree. As expected, for increasing K, the MF estimations of the effective currents (solid lines) converge to the asymptotic values $I_a^{(e)}\simeq$ 0.0284 and $I_a^{(i)}\simeq$ 0.4791 (dashed lines) for our choice of parameters. For the excitatory population, the asymptotic value of the effective input current is essentially zero, while for the inhibitory population it is positive. These results suggest that for the considered choice of parameters the dynamics of both populations will be balanced, since the quantities $I_a^{(e)}$ and $I_a^{(i)}$ do not diverge with K, however, at a macroscopic level, the excitatory population will be at the threshold, while the inhibitory one will be supra-threshold. For comparison, we have estimated $I_{eff}^{(\alpha )}$ also from the direct the network simulations (circles) for 16 ≤ K ≤ 16,384. These estimations disagree with the MF results already for K > 1,000. This is despite the fact that the population firing rates in the network are very well captured by the MF estimations at large K, as shown in Figures 2A,B. These large differences in the effective input currents are the effect of small discrepancies at the level of firing rates enhanced by the multiplicative factor $\sqrt{K}$ appearing in Equations (13). However, from the network simulations, we observe that the effective currents approach values smaller than the asymptotic ones $I_a^{(e)}$ and $I_a^{(i)}$ obtained from the neural mass model. In particular, despite the fact that from finite K simulations, it is difficult to extrapolate the asymptotic behaviors, it appears that $I_{eff}^{(e)}$ approaches a small negative value for K >> 1, while $I_{eff}^{(i)}$ converges to some finite positive value. In the following, we will see the effect of these different behaviors on microscopic dynamics. The origin of the reported discrepancies should be related to the presence of current fluctuations in the network that are neglected in the MF formulation.

The relevance of the current fluctuations for the network dynamics can be appreciated by estimating their amplitudes within a Poissonian approximation, as follows

$\begin{array}{l}\Delta I_{eff}^{(e)}=\sqrt{{\tau }_m\left[{\left(g_0^{(ee)}\right)}^2{R}^{(e)}+{\left(g_0^{(ei)}\right)}^2{R}^{(i)}\right]}\\ \Delta I_{eff}^{(i)}=\sqrt{{\tau }_m\left[{\left(g_0^{(ie)}\right)}^2{R}^{(e)}+{\left(g_0^{(ii)}\right)}^2{R}^{(i)}\right]}& (20)\end{array}$

These have been evaluated by assuming that each neuron receives on average K excitatory and inhibitory spike trains characterized by Poissonian statistics with average rates R^(e) and R⁽ⁱ⁾. However, we have neglected in the above estimation the variability of the in-degrees of each neuron. As shown in Figure 2D, these fluctuations are essentially identical for excitatory and inhibitory neurons and coincide with the MF results. In the limit K >> 1, they converge to the asymptotic values $\Delta I_{eff}^{(e)}\simeq$ 0.4623 and $\Delta I_{eff}^{(i)}\simeq$ 0.4593 (green dashed lines). It is evident that already for K > 1,000, the amplitudes of the fluctuations are of the same order or larger than the effective input currents. Thus, suggesting that the fluctuations have indeed a relevant role in determining the network dynamics and that one would observe Poissonian or sub-Poissonian dynamics for the neurons, whenever $I_a^{(\alpha )}$ is sub-threshold or supra-threshold (Lerchner et al., 2006).

In order to understand how the in-degree heterogeneity influences the network dynamics at a microscopic level, we examine the dynamics of active neurons in the function of their total in-degree $k_j^{(tot)}$. This is defined for excitatory (inhibitory) neurons as $k_j^{(tot)}=k_j^{(ee)}+k_j^{(ei)}$ ($k_j^{(tot)}=k_j^{(ii)}+k_j^{(ie)}$). Furthermore, a neuron is considered as active if it has fired at least once during the whole simulation time T_t + T_s = 100 s, therefore, if it has a firing rate larger than 0.01 Hz. As shown in Figures 3A,B, the probability distribution function (PDF) of active neurons is skewed toward values $k_j^{(tot)}>2K$ ($k_j^{(tot)}<2K$) for excitatory (inhibitory) neurons. These results reflect the fact that the excitatory (inhibitory) neurons with low (high) recurrent in-degrees $k_j^{(ee)}<<K$ ($k_j^{(ii)}>>K$) are driven below the threshold by the inhibitory activity, that is predominant in the network since R⁽ⁱ⁾ > R^(e), $g_0^{(ei)}>g_0^{(ee)}$, and $g_0^{(ii)}>g_0^{(ie)}$. The number of silent neurons for K > 1,024 is of the order of 6-10% for both inhibitory and excitatory populations, in agreement with experimental results for the barrel cortex of mice (O'Connor et al., 2010), where a fraction of 10% of neurons was identified as silent with a firing rate slower than 0.0083 Hz. It should be remarked that all the population averages we report include the silent neurons.

FIGURE 3

Figure 3. Asynchronous dynamics: Probability distribution functions (PDFs) of the total in-degrees $k_j^{(tot)}$ for excitatory (A) and inhibitory (B) active neurons for K = 16, 384. (C,D) Firing rates of the excitatory (inhibitory) neurons $r_j^{(e)}$ ($r_j^{(i)}$) vs. their total in-degrees $k_j^{(tot)}-2K$ symbols refer to K = 1,024 (red), K = 4,096 (blue), and K = 16,384 (green). The inset in (D) is an enlargement of the panel displaying the firing rates over the entire scale $k_j^{(tot)}-2K$. The magenta dashed lines in (C,D) represent the balanced state solution $({\overline{R}}_0^{(e)},{\overline{R}}_0^{(i)})$. (E,F) PDF of the excitatory (inhibitory) firing rates $r_j^{(e)}$ ($r_j^{(i)}$) for K = 16, 384, the solid (dashed) line refers to the MF results ${\overline{R}}^{(x)}$ (${\overline{R}}_0^{(x)}$) with x = e, i. The red (blue) solid line refers to a log-normal fit to the excitatory (inhibitory) PDF with mean 8.8 Hz (17.5 Hz) and SD of 3.8 Hz (2.3 Hz). The parameters are the same as in Figure 1, the firing rates have been estimated by simulating the networks for a total time T_s = 60 s, after discarding a transient T_t = 40 s.

Let us now examine how the firing rates of active neurons will modify by increasing the value of the median in-degree K. The single neuron firing rates as a function of their total in-degrees $k_j^{(tot)}$ are reported in Figures 3C,D for K=1,024, 4,096 and 16,384. A common characteristic is that the bulk neurons, those with $k_j^{(tot)}\simeq 2K$, tend to approach the firing rate values $({\overline{R}}_0^{(e)},{\overline{R}}_0^{(i)})$ (magenta dashed lines) corresponding to the expected solutions for a balanced network in the limit N >> K → ∞ (van Vreeswijk, 1996). This is confirmed by the analysis of their coefficient of variations cv_j, whose values are of order one, as expected for fluctuation driven dynamics. On the other hand, the outlier neurons, i.e., those with $k_j^{(tot)}$ far from 2K, are all characterized by low values of the coefficient of variation cv_j indicating a mean driven dynamics. However, there is a striking difference between excitatory and inhibitory neurons. For the excitatory ones, we observe that the firing rates of the outliers with $k_j^{(tot)}>>2K$ decrease for increasing K, while for the inhibitory population the increase of K leads to the emergence of outliers at $k_j^{(tot)}<<2K$ with higher and higher firing rates (refer to the inset in Figure 3D). This difference can be explained by the different values measured for $I_{eff}^{(e)}$ and $I_{eff}^{(i)}$ in the network (refer to Figure 2C). The increase of K leads for the excitatory (inhibitory) population to the emergence of neurons with very large $k_j^{(ee)}>>K$ (very small $k_j^{(ii)}<<K$) whose dynamics should be supra-threshold. However, this is compensated in the excitatory case by the rapid drop of $I_{eff}^{(e)}$ toward zero or negative values, while for the inhibitory population $I_{eff}^{(i)}$ remains positive even at the largest K we have examined.

These outliers seem to have a negligible influence on the population dynamics, as suggested by the fact that the mean firing rates are reasonably well approximated by the balanced solutions ${\overline{R}}_0^{(e)}$ and ${\overline{R}}_0^{(i)}$ and as also confirmed by examining the PDFs of the firing rates for K = 16,384. As shown in Figures 3E,F, the excitatory (inhibitory) PDF can be well fitted by a log-normal distribution with a mean 8.8 Hz (17.5 Hz) and SD of 3.8 Hz (2.3 Hz). This is considered a clear indication that the network dynamics is fluctuation driven (Roxin et al., 2011) as confirmed by recent investigations in the hippocampus and the cortex (Wohrer et al., 2013; Buzsáki and Mizuseki, 2014; Mongillo et al., 2018), as well as in the spinal motor networks (Petersen and Berg, 2016). However, the relative widths of our distributions are narrower than those reported in Mongillo et al. (2018). This difference can find an explanation in the theoretical analysis reported in Roxin et al. (2011), where the authors have shown that quite counter intuitively a wider distribution of the synaptic heterogeneities can lead to a narrower distribution of the firing rates. Indeed, in this study, we consider Lorentzian distributed in-degrees, while in Mongillo et al. (2018) Erdös-Renyi networks have been analyzed. As a further aspect, we have estimated the number of inhibitory neurons firing faster than a certain threshold ν_th, this number does not depend on the median in-degree for sufficiently large K > 5,000, however, it grows proportionally to N. In the considered cases, the fraction of these neurons is ≃1% for ν_th = 50 Hz.

From this analysis, we can conclude that at any finite K and for finite observation times, we have at a macroscopic scale an essentially balanced regime sustained by the bulk of active neurons, whose dynamics are fluctuation-driven. Furthermore, we also have a large body of silent neurons and a small fraction of mean driven outliers. These should be considered as typical features of finite heterogeneous neural circuits as shown in various experiments (O'Connor et al., 2010; Landau et al., 2016). Moreover, in the present case, we report quite different behaviors for outliers whose macroscopic effective input currents are supra- or sub-threshold.

3.3. Collective Oscillations

We will now characterize the different types of COs observable by first following a route to coherent chaos for the E-I balanced network and successively we will examine how oscillations exhibiting an abnormal level of synchronization, somehow similar to those observable during an ictal state in the brain (Lehnertz et al., 2009), can emerge in our system. Furthermore, we will consider the phenomenon of quasi-periodicity and frequency locking occurring for fluctuation driven oscillations. As the last issue, the scaling of the frequencies and amplitudes of COs with the in-degree and as a function of the external DC current is reported.

3.3.1. A Period Doubling Route to Coherent Chaos

As a first case, we will follow the path in the parameter space denoted as a dashed black line in Figure 1A. In particular, in order to characterize the different dynamical regimes, we have estimated the LS {λ_i} associated the MF equations. As shown in Figure 4, this analysis has allowed us to identify a period doubling cascade toward a chaotic region, characterized by periodic and chaotic windows. In particular, we observe a focus region (II) for $0.0015<I_0^{(e)}<50.6105$, the focus loses stability via a super-critical HB at $I_0^e\simeq 0.0015$ giving rise to COs. One observes a period doubling cascade [regime (V)] taking place in the interval $I_0^{(e)}\in \left[0.00006177;0.00047297\right]$ followed by a regime of COs at lower values of $I_0^{(e)}$. The chaotic dynamics refer to the MF evolution, and it can be, therefore, definitely identified as collective chaos (Nakagawa and Kuramoto, 1993; Shibata and Kaneko, 1998; Olmi et al., 2011). A peculiar aspect of this period doubling cascade is that the chaotic dynamics remain always confined in four distinct regions without merging in a unique interval as it happens e.g., for the logistic map at the Ulam point (Ott, 2002). This is due to the fact that the population dynamics display period four oscillations characterized by four successive bursts, whose amplitudes (measured by $R_{max}^{(e)}$) varies chaotically but each one remains restricted in an interval not overlapping with the other ones.

FIGURE 4

Figure 4. Coherent chaos. (a,b) First (red) λ₁ and second λ₂ (blue) (LEs) for the MF vs. the DC current $I_0^{(e)}$ for the parameter cut corresponding to the dashed black line in Figure 1A. The dashed vertical lines in (a) indicate a super-critical Hopf bifurcation (HB) from a stable focus to periodic COs and the region of the period doubling (PD) cascade. The symbols denote three different types of MF solutions: namely, stable focus (green triangle); periodic oscillations (blue square) and chaotic oscillations (red circle). (c,d) Bifurcation diagrams for the same region obtained by reporting the maximal value of the instantaneous firing rate R^(e) measured from MF simulations. The parameters are the same as in Figure 1, other parameters set as ${\Delta }_0^{(ii)}=0.3$, ${\Delta }_0^{(ee)}=2.0$, K = 1,000.

Let us now examine the network dynamics for the three peculiar MF solutions indicated in Figure 4a corresponding to a stable focus (II) characterized by LE (λ₁ = λ₂ = −0.0299, λ₃ = λ₄ = −0.101) for $I_0^{(e)}=0.006$ (green triangle), to a stable oscillation (IV) with (λ₁ = 0.0, λ₂ = −0.0343, λ₃ = −0.0555, λ₄ = −0.1732) for $I_0^{(e)}=0.0009$ (blue square), and to collective chaos (v) with (λ₁ = 0.0033, λ₂ = 0.0, λ₃ = −0.0809, λ₄ = −0.1855) for $I_0^{(e)}=0.00021$ (red circle). As shown in Figure 5, for all these three cases, the network dynamics is always characterized by oscillations: namely, O_P for the regimes (IV) and (V) and fluctuation induced O_F for to the stable MF focus.

FIGURE 5

Figure 5. Different types of collective oscillations(COs). Row (a) refers to the chaotic state observable for $I_0^{(e)}=0.00021$ in the MF denoted by a red circle in Figure 4a; row (b) to the oscillatory state of the MF observable for $I_0^{(e)}=0.0009$ denoted by a blue square in Figure 4a; row (c) to the stable focus for the MF observable for $I_0^{(e)}=0.006$ denoted by a green triangle in Figure 4a. The first column displays the population firing rates vs. time obtained from the network dynamics, the second, the corresponding MF attractors in the planes identified by (R^(e), V^(e)) and (R⁽ⁱ⁾, V⁽ⁱ⁾), the third, the raster plots, and the fourth, the PDFs of the excitatory firing rates $r_j^{(e)}$. Red (blue) color refers to excitatory (inhibitory) populations, the solid vertical lines in column 4 to the mean firing rate and the blue solid line to a fit to a log-normal distribution. Parameters as in Figure 2, apart from ${\Delta }_0^{(ii)}=0.3$, ${\Delta }_0^{(ee)}=2.0$, K = 1,000. For the estimation of the firing rates we employed N^(e) = 40,000 and N⁽ⁱ⁾ = 10,000, while for the raster plots, N^(e) = 10,000 and N⁽ⁱ⁾ = 2,500. The total integration time has been of 120 s after discarding a transient of 80 s.

A typical feature of the O_P oscillations is that the excitatory neurons start to fire followed by the inhibitory ones, furthermore, the peak of activity of the excitatory population usually precedes that of the inhibitory neurons of a time interval Δt. Then the inhibitory burst silences the excitatory population for the time needed to recover toward the firing threshold. This recovering time sets the frequency ν_CO of the COs. In our set-up, the excitatory bursts are wider than the inhibitory ones due to the fact that ${\Delta }_0^{(ee)}>{\Delta }_0^{(ii)}$. All these features are quite evident from the population firing rates shown in Figures 5a1,b1 and the raster plots in panels Figures 5a3,b3. These are typical characteristics of a PING-like mechanism reported for the generation of γ oscillations in the cortex (Tiesinga and Sejnowski, 2009), despite the fact that the frequencies of the COs shown in panels (a) and (d) are of the order of few Hz. Fluctuation driven oscillations O_F emerging in the network are radically different, as shown in Figure 5c1, in this case, the excitatory and inhibitory populations deliver almost simultaneous bursts. Further differences among O_P and O_F oscillations can be identified at the level of single neuron activity. These can be appreciated by considering the PDFs of the excitatory firing rates $r_j^{(e)}$ reported in the fourth column of Figure 5. As shown in Figure 5c4 these firing rates are log-normally distributed for O_F oscillations, thus, confirming their fluctuation driven origin (Roxin et al., 2011; Petersen and Berg, 2016). On the other hand, for O_P oscillations, we observe with respect to a log-normal distribution an excess of high firing neurons and a lack of low firing ones (refer to Figures 5a4,b4). This seems to indicate the presence of a larger number of mean driven excitatory neurons. Indeed this is the case, for $I_0^{(e)}=$ 0.00021 and $I_0^{(e)}=$ 0.0009, the percentage of active excitatory neurons driven by average effective currents supra-threshold $i_{eff,j}^{(e)}$ is ≃ 1.7 − 1.2%, while for $I_0^{(e)}=$0.006, it drops to ≃0.6%. The percentage of active inhibitory neurons on average supra-threshold is quite limited in both cases being of the order of 0.25–0.13%. Another interesting feature distinguishing the two kinds of oscillations is the fact that for O_P, the excitatory supra-threshold neurons have a firing rate $r_j^{(e)}>{\nu }_{CO}$ and that the few neurons with firing rates locked to ν_CO are on average exactly balanced, i.e., they have $i_{eff,j}^{(e)}\simeq 0$. The situation is different for the O_F oscillations, where we observe a group of sub-threshold excitatory and inhibitory neurons firing locked with the population bursts. In both cases, most parts of neurons are definitely sub-threshold firing at frequencies smaller than ν_CO, as expected for an E-I balanced network displaying fast network oscillations associated with irregular neural discharges (Brunel and Wang, 2003).

In order to understand the different mechanisms at the basis of O_P and O_F oscillations, let us examine how the delay Δt between excitatory and inhibitory bursts, observed for O_P oscillations, modifies as a function of the membrane time constant of the inhibitory population ${\tau }_m^{(i)}$. An increase of ${\tau }_m^{(i)}$ of ≃5 ms has the effect of reducing the delay of almost a factor six from Δt ≃ 28 ms to Δt ≃ 5 ms, as shown in Figure 6A. The increase of ${\tau }_m^{(i)}$ leads to an enhanced inhibitory action since the integration of the inhibitory membrane potentials occurs on longer time scales, and this promotes a higher activity of the inhibitory population. Indeed, this is confirmed from the drop of the effective input currents from an almost balanced situation where the average $I_{eff}^{(e)}$ and $I_{eff}^{(i)}$ are almost zero to a situation where they are definitely negative (refer to Figure 6B). Thus, for increasing ${\tau }_m^{(i)}$, the percentage of neurons below threshold also increases and as a consequence the dynamics become more and more noise driven, as testified by the increase of the current fluctuations $\Delta I_{eff}^{(e,i)}$ as shown in Figure 6C. In summary, the delay is due to the fact that despite the effective inhibitory and excitatory currents are essentially equal, as shown in Figure 6B, the wider distribution of the excitatory in-degrees promotes the presence of excitatory neurons supra-threshold that are the ones igniting the excitatory burst before the inhibitory one. The delay Δt decreases whenever the number of these supra-threshold neurons decreases, and it will vanish when the dynamics will become essentially fluctuation driven as in the case of O_F oscillations.

FIGURE 6

Figure 6. Pyramidal-interneuron gamma (PING)-like O_P COs. (A) Firing delays Δt between the excitatory population peak and the inhibitory one vs. ${\tau }_m^{(i)}$. Effective mean input currents (Equation 13) (B) and current fluctuations (Equation 20) (C) vs. ${\tau }_m^{(i)}$, the excitatory (inhibitory) population are denoted by red (blue) circles. All the data reported in this study refer to MF simulations. The parameters are $I_0^{(e)}=0.0009$, ${\Delta }_0^{(ii)}=0.3$, ${\Delta }_0^{(ee)}=2.0$, K = 1,000, and ${\tau }_m^{(e)}=20$ ms.

3.3.2. From Fluctuation Driven to Abnormally Synchronized Oscillations

As the second range of parameters, we consider the cut in the parameter plane shown in Figure 1A as a purple dashed line. For these parameters, we report in Figures 7a,b the average in time of the excitatory and inhibitory population rate as a function of the excitatory DC current $I_0^{(e)}$. In particular, we compare network simulations (red and blue circles) with the MF results (red and blue lines). These predict a stable focus (solid lines) up to $I_0^{(e)}=74.1709$, where a sub-critical HB destabilizes such solution giving rise to an unstable focus (dashed lines). In panel (a), and (b), we have also reported as green dot-dashed lines the extrema of R^(e) and R⁽ⁱ⁾ corresponding to the unstable oscillations emerging at the HB. For currents below the HB, we observe a good agreement among the average network activity and the MF results.

FIGURE 7

Figure 7. From fluctuation driven to abnormally synchronized oscillations. Firing rates R^(e) (a) and R⁽ⁱ⁾ (b) as a function of $I_0^{(e)}$ for E-I network (circles) and neural mass model (lines) for the parameter cut corresponding to the dashed purple line in Figure 1A. For the neural mass model: solid (dashed) line shows stable (unstable) focus solution ${\overline{R}}^{(e)}$ and ${\overline{R}}^{(i)}$; green dot-dashed lines refer to the extrema of R^(e)(R⁽ⁱ⁾) for the unstable LC present in region (II). The unstable LC emerges at the sub-critical HB for $I_0^{(e)}=74.1709$ separating region (II) from (I), where the focus becomes unstable. Raster plots and PDFs of the excitatory firing rates $r_j^{(e)}$ are reported for specific cases: namely, $I_0^{(e)}=0.128$ (c1,c2), $I_0^{(e)}=1.024$ (d1,d2), and $I_0^{(e)}=100$ (e1,e2). The solid vertical lines in (c2,d2,e2) refer to the mean firing rate. Parameters as in Figure 1, other parameters are set as ${\Delta }_0^{(ii)}=0.3$, ${\Delta }_0^{(ee)}=1.58$, K = 1,000 N^(e) = 10,000, and N⁽ⁱ⁾ = 2,500.

In particular, below the HB, while the MF predicts only the existence of a stable focus, the network dynamics reveals quite interesting features. As shown in Figure 7d1, the system dynamics is indeed asynchronous for intermediate current values, in this study, $I_0^{(e)}=1.024$, however, at lower currents, we observe fluctuation driven oscillations O_F as evident from the raster plot displayed in Figure 7c1 for $I_0^{(e)}=0.128$. As shown in Figures 7c2,d2, both these regimes are characterized by log-normal distributions of the firing rates, thus, indicating that the dynamics are fluctuation driven.

As reported in Montbrió et al. (2015), when the network dynamics become strongly synchronous (as expected for very high excitatory DC external current), the MF formulation fails since the population rates predicted within the MF formulation diverge. However, as shown in Figures 7e1,e2, due to finite size effects, we observe in the network a strong synchronous COs of type O_S corresponding to the MF region (I) where the MF model predicts no stable solution. These abnormally synchronized oscillations are also characterized by a quite fast frequency of oscillation ν_CO ≃ 800 − 1,000 Hz. Furthermore, similarly to the O_P oscillations, they emerge due to a PING-like mechanism. This is evident from the raster plot in Figure 7e1, where excitatory neurons fire almost synchronously followed, after an extremely short delay, by the inhibitory ones whose activity silence all the network until the next excitatory burst. Quite astonishingly, the mean population rates measured in the network are reasonably well captured by the MF solutions associated with the unstable focus even beyond the HB, despite the network is now displaying COs (as shown in Figures 7a,b).

The emergence of COs in the network can be characterized in terms of the coherence indicator ρ (Equation 4) for the whole population of neurons. This indicator is reported in Figure 8A as a function of $I_0^{(e)}$ for the same parameters previously discussed in Figure 7 and for two different values of the median in-degree : K = 100 (red circles) and K = 4,000 (blue circles). For both values of K, we observe an almost discontinuous transition in the value of the coherence indicator at the sub-critical HB from $\rho \simeq 1/\sqrt{N}$, expected for an asynchronous dynamics, to values ρ≃1 corresponding to full synchronization. This discontinuous transition leads to the emergence of abnormally synchronized oscillations O_S in the network. Moreover, at sufficiently high in-degrees, we observe the emergence of a new coherent state for low DC currents $I_0^{(e)}<$ 1.024 characterized by a finite value of the coherence indicator, namely, ρ≃0.3. The origin of these oscillations can be better understood by examining the coefficient of variation CV averaged over the whole population, this is reported in Figure 8C for the same interval of excitatory DC current and the same in-degrees as in Figure 8A. It is evident that the CV assumes finite values only for small input currents, namely $I_0^{(e)}<1.024$, indicating the presence of not negligible fluctuations in the network dynamics. Furthermore, by increasing K, these fluctuations, as measured by the CV, increase as expected for a balanced network. This analysis suggests that these oscillations cannot exist in absence of fluctuations in the network, and therefore, they are of the O_F type. Furthermore, the network should be sufficiently connected in order to sustain these COs, as one can understand from Figures 8B,D, where ρ and CV are reported as a function of K for three different values of $I_0^{(e)}$. Indeed, for these parameter values, no O_F oscillation is observable for K < 400, even in presence of finite values of the CV.

FIGURE 8

Figure 8. From fluctuation driven to abnormally synchronized oscillations. Coherence indicator ρ (Equation 4) for the whole network of excitatory and inhibitory neurons vs. the excitatory DC current $I_0^{(e)}$ (A) and the median in-degree K (C). Coefficient of variation CV for the whole network vs. $I_0^{(e)}$ (B) and K (D). In (A,C), the symbols refer to different values of the median in-degree:namely, K = 100 (red circles) and K = 4,000 (blue circles). In (B,D), the symbols refer to different excitatory DC currents: namely, $I_0^{(e)}=0.01$ (green circles), $I_0^{(e)}=0.1$ (purple circles), and $I_0^{(e)}=1.0$ (orange circles). Parameters as in Figure 1, other parameters ${\Delta }_0^{(ii)}=0.3$, ${\Delta }_0^{(ee)}=1.58$, N^(e) = 40,000, and N⁽ⁱ⁾ = 10,000.

As previously discussed in di Volo and Torcini (2018), the balance between excitation and inhibition generates endogenous fluctuations that modify the collective dynamics with respect to that predicted by the MF model, where the heterogeneity of the input currents, due to distributed in-degrees, is taken in account only as a quenched form of disorder and not as a dynamical source of the noise. However, also from this simplified MF formulation, one can obtain relevant information on the O_F oscillations, indeed as we will see in the next sub-section, the relaxation frequencies toward the stable MF focus represent a good estimation of the oscillation frequencies measured in the network. This suggests that the fluctuations present at the network level can sustain COs by continuously exciting the focus observed in the effective MF model with quenched disorder.

3.3.3. Fluctuation Driven Oscillations: From Quasi-Periodicity to Frequency Locking

As announced, this sub-section will be devoted to the characterization of the fluctuation driven oscillations O_F emerging in the region (II) reported in Figure 1. As the MF is now characterized by a stable focus with two couples of complex conjugate eigenvalues, there are two frequencies that can be excited by the irregular firing of neurons. Accordingly, as reported in di Volo and Torcini (2018), we expect the collective dynamics to be characterized by quasi-periodic dynamics with two (incommensurable) frequencies. These frequencies can be estimated by computing the power spectrum S(ν) of global quantities, e.g., the mean membrane potential V(t). In the case of periodic dynamics, S(ν) is characterized by one main peak in correspondence of the CO frequency and minor peaks at its harmonics, while in the quasi-periodic case, the power spectrum shows peaks located at the two fundamental frequencies and all their linear combinations. Indeed, as shown in Figure 9A, the power spectrum exhibit several peaks over a continuous profile and the peak frequencies can be obtained as a linear combination of two fundamental frequencies (ν₁, ν₂). As already mentioned, the noisy background is due to the fluctuations present in the balanced network. It is evident from Figure 9B, that these two fundamental frequencies are well reproduced by the two relaxation frequencies ${\nu }_1^R$ and ${\nu }_2^R$ toward the MF focus, in particular for $I_0^{(e)}\ge 0.256$. At smaller currents, while the first frequency is well reproduced by ${\nu }_1^R$, the second one is under-estimated by ${\nu }_2^R$. This is due to the phenomenon of frequency locking among the two collective rhythms present in the system: when the two frequencies become commensurable, we observe a common periodic CO. The locking order can be estimated by plotting the ratio between the two frequencies, indeed for low currents and K = 8, 192, the ratio is almost constant and equal to four denoting a 1:4 frequency locking (see Figure 9C). Furthermore, by fixing $I_0^{(e)}=0.128$ and by varying K the ratio ν₁/ν₂ can display different locked states, passing from locking of type 1:2 at low K, to 1:4 at larger values, as shown in the inset of Figure 9C.

FIGURE 9

Figure 9. From quasi-periodicity to frequency locking. (A) Power spectra S(ν) of the mean membrane potential obtained from network simulations. (B) The two fundamental frequencies ν₁(ν₂) vs. $I_0^{(e)}$. (C) Frequency ratio ν₁/ν₂ vs. $I_0^{(e)}$, in the inset ν₁/ν₂ is shown vs. K. (D) Coherence parameter ρ vs. $I_0^{(e)}$, in the inset the corresponding CV is reported. In (B,C), the symbols (solid lines) refer to ν₁ and ν₂ as obtained from the peaks of the power spectra S(ν) for V(t) obtained from the network dynamics (to the two relaxation frequencies ${\nu }_1^R$ and ${\nu }_2^R$ associated to the stable focus solution for the MF). Parameters as in Figure 1, other parameters are set as ${\Delta }_0^{(ii)}=0.3$, ${\Delta }_0^{(ee)}=1.58,N^{(e)}=$80,000, N⁽ⁱ⁾ = 20,000, K = 8,192, and $I_0^{(e)}=$ 0.128 in the inset of (C).

As evident from Figures 9B,C, the locking phenomenon arises only in the network simulations and is not captured by the MF model. Furthermore, frequency locking occurs at low currents $I_0^{(e)}<0.1$ where the dynamics of the neurons are driven by the intrinsic current fluctuations present in the network but not in the MF. Indeed for low DC currents the level of synchronization within the populations measured by ρ decreases with $I_0^{(e)}$, while the CV increases (as shown in Figure 9D). These features suggest that this phenomenon is somehow similar to what was reported in Meng and Riecke (2018) for two coupled inhibitory neural populations subject to external uncorrelated noise. Meng and Riecke (2018) observed an increase of the locking region among collective rhythms by increasing the amplitude of the additive noise terms, this joined to a counter-intuitive decrease of the level of synchronization among the neurons within each population. However, in Meng and Riecke (2018), the neurons are subject to independent external noise sources, while in our case, the sources of fluctuations are intrinsic to the system and induced by the structural heterogeneity. Due to the network sparseness, the current fluctuations experienced by each neuron can be assumed to be indeed uncorrelated (Brunel and Hakim, 1999). Therefore, we are facing a new phenomenon that we can identify as a frequency locking of collective rhythms promoted by self-induced uncorrelated fluctuations. Indeed, the locking disappears for increasing external DC currents $I_0^{(e)}>0.1$, when the coherence parameter ρ displays an abrupt jump toward higher values and the CV≃0, thus, indicating that in this regime, the neuron dynamics becomes essentially mean driven and highly synchronized.

3.3.4. Features of COs for Large In-degrees and DC Currents

The dynamics of balanced networks are usually characterized in the limit N >> K >> 1 by the emergence of a self-sustained asynchronous regime. However, LC solutions have been already reported for balanced networks in the seminal article van Vreeswijk and Sompolinsky (1996). These solutions can be either unbalanced or balanced, however, in this latter case, they were characterized by oscillations of vanishing small amplitude. van Vreeswijk and Sompolinsky (1996) have shown that balanced COs are not observable in their model in the limit N >> K → ∞ but only for finite K. Therefore, it is important to address in our case if COs can still be observable in the limit N >> K >> 1. Thus, in the following, we will investigate the dependence of COs features on the median in-degree K and the external DC currents.

Let us first consider fluctuation driven O_F oscillations, in this case, we have an analytical prediction (Equation 19) for the scaling of the fundamental frequencies ${\nu }_k^R$ associated with the relaxation toward the macroscopic focus, which should grow proportionally to $\sqrt{I^{(e)}}$. As shown in Figures 10A,C, indeed this scaling is observable for sufficiently large K and $I_0^{(e)}$. It is also evident the extreme good agreement between results obtained from the network simulations and the theoretical predictions (Equation 19), at least for the values of K reachable with our simulations. Furthermore, the frequencies of COs cover an extremely large range of values from few Hz to KHz, and this range of frequencies can be spanned by varying either K or the external DC current $I_0^{(e)}$ as shown in Figures 10A,C.

FIGURE 10

Figure 10. Frequencies and amplitudes of O_F oscillations. The two fundamental frequencies ν₁ and ν₂ vs. $I_0^{(e)}$ (A) and K (C) and the average firing rates vs. $I_0^{(e)}$ (B) and K (D) for the excitatory (red) and inhibitory (blue) populations. In the inset in (C), the effective mean input currents $I_{eff}^{(e)}$ ($I_{eff}^{(i)}$) of the excitatory (inhibitory) population are shown vs. K. The dashed line in (A,C) corresponds to a power law-scaling $\propto {I_0^{(e)}}^{1/2}$ (∝K^1/4) for the frequencies of the COs. The solid red (blue) line in (B,D) denotes the asymptotic MF result ${\overline{R}}^{(e)}$ (${\overline{R}}^{(i)}$). Network (MF) simulations are denoted as stars (circles). The MF data refer to the stable focus, in particular, in (A,C), these are the two relaxation frequencies ${\nu }_1^R$ and ${\nu }_2^R$. Parameters as in Figure 1, other parameters: (A,B) K = 1, 000, ${\Delta }_0^{(ee)}=1.58$, ${\Delta }_0^{(ii)}=0.3$; (C,D) $I_0^{(e)}=0.001$, ${\Delta }_0^{(ee)}=1.3$, ${\Delta }_0^{(ii)}=0.3$; for the network simulations, we employed N^(e) = 80,000 and N⁽ⁱ⁾ = 20,000.

To better characterize these regimes, we have also evaluated the average firing rates R^(e) and R⁽ⁱ⁾. These quantities are displayed for O_F oscillations in Figures 10B,D as a function of $I_0^{(e)}$ and K, respectively. From the network simulations (stars), we observe that R^(e) and R⁽ⁱ⁾ grow with $I_0^{(e)}$, and they are astonishingly quite well reproduced by the MF data (circles) for sufficiently large DC currents, despite the MF results refer to a stable focus and not to COs. Instead, at smaller currents (namely, $I_0^{(e)}=0.001$), the network data overestimates the MF results and the excitatory and inhibitory firing rates for K >> 1 seem to converge to a common constant value larger than those corresponding to the asynchronous regimes. For sufficiently large K, due to the prevalence of inhibition over excitation in the present model, we expect that the system will be sub-threshold, since the average excitatory and inhibitory firing rates are essentially coincident. Indeed this is confirmed by the analysis of the mean effective input currents $I_{eff}^{(e)}$ and $I_{eff}^{(i)}$ shown in the inset of Figure 10C. While for the MF focus, the dynamics appear as almost exactly balanced for all the considered median in-degree K since $I_{eff}^{(e)}\simeq I_{eff}^{(i)}\simeq 0$, for the network dynamics $I_{eff}^{(e)}$ and $I_{eff}^{(i)}$ are definitely negative for K > 1, 000. This does not prevent the emergence of COs driven by fluctuations at large K, as indeed observed.

These results seem to indicate that for N >> K → ∞, the network will not converge in this case toward a balanced regime characterized by constant effective input currents. On the contrary from our analysis, it emerges that the system will become more and more sub-threshold for increasing K > 1,000. However, the system always exhibits fluctuation driven dynamics, since we measured CV≃ 0.6-0.8 at least in the range K≃100−10000 accessible to network simulations.

Let us now examine the O_P oscillations. As shown in Figures 11A,C, the frequencies ν_CO as estimated from the MF model (open circles) reveal an almost perfect increase proportional to $\sqrt{I^{(e)}}$ analogous to the one reported for O_F oscillations. The data obtained from network simulations (stars) converge toward the MF results for sufficiently large K and $I_0^{(e)}$.

FIGURE 11

Figure 11. Frequencies and amplitudes of O_P oscillations. COs' frequency ν_CO vs. $I_0^{(e)}$ (A) and K (C) and mean firing rates vs. $I_0^{(e)}$ (B) and K (D) for the excitatory (red) and inhibitory (blue) populations. The dashed line in (A,C) corresponds to a power law-scaling $\propto {I_0^{(e)}}^{1/2}$ (∝K^1/4) for the frequencies. In the inset in (c), the effective mean input currents $I_{eff}^{(e)}$ ($I_{eff}^{(i)}$) of the excitatory (inhibitory) population are shown vs. K. The solid red (blue) line in (B,D) denotes the asymptotic MF result ${\overline{R}}^{(e)}$ (${\overline{R}}^{(i)}$). The data obtained from network (MF) simulations are denoted as stars (circles). The data reported in (A–D) refer to the open circles in Figures 1A,B, respectively. For network simulations, we employed N^(e) = 80,000 and N⁽ⁱ⁾ = 20,000.

The mean firing rates R^(e) and R⁽ⁱ⁾ grow with $I_0^{(e)}$ for fixed K and appear to converge toward a constant value for sufficiently large K for fixed $I_0^{(e)}$, refer to Figures 11B,D. Moreover, the network simulations (stars) approach the MF results (open circles) at large DC currents and median in-degrees. However, while in the MF the asymptotic values of R^(e) and R⁽ⁱ⁾ remain distinct even at large K, these seem to become identical in the network simulations. This reflects in the fact that while the MF is perfectly balanced in the whole range of examined in-degrees, since $I_{eff}^{(e)}\simeq I_{eff}^{(i)}\simeq 0$, the network simulations reveal almost balanced effective input currents up to K≃ 1,000 and above such median in-degree a prevalence of the inhibitory drive (refer to inset of Figure 11C).

For both kinds of COs, we observe that while ν_CO diverges with K, the mean firing rates approach a constant value, thus, suggesting that the percentage of neurons participating in each population burst should vanish in the limit K → ∞. This result indicates that COs will finally disappear, however, more refined analyses are needed to derive the asymptotic behavior of the system in the large K limit, ref to di Volo et al. (2021) for a detailed discussion of this aspect for purely inhibitory networks.

4. Discussion

We have extensively characterized the macroscopic regimes emerging in a sparse balanced E-I network made of spiking QIF neurons with Lorentzian distributed in-degrees. The considered neuronal model joined to the peculiar choice of the distribution has allowed us to derive an exact low dimensional neural mass model describing the MF dynamics of the network in terms of the mean membrane potentials and of the population rates of the two populations (Montbrió et al., 2015; di Volo and Torcini, 2018). The low-dimensionality of the MF equations enabled us to study analytically the stationary solutions and their stability as well as to obtain the bifurcation diagrams associated with the model and to identify the possible macroscopic states.

4.1. Asynchronous Regimes

The stationary solutions of the MF correspond to the asynchronous regime, which is the regime usually analyzed in the context of balanced dynamics (van Vreeswijk and Sompolinsky, 1996; Renart et al., 2010; Litwin-Kumar and Doiron, 2012). In the present case, we have analytically obtained the stationary solutions for the mean membrane potentials and average firing rates for Lorentzian distributed in-degrees for any finite value of the median K and an HWHM scaling as ${\Delta }_k^{(\alpha \alpha )}={\Delta }_0^{(\alpha \alpha )}{(K)}^{\eta }$ with η = 1/2. The MF estimations for the population firing rates are pretty well reproduced by the network simulations in the examined range of in-degrees K. Furthermore, from the analytic expression of the stationary firing rates (Equation 8), it is evident that for K > > 1, the asymptotic rates would not depend on the structural heterogeneity and correspond to those usually found for balanced homogeneous or Erdös-Renyi networks (van Vreeswijk and Sompolinsky, 1996; Monteforte and Wolf, 2010). This is due to the fact that the ratio ${\left({\Delta }_k^{(\alpha \alpha )}\right)}^2/K$ remains constant for K → ∞. The final scenario will depend on the scaling exponent η, in particular, by assuming η = 3/4, the asymptotic firing rates ${\overline{R}}_0^{(\alpha )}$ will explicitly depend on the parameters ${\Delta }_0^{(\alpha \alpha )}$ controlling the structural heterogeneity. Whenever η > 3/4, the balanced state breaks down, and we face a situation similar to those investigated in Landau et al. (2016) and Pyle and Rosenbaum (2016)¹.

However, despite the system approaching a balanced state, as testified by the fact that the effective input currents converge to finite values $I_a^{(\alpha )}$, and the current fluctuations stay finite for K → ∞, the balanced regime is not necessarily a sub-threshold one. Indeed, we have observed that we can have either sub-threshold or supra-threshold situations depending on the model parameters in agreement with the results previously reported in Lerchner et al. (2006). Moreover, the excitatory and inhibitory populations can achieve balanced regimes characterized by different asymptotic dynamics, where $I_a^{(i)}$ and $I_a^{(e)}$ have opposite signs.

While at a macroscopic level, the population activity for N >> K >> 1 approach is essentially that of a homogeneous balanced system, as shown in Figures 2A,B, the structural heterogeneity has a large influence on the single neuron dynamics, at least at finite K and finite investigation times. In particular, in analogy with experiments (Gentet et al., 2010; Mongillo et al., 2018), we considered a situation where the inhibitory drive prevails on the excitatory one. In this condition, microscopically the neural populations split into three groups: silent neurons, definitely sub-threshold; bulk neurons, which are fluctuation driven; and mean driven outlier neurons. In particular, excitatory (inhibitory) neurons with low (high) intra-population in-degrees are silenced due to the prevalence of synaptic inhibition. The silent neurons represent 6-10% of the whole population in agreement with experimental results for the mice cortex (O'Connor et al., 2010). Bulk neurons have in-degrees in the proximity of the median, and their firing rates approach the MF solution ${\overline{R}}_0^{(\alpha )}$ for increasing K. Outlier neurons represent a minority group almost disconnected from their own population, whose asymptotic behavior for K >> 1 is controlled by the sign of the effective mean input current.

4.2. Coherent Dynamics

The emergence of COs is observable in this balanced network whenever the level of heterogeneity in the inhibitory population is not too large, thus, suggesting that the coherence among inhibitory neurons is fundamental to support collective rhythms (Whittington et al., 2000). Indeed we observed two main mechanisms leading to COs: one that can be identified as PING-like and another one as fluctuation driven. The PING-like mechanism is present whenever the excitatory neurons are able to deliver an almost synchronous excitatory volley that in turn elicits a delayed inhibitory one. The period of the COs is determined by the recovery time of the excitatory neurons from the stimulus received from the inhibitory population. This mechanism is characterized by a delay between the firing of the pyramidal cells and the interneuronal burst as reported also in many experiments (Buzsáki and Wang, 2012). We have shown that this delay tends to vanish when the inhibitory action increases leading the system from a balanced situation to a definitely sub-threshold condition where the neural activity is completely controlled by fluctuations. In this latter case, the excitatory and inhibitory neurons fire almost simultaneously driven by the current fluctuations. These transform the relaxation dynamics toward a stable focus, observable in the MF, to sustained COs via a mechanism previously reported for inhibitory networks (di Volo and Torcini, 2018; Bi et al., 2020).

The PING-like COs undergo period doubling cascades by varying K and/or $I_0^{(e)}$ finally leading to collective chaos (Nakagawa and Kuramoto, 1993; Shibata and Kaneko, 1998). The nature of this chaotic behavior is definitely macroscopic since it is captured by the neural mass model obtained within the MF formulation, as shown by analyzing the corresponding LS. This kind of chaos implies irregular temporal fluctuations joined to coherence at the spatial level over a large part of the network resembling coherent fluctuations observed across spatial scales in the neocortex (Smith and Kohn, 2008; Volgushev et al., 2011; Okun et al., 2012; Achermann et al., 2016). Collective (or coherent) chaos has been previously shown to be a ubiquitous feature for balanced random spiking neural networks massively coupled, where K is proportional to N (Politi et al., 2018; Ullner et al., 2018). In this study, we have generalized such results to balanced random networks with sparse connectivity, where K is independent by N. Recently, it has been claimed that the presence of structured feed-forward connectivity in a random network is needed to observe coherent chaos (Landau and Sompolinsky, 2018). However, as evident from our results and those reported in Ullner et al. (2018) and Politi et al. (2018), coherent chaos can naturally emerge in a recurrent neural network in absence of any structured connectivity introduced ad hoc to promote collective behaviors. Furthermore, we have shown that collective chaos can emerge in random balanced networks with instantaneous synapses and the absence of any delay, refer to Ullner et al. (2018).

Fluctuation driven COs are usually observable in our system as quasi-periodic collective motions characterized by two incommensurate frequencies. However, whenever the current fluctuations become sufficiently strong, the two frequencies can lock and give rise to a collective periodic motion. Furthermore, the locking region is characterized by a low level of synchrony in the network. These results resemble those reported in Meng and Riecke (2018) for two interconnected inhibitory neural networks subject to external uncorrelated noise. In particular, the authors have shown that uncorrelated noise sources enhance synchronization and frequency locking among the COs displayed by the two networks, despite the reduced synchrony among neurons within each network. At variance with Meng and Riecke (2018), in our case, the noise sources are intrinsic to the neural dynamics, but they can be as well considered as uncorrelated due to the sparseness in the connections (Brunel and Hakim, 1999; Brunel, 2000). Therefore, we are reporting a new example of frequency locking among collective rhythms promoted by self-induced uncorrelated fluctuations.

According to analytical arguments, the frequencies of the COs grow proportionally to the square root of the excitatory DC current. This on one side allows simply by varying the parameters $I_0^{(e)}$ or K, to cover with our model a broad range of COs' frequencies analogous to those found experimentally in the cortex (Chen et al., 2017). On another side, it implies that the frequencies of COs diverge as K^1/4, while the average firing rates seem to converge to a common value for sufficiently large K. These results seem to indicate that for large K, the network will become more and more unbalanced, with a prevalence of inhibition, while the amplitude of COs will tend to vanish. However, this analyses is not conclusive and more detailed analysis are required to capture the asymptotic behavior of the system in the limit N >> K >> 1.

4.3. Future Developments

The examined neural mass model has been derived by taking into account the random fluctuations due to the sparseness in the network connectivity only as a quenched disorder affecting the distribution of the effective synaptic couplings (Montbrió et al., 2015; di Volo and Torcini, 2018). The current fluctuations can be correctly incorporated in an MF formulation by developing a Fokker-Planck formalism for the problem, however, this will give rise to high (infinite) dimensional MF models (Brunel and Hakim, 1999; Brunel, 2000). We are currently developing reduction formalisms for the Fokker-Planck equation to obtain low dimensional neural mass models which will include the intrinsic current fluctuations (di Volo et al., 2021; Goldobin et al., 2021).

Relevant topics to investigate in the future to assess the generality of the reported results are their dependence on the chosen spiking neuron model and network architecture. In particular, for random networks, it is important to understand the role played by the distribution of the in-degrees, this is also in view of the recent findings reported in Klinshov et al. (2021).

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

HB and MV performed the simulations and data analysis. MV and AT were responsible for the state-of-the-art review and the article write-up. All the authors conceived and planned the research.

Funding

AT received financial support by the Excellence Initiative I-Site Paris Seine (grant no. ANR-16-IDEX-008) (together with HB), by the Labex MME-DII (grant no. ANR-11-LBX-0023-01), by the ANR Project ERMUNDY (grant no. ANR-18-CE37-0014) (together with MV), and all part of the French programme Investissements d'Avenir.

Conflict of Interest

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.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors acknowledge extremely useful discussions with D. G. Goldobin, G. Mongillo, E. Montbrió, S. Olmi, and A. Politi.

Footnotes

1. ^In such cases, balance has been recovered either by rewiring the post-synaptic connections (Pyle and Rosenbaum, 2016) or by introducing some sort of homeostatic plasticity or of spike-frequency adaptation (Landau et al., 2016).

References

Achermann, P., Rusterholz, T., Dürr, R., König, T., and Tarokh, L. (2016). Global field synchronization reveals rapid eye movement sleep as most synchronized brain state in the human eeg. R. Soc. Open Sci. 3:160201. doi: 10.1098/rsos.160201

PubMed Abstract | CrossRef Full Text | Google Scholar

Angulo-Garcia, D., Luccioli, S., Olmi, S., and Torcini, A. (2017). Death and rebirth of neural activity in sparse inhibitory networks. New J. Phys. 19:053011. doi: 10.1088/1367-2630/aa69ff

CrossRef Full Text | Google Scholar

Atallah, B. V., and Scanziani, M. (2009). Instantaneous modulation of gamma oscillation frequency by balancing excitation with inhibition. Neuron 62, 566–577. doi: 10.1016/j.neuron.2009.04.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Barral, J., and Reyes, A. D. (2016). Synaptic scaling rule preserves excitatory-inhibitory balance and salient neuronal network dynamics. Nat. Neurosci. 19:1690. doi: 10.1038/nn.4415

PubMed Abstract | CrossRef Full Text | Google Scholar

Benettin, G., Galgani, L., Giorgilli, A., and Strelcyn, J.-M. (1980). Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. part 1: Theory. Meccanica 15, 9–20. doi: 10.1007/BF02128236

CrossRef Full Text | Google Scholar

Berg, R. W., Alaburda, A., and Hounsgaard, J. (2007). Balanced inhibition and excitation drive spike activity in spinal half-centers. Science 315, 390–393. doi: 10.1126/science.1134960

PubMed Abstract | CrossRef Full Text | Google Scholar

Berg, R. W., Willumsen, A., and Lindén, H. (2019). When networks walk a fine line: balance of excitation and inhibition in spinal motor circuits. Curr. Opin. Physiol. 8, 76–83. doi: 10.1016/j.cophys.2019.01.006

CrossRef Full Text | Google Scholar

Bi, H., Segneri, M., di Volo, M., and Torcini, A. (2020). Coexistence of fast and slow gamma oscillations in one population of inhibitory spiking neurons. Phys. Rev. Res. 2:013042. doi: 10.1103/PhysRevResearch.2.013042

CrossRef Full Text | Google Scholar

Brunel, N. (2000). Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J. Comput. Neurosci. 8, 183–208. doi: 10.1023/A:1008925309027

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunel, N., and Hakim, V. (1999). Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. Neural Comput. 11, 1621–1671. doi: 10.1162/089976699300016179

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunel, N., and Wang, X.-J. (2003). What determines the frequency of fast network oscillations with irregular neural discharges? i. synaptic dynamics and excitation-inhibition balance. J. Neurophysiol. 90, 415–430. doi: 10.1152/jn.01095.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruno, R. M., and Sakmann, B. (2006). Cortex is driven by weak but synchronously active thalamocortical synapses. Science 312, 1622–1627. doi: 10.1126/science.1124593

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., and Mizuseki, K. (2014). The log-dynamic brain: how skewed distributions affect network operations. Nat. Rev. Neurosci. 15, 264–278. doi: 10.1038/nrn3687

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., and Wang, X.-J. (2012). Mechanisms of gamma oscillations. Annu. Rev. Neurosci. 35:203–225. doi: 10.1146/annurev-neuro-062111-150444

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, G., Zhang, Y., Li, X., Zhao, X., Ye, Q., Lin, Y., et al. (2017). Distinct inhibitory circuits orchestrate cortical beta and gamma band oscillations. Neuron 96, 1403–1418. doi: 10.1016/j.neuron.2017.11.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Dehghani, N., Peyrache, A., Telenczuk, B., Le Van Quyen, M., Halgren, E., Cash, S. S., et al. (2016). Dynamic balance of excitation and inhibition in human and monkey neocortex. Sci. Rep. 6:23176. doi: 10.1038/srep23176

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A., and Paré, D. (1999). Impact of network activity on the integrative properties of neocortical pyramidal neurons in vivo. J. Neurophysiol. 81, 1531–1547. doi: 10.1152/jn.1999.81.4.1531

PubMed Abstract | CrossRef Full Text | Google Scholar

di Volo, M., Segneri, M., Goldobin, D. S., Politi, A., and Torcini, A. (2021). Coherent oscillations in balanced neural networks driven by endogenous fluctuations. bioRxiv [Preprint]. doi: 10.1101/2021.10.18.464823

CrossRef Full Text | Google Scholar

di Volo, M., and Torcini, A. (2018). Transition from asynchronous to oscillatory dynamics in balanced spiking networks with instantaneous synapses. Phys. Rev. Lett. 121:128301. doi: 10.1103/PhysRevLett.121.128301

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, B. (2007). XPPAUT. Scholarpedia 2:1399. revision #136177. doi: 10.4249/scholarpedia.1399

CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Kopell, N. (1986). Parabolic bursting in an excitable system coupled with a slow oscillation. SIAM J. Appl. Math. 46, 233–253. doi: 10.1137/0146017

CrossRef Full Text | Google Scholar

Gentet, L. J., Avermann, M., Matyas, F., Staiger, J. F., and Petersen, C. C. (2010). Membrane potential dynamics of gabaergic neurons in the barrel cortex of behaving mice. Neuron 65, 422–435. doi: 10.1016/j.neuron.2010.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldobin, D. S., Di Volo, M., and Torcini, A. (2021). A reduction methodology for fluctuation driven population dynamics. Phys. Rev. Lett. 127:038301. doi: 10.1103/PhysRevLett.127.038301

PubMed Abstract | CrossRef Full Text | Google Scholar

Golomb, D. (2007). Neuronal synchrony measures. Scholarpedia 2, 1347. doi: 10.4249/scholarpedia.1347

CrossRef Full Text | Google Scholar

Haider, B., Duque, A., Hasenstaub, A. R., and McCormick, D. A. (2006). Neocortical network activity in vivo is generated through a dynamic balance of excitation and inhibition. J. Neurosci. 26, 4535–4545. doi: 10.1523/JNEUROSCI.5297-05.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Isaacson, J. S., and Scanziani, M. (2011). How inhibition shapes cortical activity. Neuron 72, 231–243. doi: 10.1016/j.neuron.2011.09.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Kadmon, J., and Sompolinsky, H. (2015). Transition to chaos in random neuronal networks. Phys. Rev. X 5:041030. doi: 10.1103/PhysRevX.5.041030

CrossRef Full Text | Google Scholar

Klinshov, V., Kirillov, S., and Nekorkin, V. (2021). Reduction of the collective dynamics of neural populations with realistic forms of heterogeneity. Phys. Rev. E 103:L040302. doi: 10.1103/PhysRevE.103.L040302

PubMed Abstract | CrossRef Full Text | Google Scholar

Landau, I. D., Egger, R., Dercksen, V. J., Oberlaender, M., and Sompolinsky, H. (2016). The impact of structural heterogeneity on excitation-inhibition balance in cortical networks. Neuron 92, 1106–1121. doi: 10.1016/j.neuron.2016.10.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Landau, I. D., and Sompolinsky, H. (2018). Coherent chaos in a recurrent neural network with structured connectivity. PLoS Comput. Biol. 14:e1006309. doi: 10.1371/journal.pcbi.1006309

PubMed Abstract | CrossRef Full Text | Google Scholar

Le Van Quyen, M., Muller, L. E., Telenczuk, B., Halgren, E., Cash, S., Hatsopoulos, N. G., et al. (2016). High-frequency oscillations in human and monkey neocortex during the wake-sleep cycle. Proc. Natl Acad. Sci. U.S.A. 113, 9363–9368. doi: 10.1073/pnas.1523583113

PubMed Abstract | CrossRef Full Text | Google Scholar

Lefort, S., Tomm, C., Sarria, J.-C. F., and Petersen, C. C. (2009). The excitatory neuronal network of the {C2} barrel column in mouse primary somatosensory cortex. Neuron 61, 301–316. doi: 10.1016/j.neuron.2008.12.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Lehnertz, K., Bialonski, S., Horstmann, M.-T., Krug, D., Rothkegel, A., Staniek, M., et al. (2009). Synchronization phenomena in human epileptic brain networks. J. Neurosci. Methods 183, 42–48. doi: 10.1016/j.jneumeth.2009.05.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Lerchner, A., Ursta, C., Hertz, J., Ahmadi, M., Ruffiot, P., and Enemark, S. (2006). Response variability in balanced cortical networks. Neural Comput. 18, 634–659. doi: 10.1162/neco.2006.18.3.634

PubMed Abstract | CrossRef Full Text | Google Scholar

Litwin-Kumar, A., and Doiron, B. (2012). Slow dynamics and high variability in balanced cortical networks with clustered connections. Nat. Neurosci 15, 1498–1505. doi: 10.1038/nn.3220

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, J. H., and Riecke, H. (2018). Synchronization by uncorrelated noise: interacting rhythms in interconnected oscillator networks. Sci. Rep. 8, 1–14. doi: 10.1038/s41598-018-24670-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Mongillo, G., Rumpel, S., and Loewenstein, Y. (2018). Inhibitory connectivity defines the realm of excitatory plasticity. Nat. Neurosci. 21, 1463–1470. doi: 10.1038/s41593-018-0226-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Montbrió, E., Pazó, D., and Roxin, A. (2015). Macroscopic description for networks of spiking neurons. Phys. Rev. X 5:021028. doi: 10.1103/PhysRevX.5.021028

CrossRef Full Text | Google Scholar

Monteforte, M., and Wolf, F. (2010). Dynamical entropy production in spiking neuron networks in the balanced state. Phys. Rev. Lett. 105:268104. doi: 10.1103/PhysRevLett.105.268104

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakagawa, N., and Kuramoto, Y. (1993). Collective chaos in a population of globally coupled oscillators. Progr. Theor. Phys. 89, 313–323. doi: 10.1143/ptp/89.2.313

CrossRef Full Text | Google Scholar

O'Connor, D. H., Peron, S. P., Huber, D., and Svoboda, K. (2010). Neural activity in barrel cortex underlying vibrissa-based object localization in mice. Neuron 67, 1048–1061. doi: 10.1016/j.neuron.2010.08.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Okun, M., and Lampl, I. (2008). Instantaneous correlation of excitation and inhibition during ongoing and sensory-evoked activities. Nat. Neurosci. 11, 535. doi: 10.1038/nn.2105

PubMed Abstract | CrossRef Full Text | Google Scholar

Okun, M., Yger, P., Marguet, S. L., Gerard-Mercier, F., Benucci, A., Katzner, S., et al. (2012). Population rate dynamics and multineuron firing patterns in sensory cortex. J. Neurosci. 32, 17108–17119. doi: 10.1523/JNEUROSCI.1831-12.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Olmi, S., Politi, A., and Torcini, A. (2011). Collective chaos in pulse-coupled neural networks. EPL 92, 60007. doi: 10.1209/0295-5075/92/60007

CrossRef Full Text | Google Scholar

Ostojic, S. (2014). Two types of asynchronous activity in networks of excitatory and inhibitory spiking neurons. Nat. Neurosci. 17, 594–600. doi: 10.1038/nn.3658

PubMed Abstract | CrossRef Full Text | Google Scholar

Ott, E. (2002). Chaos in Dynamical Systems. Cambridge: Cambridge University Press.

Google Scholar

Ott, E., and Antonsen, T. M. (2008). Low dimensional behavior of large systems of globally coupled oscillators. Chaos 18, 037113. doi: 10.1063/1.2930766

PubMed Abstract | CrossRef Full Text | Google Scholar

Petersen, P. C., and Berg, R. W. (2016). Lognormal firing rate distribution reveals prominent fluctuation-driven regime in spinal motor networks. Elife 5:e18805. doi: 10.7554/eLife.18805

PubMed Abstract | CrossRef Full Text | Google Scholar

Pikovsky, A., and Politi, A. (2016). Lyapunov Exponents: A Tool to Explore Complex Dynamics. Cambridge: Cambridge University Press.

Google Scholar

Politi, A., Ullner, E., and Torcini, A. (2018). Collective irregular dynamics in balanced networks of leaky integrate-and-fire neurons. Eur. Phys. J. Special Top. 227, 1185–1204. doi: 10.1140/epjst/e2018-00079-7

CrossRef Full Text | Google Scholar

Pyle, R., and Rosenbaum, R. (2016). Highly connected neurons spike less frequently in balanced networks. Phys. Rev. E 93:040302. doi: 10.1103/PhysRevE.93.040302

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramirez, J.-M., and Baertsch, N. A. (2018). The dynamic basis of respiratory rhythm generation: one breath at a time. Annu. Rev. Neurosci. 41,475–499. doi: 10.1146/annurev-neuro-080317-061756

PubMed Abstract | CrossRef Full Text | Google Scholar

Renart, A., de la Rocha, J., Bartho, P., Hollender, L., Parga, N., Reyes, A., et al. (2010). The asynchronous state in cortical circuits. Science 327, 587–590. doi: 10.1126/science.1179850

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenbaum, R., and Doiron, B. (2014). Balanced networks of spiking neurons with spatially dependent recurrent connections. Phys. Rev. X 4:021039. doi: 10.1103/PhysRevX.4.021039

CrossRef Full Text | Google Scholar

Roxin, A., Brunel, N., Hansel, D., Mongillo, G., and van Vreeswijk, C. (2011). On the distribution of firing rates in networks of cortical neurons. J. Neurosci. 31, 16217–16226. doi: 10.1523/JNEUROSCI.1677-11.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Shadlen, M. N., and Newsome, W. T. (1994). Noise, neural codes and cortical organization. Curr. Opin Neurobiol. 4, 569–579. doi: 10.1016/0959-4388(94)90059-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Shadlen, M. N., and Newsome, W. T. (1998). The variable discharge of cortical neurons: implications for connectivity, computation, and information coding. J. Neurosci. 18, 3870–3896. doi: 10.1523/JNEUROSCI.18-10-03870.1998

PubMed Abstract | CrossRef Full Text | Google Scholar

Shibata, T., and Kaneko, K. (1998). Collective chaos. Phys. Rev. Lett. 81:4116. doi: 10.1103/PhysRevLett.81.4116

CrossRef Full Text | Google Scholar

Shu, Y., Hasenstaub, A., and McCormick, D. A. (2003). Turning on and off recurrent balanced cortical activity. Nature 423, 288–293. doi: 10.1038/nature01616

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, M. A., and Kohn, A. (2008). Spatial and temporal scales of neuronal correlation in primary visual cortex. J. Neurosci. 28, 12591–12603. doi: 10.1523/JNEUROSCI.2929-08.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Softky, W. R., and Koch, C. (1993). The highly irregular firing of cortical cells is inconsistent with temporal integration of random epsps. J. Neurosci. 13, 334–350. doi: 10.1523/JNEUROSCI.13-01-00334.1993

PubMed Abstract | CrossRef Full Text | Google Scholar

Tiesinga, P., and Sejnowski, T. J. (2009). Cortical enlightenment: are attentional gamma oscillations driven by ing or ping? Neuron 63, 727–732. doi: 10.1016/j.neuron.2009.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Ullner, E., Politi, A., and Torcini, A. (2018). Ubiquity of collective irregular dynamics in balanced networks of spiking neurons. Chaos 28:081106. doi: 10.1063/1.5049902

PubMed Abstract | CrossRef Full Text | Google Scholar

Ullner, E., Politi, A., and Torcini, A. (2020). Quantitative and qualitative analysis of asynchronous neural activity. Phys. Rev. Res. 2:023103. doi: 10.1103/PhysRevResearch.2.023103

CrossRef Full Text | Google Scholar

van Vreeswijk, C. (1996). Partial synchronization in populations of pulse-coupled oscillators. Phys. Rev. E 54:5522. doi: 10.1103/PhysRevE.54.5522

PubMed Abstract | CrossRef Full Text | Google Scholar

van Vreeswijk, C., and Sompolinsky, H. (1996). Chaos in neuronal networks with balanced excitatory and inhibitory activity. Science 274, 1724–1726. doi: 10.1126/science.274.5293.1724

PubMed Abstract | CrossRef Full Text | Google Scholar

Volgushev, M., Chauvette, S., and Timofeev, I. (2011). Long-range correlation of the membrane potential in neocortical neurons during slow oscillation. Prog Brain Res. 193, 181–199. doi: 10.1016/B978-0-444-53839-0.00012-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Whittington, M. A., Cunningham, M. O., LeBeau, F. E., Racca, C., and Traub, R. D. (2011). Multiple origins of the cortical gamma rhythm. Dev. Neurobiol. 71, 92–106. doi: 10.1002/dneu.20814

PubMed Abstract | CrossRef Full Text | Google Scholar

Whittington, M. A., Traub, R. D., Kopell, N., Ermentrout, B., and Buhl, E. H. (2000). Inhibition-based rhythms: experimental and mathematical observations on network dynamics. Int. J. Psychophysiol. 38, 315–336. doi: 10.1016/S0167-8760(00)00173-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Wohrer, A., Humphries, M. D., and Machens, C. K. (2013). Population-wide distributions of neural activity during perceptual decision-making. Progr. Neurobiol. 103, 156–193. doi: 10.1016/j.pneurobio.2012.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar