Scale free avalanches in excitatory-inhibitory populations of spiking neurons with conductance based synaptic currents

We investigate spontaneous critical dynamics of excitatory and inhibitory (EI) sparsely connected populations of spiking leaky integrate-and-fire neurons with conductance-based synapses. We use a bottom-up approach to derive a single neuron gain function and a linear Poisson neuron approximation which we use to study mean-field dynamics of the EI population and its bifurcations. In the low firing rate regime, the quiescent state loses stability due to saddle-node or Hopf bifurcations. In particular, at the Bogdanov-Takens (BT) bifurcation point which is the intersection of the Hopf bifurcation and the saddle-node bifurcation lines of the 2D dynamical system, the network shows avalanche dynamics with power-law avalanche size and duration distributions. This matches the characteristics of low firing spontaneous activity in the cortex. By linearizing gain functions and excitatory and inhibitory nullclines, we can approximate the location of the BT bifurcation point. This point in the control parameter phase space corresponds to the internal balance of excitation and inhibition and a slight excess of external excitatory input to the excitatory population. Due to the tight balance of average excitation and inhibition currents, the firing of the individual cells is fluctuation-driven. Around the BT point, the spiking of neurons is a Poisson process and the population average membrane potential of neurons is approximately at the middle of the operating interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[V_{Rest}, V_{th}]$$\end{document}[VRest,Vth]. Moreover, the EI network is close to both oscillatory and active-inactive phase transition regimes. Supplementary Information The online version contains supplementary material available at 10.1007/s10827-022-00838-4.


Diffusion Approximation
The Poisson input in the limit of a high firing rate and small synaptic weights can be approximated by a diffusion process. Suppose, in the time interval [t, t+ dt], N (t, t + dt) excitatory spikes arrive at the cell each with synaptic strength w e . As the spike arrival is a Poisson process with the rate λ the distribution of N (t, t+dt) is Poisson and all the cumulants of the random variable N are equal to λdt. This leads to the following cumulant for I(t, t + dt) = w e N (t, t + dt) : κ 1 = I t,t+dt dt = w e λdt κ 2 =V ar(I t,t+dt ) = w 2 e λdt κ 3 =w 3 e λdt Higher cumulants can be ignored if we assume that w 3 e λ goes to zero in the limit of a high number of afferent inputs. In theory, this can be achieved by assuming weights to scale as w e = W √ k where k is the number of presynaptic neurons. In this case, λ ∼ O(k) and the average excitatory and inhibitory currents are each of order O( √ k), the variance of the current is of O(1) and higher cumulants vanish in the limit of large k. In this case, one can take I(t, t + dt) as a Gaussian random variable with mean and variance given by the above equation. It can also be written as where W t is a Wiener process.
In the conductance based model the input to the cell causes a change in the conductance. As the input is stochastic the conductance is also a stochastic variable which can be written as The second term is the integral of a Wiener process with exponential kernel. For a stochastic process Y t = t − inf f (s)dW s , we can easily verify that: If we consider a stationary and homogeneous Poisson process as the input then g(t) will also attain a stationary probability distribution. Using the above equation, mean and variance of g(t) reach the limits g =τ g 0 w e λ V ar(g) = τ g 2 0 w 2 e λ 2 (5) With the above scaling of the weights, higher cumulants vanish and one can assume that g(t) reaches a stationary Gaussian probability distribution with mean and variance as above. The covariance of this process can be derived by direct multiplication and averaging over the noise terms of g(t) and g(t ) from Equation 5 to reach the following equation when t > t: The same procedure applies to inhibitory currents. Figure 1 shows statistics of conductance g(t) and the validity of diffusion approximation. If we assume that the synaptic time scale τ is very small in comparison to the membrane potential time scale, then one can ignore the cross correlation and consider Y t = g(t) − g as Gaussian white noise. In this limit, This leads to a stochastic differential equation for the membrane potential evolution in the conductance-based model, when v(t) < V th : Here ξ exc (t) and ξ inh (t) are purely random Gaussian processes: The first line of Equation 8 is the deterministic evolution of the potential. When the fixed point of the deterministic term, v inf det = V st = a b as defined by Equation 8 is greater than V th , the effect of fluctuations is marginal and the firing of the neuron is governed by the drift term. However, when V st is below the threshold, fluctuations in the input can result in firing of the neuron.
In the following, we calculate the variance and the mean of the membrane potential and the output firing rate of the neuron in the Gaussian approximation. The Fokker-Planck equation corresponding to Eq.8 in the Itô interpretation is with the boundary condition p(V th , t) = 0. The density current at v = V th is equivalent to the firing rate, and this current is fed back to the equation at v = V r resulting in a discontinuity of the membrane potential derivative.
The stationary probability distribution and firing rate are obtained by solving the equation for V r < v < V th . Together with the normalization requirement V th −∞ p(v)dv + r 0 * t ref = 1 one can solve equation 10. numerically for both the stationary probability distribution and the stationary firing rate.
When V st is sufficiently smaller than V th , i.e., in the low firing rate regime, we can ignore the non-linearity caused by the threshold and write down the evolution of the mean and variance of the membrane potential as follows (See Fig.2) This leads to the stationary value for average and variance of the membrane voltage In the low firing regime, the stationary probability distribution can be approximated as (see Fig.3): The stationary firing rate is derived from equation 10 by plugging in the Gaussian approximation for the stationary potential probability density P (V, t → ∞) = N ( V , σ V (st) ): We can improve the approximation for the potential distribution and the firing rate by considering the autocorrelation in the conductance. When τ is not negligible but small, we can use the τ expansion method to account for first-order corrections to the Fokker-Planck equation (section 1.1.1).

Tau expansion
The equation where η(t) is colored Gaussian noise with correlation corresponds to the following Fokker-Planck equation derived by expansion with respect to τ : In our case H e = τ 2 g 2 0 w 2 e λ e 2 = D e 2 , resulting in: where shows that these corrections lead to a better approximation of the stationary membrane potential variance in a low firing rate regime. The stationary variance of the membrane potential decreases at higher values of τ syn ( Figures 4A and 4B) . Higher rates in excitatory and inhibitory input would lead to lower stationary variance and lower rates(section 1.1.2).

Neuron response in higher values of input rates and synaptic time constant
As it can be seen from Fig.2 , the stationary standard deviation of the membrane potential in the case of Poisson input does not show high sensitivity to input rates when the stationary mean potential is fairly away from the threshold. In Fig.4, increasing the excitatory input rate by 60% causes a 2% increase of σ(V ) st on average for different values of v st . This can be predicted from   Fig. 4: (A1-A2) Stationary variance of the membrane potential and the output firing rates, when the excitatory rate (x-axis) and inhibitory rates are balanced so that the average membrane potential is at −57mv for three different values of τ syn = [1 (Blue), 3 (Red), 5 (Black)] ms. (B1-B2) Same with average membrane potential tuned at the threshold value −50mv. (C) Stationary potential variance for a neuron receiving two different excitatory input rates, 1000Hz(red) and 1600Hz (black), and corresponding inhibitory input, which places the average potential at a specific value shown in the x-axis. Solid lines are simulation results, dashed lines correspond to low firing regimes with tau approximation (Equation 13), and fine dashed lines show Gaussian approximation. (τ syn = 2.5ms) equation 11 as both the input noise and drift term in the numerator and the denominator depend linearly on input rates. Suppose v st is fixed for a set of inhibitory and excitatory rates, which means there is a linear relation of the form ρ E = κρ I + c originating from the condition on the fixed stationary average membrane potential. In the limit of high rates, the stationary variance approaches a constant value: when α γ > β η , the variance of the stationary membrane potential (and hence also the output rate) increases by proportional increase of both inhibitory and excitatory rates and reaches a constant value α γ . This condition translates into We remind the reader that the Gaussian approximation is only legitimate in the case of small τ and for a low firing rate regime. We want to consider cases in which these two conditions are not satisfied. Firstly, at higher values of τ syn and at the stationary values of the average membrane potential lower than the threshold (low firing regime), there is an inversion of the mentioned scenario. In this case, at sufficiently high values of the input rates, conditioned on constant average membrane potential, the variance of the membrane potential and accordingly the output rate decrease (Fig.4A). This is due to the filtering effect of the input by the gradual decay of the synaptic conductances. Using the method of τ -expansion of section 1.1.1, we will consider autocorrelation in synaptic conductance and partially take into account this effect. The filtering of the high-frequency signal by the slow conductances is a mechanism of gain control. In general, the decay and the rise time of the inhibitory synapses are longer than the excitatory ones which would highlight the inhibitory input as its overall strength increases through temporal persistence. In addition, the voltage-dependent inhibitory current is higher at higher values of the membrane potential. The balanced average membrane potential is somewhere in the mid-range. Longer synaptic decay time constant, higher synaptic strength, delay, and potential dependence of the inhibitory synapses increase the overall inhibition strength and compensate for the smaller number of them in comparison with excitatory synapses. It should be noted that output rates on the constant voltage level line with balanced input at V = −57mv vary linearly with the input rates at moderate rate values corresponding to the low firing rate regime (Fig.4A).
On the other hand, when the stationary average membrane potential is located right at the threshold value, in conflict with the low firing regime assumption, equation 14 does not hold and higher rates of balanced input lead to a higher output rate independent of the value of τ syn (Fig.4B). Moreover, the output rate varies like √ I in . For analyzing the low firing rate, we could linearize the output rate around the midpoint of the neuron potential range, i.e., at V = −57mv.

Condition for Poisson output
In the model of conductance based synaptic currents, the output spike train is Poisson in the regime of fluctuation driven spiking when the stationary average potential is away from the threshold. To analyze the condition for Poisson firing, we investigate the case where the stationary membrane potential is below firing threshold. When excitatory and inhibitory currents are matched in a way that the stationary membrane potential V st is close to the firing threshold V th , the approximations in the previous section are no longer valid. In the limit of a high number of affarents the drift term would set the average membrane potential to its stationary value in a short time. Here we calculate the mean firing time and variance of it based on the approximation that fluctuations in the input close to the stationary average potential are weakly dependent on the voltage level.
Therefore, our problem is reduced to the well known problem of the first passage time of a Brownian particle evolving as dx dt = −kx + ξ(t) to reach the threshold a. ξ is white noise with variance σ. We want to approximate results for first passage time moments when the stationary membrane potential is close to the threshold. The goal here is to obtain approximate analytical results for firing rate and CV of the spike time interval to identify conditions for output Poisson firing and linearization of the output rate.
We apply the Siegert formula (Siegert (1951)) for the first passage time in our case. Let us take . Thus, the random variable x evolves as dx dt = −bx + σ(x)ξ(t). where ξ is a white noise with unit variance and For the case of Poisson current to the cell, if the mean value is close to the threshold, the magnitude of fluctuations does not vary much in the interval [V rest , V th ] , therefore in the following we neglect dependence of σ on x.

First passage time statistics
To determine first passage time(FPT) statistics let us first review the recursive formula for FPT moments as discussed in Siegert (1951). Suppose the stochastic process X t with conditional probability density P (x, t | x 0 ) satisfies the Fokker-Planck equation with initial and boundary conditions P (x, t 0 | x 0 ) = δ(x − x 0 ) and P (±∞, t | x 0 ) = 0, respectively. The first passage time probability density ρ(Θ | t, x 0 ) of the stochastic process following the above Fokker-Planck equation satisfies Based on this equation we can write a recursion formula for moments of FPT with diffusion strength B(y) and drift term A(y): with t 0 = 1, and W (k) is the stationary probability distribution In particular, for the first moment we have For the case of Langevin equation in section 1.2, writing x th and using equations 17 and 18, the average first passage time can be written as where erf (.) is the Gauss error function.
Taking y = x th b σ 2 and θ = Θ b σ 2 and using the following series expansion for small y y 0 e z 2 dz = y + y 3 3 + y 5 10 ...
we arrive at the following approximation for the mean passage time, It can be seen from equation 14 that the factor b σ 2 in large balanced rates of excitatory and inhibitory inputs asymptotically approaches a constant value. Therefore, we can approximate κ to be very weakly dependent on input rates and take it as a constant factor.
we can write down the average rate of firing, ρ = 1 , in the case that the average stationary potential is close to threshold as In Figure 5, we have plotted this rate approximation which causes only a very small error when the average potential is near threshold. If v = a b is constant for balanced inhibitory and excitatory input rates, then the output rate of the neuron would be linearly proportional to the input rates via the factor b. Moreover the output rate would decrease as 1 x th with the increase of the distance of the stationary average potential from the threshold. Next, we want to investigate the variance of first passage time. From the recursion formula 16 we have Fig. 5: Firing rate approximation by Gaussian assumption of Eq.12 in section 1 (blue) and by near-threshold high firing assumption(red) of Eq.21 are compared with simulation results (black curve). Here, we fix the inhibitory input and select a value of the excitatory input rate that leads to a specific mean stationary membrane potential shown on the x-axis.
After some straightforward calculations we arrive at where functions φ and ψ are multi-variable integrals containing powers of erf and e x 2 in the integrand with series expansion : Rewriting the first term in brackets in terms of first passage time and considering just terms linear in x th where the constant C, coming from the integral, is negative. From the recursion formula (Eq.16)and Eq. 22,we end up with the following expression for the CV of the time interval between spikes: where C is a negative constant (see Eq.22 ) and monotonically goes to zero as x th increases. In the limit of large x th , the second and the third term both go to zero and CV approaches 1. However, in the near threshold approximation the maximum of the third term occurs where CV is approaching 1. Expanding in powers of x th , we arrive at As shown in section 1.1.2 (Eq.14), σ √ b reaches a constant value for high input rates. This can be used to determine the value of V st that leads to maximal CV. At V P := V opt st ≈ −0.56mv, the CV for different input rates has a maximum independently of the rate values.
2 Bifurcation analysis of a network with Sigmoid gain function

Location of the BT point
Here we want to approximate local bifurcation lines in the vontrol parameter space of a neuronal population with sigmoid gain function. At equilibrium, the population rates satisfy: and the equation for the inhibitory nullcline can be written as : The condition on zero trace T r(J) = J 11 + J 22 = 0 of the Jacobian (Eq.31 main text) parameterized by the inhibitory nullcline curve(Eq.26) determines the value for w EE at which a Hopf bifurcation can occur. Ignoring terms related to α, which are relatively small, from equating the trace to zero, we have . For each (ρ E , ρ I ) point on the inhibitory nullcline (Eq.26), the above equation gives a value of w EE which sets the trace of the Jacobian to zero at this point. The second equation in Eq.25 which corresponds to the excitatory nullcline determines ρ E Ext parameterized by ρ Inh . Next, we should check the positivity of the determinant to sketch the Hopf bifurcation line in the w EE − ρ E Ext plane. Neglecting non-linearities caused by α, the determinant of the Jacobian conditioned on zero trace is At extremely low values of the rates (near zero), the determinant is negative because of the −1 in the above formula. Conditioned on a sufficient amount of inhibitory feedback strength, which is proportional to |c IE c EI |, the determinant becomes positive at some point and the low fixed point loses stability through a Hopf bifurcation. A point where both determinant and trace of J are zero, is called a Bogdanov-Takens (BT) bifurcation point.
Inserting ρ E (1 − ρ E ρ max ) from det(J) |T r(J)=0 = 0 into the denominator of the formula for w H EE and introducing the parameter γ = ρ I (1 − ρ I ρ max ) , at the BT point in the low rate regime: When | c II γ | is at a moderate value, i.e. sufficiently greater than one : If we take number of connections to satisfy On the semi-linear part of the inhibitory nullcline we have an approximate relation between rates of the form The function γ(ρ I ) has a maximum at ρ max 2 . With this in mind, the condition for positive determinant at a potential Hopf bifurcation fixed point at lower rates of the linear regime is Ascending on the inhibitory nullcline, we reach the nonlinear high branch of the curve, where a linear relation between rates no longer holds and the second derivative of ρ * E (ρ I ) will increase. In addition, γ(ρ inh ) decreases towards zero. Taking into account these two facts, on the high branch det(J) L |T r(J)=0 decreases and passes through zero at another Bogdanov Takens point at high rate values. If inhibitory feedback is not strong enough, the conditions for Hopf bifurcation are not satisfied.
To sketch the saddle-node bifurcation line we should look at solutions to det(J) = 0. Inserting ρ E (ρ I ) from Equation 26 into det(J), for each point on the inhibitory nullcline, there exists some w EE for which det(J) = 0. The only condition to check is w EE > 0. Again the condition that the excitatory nullcline intersects the inhibitory one at the fixed point determines ρ Ext . Along the semi-linear section of the nullcline, the condition det(J) = 0 translates into the alignment of the slopes of the linearized nullclines. Therefore, along this section w EE varies very little.

Normal form of BT at the low firing rate regime
Consider the system in the vicinity of the BT bifurcation point Suppose that at µ = 0, the system has a fixed point at z 0 with a Jacobian with a zero eigenvalue of multiplicity two. At the BT point, there exist two generalized eigenvectors q 0 and q 1 such that Also for J T we select vectors p 0,1 by J T p 1 = 0, J T p 0 = p 1 with normalization p 0 , q 0 = p 1 , q 1 = 1 p 0 , q 1 = p 1 , q 0 = 0 By linear change of coordinates with transformation matrix T = (q 0 , q 1 ) , i.e., x = T z , our system can be written as Introducing a direction-preserving time reparametrization and smooth invertible parameter changes, we can tranform the system to the normal form: y 1 dτ = y 2 y 2 dτ = 1 + 2 y 1 + a 2 y 2 1 + b 2 y 1 y 2 + O(y 1 y 2 3 ) 1,2 (µ) are transformed bifurcation parameters , a 2 = g xx /2 and b 2 = g xy +f xx and dt = (1 + Θx 1 )dτ . The y coordinates relate to the original z = T −1 x coordinates via : u 1 = x 1 , y 1 = u 1 u 2 =ẋ 1 , y 2 = u 2 + Θu 1 u 2 (36) where Θ = g yy + 2f xy . Fixed points of the normal form of Equation 35 are (y 1 , y 2 ) = (± ( − 1 a 2 , 0). Taking a 2 > 0, when 1 < 0, there exist two fixed points, with Jacobian   0 1 ±2 − 1 a 2 2 ± b 2 − 1 a 2   y + 0 is a saddle in 1 < 0 for all 2 , while y − 0 is a sink for 2 < b 2 − 1 a 2 and a source for 2 > b 2 − 1 a 2 . When b 2 > 0, then the line 2 = b 2 − 1 a 2 is a subcritical Hopf bifurcation and when b 2 > 0 the same line is a supercritical Hopf bifurcation. To summarize, defining σ = sgn(a 2 * b 2 ), if σ is negative then a stable limit cycle appears and the Hopf bifurcation is supercritical (Fig.??B1), but if σ is positive, we have a sub-critical Hopf bifurcation (Fig.??B2). As shown in Figure ??, near the BT point apart from local bifurcations, i.e., Hopf and saddle-node, there is saddle-node separatrix loop bifurcation which annihilates the stable or unstable limit cycle that is produced by a super-or sub-critical Hopf bifurcation, respectively. By writing the normal form we can analyze the type of BT from explicit linearization. For the case of J BT in equation ??, generalized eigenvectors are q 0 = 1 α/β , q 1 = 1 (α − 1)/β , p 1 = (1, −β/α) and p 0 = (1/α − 1, β/α). Therefore, new parameterized coordinates are x 1 = E/β − α/β(E − I) and x 2 = E − I. The normal form parameters a 2 and b 2 are a 2 = 1 2 < p 1 , F (q 0 , q 0 ) > b 2 =< p 0 , F (q 0 , q 0 ) > + < p 1 , F (q 0 , q 1 ) > where F 1 (q 0 , q i ) = l,k f lk q l 0 q k i and F 2 (q 0 , q i ) = l,k g lk q l 0 q k i . Using the logistic gain function, second derivatives in E − I coordinates can be written as g EE ∝ W 2 IE g I (1 − g I )(1 − 2g I ) g EI ∝ −|W II ||W IE |g I (1 − g I )(1 − 2g I ) where g EI < 0, f EE > 0 and g EE > 0. After straightforward but lengthy calculation we can confirm sgn(σ) = sgn(a 2 )sgn(b 2 ) < 0 using |W II |g I (1 − g I ) = W EE f E (1 − f E ) and β α = c EI c EE = c II c IE at the lower BT point.