Hopf Bifurcation in Mean Field Explains Critical Avalanches in Excitation-Inhibition Balanced Neuronal Networks: A Mechanism for Multiscale Variability
1 Hopf Bifurcation in Mean Field Explains Critical Avalanches in Excitation-Inhibition Balanced Neuronal Networks: A Mechanism for Multiscale Variability
Junhao Liang , Tianshou Zhou and Changsong Zhou Department of Physics, Centre for Nonlinear Studies and Beijing-Hong Kong-Singapore Joint Centre for Nonlinear and Complex Systems (Hong Kong), Institute of Computational and Theoretical Studies, Hong Kong Baptist University, Kowloon Tong, Hong Kong School of Mathematics and
Guangdong Key Laboratory of Computational Mathematics, Sun Yat-sen University, Guangzhou, P. R. China Department of Physics, Zhejiang University, 38 Zheda Road, Hangzhou, China * Correspondence: [email protected]
Keywords: Excitation-Inhibition balance, Neuronal avalanches, Criticality, Cortical dynamics, Mean-field theory
Abstract
Cortical neural circuits display highly irregular spiking in individual neurons but variably sized collective firing, oscillations and critical avalanches at the population level, all of which have functional importance for information processing. Theoretically, the balance of excitation and inhibition inputs is thought to account for spiking irregularity and critical avalanches may originate from an underlying phase transition. However, the theoretical reconciliation of these multilevel dynamic aspects in neural circuits remains an open question. Herein, we study excitation-inhibition (E-I) balanced neuronal network with biologically realistic synaptic kinetics. It can maintain irregular spiking dynamics with different levels of synchrony and critical avalanches emerge near the synchronous transition point. We propose a novel semi-analytical mean-field theory to derive the field equations governing the network macroscopic dynamics. It reveals that the E-I balanced state of the network manifesting irregular individual spiking is characterized by a macroscopic stable state, which can be either a fixed point or a periodic motion and the transition is predicted by a Hopf bifurcation in the macroscopic field. Furthermore, by analyzing public data, we find the coexistence of irregular spiking and critical avalanches in the spontaneous spiking activities of mouse cortical slice in vitro , indicating the universality of the observed phenomena. Our theory unveils the mechanism that permits complex neural activities in different spatiotemporal scales to coexist and elucidates a possible origin of the criticality of neural systems. It also provides a novel tool for analyzing the macroscopic dynamics of E-I balanced networks and its relationship to the microscopic counterparts, which can be useful for large-scale modeling and computation of cortical dynamics. a r X i v : . [ c ond - m a t . d i s - nn ] J a n
1. Introduction
The mammal brain consists of tens of billions of neurons, which process information and communicate through electrophysiological action potentials, also known as spikes. This large number of neurons exhibit diverse spiking behaviors across broad ranges of spatial and temporal scales (Okun et al., 2015; Stringer et al., 2019a, 2019b). Understanding the origin and dynamic mechanism of this complexity is crucial for the advancement of neurobiology, the development of therapies for brain diseases and the future design of brain-inspired intelligent systems. Two striking features can be simultaneously observed at different levels of cortical neuronal systems: 1) irregularity in spiking times, indicated by seemingly random spiking time that resembles Poisson process (Softky and Koch, 1993; Holt et al., 1996) of individual neurons; 2) variability in population firing rates, manifested in widely observed collective neural activities such as population oscillations (Brunel and Wang, 2003; Herrmann et al., 2004) and critical neural avalanches (Beggs and Plenz, 2003; Gireesh and Plenz, 2008; Friedman et al., 2012; Bellay et al., 2015), etc. Biologically, the spiking irregularity has been proposed to originate from the balance between excitation (E) and inhibition (I) inputs so that spiking of neurons is driven by fluctuations (Shu et al., 2003; Okun and Lampl, 2008; Xue et al., 2014), and has been associated with functional advantages in efficient coding and information processing (DenΓ¨ve and Machens, 2016). The emergence of collective cortical activities originates in the fact that neurons interact through recurrent networks (Abeles, 1991), in which dynamic activities can reverberate. As a result, dynamic correlations arise from structural correlations. In particular, even weak pairwise correlation is sufficient to induce strongly correlated collective network activities (Schneidman et al., 2006). Collective neural activities can emerge with different amplitudes and are often organized as critical avalanches with various sizes (Beggs and Plenz, 2003; Gireesh and Plenz, 2008; Friedman et al., 2012; Bellay et al., 2015; Fontenele et al., 2019). These avalanches are cascades of activity bursts in neuronal networks. At criticality, the size and duration of avalanches are approximately distributed according to power-laws, with critical exponents satisfying the scaling relation (Friedman et al., 2012; Fontenele et al., 2019). Avalanches in the critical state can maximize the informational complexity and variability, and are thought to have functional advantages in information processing (Kinouchi and Copelli, 2006; Shew et al., 2009, 2011). Traditional mean-field theory of E-I balanced networks (Van Vreeswijk and Sompolinsky, 1996, 1998; Renart et al., 2010) with binary neuron and instantaneous synapse explains the spiking irregularity of individual neurons. However, it fails to account for collective neural activities, because it predicts an asynchronous dynamic state with vanishing correlation in unstructured (i.e. random topology) networks. Such vanishing correlation arises due to sparse network connectivity (Van Vreeswijk and Sompolinsky, 1998) or shared excitatory and inhibitory inputs cancelling correlation in dense recurrent networks (Renart et al., 2010). In terms of rate coding, the asynchronous state is not efficient for information processing, as the population firing rate only exhibits a linear response to the input rate (Van Vreeswijk and Sompolinsky, 1996) but without firing rate variability on faster time scales. In this case, the whole network acts as a rate unit for computation. The traditional E-I balanced theory can be generalized in two directions, which are biologically more plausible: structured networks and synaptic kinetics. Firstly, heterogeneous neural network structures (Landau et al., 2016) can induce firing rate variability. For example, clustered network structures (Litwin-Kumar and Doiron, 2012) can induce slow firing rate oscillations and show stimulus-induced variability reductions. Hierarchical modular networks (Wang et al., 2011) can support self-sustained firing rate oscillations across different levels. Spatial networks involving a distance-dependent coupling rule can unveil the distance-tuned correlation relation (Rosenbaum et al., 2017; Darshan et al., 2018) and the emergence of propagating waves (Keane and Gong, 2015; Keane et al., 2018; Gu et al., 2019; Huang et al., 2019) observed in experiments. Secondly, even in unstructured networks, network firing rate oscillations can be induced by realistic synaptic filtering kinetics (Brunel and Wang, 2003; Yang et al., 2017). Such oscillations typically occur in cases where the synaptic decay time scales of inhibition are slower than excitation, which is actually a biologically plausible situation if the synaptic receptors under consideration are AMPA for excitation and GABA for inhibition (Salin and Prince, 1996; Zhou and Hablitz, 1998). More importantly, network oscillation can be sparsely participated by subgroups of neurons, thus preserving the irregular spiking feature of individual neurons (Brunel, 2000; Brunel and Hakim, 2008). Note, however, that theoretical analysis of the macroscopic dynamics of E-I neural network with synaptic kinetics is very difficult. Existing theory is very limited, e.g., by requiring very special assumptions, such as the Lorenz distribution of certain parameters (Dumont and Gutkin, 2019). Critical avalanches can also rise from neural dynamics under unstructured network topology (Beggs and Plenz, 2003; Kinouchi and Copelli, 2006), while its dynamic origin is still controversial. Previous theories have suggested that critical neural avalanches may arise at the edge of a phase transition. Early studies indicated that it may occur between a quiescent and active phase from critical branching processes (Beggs and Plenz, 2003; Haldeman and Beggs, 2005), while later experimental (Fontenele et al., 2019) and theoretical (di Santo et al., 2018; Dalla Porta and Copelli, 2019) studies also proposed that it may occur near the onset of synchrony. Mechanisms other than criticality that generate avalanches with power-law distribution have also been proposed (Martinello et al., 2017; Touboul and Destexhe, 2017; Wilting and Priesemann, 2019b). The emergence of critical avalanches has also been proposed to be closely related to the maintenance of E-I balance (Lombardi et al., 2012; Poil et al., 2012; Yang et al., 2012). Nevertheless, the exact relationship between neural criticality and E-I balance remains poorly understood. Here, we try to address the above important open questions. In particular, how E-I balance induced irregular spiking reconciles with collective neural activities? E-I neural networks can organize into a βsparse synchronyβ state (Brunel, 2000; Brunel and Hakim, 2008) where neurons are remained fluctuation-driven to spike irregularly whereas firing rate oscillation emerges in population level. However, the relationship between E-I balance and sparse synchrony is less clear. Most importantly, the mechanism by which E-I balance induced irregular spiking can coexist with critical neural activities (Bellay et al., 2015) in recurrent neural circuit remains unclear. In this work, we first re-examine the dynamics of integrate-and-fire (IF) E-I neuronal network with realistic synaptic kinetics, which can manifest individually irregular spiking activities with different synchronous levels. The network firing rate dynamics can be effectively captured by a set of macroscopic field equations derived by a novel semi-analytical mean-field theory. An advantage of our theory is that it is simple and does not require special properties of the model. The synchronous transition point where network firing rate oscillation emerges is predicted by a Hopf bifurcation in the field equations. We find that critical microscopic avalanche dynamics emerges near the onset of synchronization, with critical exponents approximately satisfying the scaling relations, which manifests the hallmark of criticality (Sethna et al., 2001). The mechanism of critical avalanches could be understood as demographic noise-driven random walks near a macroscopic bifurcation point. On this basis, we propose that the E-I balanced state in the microscopic spiking network corresponds to a stable macroscopic state in the field equations. The asynchronous state, consistent with the traditional theory, corresponds to a stable fixed point, which can be destabilized through a Hopf bifurcation, giving rise to a stable limit circle, corresponding to network firing rate oscillation. As such, the E-I balanced state can incorporate network oscillation with different synchronous levels, which accounts for the coexistence of variability in both individual and population scales. Finally, we empirically verify the coexistence of irregular spiking and collective critical avalanches in the public experimental data of spontaneous spiking activities recorded in mouse somatosensory cortex in vitro (Ito et al., 2016). Scaling relations similar to the network model are also found to hold in these critical data sets. Our own analysis further indicates the universality of the observed phenomena in model networks. The theory proposed here explains how collective neural activities coexist with irregular neuron spiking and reveals a possible origin of criticality in neural systems. It also serves as a novel tool to study the dynamics of IF networks with biologically realistic synaptic filtering kinetics, and thus has useful application in large-scale modeling of brain networks.
2. Materials and Methods
We study a leaky IF spiking neuronal circuit. Neurons are coupled by a random network with density π and size = π πΈ + π πΌ , which consists of π πΈ excitatory (E) neurons and π πΌ inhibitory (I) neurons. Thus, each neuron in the network has on average π πΈ = ππ πΈ E neighbors and π πΌ = ππ πΌ I neighbors. Each neuron also receives π π excitatory inputs modelled by independent Poisson processes with frequency π π , mimicking external inputs for the circuit under consideration. We set π = 0.2 , π πΈ : π πΌ = 4: 1 , π π = π πΈ and the network size is π = 10 , unless otherwise specified. The sub-threshold membrane potential of neuron π at time π‘ , denoted as π π (π‘) , is governed by ππ π ππ‘ = π πΌ (π π ) + π½ πΌπ β πΉ πΈ β π π (π‘ β π ππΈ ) πβπ ππ + π½ πΌπΈ β πΉ πΈ β π π (π‘ β π ππΈ ) πβπ ππΈ (1) +π½ πΌπΌ β πΉ πΈ β π π (π‘ β π ππΌ ) πβπ ππΌ , where, π π (π‘) is the membrane potential of neuron π (belonging to type πΌ = πΈ, πΌ ) at time π‘ . π ππΌ represents the πΌ neighbors of neuron π . The input sources for neurons include excitatory inputs from external neurons (population π ), inputs from recurrent excitatory neurons (population πΈ ) and inputs from recurrent inhibitory neurons (population πΌ ). The first term of Equation (1) describes the leaky current π πΌ (π π ) = (π πππ π‘πΌ β π π )/π πΌ , which has the effect to drive the membrane potential back to the leaky potentials, which are set to be π πππ π‘πΈ = π πππ π‘πΌ = β70 ππ . The membrane time constants are set as π πΈ = 20 ππ , π πΌ = 10 ππ for E and I neurons, respectively. The second to fourth terms of Equation (1) are the external, excitatory recurrent and inhibitory recurrent currents, respectively. Input currents are the summations of the filtered pulse trains. Here, π π (π‘) = β πΏ(π‘ β π‘ ππ ) π denotes the spike train of the j -th neuron. The excitatory and inhibitory synapses have latency period (delay) π ππΈ and π ππΌ respectively. For the numerical results presented, we consider π ππΈ = π ππΌ = 0 for simplicity (i.e. no transmission delay), which is a reasonable approximation for local circuits. The synaptic filter is modelled as a bi-exponential function, i.e. Equation (2). πΉ πΌ (π‘) = ππΌ βπ π [exp (β π‘π ππΌ ) β exp (β π‘π π )] , π‘ β₯ 0 . (2) In Equation (2), we set the synaptic rise time π π = 0.5 ππ for both E and I neurons, while the synaptic decay times π ππΈ , π ππΌ depend on the type of presynaptic neuron. We set π ππΈ = 2 ππ and let π ππΌ change from to study the effect of different E and I synaptic filtering time scales. Hence, in our study here, π ππΌ serves as a control parameter to induce the dynamical transition. Biologically, the inhibition decay time π ππΌ depends on the constitution of synaptic receptors (Salin and Prince, 1996; Zhou and Hablitz, 1998), and can also be changed by chemicals such as narcotics (Brown et al., 2010). We point out that suitable changes of other model parameters such as synaptic strength, network connection density, etc. may also induce similar dynamical transitions we are going to study below. The integration dynamics is as follows. When the membrane potential reaches the threshold π π‘β = β50 ππ , a spike is emitted and the membrane potential is reset to π πππ ππ‘ = β60 ππ . Then, synaptic integration is halted for 2 ms for E neurons and 1 ms for I neurons, modelling the refractory periods in real neurons. Synaptic weights are set as π½ πΈπ =0.45 ππ , π½ πΌπ = 0.72 ππ , π½ πΈπΈ = 0.36 ππ , π½ πΌπΈ = 0.72 ππ , π½ πΈπΌ = β0.81 ππ and π½ πΌπΌ =β1.44 ππ , which will satisfy the balanced condition. Network dynamics are simulated by a modified second-order Runge-Kutta scheme (Shelley and Tao, 2001) with a time step of ππ‘ =0.05 ππ . For each parameter, the network is simulated for
16 s with the first discarded to avoid a transient effect. The statistical indexes are then computed by averaging the results of 15 trials with randomly distributed initial membrane potentials. The dynamics we considered here are current based. The case of conductance-based dynamics where the input of a neuron depends on its membrane potential is further studied in Supplementary Material:
Appendix 1. Model Extensions. . We will develop a semi-analytical mean-field theory to approximate the average dynamics of the network by noting the fact that 1) the network topology is homogeneous and 2) the number of neuron is large. The following mean-field derivation holds for large enough and the network size is explicitly included in the field model equations (see Equations (11), (12) below). We denote the average membrane potential of the network as π πΌ = β©π π βͺ πβπΌ β β©π π βͺ πΌ , πΌ = πΈ, πΌ . The goal is to derive a set of self-consistent field equations governing the temporal dynamics of π πΌ . First, by noting that the convolution πΉ πΌ β π π (π‘ β π ππΌ ) obeys the equation (π ππΌ πππ‘ + 1) (π π πππ‘ + 1) [πΉ πΌ β π π (π‘ β π ππΌ )] = β πΏ(π‘ β π‘ ππ β π ππΌ ) π , (3) then we have (π ππΌ πππ‘ + 1) (π π πππ‘ + 1) β©β πΉ πΌ β π π (π‘ β π ππΌ ) πβπ ππΌ βͺ πβπΈ,πΌ = β©β β πΏ(π‘ β π‘ ππ β π ππΌ ) ππβπ ππΌ βͺ πβπΈ,πΌ . (4) Under mean-field approximation, each neuron is the same in terms of their neighbors, so that β©β β πΏ(π‘ β π‘ ππ β π ππΌ ) ππβπ ππΌ βͺ πβπΈ,πΌ β π πΌ π πΌ (π‘ β π ππΌ ) , (5) where π πΌ is the average number of neighbors of a neuron in the network and π πΌ (π‘) is the mean firing rate of πΌ type neurons, defined as π πΌ (π‘) = lim βπ‘β0 1βπ‘ β« β©β πΏ(π β π‘ ππ ) π βͺ πβπΌ ππ π‘+βπ‘π‘ . (6) In the standard definition of firing rate (Dayan and Abbott, 2001) of a neuron, the average in Equation (6) is taken over different simulation trials. Since ergodicity and the network homogeneity, neurons within a population should have the same firing rate and it can be computed through Equation (6) (i.e. population averages can be thought as sample averages) when the network is large enough in the stationary state. Formally, for measuring the firing rate from data, the time interval βπ‘ in Equation (6) has to be finite. Here we choose βπ‘ = 1 ππ (shorter than the refractory period) so that a neuron can at most have one spike between π‘ and π‘ + βπ‘ . Then, by the definition of πΏ function, π πΌ (π‘) represents the proportion of πΌ type neurons that spike between π‘ and π‘ + βπ‘ as well as the mean firing rate of πΌ type neurons at time π‘ with unit per ms . In previous analysis framework of IF neurons through continuous stochastic process theory (Burkitt, 2006), the membrane potential π π of neuron π cannot cross the spiking threshold ( π π is restricted to (ββ, π π‘β ) with π π‘β being an absorbing boundary). This is a theoretical artefact, contrary to the true neurophysiology. Furthermore, in numerical integration, the resetting is achieved by finding those neurons whose membrane potential increases over the spiking threshold in each numerical step (Shelley and Tao, 2001). This inspires us to naturally consider that a neuron π should have a spike at time π‘ if π π (π‘) > π π‘β . Formally, we can consider π π (π‘) as the membrane potential of neuron π at time π‘ before the resetting rule in each numerical step, then π πΌ (π‘) = β©π»(π π β π π‘β )βͺ πβπΌ , (7) N N ο‘ where π» is the Heaviside function π»(π₯) = {1, π₯ β₯ 00, π₯ < 0 . Equation (7) explicitly builds the link that the population firing rate is the proportion of the neurons whose membrane potential is above the spiking threshold. As a preliminary approximation, we assume π π (π‘) obey a Gaussian distribution π πΌ (π) with time-dependent mean π πΌ (π‘) and time-independent variance π πΌ2 . We will verify that this assumption is plausible in the network and dynamic regimes we studied, referring to Figure 1D below. Then, π πΌ (π‘) = β©π»(π π β π π‘β )βͺ πβπΌ = β« π πΌ (π)ππ βπ π‘β = β erf ( π π‘β βπ πΌ β2π πΌ ) , (8) where erf is the error function πππ(π₯) = β« π βπ‘ ππ‘ π₯0 . Although there is no elementary expression for the error function, it can be approximated by elementary functions. For example, a good approximation that can keep the first and second moments is πππ(π₯) β tanh ( ππ₯β6 ) . Under this approximation, we have π πΌ (π‘) = β tanh ( πβ6 π π‘β βπ πΌ β2π πΌ ) = ππ‘ββππΌππΌ πβ3 ) . (9) Here, the standard deviation of the voltage, π πΌ , acts as an effective parameter to construct the voltage-dependent temporal firing rate. Note that this approximation scheme basing only on the first-order statistics neglects several factors that affect the accurate firing rate, including higher order statistics, noise correlation and refractory time. Thus, it does not have an analytical form and should be estimated numerically, from the steady-state mean voltage π πΌπ π and mean firing rate π πΌπ π and Equation (9) that π πΌ = π π‘β βπ πΌπ π ln [(π πΌπ π ) β1 β1] πβ3 . (10) The sigmoid transfer function Equation (9) is the intrinsic nonlinear property that induces oscillation transition in the field model. Note that neural field models of Wilson-Cowan type (Wilson and Cowan, 1972) would also contain a presumed sigmoid transfer function. Field models of this type can also qualitatively reproduce some dynamic features of E-I neuronal networks. Here we explicitly construct the sigmoid function from the microscopic spiking network. Thus, the quality of the scheme depends on suitable choices of effective parameters π πΌ and once π πΌ are chosen as suitable values, our field equations can predict the dynamics of the E-I network quantitatively (see Supplementary Material: Appendix 2. Sensitivity of the Critical Points on the Effective Parameters, where we study the sensitivity of the oscillation transition with respect to the effective parameters π πΌ ). Denote Ξ¦ πΌ (π‘) = β©β πΉ πΌ β π π (π‘ β π ππΌ ) πβπ ππΌ βͺ πΈ,πΌ as the averaged synaptic time course of πΌ inputs received by a neuron and from Equations (4), (5), (9) it will obey (π ππΌ πππ‘ + 1) (π π πππ‘ + 1) Ξ¦ πΌ = π πΌ [1+exp( ππ‘ββππΌ(π‘βπππΌ)ππΌ πβ3 )] , πΌ = πΈ, πΌ . (11) Note that each neuron receives π π independent Poisson spike trains externally with rate π π . Thus, the input of each neuron has a variance π π π π . If we do not consider its filtering effect since the fast decay time of excitatory synapse, by diffusion approximation we know the external input β πΉ πΈ β π π (π‘ β π ππΈ ) πβπ ππ of each neuron can be approximated by π π π π + βπ π π π π π (π‘) , where π π (π‘) is a standard Gaussian white noise (GWN) with zero mean and unit variance. Since {π π (π‘)} π are independent GWNs, β©βπ π π π π π (π‘)βͺ πΌ can be equivalently approximated by βπ π π π /π πΌ π πΌ (π‘) , where π πΌ (π‘) is another standard GWN. Thus, taking the average β©. βͺ πΌ of the original Equation (1) and note that in the leaky IF model, π πΌ is linear, we arrive at ππ πΌ ππ‘ = π πΌ (π πΌ ) + π½ πΌπ (π π π π + β π π π π π πΌ π πΌ (π‘)) + π½ πΌπΈ Ξ¦ πΈ + π½ πΌπΌ Ξ¦ πΌ , πΌ = πΈ, πΌ . (12) In the field Equation (12), π πΈ (π‘) and π πΌ (π‘) are two independent GWNs. We find that this approximation is independent of the nature of noise in the spiking network model. In the network model Equation (1), the nature of noise from external inputs is synaptic-filtered Poisson shot noise. We further examine the case where external inputs in Equation (1) are with GWNs and the case of constant external input (i.e. no noise) and find that in such cases Equation (12) can still well approximate the macroscopic dynamics of the network (for constant external input the network dynamics is not stochastic, but the spiking activity still appears irregular due to the chaotic nature of the network (Van Vreeswijk and Sompolinsky, 1996, 1998)). Generally speaking, the mean-field theory only holds when the system size is infinity. The incorporation of noise into the field model can smooth out the systematic errors, compensate the finite size effect and make it closer to the true rate dynamics statistically. Thus, for numerical simulation of the field equations we will keep the noise terms in Equation (12) whereas the deterministic counterpart would be used for stability analysis. In summary, we have proposed a novel technique to derive a set of self-consistent field Equations (11), (12), to approximate the average dynamics of the original spiking network Equations (1), (2). The deterministic steady-state (fixed point) of the field model can be found from Equations (11), (12) by letting πππ‘ = 0 and assuming π πΌ (π‘) = 0 , resulting in algebraic equations π πΌ (π πΌ ) + π½ πΌπ π π π π + π½ πΌπΈ π πΈ π πΈ + π½ πΌπΌ π πΌ π πΌ = 0 , πΌ = πΈ, πΌ . (13) where π πΈ , π πΌ are given by Equation (9). In our case, the synaptic strengths and external input rates are the major parameters determining the value of the fixed point, while synaptic decay time would affect its stability. This is because the value of steady-state does not depend on π ππΌ , π π , which can also be seen from Equation (2) that the synaptic filter is normalized ( β« πΉ πΌ β πΏ(π‘) β0 = 1 independent of π ππΌ , π π ) so that synaptic rise and decay times would not affect the time-averaged firing rate. Thus, the scheme here cannot capture the nontrivial effect of synaptic filtering on affecting the firing rate (e.g. see the formula given by Fourcaud and Brunel (Fourcaud and Brunel, 2002)). Note that in general settings of balanced networks with dense and strong coupling (Renart et al., 2010), the quantities π π , π πΈ , π πΌ are of order π(π) but the synaptic weights are of order
π(π β1/2 ) . When π is large enough, the first term in Equation (13) can be neglected and it reduces to π½ πΌπ π π π π + π½ πΌπΈ π πΈ π πΈ + π½ πΌπΌ π πΌ π πΌ = 0 , πΌ = πΈ, πΌ , (14) which is a set of linear equations to solve the steady firing rate π πΈ , π πΌ . To guarantee a unique positive solution in this case, the sequence { π½ πΈπ π½ πΌπ , π½ πΈπΌ π½ πΌπΌ , π½ πΈπΈ π½ πΌπΈ } should be in ascending or descending order, which is the so-called balanced condition in the traditional theory (Van Vreeswijk and Sompolinsky, 1996, 1998; Renart et al., 2010). Thus, the theory here is a generalization of the traditional theory of balanced network. In the traditional theory of asynchronous dynamics, the E-I balanced state can be considered as the existence of a stable fixed point of Equations (11), (12). Now we can consider how the stability of this fixed point can be changed with the aid of this dynamic form. For the case without synaptic delay (i.e. π ππΈ = π ππΌ = 0 ), the field equations are ordinary differential equations. By taking π =(π πΈ , π πΌ , Ξ¦ πΈ , πΞ¦ πΈ ππ‘ , Ξ¦ πΌ , πΞ¦ πΌ ππ‘ ) T , the field model Equations (11), (12) can be written in the first-order form ππππ‘ = πΉ(π) without considering the noise. The Jacobian matrix at the steady state is π½ = ( β πΈ
0 π½ πΈπΈ
0 π½ πΈπΌ
00 β πΌ π½ πΌπΈ
0 π½ πΌπΌ π πΈ π πΈβ² (π πΈ )π ππΈ π π
0 β ππΈ π π β ππΈ β π
0 0 0 0 0 0 0 1 0 π πΌ π πΌβ² (π πΌ )π ππΌ π π
0 0 β ππΌ π π β ππΌ β π ) , (15) with π πΌβ² (π πΌ ) = πexp [(π π‘β βπ πΌ )π/(β3π πΌ )]β3π πΌ (1+exp[ (ππ‘ββππΌ)πβ3ππΌ ]) estimated at the steady-state value of π πΌ given by Equation (13). The eigenvalues of π½ can determine the stability of the steady state. Note that many models use single exponential function as the synaptic filter, i.e. π π = 0 , and in this case the dynamic form becomes 4-dimensional with π = (π πΈ , π πΌ , Ξ¦ πΈ , Ξ¦ πΌ ) T . For models without considering synaptic filtering effect (that is, the case of instantaneous synapse where π π = π ππΌ = 0 ), the dynamic form becomes 2-dimensional, which can be considered as the dynamic form of traditional theory without synaptic kinetics. In these cases, the stability analysis can be performed in a similar way. In the presence of synaptic transmission delay, the field equations would become delay differential equations. In this case, stability analysis would in general become more difficult. For statistically analyzing neural dynamics, the neuron spike train series have to be constructed from the model simulation data or the experimental data as follows. The time axis is first divided into consecutive time windows with sizes Ξπ‘ ms. The number of spikes of neuron π are then counted in each window to obtain a discrete sequence π π (π‘) , which is designated as the spike count series of neuron π with time windows Ξπ‘ . Alternatively, the number of spikes of the whole neuron population can be counted in each window. This constructs the population spike count series π πΌ (π‘) for the E and I populations, respectively. Furthermore, π πΌ (π‘) = π πΌ (π‘)π πΌ Ξπ‘ is the population averaged firing rate series and the power spectrum of the population firing rate series can indicate the collective oscillation property. For computing certain quantities, such as the correlation between neurons, the spike count series is filtered by a square kernel πΎ π (π‘) = { , π‘ β [βπ, 0]0, ππ‘βππ with length π and the ensuing filtered spike train is defined as πΜ π (π‘) = πΎ π β π π = β π π (π‘ β π ) πβ1π =0 . The spiking time irregularity of a neuron can be quantified using the CV (coefficient of variance), which is defined as the standard deviation of the neuron ISIs (inter-spike intervals) over its mean. Totally regular activities have CV values of 0, while Poisson processes have CV values of 1. A higher CV value indicates larger irregularity. We measure the firing rate variability on the individual neuron level and the population level by the relevant Fano Factor (FF). The FF of neuron π is defined as πΉπΉ π = π£ππ(π π (π‘))/β©π π (π‘)βͺ . The average and the variance here are taken across the spike count series. Similarly, population FF is defined as πΉπΉ πΌ = π£ππ(π πΌ (π‘))/β©π πΌ (π‘)βͺ . Note that FF depends on the time window size Ξπ‘ to construct the spike count series (we will use Ξπ‘ =50 ππ ). The population FF can also be computed using the macroscopic field equations, as the field equations predicts the population firing rate, which can be transferred to population spike counts (by multiplying the number of neurons in the population).
The synchrony of the network can be characterized from two aspects: the cross-correlation of spiking times and the coherence of the membrane potential (Golomb, 2007). The former quantifies the coherence of threshold events whereas the later quantifies the coherence of the subthreshold dynamics. We employ the commonly used Pearson correlation coefficient (PCC) to quantify the synchrony of the spiking time. The spike count series of neuron π is first constructed with time window Ξπ‘ = and then filtered by a square kernel with length π = 50 ππ . The PCC between neuron π and π is defined as π ππ = πππ£(π π (π‘),π π (π‘))βπ£ππ(π π (π‘))π£ππ(π π (π‘)) . The details such as filtering in calculating PCC would affect its exact value (Cohen and Kohn, 2011), but not the qualitative change. The index β©π ππ βͺ π,πβπΈ , PCC averaged over all excitatory neuron pairs, is used to quantify the network synchrony degree of threshold events. The voltage series of each neuron is constructed for each millisecond that π π (π) βΆ= π π (π‘)| π‘=π ππ . Voltage coherence is defined as = βπ πΌ2 /β©π π2 βͺ πΌ , where π πΌ2 = β©π πΌ2 βͺ π‘ β β©π πΌ βͺ π‘2 is the variance of the mean voltage π πΌ = β©π π βͺ πβπΌ and π π2 = β©π π2 βͺ π‘ β β©π π βͺ π‘2 is the variance of the voltage of neuron π . The voltage coherence of the excitatory population is used to quantify the coherence of the subthreshold dynamics. We further use the CV of the excitatory population firing rate series constructed with Ξπ‘ = 1 ππ to quantify temporal firing rate variability at short timescales, which is another way of indicating the network synchrony, as stronger synchrony indicates larger population firing rate variability at short timescales.
We measure the neuronal avalanches of the excitatory neuron population from its population spike count series π πΈ (π‘) constructed with window (bin) size Ξπ‘ . An avalanche is defined as a sequence of consecutive non-empty bins, separated by empty bins (with no spiking inside). The size π of an avalanche is defined as the total number of spikes within the period and the duration π is defined as the number of time bins it contains. To compare different data sets in a unified standard, window size is chosen as the average ISI of the merged spiking train (constructed by merging the spike trains of all neurons). Thus, it depends on the mean firing rate of neurons in different data sets. This choice has been described as the βoptimalβ window size to measure avalanches (Beggs and Plenz, 2003), as excessively small or large windows would lead to systematic bias. A further advantage of this choice is that the measured size and duration would approximately lie in the same scale ranges, allowing better comparison. Note that the definition of avalanche here is the time-binning (non-causal) avalanche, which corresponds to the case of experimental measurement (Beggs and Plenz, 2003; Bellay et al., 2015), but is different from the causal avalanche (Williams-Garcia et al., 2017) studied by many physical models rooted in critical branching processes. Avalanche size obeying a power-law distribution is a typical sign of scale-invariance, which is a statistical property. We use a simple method to quantify how far the avalanche size distribution deviates from a power-law distribution. The avalanche size frequency distribution histogram (π, π(π)) , is first obtained with 80 plotting-bins from minimum to maximum size. The least squares method is then used to find the best-fit-line in doubly logarithmic coordinates, such that β [πππ(π) β (π + π πππ)] is minimised. After obtaining the best-fit coefficients (π , π ) , the fitted frequency distribution values were estimated as π πππ‘ (π) = 10 π +π πππ . Finally, the normalised distance, defined as π· = β π|π(π) β π πππ‘ (π)| π / β ππ(π) π , the average size difference per avalanche between the actual and fitting frequency distribution normalized by the mean avalanche size, is used to measure the distance to the best fitting power-law distribution. This distance acts as an index to indicate the degree of criticality for comparison. The maximum likelihood estimation (MLE) method using the NCC toolbox (Marshall et al., 2016) is used to estimate the critical exponents. The toolbox provides a doubly truncated algorithm based on the MLE to find the range that passed the truncation based KS statistics test (Marshall et al., 2016). This truncation scheme can avoid the noise (in the small avalanche size range) and finite size effect (in the large avalanche size range) interruption in estimating the critical exponents. We find the largest truncated range that can pass the KS-based test with a p value larger than 0.1 and boarder than one-third of the whole range on the logarithm scale. This means that the data can produce a KS-statistics value that is less than the values generated by at least 10% of the power-law models in the truncated range. The estimated slopes within the truncated ranges in the avalanche size and duration distributions define the critical exponents π(π)~π βπ and π(π)~π βπΌ . A third exponent is defined as β©πβͺ(π)~π , where β©πβͺ(π) is the average size of avalanches with the same duration π . This exponent is fitted using a weighted least squares method (Marshall et al., 2016) to those avalanches that fall into the truncated duration range for estimating πΌ . We used the public experimental data measuring the neuronal spiking activity in mouse somatosensory cortex cultures in vitro (Ito et al., 2016). A total of 25 data sets were used and the length of each record is one hour, with the exception of Set 19, which was 48 minutes long. The recordings were performed on organo-typical slice cultures after two to four weeks in vitro , without stimulation (Ito et al., 2014). Spiking times were sorted with a PCA-based algorithm (Litke et al., 2004) to locate the signals recorded with a large and dense multi-electrode array of 512 electrodes.
Under in vitro conditions, spiking in the culture clearly shows up-down state transition. Active spiking periods (up-state) and silent periods (down-state) alternated slowly with a frequency of circa
Ξπ‘ = 10 ππ and then filtered by a square kernel with length
π = 100 ππ . Then, the up-state is defined as the time periods that the population firing rates are higher than 30% of the maximum rate of the given dataset and last longer than 1s. We further examine the power spectral density of the population firing rate series (before further filtering) in the up-state, the distribution of firing rates of the neurons, CV of ISIs and the avalanches in the up-states. The time bin for measuring the avalanche was the mean ISI, averaged through the up-state, of the merged spike train. To test critical properties, we first determine whether the size and duration distribution is close to a power-law and then estimate its critical exponents. As in our analysis of modelling data, the doubly truncated and statistical test algorithm from the NCC toolbox (Marshall et al., 2016) is used. We accept the power-law distribution of a data set if the following two conditions can be jointly satisfied: 1) the truncated range has to be boarder than one-third of the whole range in the logarithm scale and; 2) the data in the truncated range can pass the KS-based test with p value larger than 0.1. For the cases deemed to be power-law distribution, the truncated ranges, estimated exponents and p values in the KS-based test are shown. For the data sets that have power-law avalanche distributions in both size and duration, we further compare the critical exponents and the value πΌβ1πβ1 to see whether scaling relation (Equation (16) below) holds. To test whether power-law distribution is the intrinsic structure of the data, we also analyze the surrogated data constructed by randomly shuffling the ISIs of neurons in each up-state period. This random shuffling can destroy the intrinsic temporal correlation structures of the avalanches.
3. Results
To begin, we first intuitively illustrate the non-trivial effect of synaptic filtering (
Figure 1A ) in shaping the network dynamics. For fast inhibition decay time, the network spiking dynamics is asynchronous (Renart et al., 2010), as shown in
Figure 1B . Such a βstrict balanceβ is the case of traditional E-I balance theory (Van Vreeswijk and Sompolinsky, 1996, 1998) with inhibition domination, where network inhibitory feedback can cancel the excitatory current spill on a fast time scale (
Figure 1H ). This scenario of strict balance breaks down when inhibition becomes slow, which induces an effective delay in the inhibitory cancellation (
Figure 1I ). Such a delay before the cancellation would induce a window during which excitatory current spills over the network, resulting in the collective spiking shown in
Figure 1C . The strength of the collective activity depends on the length of this excitation-dominant window. As inhibition becomes dominant again after a delay, a temporal quiescent episode ensues after the collective spiking. Thus, slow synaptic inhibition induces network oscillation, as was also shown in previous studies (Brunel and Wang, 2003; Yang et al., 2017; Huang et al., 2019). In this βloose balanceβ scenario, the network maintains balance on a slower time scale. In this case, as excitation and inhibition dominate alternately, the strict balance is temporally broken on a fast time scale. Note that in the loose balance state, the delay in the inhibitory cancellation is merely a few milliseconds (
Figure 1I ), consistent with the inhibition tracking delays observed in experiments (Okun and Lampl, 2008). Alternatively, these different collective behaviors can be seen from the network population firing rate dynamics. The network firing rate in the strict balance, asynchronous state is almost steady (
Figure 1E ), whereas the loose balance induces firing rate oscillation (
Figure 1F ), with fast firing rate variability at the population level. The emergence of oscillation implies a regularity of the population dynamics of the network. However, it is important to stress that the spiking of the neurons is still highly irregular. There are at least two reasons for this irregularity. First, the population oscillation is quasi-periodic rather than periodic due to the stochastic nature of the network dynamics. Second, and more importantly, each collective burst is randomly participated in by partial neurons, manifested by a firing rate much smaller than the oscillation frequency, referring to Figure 2E . This can also be understood from the fact that the firing rate is still low at the peaks of spiking activity corresponding to spiking synchrony (
Figure 1F ). As such, the dynamical phenomenon here with loose balance was also referred to as βsparse synchronyβ (Brunel and Hakim, 2008). The histogram counts of the coefficient of variance (CV) of the inter-spike intervals (ISIs) of neurons (see Materials and Methods) in the network are shown in
Figure 1G . In both asynchronous and synchronous states, the overall value of the CV is around 1, indicating that the spiking irregularity approximately resembles the Poisson process (Holt et al., 1996). As our mean-field theory assumes that the membrane potentials of the neurons in the network approximate Gaussian, we further analyze the distribution of the membrane potentials. As shown in
Figure 1D , in both fast and slow inhibition cases the distribution is not totally Gaussian and instead right skewed, a feature of the finite threshold IF dynamic (Brunel, 2000; Keane and Gong, 2015). However, the skewness (defined as the third central moment divided by the cubic of the standard deviation) and kurtosis (defined as the fourth central moment divided by the quartic of the standard deviation) measuring the deviation from Gaussian are (-0.67, 3.12) and (-0.66, 3.3) in the fast and slow inhibition cases respectively, which are not far from (0, 3) for the case of Gaussian distribution. Thus, our theory that assumes a Gaussian distribution is still effective, as will be shown later. Figure 1.
Synchronization and network oscillation induced by slow inhibition in balanced networks . Fast inhibition induces strict balance with asynchronous spiking and an almost steady network firing rate. Slow inhibition results in loose balance, with synchronous grouped spiking and network oscillation. Individual neurons spike irregularly in both cases. (A)
Normalized E/I post-synaptic response when receiving a pre-synaptic spike. (B, C)
Raster plot of the spiking time. Each blue/red point corresponds to a spike of the E/I neuron. (D)
The distribution of membrane potential. (E, F)
The corresponding firing rate of E population (blue) and the firing rate predicted by the field equations (green). (H, I)
The blue/red curves represent the average input E/I current of a neuron in the network and the black curves represent the difference between them. (G)
The pdf. of the CV of ISI of E neurons. The parameters are set as π π = 5 π»π§ ; π ππΌ = 1 ππ in the fast inhibition case ( B, E, H and left panels of
D, G ) and π ππΌ = 3.5 ππ in the slow inhibition case ( C, F, I and right panels of
D, G ). The network dynamic transition induced by a looser E-I balance can be characterized by the dynamics of population firing rate (
Figure 1E, F ). However, it is theoretically challenging to analyze the population dynamics of IF networks with synaptic kinetics. Here, we propose a novel mean-field approximation theory to derive the macroscopic dynamic equations of an IF network, i.e. Equations (11), (12). This technique for deriving the macroscopic field equations is highly generalizable. Extensions to the cases of time-varying inputs and conductance-based dynamics are presented in Supplementary Material:
Appendix 1. Model Extensions. For numerically estimated π πΌ (i.e. use Equation (10) to compute π πΌ under different external input strength π π ), the mean firing rate of the network can be correctly estimated from the field Equations (11), (12) ( Figure 2A blue). We further consider whether the field equations can predict the firing rate with fixed parameters π πΌ . Following the derivation of Equation (12), if we simply assume ππ π ππ‘ β ππππ π‘ β π π π πΌ + π½ πΌπ βπ π π π π π (π‘) , then π π is Gaussian distributed with standard deviation π πΌ = βπ½ πΌπ2 π π π π π πΌ . Figure 2A green shows the corresponding results of fixed π πΌ =βπ½ πΌπ2 π π π π π πΌ for π π = 5 π»π§ . It cannot predict the exact firing rate but can correctly predict the linear response to external input, a property of asynchronous balanced network (Van Vreeswijk and Sompolinsky, 1996). The dynamic difference between the asynchronous state and synchronous state can already be predicted by the field equations in terms of the population firing rate dynamics, as shown in Figure 1E, F . The mechanism is explained by a Hopf bifurcation in the field equations through stability analysis (see Materials and Methods for further details), as shown in
Figure 2B . In the case of fast inhibition, the fixed point of the field model is generally a stable focus, whose Jacobian has complex eigenvalues with negative real parts. In this case, the network firing rate only fluctuates mildly due to noise perturbations. When π ππΌ increases, the fixed point will lose its stability through a supercritical Hopf bifurcation, as indicated by a pair of its dominant complex conjugate eigenvalues π = πΌ + ππ crossing the imaginary axis. The Hopf bifurcation predicts that the stable fixed point will give way to a stable periodic solution, whose amplitude grows from zero. The frequency of the periodic motion can be estimated as π/2π in the linear order. This prediction is approximately equal to the numerically measured peak frequency of the network excitatory firing rate oscillation near the critical bifurcation point, as shown by the blue circles in Figure 2B . Previous studies have shown that E-I networks can undergo a transition to oscillation through perturbation analysis by assuming the form of the network steady-state firing rate solution (Brunel and Wang, 2003; Brunel and Hakim, 2008) and the phenomenon has been associated with a Hopf bifurcation of rate equations with effective transmission delay (Brunel and Hakim, 2008). Here, our mean-field approach with the aid of macroscopic field equations derived from microscopic neuronal network straightforwardly reveals the Hopf bifurcation mechanism during this transition to sparse synchrony state. As shown in
Figure 2C , the network synchrony increases dramatically after π ππΌ crosses the Hopf bifurcation point. Figure 2D shows that the network temporal firing rate variability (see Materials and Methods for further details) increases conspicuously during this transition, which is also qualitatively predicted by the field equations. The CV of ISIs averaged over the excitatory population is shown in
Figure 2E . During the onset of collective oscillation, the averaged CV of ISIs first slightly decreases and then increases (as the strong bursting oscillation activity develops at around π ππΌ β 4 ππ , referring to Figure 3A bottom panel) while its overall value is around 1. The coexistence of irregular spiking and collective oscillation can be understood from the ratio between the mean firing rate of neurons and peak frequency of network oscillation. If each oscillation is participated by all the neurons, then the ratio should approach one. However, as can be seen from
Figure 2E , the synchrony is sparse (Brunel and Hakim, 2008) in that this ratio is nearly 0.1 after oscillation onset, implying that each oscillation is randomly participated by only 10% of the neurons. Equivalently, each neuron only spikes sparsely and randomly participates in about 10% of the collective oscillations, giving rise to high variability of the ISI. As π ππΌ further increases to bursting onset, this ratio rapidly increases to approaching one, implying that almost all the neurons participate in each oscillation in the bursting state, where the inhibitory feedback is too slow such that the excitatory current can spill over the whole network in the excitatory dominant period. In addition to spiking time variability, we further examine the fluctuations of firing rate in a short time scale (time window
50 ππ ) by Fano Factor (FF) (Teich et al., 1997). As shown in
Figure 2F , the change of FF of the individual neurons (see Materials and Methods for further details) is similar to that of CV (
Figure 2E ), indicating that the individual firing rate fluctuation does not increase with the onset of collective synchrony. In contrast, the change of population FF (
Figure 2F , and see Materials and Methods for details) is similar to the spiking correlation (
Figure 2C ) and temporal variability of firing rate (
Figure 2D ). This trend was further qualitatively confirmed by the field model prediction, as shown in
Figure 2F . Thus, it appeared that the observed collective oscillation is a network property rather than an individual neuron property. Hence, the network collective oscillation activity is compatible with individual neuron irregular spiking. The latter property is the prominent feature of E-I balance (Okun and Lampl, 2009). In summary, strict balance with asynchronous network spiking is predicted by a stable focus whereas loose balance with sparsely synchronous network spiking is predicted by a stable limit circle in the field equations. Both strict and loose balance conditions support the irregular spiking of individual neurons.
Figure 2.
The mean-field theory prediction of the dynamic transition of the spiking network. (A)
The excitatory firing rate in the asynchronous steady-state ( π ππΌ = 1 ππ ) under different strengths of external input and different network sizes. Blue symbols, obtained by the field equations with numerically estimated π πΌ under different input rates π π , fits well with the network simulations (blue dashed lines). Green symbol-lines are obtained by the field equations fixed π πΌ =βπ½ πΌπ2 π π π π π πΌ for π π = 5 π»π§ . Synaptic strengths π½ πΌπ½ are multiplied by 1.41 and 1.15 in the cases of π = 5000 and
π = 7500 , respectively to maintain the usual scale π½ πΌπ½ ~π β1/2 . (B) Deterministic field equations predict that a Hopf bifurcation occurs with the increase of inhibitory decay time π ππΌ at a critical value around π ππΌ β 3 ππ (indicated by a vertical dashed line). The real part and imaginary part divided by of the dominant eigenvalue are denoted by the solid and dashed lines, respectively. Blue circles are peak frequencies of the excitatory firing rate oscillation from network simulation. (C) The PCC and coherence index show the emergence of network oscillation as π ππΌ increases. (D) The CV of population firing rate increases with the increase of π ππΌ . (E) The average CV of ISIs and the ratio between mean firing rate and peak frequency of network oscillation. (F)
The average FF of individual neurons (left axis), and the FF of the network population spiking counts in network and field model (right axis). Time windows for measuring FFs are
50 ππ . Averages are taken across the excitatory population for the measurements in (C-F) . The external input is set as π π = 5 π»π§ in (B-F) . For our mean-field theory of balanced networks, the Hopf bifurcation predicts that a global oscillation emerges in the field equations with oscillation amplitude growing continuously from zero. This picture is similar to a second-order phase transition in statistical physics when oscillation amplitude is taken as the order parameter. Thus, we measure the spiking avalanches (see Materials and Methods for further details) to examine whether this kind of criticality can result in scale-invariant spiking behaviour, i.e. critical avalanches, as observed in experiments (Beggs and Plenz, 2003; Friedman et al., 2012; Bellay et al., 2015; Fontenele et al., 2019).
Figure 3B illustrates the time course of an avalanche.
To more clearly compare the feature of avalanche dynamics, raster plots of the spiking time of excitatory neurons in three typical cases are shown in
Figure 3A . These plots illustrate the asynchronous state (upper panel), the onset of synchrony and oscillation (middle panel) and the developed collective oscillation (lower panel). Power spectrum density analysis of network firing rate is shown in
Figure 3D . The asynchronous state appears without collective oscillation (frequency peak) and the avalanche size and duration distributions are exponential-like, as shown by the blue curves in
Figure 3E, F . It is thus subcritical. Alternatively, the asynchronous state can be considered as very noisy population oscillations that are participated by neurons very sparsely. The onset of synchrony is usually accompanied by a frequency peak of fast oscillation situated within the gamma bands, which is thought to have functional importance in various cognitive processes (Herrmann et al., 2004). This fast oscillation is temporally organized as scale-invariant avalanches, with power-law-like size and duration distributions, as shown by the green curves in
Figure 3E, F and it is thus a critical state. It can be seen from
Figure 2C that the average pairwise correlation is still low in the critical state, since the avalanches at this stage are only randomly participated by a small portion of neurons (
Figure 2E ), which reconciles the coexistence of weak pairwise correlation and strong clustered spiking patterns (Schneidman et al., 2006). The collective oscillation state has more slow frequency peaks and the avalanche size and duration distributions have heavier tails (corresponding to the red curves in
Figure 3E, F ) compared with the critical state, which are features of a supercritical state. This is because avalanches from collective oscillations with specific harmony can produce typical scales, giving rise to heavy tail in the avalanche size or duration distribution. We examine the distance π· between the avalanche size distribution and its fitted power-law distribution (see Materials and Methods for further details).
Figure 3D shows that the network was closest to criticality when poised near the Hopf bifurcation point predicted by the mean-field theory. The existence of power-law distribution is only partial evidence of criticality, as other mechanisms could generate power-law distribution (Yadav et al., 2017). We further examined the scaling relation (Sethna et al., 2001) in the critical state (see
Materials and Methods for further details). The estimated slopes within the truncated ranges in the avalanche size and duration distributions define the critical exponents
π(π)~π βπ and π(π)~π βπΌ . A third exponent is defined as β©πβͺ(π)~π , where β©πβͺ(π) is the average size of avalanches with the same duration π . The third scaling feature can be found in both the subcritical and critical cases ( Figure 3G ), in accordance with previous experimental findings (Friedman et al., 2012). We find that the scaling relation (Sethna et al., 2001) πΌβ1πβ1 = (16) approximately holds at the critical state. The exact value of critical exponents may depend on the details of the system. To demonstrate this, we slightly vary the input strength π π from 4 to 8 Hz and the network still shows significant critical properties at the critical value π ππΌ = 3 ππ . As can be seen in Figure 3H , different exponents keeping the scaling relation Equation (16) can be found by varying the input strength π π . This phenomenon, together with the approximate linear relationship between exponents π and πΌ (inset of Figure 3H ), is in accordance with the results measured in vivo in the primary visual cortexes (V1) of various animals (Fontenele et al., 2019) across a wide range of neural activity states, as well as our experimental data analysis later. The scaling relation expressed by Equation (16) provides additional evidence that the avalanches in the microscopic network occurring near the bifurcation point of the mean-field model possess the properties of criticality. Note that an avalanche can be understood as a temporal propagation of spiking activity in a network. These temporal propagations occur when the excitatory current temporally spill over the network and the avalanche size depends on the strength of the propagations, referring to
Figure 1H, I . In the strict balanced state where excitation is cancelled by inhibition very fast, only small size avalanches can occur and the avalanche distribution is thus subcritical. On the contrary, in the loose balanced state with sparse network oscillations, large avalanches with typical scales can be induced by the temporal domination of excitation and the avalanche distribution can become supercritical. Only near the transition point where the macroscopic dynamic is also noise-driven, avalanches occur with all scales and the avalanche distribution thus can be scale-free. More specifically, on the macroscopic scale, the dynamical process of avalanche corresponds to a noise-induced excursion of the population firing rate. As our measurement of avalanches in the network by individual spiking times requires a fine time scale, and the information in this fine scale is averaged out in the field model, thus, information describing small avalanches vanishes in the field model, which only predicts the global firing rate dynamics.
Although it is still difficult to directly link the microscopic avalanches dynamics to the macroscopic firing rate dynamics, the scale-invariance property of criticality inspires that the properties of burst activities (avalanche) in the macroscopic scale can shed light on the origin of power-law scaling in the microscopic network. To consider the avalanche in the macroscopic scale, one should inspect the fluctuation behavior of a macroscopic signal π₯(π‘) of the network, such as the mean firing rate, etc. An avalanche is a process that starts growing at π₯(0) = π₯ π‘β + π (π β 0 + ) and at the first time it goes back to π₯(π) =π₯ π‘β at time π‘ = π . Here, π₯ π‘β is a threshold above which the avalanche is defined. Such an approach is often used to measure neural avalanche in macroscopic signal such as magnetoencephalography (MEG) (Shriki et al., 2013). The quantity π turns out to be the first-passage time (FPT) (back to π₯ π‘β ) of this process and it defines the duration of an avalanche. The area π = β« π₯(π‘)ππ‘ π0 in this process defines the size of the avalanche. The scale-free behavior at the synchronous transition point may be understood as the general feature of dynamical systems near a bifurcation point when subjected to demographic noise ππ₯ππ‘ ~βπ₯π(π‘) (di Santo et al., 2017). This is because in the critical state of the E-I network, neurons are subjected to Poisson-like noise input with very weak correlation. Thus, according to a Gaussian approximation given by the central limit theorem, the overall fluctuation of population activity density scales with its square root as given by the central limit theorem (Benayoun et al., 2010; di Santo et al., 2018). Noise-driven random walker theory predicts that power-law distribution of the avalanche is a general feature of a dynamical system subjected to such kind of noise when near the Hopf bifurcation (actually, general bifurcations) point. Specifically, although the dynamic form of this macroscopic signal may not be explicit in the field equations, we can heuristically consider a situation that π(π‘) = π₯ β π₯ π‘β obeys an intrinsic dynamic as the normal form of the amplitude dynamics of Hopf bifurcation but driven by noisy force modeled as GWN. The Langevin equation it obeys is ππππ‘ = ππ β π + π(π, π‘) . (19) The first part of Equation (19) is the normal form of the oscillation amplitude dynamics of a supercritical Hopf bifurcation (Marsden and McCracken, 2012), where periodic motion arises when π increases across the bifurcation point π = 0 . The second part π(π, π‘) is the noisy driving force, where the fluctuation scales with square root of the activity. Thus, it has the form π(π, π‘) = β +βππ(π‘) , where β is the mean bias, including the effect from recurrent excitatory, recurrent inhibitory and external inputs. π(π‘) is a standard GWN. The fact is that avalanche dynamics given by the first-passage process of Equation (19) can be mapped to the case of random walks in logarithmic potential (di Santo et al., 2017) by a scaling analysis (details in (di Santo et al., 2017)). Under this approach, the FPT distribution of the avalanche process can be solved by absorbing boundary approach in an analytical way, resulting in π(π = π‘) = (β2π) Ξ( ββ) π‘ exp (β π‘ ) . (20) This explains the power-law distribution relation π(π)~π βπΌ with πΌ = β 2β and other exponents can be obtained by scaling argument (di Santo et al., 2017), although the relation between the avalanche critical exponents in the macroscopic scale and the microscopic scale needs further exploration. In all, irregular microscopic spiking leads to macroscopic fluctuation, which becomes the dominant effect that shapes the dynamic when near the macroscopic bifurcation point. The scale-free avalanche dynamics in the microscopic spiking network at the critical state may be understood as the scaling features of the first-passage dynamics near the Hopf bifurcation point in the macroscopic field model. Figure 3.
Critical dynamics near the onset of collective oscillation. (A)
The raster plots of the spiking time of the excitatory neurons at π ππΌ = 1.2 ππ (upper panel), π ππΌ = 3 ππ (middle panel) and π ππΌ = 4.3 ππ (lower panel). (B) An example of the time course of an avalanche. (C)
The distance between the avalanche size distribution and its best-fitting power-law distribution. (D)
The power spectrum density of the network excitatory firing rate. (E)
The probability density distribution of the avalanche size. (F)
The probability density distribution of the avalanche duration. (G)
The mean avalanche size with respect to a given avalanche duration. For (D-G) , the blue, green and red curves correspond to the cases of π ππΌ = 1.2, 3, 4.3 ππ , respectively. (H) Scaling relation Equation (16) approximately holds for critical states ( π ππΌ = 3 ππ ) with different input strengths , where the dashed line represents Equation (16) exactly holding. The inset shows an approximate linear relationship between exponents π and πΌ . Here, π π = 5 π»π§ is used in (A-D) and the example of (E-G) is π π = 8 π»π§ . In the following, we further analyze public experimental data measuring the in vitro neuronal spiking activity of mouse somatosensory cortex cultures (Ito et al., 2016) (25 data sets in total) to seek for the phenomena in the E-I balanced model. For the experimental data, spiking in the culture clearly shows up-down state transition. Active spiking periods (up-state) and silent periods (down-state) alternate slowly with a frequency of circa o Q density (PSD) in Figure 4B . On the contrary, spikes were too few in the down-state, thus up-states are more likely to represent normal neural activities and closer to the dynamic regime we have studied in the model. Furthermore, the neurons recorded in data exhibit heterogeneity (Ito et al., 2014), i.e. board distribution of individual firing rate, which is different from our modeling results with homogeneous random network topology. However, the spiking of neurons in up-state is more homogeneous, which can be seen from the evener distribution of firing rates of neurons (
Figure 4A ). Taken together, we expect that our modeling results may explain partial properties of the experimental data in up-states. While most of the data sets display heavy tails in the distributions of avalanche size and duration (Supplementary Figure 3 and Supplementary Figure 4), we find that nine of the 25 data sets exhibit both size and duration distributions that correspond to power-laws according to our standards (see Materials and Methods for further details). Among them, seven sets contain critical exponents that approximately satisfy the scaling relation Equation (16) with errors < 0.1. Further details of the analysis results of the data sets are presented in Supplementary Table 1. Estimating the critical exponents of the data sets. These estimated critical exponents πΌ , π and are also close to the ranges found in our model. The statistical results of these seven data sets are shown in Figure 4C, G . The boxplots of the distribution of CV of ISI across neurons in each set are shown in
Figure 4C . Typical values of CV of neuron ISI are around 1~1.5, a range similar to the model prediction, indicating irregular spiking.
Figure 4D to F illustrates the avalanche size and duration statistics of Set 1. Power-law distribution of avalanche size and duration can be observed. For surrogated data, where the inter-spike intervals of neurons were randomly shuffled, power-law distributions cannot be maintained, indicating that the critical properties are intrinsic in the data. As shown in
Figure 4G , the scaling relationship between critical exponents (i.e. Equation (16)) approximately holds for these data sets. For these critical data sets, we also find a linear relationship between exponents π and πΌ (as seen in the inner panel of Figure 4G ) whereby πΌ = 1.36π β 0.37 , while the slope is slightly different from the value in our model with πΌ = 1.21π β 0.154 (inner panel,
Figure 3H ) and a previous study (Fontenele et al., 2019) with πΌ = 1.28π β 0.28 . These results, obtained from experimental data, indicate the phenomena of the model that irregular neuron spiking can collectively organize as scale-invariant critical avalanches are widespread. Figure 4.
Coexistence of irregular spiking and critical avalanches in experimental data sets. (A)
The frequency distribution of the neurons firing rate in up and down states. (B)
The power spectral density of the population firing rate in up and down states. In (A, B) , the examples are shown by an up and a down period of data Set 1. (C)
Boxplots of the CV distribution across neurons in different data sets. (D-F) : critical avalanche properties of data Set 1. (D)
Probability density distribution of the avalanche size. (E)
Probability density distribution of the avalanche duration. (F)
Mean avalanche size with respect to the given avalanche duration. The power-law distributions of size and duration are destroyed by surrogated data with shuffled ISI (denoted by blue curves). (G)
Scaling relations between critical exponents of different data sets, similar to
Figure 3H .
4. Discussion
In summary, we have shown that E-I balanced IF networks with synaptic filtering kinetics can reconcile the coexistence of irregular spiking and collective critical avalanches near a synchronous transition point. The mechanism is unveiled by the macroscopic field equations derived by a novel mean-field theory which effectively capture the network dynamics. Finally, we further show that the phenomenon can be widely observed in the spontaneous spiking activities recorded in vitro in mouse somatosensory cortex cultures. Traditionally, E-I balance (Okun and Lampl, 2009) has been deemed to be the origin of neuron spiking irregularity. Classical mean-field theory (Van Vreeswijk and Sompolinsky, 1996, 1998; Renart et al., 2010) predicts that irregular spiking is an asynchronous state with low spiking time correlation. However, input synchrony (Stevens and Zador, 1998) has also been shown to contribute to the property of irregular spiking. This varied understanding renders a question on how to effectively characterize the dynamic of the E-I balance state. Here, we study E-I networks that incorporates synaptic kinetics and the critical transition from strict to looser balance, where individual neurons continue to spike irregularly during these different dynamic states. Thus, individual irregular spiking is compatible with asynchronous (Renart et al., 2010) or sparse synchrony (Brunel and Hakim, 2008) state (i.e. synchronous inputs). Our theory shows that the balanced state can be a stable fixed point or a stable limit cycle on the macroscopic scale. The former induced by strict balance corresponds to the asynchronous state in accordance with the traditional theory, while the latter induced by loose balance corresponds to collective network oscillation, with critical dynamics during the transition between them. The critical dynamics where neurons are weakly synchronous also provides an explanation of how weak pairwise correlation can induce abundant collective behavior (Schneidman et al., 2006). Our study thus gives an effective characterization that the E-I balanced state can be a stable state of different characteristics (fixed point or limit cycle) in the macroscopic scale.
The criticality of neural systems has been long debated (Wilting and Priesemann, 2019a). Previous experimental and theoretical studies (Beggs and Plenz, 2003; Lombardi et al., 2012; Poil et al., 2012; Yang et al., 2012) have suggested that critical avalanches exist in the E-I balanced state, while many simplified models use unrealistic neural dynamics (such as spreading processes on networks) and the definition of E-I balance is usually ambiguous. Our study of IF network with realistic synaptic kinetics considers realistic irregular asynchronous or synchronous spiking in E-I balanced state, and reveals that that critical avalanches exist near the synchronous transition point between these states. Compared with well-studied simplified physical models, our discovery offers a more biologically plausible explanation of the origin of scale-free dynamics in neural systems, i.e. the dynamics are caused by a potential phase transition indicated by a Hopf bifurcation on the macroscopic scale of E-I balanced neuronal networks. Our theory is thus consistent with the understanding that criticality occurs at the edge of a synchronous transition (di Santo et al., 2018; Fontenele et al., 2019) with intermediate levels of spiking variability (Fontenele et al., 2019), and that critical avalanches can temporally organize as collective oscillations (Gireesh and Plenz, 2008). Significantly, in our model, the synaptic transmission is not instantaneous due to the filtering effect, which renders difficulty in the distinction of the spatially causal relation (Martinello et al., 2017) between successive spiking. The time binning avalanches we measure here are temporally causal (as in experimental measurements) but not necessarily spatially causal in the network. Thus, our estimated critical exponents do not agree with the spatially causal avalanches produced by critical branching processes (Haldeman and Beggs, 2005) (i.e. directed percolation (DP) class), which predicts the classical results π = 1.5 , πΌ = 2 , = 2 . Indeed, critical exponents different from the DP class have been found in previously reported experiments (Palva et al., 2013; Fontenele et al., 2019). Our own analysis of the experimental data ( Figure 4G ) also confirmed the variation of the exponents while maintaining the scaling relationship Equation (16). Furthermore, the linear relationship between exponents π and πΌ , also found in recent studies (Dalla Porta and Copelli, with synaptic kinetics In traditional theory, the firing rate of an IF neuron can be estimated using statistics relating to the first-passage process of membrane potential with diffusion approximation (Burkitt, 2006). Under the assumption of inputs with Gaussian white noise (without auto-correlation), the firing rate statistics can be derived from the first-passage time theory of the Ornstein-Uhlenbeck (OU) process (Ricciardi and Sacerdote, 1979). However, the effect of synaptic filtering will induce memory effects so that the input noise is no longer white, which prevents a complete analytic treatment of the system (for fast decay synapses, the firing rate has been derived asymptotically through a singular perturbation method) (Fourcaud and Brunel, 2002). Self-consistent relation of steady-state firing rate of a network (Shriki et al., 2003; Renart et al., 2004) can be derived by a presumed f-I curve, which expresses the relationship between input current and output rate. A more challenging but also more useful issue that has attracted much recent attention is to find macroscopic transient dynamic descriptions of the neuronal networks. Schaffer et al. (Schaffer et al., 2013) derived the complex firing rate equations of IF networks through the eigenfunction expansion of the Fokker-Planck equation under diffusion approximation.
Deriving rate equations of adaptive nonlinear IF networks has also been studied (Augustin et al., 2017) under some effective approximation of the Fokker-Planck equation. Montbrio et al. (MontbriΓ³ et al., 2015) derived the rate equations for quadratic IF networks using the Lorentzian ansatz. This approach has been generalized to including gap junctions (Laing, 2015) and synaptic filtering kinetics (Dumont and Gutkin, 2019). Schwalger et al (Schwalger et al., 2017) developed a method to derive the stochastic rate equations of adaptive IF networks based on mean-field approximation of the renewal equation. In general, analytical theory works for specific conditions and thus is highly specific and with strong complexity. Furthermore, most proposed theories failed to capture the synchronous transitions induced by the synaptic filtering effect studied in our model. Instead of being theoretically perfect, here we seek for simplicity and effectiveness. The key feature of our framework (see detailed derivations in Materials and Methods) is that the mean-field equations of the macroscopic dynamical variables can be closed by the voltage-dependent mean firing rate (i.e. Equation (9)). While many previous theoretical analyses of E-I network usually require the assumption of low input correlation (VanVreeswijk and Sompolinsky, 1998; Brunel, 2000), the formula here is constructed by counting the number of neuron spiking in a small time window but not required to explicitly consider the correlation of neurons, which allows it to essentially capture the sub and supra threshold microscopic dynamics of a spiking network. Our novel scheme to derive the dynamic equations governing the first-order statistics of the neural network is not completely analytical, since the derivation neglects several factors that shape the network dynamics, including higher order statistics, noise correlation and refractory time. All these neglected factors have been incorporated into the effective parameters, π πΌ , given by Equation (10). Nonetheless, this semi-analytical theory provides a simple yet useful method to analyze dynamic features related to the network firing rate, such as synchrony, criticality and response to dynamic stimuli. Importantly, the theory successfully explains the dynamic origin of the collective oscillation induced by prolonged inhibitory synaptic decay times in the biologically plausible E-I balanced network. The semi-analytical nature of the theory results in the simple form of the macroscopic field equations, which facilitates further generalization. In summary, our theory extends the traditional understanding of E-I balance, the asynchronous state to more realistic dynamics, with the different levels of synchronization that are typically observed in various experiments. It uncovers a possible origin of criticality in neural systems whereby the Hopf transition point in E-I balanced neural circuits can simultaneously reconcile irregular spiking and critical avalanches, a phenomenon that is further observed with in vitro experimental data. Our novel semi-analytical mean-field theory offers a simple yet useful and highly generalizable theoretical tool to analyze the dynamics of integrate-and-fire neuronal networks with realistic synaptic kinetics, building a foundation for the integrated study of neural dynamics on different spatial scales. Thus, it is straightforward that the mean-field analysis introduced here can be used in the study of neural dynamics across a diverse range of topics. It can be easily generalized to include other factors by assuming different macroscopic variables in the field model. Possible extensions of the analysis include studying the effects of multiple neuronal populations and synaptic receptor types, cluster or spatially extended network connections, adaptive behaviors such as short-term plasticity, etc. Thus, our work allows further exploration of the mechanism that determines the role of synaptic kinetics in working memory retrieval (Mongillo et al., 2008), self-organized critical phenomena due to plasticity (Levina et al., 2007; Millman et al., 2010) and spatially causal avalanches or waves (Keane and Gong, 2015; Keane et al., 2018; Gu et al., 2019) arising from competition between Hopf and Turing instability (Huang et al., 2019). As a theory that links microscopic neuronal spiking and macroscopic collective activity that are consistent in several aspects with real neural dynamics with regard to E-I balance and neural criticality, our theory also establish a base to model large-scale brain connectomes (Haimovici et al., 2013; Chaudhuri et al., 2015; Wang et al., 2019), to study large-scale brain networks and information processing with realistic multiscale complex dynamics. The potential applications of our theory listed here will receive the further exploration they deserve in due course. Funding
This work was supported by the Hong Kong Baptist University (HKBU) Strategic Development Fund, the Hong Kong Research Grant Council (GRF12200217), the HKBU Research Committee and Interdisciplinary Research Clusters Matching Scheme 2018/19 (RC-IRCMs/18-19/SCI01), the National Natural Science Foundation of China (Grants No. 11775314, No. 91530320 and No. 11975194) and Key-Area Research and Development Program of Guangzhou (Grants No. 202007030004) . This research was conducted using the resources of the High-Performance Computing Cluster Centre at HKBU, which receives funding from the RGC and the HKBU.
References
Abeles, M. (1991).
Corticonics: Neural circuits of the cerebral cortex . Cambridge University Press. Augustin, M., Ladenbauer, J., Baumann, F., and Obermayer, K. (2017). Low-dimensional spike rate models derived from networks of adaptive integrate-and-fire neurons: comparison and implementation.
PLoS Comput. Biol.
13, e1005545. Beggs, J. M., and Plenz, D. (2003). Neuronal avalanches in neocortical circuits.
J. Neurosci.
23, 11167β11177. Bellay, T., Klaus, A., Seshadri, S., and Plenz, D. (2015). Irregular spiking of pyramidal neurons organizes as scale-invariant neuronal avalanches in the awake state.
Elife
4, e07224. Benayoun, M., Cowan, J. D., van Drongelen, W., and Wallace, E. (2010). Avalanches in a stochastic model of spiking neurons.
PLoS Comput. Biol.
6, e1000846. Brown, E. N., Lydic, R., and Schiff, N. D. (2010). General anesthesia, sleep, and coma.
N. Engl. J. Med.
J. Comput. Neurosci.
8, 183β208. Brunel, N., and Hakim, V. (2008). Sparsely synchronized neuronal oscillations.
Chaos An Interdiscip. J. Nonlinear Sci.
18, 15113. 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. Burkitt, A. N. (2006). A review of the integrate-and-fire neuron model: I. Homogeneous synaptic input.
Biol. Cybern.
95, 1β19. Chaudhuri, R., Knoblauch, K., Gariel, M.-A., Kennedy, H., and Wang, X.-J. (2015). A large-scale circuit mechanism for hierarchical dynamical processing in the primate cortex.
Neuron
88, 419β431. Cohen, M. R., and Kohn, A. (2011). Measuring and interpreting neuronal correlations.
Nat. Neurosci.
14, 811. Dalla Porta, L., and Copelli, M. (2019). Modeling neuronal avalanches and long-range temporal correlations at the emergence of collective oscillations: Continuously varying exponents mimic M/EEG results.
PLOS Comput. Biol.
15, e1006924. Darshan, R., van Vreeswijk, C., and Hansel, D. (2018). Strength of correlations in strongly recurrent neuronal networks.
Phys. Rev. X
8, 31072. Dayan, P., and Abbott, L. F. (2001).
Theoretical neuroscience: computational and mathematical modeling of neural systems . MIT press Cambridge, MA. Denève, S., and Machens, C. K. (2016). Efficient codes and balanced networks.
Nat. Neurosci.
19, 375. di Santo, S., Villegas, P., Burioni, R., and MuΓ±oz, M. A. (2017). Simple unified view of branching process statistics: Random walks in balanced logarithmic potentials.
Phys. Rev. E
95, 32115. di Santo, S., Villegas, P., Burioni, R., and MuΓ±oz, M. A. (2018). Landau--Ginzburg theory of cortex dynamics: Scale-free avalanches emerge at the edge of synchronization.
Proc. Natl. Acad. Sci.
PLoS Comput. Biol.
15, e1007019. Fontenele, A. J., de Vasconcelos, N. A. P., Feliciano, T., Aguiar, L. A. A., Soares-Cunha, C., Coimbra, B., et al. (2019). Criticality between cortical states.
Phys. Rev. Lett.
Neural Comput.
14, 2057β2110. Friedman, N., Ito, S., Brinkman, B. A. W., Shimono, M., DeVille, R. E. L., Dahmen, K. A., et al. (2012). Universal critical dynamics in high resolution neuronal avalanche data.
Phys. Rev. Lett.
Proc. Natl. Acad. Sci.
Scholarpedia
2, 1347. Gu, Y., Qi, Y., and Gong, P. (2019). Rich-club connectivity, diverse population coupling, and dynamical activity patterns emerging from local cortical circuits.
PLoS Comput. Biol.
15, e1006902. Haimovici, A., Tagliazucchi, E., Balenzuela, P., and Chialvo, D. R. (2013). Brain organization into resting state networks emerges at criticality on a model of the human connectome.
Phys. Rev. Lett.
Phys. Rev. Lett.
94, 58101. Herrmann, C. S., Munk, M. H. J., and Engel, A. K. (2004). Cognitive functions of gamma-band activity: memory match and utilization.
Trends Cogn. Sci.
8, 347β355. Holt, G. R., Softky, W. R., Koch, C., and Douglas, R. J. (1996). Comparison of discharge variability in vitro and in vivo in cat visual cortex neurons.
J. Neurophysiol.
75, 1806β1814. Huang, C., Ruff, D. A., Pyle, R., Rosenbaum, R., Cohen, M. R., and Doiron, B. (2019). Circuit models of low-dimensional shared variability in cortical networks.
Neuron Ito, S., Yeh, F.-C., Hiolski, E., Rydygier, P., Gunning, D. E., Hottowy, P., et al. (2014). Large-scale, high-resolution multielectrode-array recording depicts functional network differences of cortical and hippocampal cultures.
PLoS One
9, e105324. Ito, S., Yeh, F.-C., Timme, N. M., Hottowy, P., Litke, A. M., and Beggs, J. M. (2016). Spontaneous spiking activity of hundreds of neurons in mouse somatosensory cortex slice cultures recorded using a dense 512 electrode array.
CRCNS. org . Keane, A., and Gong, P. (2015). Propagating waves can explain irregular neural dynamics.
J. Neurosci.
35, 1591β1605. Keane, A., Henderson, J. A., and Gong, P. (2018). Dynamical patterns underlying response properties of cortical circuits.
J. R. Soc. Interface
15, 20170960. Kinouchi, O., and Copelli, M. (2006). Optimal dynamical range of excitable networks at criticality.
Nat. Phys.
2, 348. Laing, C. R. (2015). Exact neural fields incorporating gap junctions.
SIAM J. Appl. Dyn. Syst.
14, 1899β1929. 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. Levina, A., Herrmann, J. M., and Geisel, T. (2007). Dynamical synapses causing self-organized criticality in neural networks.
Nat. Phys.
3, 857. Litke, A. M., Bezayiff, N., Chichilnisky, E. J., Cunningham, W., Dabrowski, W., Grillo, A. A., et al. (2004). What does the eye tell the brain?: Development of a system for the large-scale recording of retinal output activity.
IEEE Trans. Nucl. Sci.
51, 1434β1440. Litwin-Kumar, A., and Doiron, B. (2012). Slow dynamics and high variability in balanced cortical networks with clustered connections.
Nat. Neurosci.
15, 1498. Lombardi, F., Herrmann, H. J., Perrone-Capano, C., Plenz, D., and De Arcangelis, L. (2012). Balance between excitation and inhibition controls the temporal organization of neuronal avalanches.
Phys. Rev. Lett.
The Hopf bifurcation and its applications . Springer Science & Business Media. Marshall, N., Timme, N. M., Bennett, N., Ripp, M., Lautzenhiser, E., and Beggs, J. M. (2016). Analysis of power laws, shape collapses, and neural complexity: new techniques and matlab support via the NCC toolbox.
Front. Physiol.
7, 250. Martinello, M., Hidalgo, J., Maritan, A., Di Santo, S., Plenz, D., and MuΓ±oz, M. A. (2017). Neutral theory and scale-free neural dynamics.
Phys. Rev. X
7, 41071. Millman, D., Mihalas, S., Kirkwood, A., and Niebur, E. (2010). Self-organized criticality occurs in non-conservative neuronal networks during up states.
Nat. Phys.
6, 801. Mongillo, G., Barak, O., and Tsodyks, M. (2008). Synaptic theory of working memory.
Science (80-. ).
Phys. Rev. X
5, 21028. Mora, T., Deny, S., and Marre, O. (2015). Dynamical criticality in the collective activity of a population of retinal neurons.
Phys. Rev. Lett.
Nat. Neurosci.
11, 535. Okun, M., and Lampl, I. (2009). Balance of excitation and inhibition.
Scholarpedia
4, 7467. Okun, M., Steinmetz, N. A., Cossell, L., Iacaruso, M. F., Ko, H., BarthΓ³, P., et al. (2015). Diverse coupling of neurons to populations in sensory cortex.
Nature
Proc. Natl. Acad. Sci.
J. Neurosci.
32, 9817β9823. Renart, A., Brunel, N., and Wang, X.-J. (2004). Mean-field theory of irregularly spiking neuronal populations and working memory in recurrent cortical networks.
Comput. Neurosci. A Compr. approach , 431β490. Renart, A., De La Rocha, J., Bartho, P., Hollender, L., Parga, N., Reyes, A., et al. (2010). The asynchronous state in cortical circuits.
Science (80-. ).
Biol. Cybern.
35, 1β9. Rosenbaum, R., Smith, M. A., Kohn, A., Rubin, J. E., and Doiron, B. (2017). The spatial structure of correlated neuronal variability.
Nat. Neurosci.
20, 107. Salin, P. A., and Prince, D. A. (1996). Spontaneous GABAA receptor-mediated inhibitory currents in adult rat somatosensory cortex.
J. Neurophysiol.
75, 1573β1588. Schaffer, E. S., Ostojic, S., and Abbott, L. F. (2013). A complex-valued firing-rate model that approximates the dynamics of spiking networks.
PLoS Comput. Biol.
9, e1003301. Schneidman, E., Berry II, M. J., Segev, R., and Bialek, W. (2006). Weak pairwise correlations imply strongly correlated network states in a neural population.
Nature
PLoS Comput. Biol.
13, e1005507. Sethna, J. P., Dahmen, K. A., and Myers, C. R. (2001). Crackling noise.
Nature
J. Comput. Neurosci.
11, 111β119. Shew, W. L., Yang, H., Petermann, T., Roy, R., and Plenz, D. (2009). Neuronal avalanches imply maximum dynamic range in cortical networks at criticality.
J. Neurosci.
29, 15595β15600. Shew, W. L., Yang, H., Yu, S., Roy, R., and Plenz, D. (2011). Information capacity and transmission are maximized in balanced cortical networks with neuronal avalanches.
J. Neurosci.
31, 55β63. Shriki, O., Alstott, J., Carver, F., Holroyd, T., Henson, R. N. A., Smith, M. L., et al. (2013). Neuronal avalanches in the resting MEG of the human brain.
J. Neurosci.
33, 7079β7090. Shriki, O., Hansel, D., and Sompolinsky, H. (2003). Rate models for conductance-based cortical neuronal networks.
Neural Comput.
15, 1809β1841. Shu, Y., Hasenstaub, A., and McCormick, D. A. (2003). Turning on and off recurrent balanced cortical activity.
Nature
J. Neurosci.
13, 334β350. Stevens, C. F., and Zador, A. M. (1998). Input synchrony and the irregular firing of cortical neurons.
Nat. Neurosci.
1, 210. Stringer, C., Pachitariu, M., Steinmetz, N., Carandini, M., and Harris, K. D. (2019a). High-dimensional geometry of population responses in visual cortex.
Nature
Science (80-. ).
JOSA A
14, 529β546. TkaΔik, G., Mora, T., Marre, O., Amodei, D., Palmer, S. E., Berry, M. J., et al. (2015). Thermodynamics and signatures of criticality in a network of neurons.
Proc. Natl. Acad. Sci.
Phys. Rev. E
95, 12413. Van Vreeswijk, C., and Sompolinsky, H. (1996). Chaos in neuronal networks with balanced excitatory and inhibitory activity.
Science (80-. ).
Neural Comput.
10, 1321β1371. Wang, R., Lin, P., Liu, M., Wu, Y., Zhou, T., and Zhou, C. (2019). Hierarchical connectome modes and critical state jointly maximize human brain functional diversity.
Phys. Rev. Lett.
Front. Comput. Neurosci.
5, 30. Williams-Garcia, R. V, Beggs, J. M., and Ortiz, G. (2017). Unveiling causal activity of complex networks.
EPL (Europhysics Lett.
Biophys. J.
12, 1β24. Wilting, J., and Priesemann, V. (2019a). 25 years of criticality in neuroscience: established results, open controversies, novel concepts.
Curr. Opin. Neurobiol.
58, 105β111. Wilting, J., and Priesemann, V. (2019b). Between perfectly critical and fully irregular: a reverberating model captures and predicts cortical spike propagation.
Cereb. Cortex
29, 2759β2770. Xue, M., Atallah, B. V, and Scanziani, M. (2014). Equalizing excitation--inhibition ratios across visual cortical neurons.
Nature
Phys. Rev. E
96, 22215. Yang, D.-P., Zhou, H.-J., and Zhou, C. (2017). Co-emergence of multi-scale cortical activities of irregular firing, oscillations and avalanches achieves cost-efficient information capacity.
PLoS Comput. Biol.
13, e1005384. Yang, H., Shew, W. L., Roy, R., and Plenz, D. (2012). Maximal variability of phase synchrony in cortical networks with neuronal avalanches.
J. Neurosci.
32, 1061β1072. Zhou, F.-M., and Hablitz, J. J. (1998). AMPA receptor-mediated EPSCs in rat neocortical layer II/III interneurons have rapid kinetics.
Brain Res. upplementary Material
Appendix 1. Model Extensions
The field model Equations (11), (12) can be conveniently used to compute the transient firing rate dynamics of the network in response to the time-varying external input. For inhomogeneous Poisson external inputs with time-dependent firing rate π π (π‘) , the constant term π π in Equation (12) can just be modified to the time-dependent term π π (π‘) and the field model can be simulated directly in this way. Note that a good estimation of the effective parameter π πΈ , π πΌ may depend on π π , as the estimation result given by Equation (10) depends on π π . However, we point out that if π π (π‘) does not change very strongly, the parameters π πΈ , π πΌ can be kept constant throughout the change of π π (π‘) once they have been estimated. To demonstrate this, we set a time-varying input π π (π‘) = { π , π‘ β [0, 200)4π , π‘ β [200, 400)2π , π‘ β [400, 500)2π(1 + sin [ π100 (π‘ β 500)]) , π‘ β [500, 800] (S1) with π = 5π»π§ for each neuron as shown in Supplementary Figure S1A . Here, we do not consider the synaptic transmission delay and the synaptic decay times are set as π ππΈ = π ππΌ = 4 ππ , π π =0 ππ , where the network would be in asynchronous dynamics. Note that the external input π π (π‘) contains constant part, discontinuous jumps and continuous changes. Simulations show that the firing rate changes accordingly respondent to the change of input, as can be seen from the raster plot in Supplementary Figure S1B . At the same time, we can simulate the field model Equations (11), (12) with the same time-varying input and with fixed parameters π πΈ , π πΌ used in Figure 2A . As shown in
Supplementary Figure S1C , the field model predicts the changing trend of the average firing rate of the network. It should be noted that this simple scheme ignores some complex nature of the firing rate response properties in the presence of synaptic filtering (Moreno-Bote and Parga, 2004; Ledoux and Brunel, 2011).
The present mean-field theory can be directly generalized to conductance-based (COB) model where the postsynaptic inputs received by each neuron depend on the membrane potential of the neuron. Specifically, we further study a COB model with the dynamic equation Equation (1) replaced by ππ π ππ‘ = π πΌ (π π ) + (π πΈπππ£ β π π )[π πΌπ β πΉ πΈ β π π (π‘ β π ππΈ ) πβπ ππ + (S2) π πΌπΈ β πΉ πΈ β π π (π‘ β π ππΈ ) πβπ ππΈ ] + (π πΌπππ£ β π π )π πΌπΌ β πΉ πΌ β π π (π‘ β π ππΌ ) πβπ ππΌ Here, the reversal potential for excitatory and inhibitory synaptic currents are π πΈπππ£ = 0 ππ and ( ) o Q t
πΌπππ£ = β70 ππ respectively. The synaptic strengths of conductance are set as π πΈπ = 0.025 , π πΌπ = 0.04 , π πΈπΈ = 0.02 , π πΌπΈ = 0.04 , π πΈπΌ = 0.27 , π πΌπΌ = 0.48. Other notations, parameters and settings are the same as the current-based (CUB) case. Similar to the CUB model, the COB model shows emergence of collective oscillation induced by slow inhibition. Such a critical transition can also be predicted by our mean-field theory as a Hopf bifurcation, while the derivation of the field equation is slightly different. In the CUB case, the field equation Equation (12) can be obtained by taking the average β©. βͺ πΌ of the original equation Equation (1) under mean-field assumption. In the COB case, we still can take the average β©. βͺ πΌ of the Equation (S2) under mean-field assumption, but have to proceed with the decoupling assumption that β©π π [π πΌπ β πΉ πΈ β π π (π‘ β π ππΈ ) πβπ ππ + π πΌπΈ β πΉ πΈ β π π (π‘ β π ππΈ ) πβπ ππΈ ]βͺ πΌ β (S3) β©π π βͺ πΌ β©π πΌπ β πΉ πΈ β π π (π‘ β π ππΈ ) πβπ ππ + π πΌπΈ β πΉ πΈ β π π (π‘ β π ππΈ ) πβπ ππΈ βͺ πΌ This is based on the fact that in an E-I balanced network where neurons spike irregularly, one expects that at any given time π‘ , the correlation between the membrane potential and the recurrent E, I conductance input for different neurons is small. As such, we get the field equations ππ πΌ ππ‘ = π πΌ (π πΌ ) + (π πΈπππ£ β π πΌ ) [π πΌπ (π π π π + β π π π π π πΌ π πΌ (π‘)) + π πΌπΈ Ξ¦ πΈ ] + (S4) (π πΌπππ£ β π πΌ )π πΌπΌ Ξ¦ πΌ , πΌ = πΈ, πΌ to replace Equation (12), where Ξ¦ πΌ (π‘) = β©β πΉ πΌ β π π (π‘ β π ππΌ ) πβπ ππΌ βͺ πΈ,πΌ still obeys Equation (11). Thus, Equations (11), (S4) constitute the field equations of the COB model Equation (S2). The sigmoid relation Equation (9) can still be assumed and π πΈ , π πΌ can be estimated in a numerical way through Equation (10). A summarization and comparison between the field equations of CUB model Equation (1) and COB model Equation (S2) is as follows. CUB: { ππ πΌ ππ‘ = π πΌ (π πΌ ) + π½ πΌπ (π π π π + β π π π π π πΌ π πΌ (π‘)) + π½ πΌπΈ Ξ¦ πΈ + π½ πΌπΌ Ξ¦ πΌ (π ππΌ πππ‘ + 1) (π π πππ‘ + 1) Ξ¦ πΌ = π πΌ [1+exp( ππ‘ββππΌ(π‘βπππΌ)ππΌ πβ3 )] , πΌ = πΈ, πΌ (S5) COB: { ππ πΌ ππ‘ = π πΌ (π πΌ ) + (π πΈπππ£ β π πΌ ) [π πΌπ (π π π π + β π π π π π πΌ π πΌ (π‘)) + π πΌπΈ Ξ¦ πΈ ] + (π πΌπππ£ β π πΌ )π πΌπΌ Ξ¦ πΌ (π ππΌ πππ‘ + 1) (π π πππ‘ + 1) Ξ¦ πΌ = π πΌ [1+exp( ππ‘ββππΌ(π‘βπππΌ)ππΌ πβ3 )] , πΌ = πΈ, πΌ (S6) The calculation of the steady-state and its stability analysis at zero transmission delays can be performed in the same way as in CUB model. The qualitative results are similar to the CUB model. There is a critical value π ππΌβ such that when π ππΌ < π ππΌβ the steady-state is a stable focus, corresponding to the asynchronous strict balance state of the network. For π ππΌ > π ππΌβ , the steady-state destabilizes through a supercritical Hopf bifurcation ( Supplementary Figure S1D ), orresponding to the onset of collective oscillation in the network, as shown by the PCC in
Supplementary Figure S1E . The spiking of individual neurons are still irregular, as can be seen from the high CV of ISIs in
Supplementary Figure S1F . Furthermore, near the Hopf bifurcation point, the COB model exhibits critical properties in terms of avalanche dynamics similar to the results of CUB model. Overall, the quality of theoretical prediction in the COB case is worse than the CUB case. Indeed, the COB input would lead membrane potential more bias to a Gaussian distribution (Richardson and Gerstner, 2005), an assumption in our derivation. A complete analytical treatment of COB model is very challenging (Renart et al., 2004). However, the semi-analytical mean-field approach here constitutes an effective description of the macroscopic dynamics of E-I network, which has an advantage that it works for both CUB and COB dynamics.
Supplementary Reference
Ledoux, E., and Brunel, N. (2011). Dynamics of networks of excitatory and inhibitory neurons in response to time-dependent inputs.
Front. Comput. Neurosci.
5, 25. Moreno-Bote, R., and Parga, N. (2004). Role of synaptic filtering on the firing response of simple model neurons.
Phys. Rev. Lett.
92, 28102. Renart, A., Brunel, N., and Wang, X.-J. (2004). Mean-field theory of irregularly spiking neuronal populations and working memory in recurrent cortical networks.
Comput. Neurosci. A Compr. approach , 431β490. Richardson, M. J. E., and Gerstner, W. (2005). Synaptic shot noise and conductance fluctuations affect the membrane voltage with equal significance.
Neural Comput.
17, 923β947.
Appendix 2. Sensitivity of the Critical Points on the Effective Parameters
The mean-field scheme to derive the field equations introduces two effective parameters π πΈ , π πΌ to construct the voltage-dependent firing rate relation Equation (9) and they are the crucial quantities that determine the quality of the scheme. Thus, it is important to know how the theoretically predicted critical point π ππΌβ depends on the choice of π πΈ , π πΌ . Although the critical point in the field model can be thought as a Hopf bifurcation point, the concept of critical point is not decisive in the E-I spiking neuronal network. This can be seen from Figure 2C and
Supplementary Figure S1B which show that the spiking correlation increases in a somewhat continuous way as π ππΌ increases. Furthermore, the critical properties of avalanche shown in Figure 3 are statistical properties so that one can expect that for parameters close to the critical value these properties can still maintain in a statistically significant manner. As a first approximation, the critical point in the E-I spiking network can be defined as the parameter value where the distance of the avalanche size distribution to its best fit power law distribution, π· defined in Methods, is minimal, as shown in Figure 3C . We find that for large network size, if π πΈ , π πΌ are estimated in the numerical way through Equation (10), the critical point in the spiking network, is very close to the Hopf bifurcation point in the field model. We denote π ππΌβ as the Hopf bifurcation point in this βoptimalβ estimation of the parameters π πΈ , π πΌ using Equation (10). We compute the Hopf bifurcation point π ππΌ π»πππ (π πΈ , π πΌ ) predicted by the field model for different values of π πΈ , π πΌ and compare it to the βrealβ critical point (estimated by π ππΌβ ). The difference π ππΌ π»πππ (π πΈ , π πΌ ) β π ππΌβ of the CUB model and the COB model can be seen in Supplementary Figure S2 . From
Supplementary Figure S2 , we can see that once π πΈ , π πΌ are estimated with suitable values, such as by Equation (10), the prediction of the critical synchronous transition point is very precise. We also notice that the bifurcation point predict by the field model seems to mainly depend on the difference π πΈ β π πΌ . Once π πΈ β π πΌ lies on suitable ranges, the predicted critical point will be very close to the βrealβ one. It can also be noticed that in the CUB model the critical point is not sensitive to the values of π πΈ , π πΌ compared with the COB model, as can be seen from Supplementary Figure S2A to C that the difference π ππΌ π»πππ (π πΈ , π πΌ ) β π ππΌβ is still low for a large range of parameter values. On the contrary, the sensitivity in the COB case implies that the COB model has more complicated intrinsic dynamic nature that has to be further explored. Supplementary Figure 1. Results of the generalized models. (A-C)
Network dynamics in response to time-varying external input. (A)
The time-dependent input function π π (π‘) used in simulation. (B) Raster plot of the spiking time of neurons (only part of the N =10000 neurons are shown). The excitatory/ inhibitory neurons are indicated in blue/red. (C) Comparison of the mean firing rate of excitatory population obtained by network simulation and field model simulation. (D-F)
Mean-field theory prediction of the transition from asynchronous to synchronous state in COB model. (D)
Field equations predict that a Hopf bifurcation occurs as the increase of inhibitory decay time π ππΌ at a critical value around π ππΌ β 4.4 ππ . The real and imaginary part (divided by ) of the dominant eigenvalue are given by the solid and dashed lines, respectively. (E) The PCC index shows the emergence of network oscillation as the increase of π ππΌ across the bifurcation point. (F) The CV of ISI at different value of π ππΌ . The parameters in COB model are set as π ππΈ = π ππΌ = 0 ππ , π ππΈ = 2 ππ , π π = 0 ππ and π π = 5 π»π§ . Supplementary Figure 2.
Sensitivity of the predicted critical point on the effective parameters.
The difference value π ππΌ π»πππ (π πΈ , π πΌ ) β π ππΌβ is shown by color for current-based (CUB) model (A-C) and conductance-based (COB) model (D-F) . White dots in (A-F) indicate the positions of the numerical estimation results by Equation (10) and the corresponding estimated critical value π ππΌβ is shown on top of the plots. Black dots in (A-C) indicate the positions of the theoretical estimation results by fixed π πΈ , π πΌ used in Figure 2A in the CUB model. Parameters are π ππΈ = 2 ππ for (A, D) , π ππΈ = 3 ππ for ( B, E ), π ππΈ = 4 ππ for ( C, F ) and π π = 5 π»π§ for all the cases. Supplementary Figure 3.
Avalanche size distributions of the data sets.
According to our standard, the up-state of data Set 1,3,4,5,7,9,13,14,15,17,21,25, plotted by green color, exhibit significant power-law size distribution
π(π)~π βπ . Refer to Supplementary Table 1 for Details.
Supplementary Figure 4.
Avalanche duration distributions of the data sets.
According to our standard, the up-state of data Set 1,3,4,5,6,7,8,9,10,12,13,14,17,18,21,24, plotted by green color, exhibit significant power-law duration distribution
π(π)~π βπΌ . Refer to Supplementary Table 1 for Details.
Supplementary Table 1. Estimating the critical exponents of the data sets.
The number of neurons, length of the up-states, the time bin used in measuring the avalanches, maximum avalanche size and duration, those estimated critical exponents, data ranges after truncations and the p-values in KS test of the fitted power laws in each data set are shown. Here, the avalanche size and duration are in their original linear scale. S e t N o . N e u r o n a m o un t s U p - s t a t e l e n g t h ( m i n s ) b i n s i z e ( m s ) m a x s i z e m a x du r a t i o s i z e r a n g e Ο p v a l u e du r a t i o n r a n g e Ξ± p v a l u e / Ο Ξ½ z ( Ξ± - ) / ( Ο - ) . . [ , ] . . [ , ] . . . . . . . . [ , ] . . [ , ] . . . . . . [ , ] . . [ , ] . . . . . . [ , ] . . [ , ] . . . . . . [ , ] . . . . . [ , ] . . [ , ] . . . . . . [ , ] . . . . [ , ] . . [ , ] . . . . . . [ , ] . . . . . . . [ , ] . . . . . [ , ] . . [ , ] . . . . . . [ , ] . . [ , ] . . . . . . [ , ] . . . . . . [ , ] . . [ , ] . . . . . . [ , ] . . . . . . . . . [ , ] . . [ , ] . . . . . . . . . . [ , ] . . . . . [ , ] . ..