Self-organized criticality in a mesoscopic model of excitatory-inhibitory neuronal populations by short-term and long-term synaptic plasticity

Dynamics of an interconnected population of excitatory and inhibitory spiking neurons wandering around a Bogdanov-Takens (BT) bifurcation point can generate the observed scale-free avalanches at the population level and the highly variable spike patterns of individual neurons. These characteristics match experimental findings for spontaneous intrinsic activity in the brain. In this paper, we address the mechanisms causing the system to get and remain near this BT point. We propose an effective stochastic neural field model which captures the dynamics of the mean-field model. We show how the network tunes itself through local long-term synaptic plasticity by STDP and short-term synaptic depression to be close to this bifurcation point. The mesoscopic model that we derive matches the directed percolation model at the absorbing state phase transition.

1. Introduction

Neural networks as complex dynamical systems with many degrees of freedom varying over different time scales can be seen as a self-tuning system that attains a dynamical regime where the system can carry out its task. On the other hand, spontaneous intrinsic activity of cortical neural assemblies in absence of any information processing task can be perceived as a substrate for the neural dynamics which can give us insights into the preferred dynamical regime and the goal of self-organization processes. Dynamic and functional characteristics of spontaneous activity are connected to the structural architecture of the brain as well as the ongoing self-organization process. Experimental findings on different temporal and spatial resolutions highlight the scale-free characteristic of spontaneous activity. When recorded by coarse-grained methods like EEG and MEG, spontaneous brain activity shows nested oscillations with a power spectrum that indicates scale-free properties, i.e., P(f) ∝ 1/f^β (Linkenkaer-Hansen et al., 2001; Miller et al., 2009; Hardstone et al., 2012). Microelectrode recordings of smaller cortical regions show activity in the form of avalanches with power-law distributions of size and duration in different setups such as cultured slices of rat cortex (Beggs and Plenz, 2003), awake monkeys (Petermann et al., 2009), in cerebral cortex and hippocampus of anesthetized, asleep, and awake rats (Ribeiro et al., 2010), and the visual cortex of an anesthetized cat (Hahn et al., 2010).

One explanation for the scale-free characteristics is that cortical networks operate near a critical second-order phase transition (Chialvo, 2004; Tagliazucchi and Chialvo, 2013). On one hand, close to the edge of an active-inactive phase transition, local populations of neurons in the cortex would be in an idle state ready for processing information, but at the same time away from overactivation. On the other hand, close to the phase transition, an order parameter, a macroscopic state different from the inactive state comes into existence. The emergence of a macroscopic mode of activity acts as the coordinator of individual neurons which are enumerated in quantity and prone to various kinds of noise to produce a cooperative large-scale activity (Chialvo, 2010).

The type of the phase transition and the dynamical regime that produces the aforementioned characteristics of avalanches have been studied in different neuronal models. Active-inactive phase transition in a purely excitatory population of neurons (Levina et al., 2007, 2009) and synchronization phase transition in an excitatory population coupled with short-term depression of synaptic resources (di Santo et al., 2018) have been proposed as the origin of avalanche dynamics. The main hypothesis is that system resides near the bifurcation point of the quiescent state. In an Inh.-Exc. network, Benayoun et al. proposed a stochastic model of spiking neurons which matches the Wilson-Cowan mean field in the limit of infinite system size that shows scale-free avalanches in the balanced state in which the sum of excitation and inhibition is much larger than the net difference between them (Benayoun et al., 2010). Under symmetry conditions on weights, the Jacobian has negative eigenvalues close to zero in the balanced state. de Candia et al. (2021) showed that this stochastic model exhibits second-order phase transition with scale-free avalanches and the interaction of noise and nonlinearity is the origin of this behavior. In this model, the Poisson firing of the neurons is preassumed and symmetric synaptic connections and O(N⁻¹) scaling of weights is required for applying the linear noise approximation. Furthermore, the origin of the scale-free behavior and the bifurcation diagram of the model in a wider regime of parameters has not been studied. To investigate further the emergence of avalanches in the EI network, in Ehsani and Jost (2022), we used a bottom-up approach to obtain a single neuron's gain function subjected to fluctuating conductance-based synaptic currents and studied the linear Poisson firing regime for EI homogenous sparse population using Fokker-Planck formalism. We showed that the Bogdanov-Takens bifurcation point of the mean-field equations for dynamics of a sparse homogenous excitatory and inhibitory population of spiking neurons is the operating point of the system producing the characteristic spontaneous activity in the form of scale-free avalanches and Poisson firing of neurons. In this regime, the system is close to both the saddle-node bifurcation point at the low firing rate regime and Hopf-bifurcation of the quiescent fixed point. Tight temporal balance of inhibition and excitation at this state and Poisson firing of neurons is the origin of critical behavior. By mapping single population dynamics to a branching process, we have explained the emergence of power law exponents in the spiking neural network which coincides with critical exponents of the branching process.

To tune the system at the critical point, many modeling approaches and adaptive mechanisms have been suggested during the decades of research on the critical brain hypothesis. The self-organizing principle is mainly based on the model of neuronal dynamics and the choice of the control parameter. In the excitatory neuronal population, a SOC model that attracted much attention is introduced by Levina et al. (2007, 2009). In their model, short-term depression of excitatory synapses coupled with internal neuronal dynamics self-tunes the system at the edge of the active-inactive phase transition. In addition to self-organization by short-term depression in synapses which is also used in Peng and Beggs (2013) and di Santo et al. (2018), self-organization by other control parameters like degree of connectivity or synaptic strength has been studied (Bornholdt and Roehl, 2003; Rybarsch and Bornholdt, 2014). In Brochini et al. (2016), self-organization in stochastic spiking neuron model by short-term plasticity of the gain function instead of synaptic weights is introduced. Meisel and Gross (2009) introduced a self-organizing excitatory neural network by STDP. Scarpetta and Candia (2013) studied transient replay of stored patterns by STDP in a balanced network at a macroscopic phase transition induced by noise and also observed scale free avalanches and studied their statistics near the noise induced phase transition (Scarpetta et al., 2013).

Here, we consider the self-tuning of the system at the BT critical point of the EI interconnected network. The self-organizing parameter in our network is the balance of opposing forces resulting from the activities of inhibitory and excitatory populations, and the self-organizing mechanisms are long-term synaptic plasticity through the mechanism of Spike Timing Dependent Plasticity (STDP) and homeostatic short-term depression of the synapses. The former tunes the overall strength of excitatory and inhibitory pathways to be close to a balanced regime of these currents and the latter, which is based on the finite amount of resources in brain areas, acts as an adaptive mechanism that tunes micro populations of neurons subjected to fluctuating external inputs to attain the balance in a wider range of external input strengths. The importance of both types of synaptic plasticity in the emergence and tunning of the critical state has been also mentioned in a review article by Zeraati et al. (2021). Here, we show this phenomenon in the context of EI networks with long-term synaptic plasticity in all types of synapses which has not been studied previously.

For analytical analysis of STDP on average weight connections, we use the inhomogenous Poisson neuron assumption that has been studied in Kistler and van Hemmen (2000) and Burkitt et al. (2007). Under general conditions on inhibitory and excitatory STDP kernels, i.e., negative integral of the excitatory and positive integral of the inhibitory STDP kernels, learning results in a balanced internal state. This condition on kernels leads to stabilization of rates as also discussed in Kempter et al. (2001). We use the model of Tsodyks and Markram (1997) for short-term depression of the excitatory synapses which has been studied vastly (refer to for example Kistler and van Hemmen, 1999).

Using the Poisson firing assumption, we propose a microscopic Markovian model which captures the internal fluctuations in the network due to the finite size and matches the macroscopic mean-field equation by coarse-graining. Near the critical point, a phenomenological mesoscopic model for excitatory and inhibitory fields of activity is possible due to the time scale separation of slowly changing variables and fast degrees of freedom. We will show that the mesoscopic model corresponding to the neural field model near the local Bogdanov-Takens bifurcation point matches the Langevin description of the directed percolation process.

2. Materials and methods

2.1. Neuron model

We use an integrate and fire neuron model in which the change in the membrane voltage of the neuron receiving time dependent synaptic current i(t) follows

$\begin{array}{ll}C\frac{dv(t)}{dt}=g_{Leak}(v_{Leak}-v(t))+i(t)& (1)\end{array}$

for v(t) < v_th. When the membrane voltage reaches v_th = −50mv, the neuron spikes and immediately its membrane voltage resets to v_rest which is equal to v_Leak = −65mv.

In the following, we want to concentrate on a model with just one type of inhibitory and one type of excitatory synapses, which can be seen as the average effect of the two types of synapses. We can write the synaptic inhibitory and excitatory current as

$\begin{array}{ll}\begin{array}{c}\hfill i(t)=g_{inh}(t)*(V_{Rinh}-v(t))+g_{exc}(t)*(V_{Rexc}-v(t))\end{array}& (2)\end{array}$

V_Rinh and V_Rexc are the reverse potentials of excitatory and inhibitory ion channels and based on experimental studies we choose values of −80 and 0mv for them respectively. g_Inh(t) and g_Exc(t) are the conductances of inhibitory and excitatory ion channels. These conductances are changed by the inhibitory and excitatory input to the cell. Each spike of a presynaptic inhibitory or excitatory neuron j to a postsynaptic neuron k that is received by k at time t₀ will change the inhibitory or excitatory ion channel conductance of the postsynaptic neuron for t > t₀ according to:

$\begin{array}{l}{g}_{Inh}^k(t)=w_{kj}*g_0^{inh}*exp(-\frac{t-t_0}{{\tau }_{syn}^{inh}})\\ g_{Exc}^k(t)=w_{kj}*g_0^{exc}*exp(-\frac{t-t_0}{{\tau }_{syn}^{exc}})& (3)\end{array}$

Here, we assume that the rise time of synaptic conductances is very small compared to other time scales in the model, and therefore, we modeled the synaptic current by a decay term with synaptic decay time constant τ_syn which we assume to be the same value of 5ms for both inhibitory and excitatory synapses.

2.2. Network architecture

In the remainder of this study, in the simulation, we consider a population of $N_{Exc}=2*10^4$ and N_Inh = 0.25 * N_Exc inhibitory spiking neurons with conductance-based currents introduced in this section. Each excitatory neuron in the population is randomly connected to $k_{EE}=\frac{N_{Exc}}{100}=200$ excitatory and $k_{EI}=\frac{k_{EE}}{4}$ inhibitory neurons and each inhibitory neuron is connected to k_IE = k_EE and $k_{II}=\frac{k_{EE}}{4}$ excitatory and inhibitory neurons, respectively. The weights of excitatory synaptic connections are in a range that 10 − 20 synchronous excitatory spikes suffice to depolarize the target neuron to the level of its firing threshold when it is initially at rest at the time of input arrival. Weights are being drawn from a log-normal probability density with low variance. Therefore, approximately $O(\sqrt{k_{EE}})$ spikes are adequate for firing. Assuming homogeneity in the population as we have discussed in the introduction we can build a mean-field equation for the excitatory and inhibitory population in this sparse network, assuming each neuron receives input with the same statistics.

2.3. Avalanche regime of activity as desired operating point of the system

In Ehsani and Jost (2022), we have investigated the dynamics of excitatory and inhibitory (EI) sparsely connected populations of spiking leaky integrate neurons with conductance-based synapses. We have seen that close to the Bogdanov-Takens bifurcation point of the mean field equation, the output firing of the population is in the form of avalanches with scale free size and duration distribution. This matches the characteristics of low firing spontaneous activity in the cortex. By linearizing gain functions and excitatory and inhibitory nullclines, we approximated 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 [V_rest, V_th]. Moreover, the EI network is close to both oscillatory and active-inactive phase transition regimes.

At equilibrium, population rates satisfy a system of fixed point equations of the form:

$\begin{array}{l}{\rho }_I=g^I({\rho }_I,c_{EI}{\rho }_E+c_{II}{\rho }_I+d{\rho }_{Ext}^I)-z_0\\ {\rho }_E=g^E({\rho }_I,c_{EE}{\rho }_E+c_{EI}{\rho }_I+d{\rho }_{Ext}^E)-z_0& (4)\end{array}$

where c_xy = ck_xyw_xy(V_Ry − 〈V_x〉). k_xy is the average number of synaptic connections between neurons of population y to neurons in population x with an average strength of w_xy. z₀ is a constant that depends on V_rest, V_th, the maximal rates and the SD of the input. V_Ry is the reverse potential level of a neuron of type y and 〈V_x〉 is the average potential level of neurons in population x that can be written in fluctuation driven the firing regime as:

$\begin{array}{ll}<V_x>=\frac{g_L{V}_L+g_{exc}^0{w}_{xE}{\rho }_E\tau V_{Rexc}+g_{inh}^0{w}_{xI}{\rho }_I\tau V_{Rinh}}{g_L+g_{exc}^0{w}_{xE}{\rho }_E\tau +g_{inh}^0{w}_{xI}{\rho }_I\tau }& (5)\end{array}$

τ is the synaptic current decay time constant, g_L, g_exc and g_inh are the baseline conductances of leaky, excitatory and inhibitory ion gates, respectively. We took w_EE and ${\rho }_{Ext}^E$ as control parameters and analyzed solutions to Equation (4). By substituting nonlinear gain functions with their corresponding linearization in the Poisson firing regime, we showed that the low BT point is located close to the matching condition for the y-intercept and the slopes of the linearized nullclines, which are written as:

$\begin{array}{l}{c}_{EE}^{*}=\frac{c_{EI}{c}_{EI}}{c_{II}}\\ {\rho }_{Ext}^{E^{*}}=\frac{c_{EE}^{*}}{c_{IE}}({\rho }_{Ext}^I-d)+d& (6)\end{array}$

where d is a constant equal to $\frac{g_L(V_{rest}-V_{th})}{\tau *g_{exc}^0*(V_{th}-V_{R_{exc}})}$.

Figure 1 shows the location of the BT point on the local bifurcation diagram and matching condition of Equation 6. Nullclines of Equation 4 near the BT point are depicted in Figure 2. Close to the BT point, the volume of the basin of attraction of the quiescent state is small, and internal noise can make the system escape from it. Activity grows and decays back to the quiescent state along heteroclinic orbits connecting the two saddle points in the low firing rate regime which coincides with the slow manifold of the fixed points. Along this slow manifold, there is a tight temporal balance of excitation and inhibition in the forms of avalanches of highly variable sizes. This results in balanced currents at each cell leading to Poisson firing of individual neurons. At the population level, the temporal balance of excitation and inhibition cancels out the average current to the cells while the fluctuation in the current which is proportional to the average rates leads to a branching factor close to unity for both excitatory and inhibitory populations. The mapping to the branching process on a single EI network explains the emergence of power law exponents in the spiking neural network which agrees with the mean field branching process.

FIGURE 1

Figure 1. Zoom in on the local bifurcation diagram at low firing rates and the corresponding regimes of phase space with different numbers of fixed points. The dashed line is the condition on the equal slope of linearized nullclines and the semi-dashed line is the condition on equal y-intercepts. BT point (black dot) is close to the intersection of these lines. In the labeling of regions (Q) denotes quiescent state fixed point, (L) is the fixed point in low firing rate, (M) is the fixed point in the linear section, and (H) is the high firing fixed point.

FIGURE 2

Figure 2. Nullclines' arrangements near the BT point in the local bifurcation diagram. The black dashed line is saddle-separatrix loop bifurcation and the blue dotted-dashed is saddle-node on limit cycle (SNLC) bifurcation line.

2.4. Synaptic plasticity

We will analyze and simulate a network in which neurons will adapt their connections according to the Spike-timing dependent plasticity (STDP) paradigm (Gerstner et al., 1996), which provides a foundation for temporal coding.

In STDP, the weight of a connection is modified depending on the time interval between pairs of pre- and post-synaptic spikes. For every pair, the weight of the synapse changes according to the equations

$\begin{array}{ll}\Delta w(\Delta t)=\{\begin{array}{cc}{f}_{+}(w)K_{+}(\Delta t)& \text{if }\Delta t\ge 0\\ -f_{-}(w)K_{-}(\Delta t)& \text{if }\Delta t<0\end{array}& (7)\end{array}$

where Δt = t_post − t_pre is the time difference between the postsynaptic spike and the presynaptic one. The functions f₊ and f₋ model the dependence of the weight change on the current value of the synaptic weights. K₊ and K₋, called STDP kernels, usually are decaying functions of time which reflects the fact that closer pre- and post-synaptic spikes generate stronger weight changes. Usually, we model the kernels by a single exponential such as $K_{+}=A_{+}{e}^{-\frac{|\Delta t|}{{\tau }_{s_{+}}}}$ and $K_{-}=A_{-}{e}^{-\frac{|\Delta t|}{{\tau }_{s_{-}}}}$. As it is evident from equation (7) when the postsynaptic neuron fires after the presynaptic neuron, the strength of the connection increases and it decreases for the opposite temporal order. We assume the same type of the STDP rule for both inhibitory and excitatory connections although with different kernels. In the following, we suppose that the dependence of STDP on the synaptic weight is negligible and therefore replace the functions f₊ and f₋ by a constant which is then absorbed into the kernels. In this case, we have to assume a saturation level for the maximum strength of the synapses, $w_{max}^E$ and $w_{max}^I$.

We use the model of Tsodyks and Markram (1997) for short-term depression of the excitatory synapses reduces the outgoing synaptic efficacy of excitatory synapses to an excitatory neuron in case of a high rate of presynaptic activity. To model the STP effect, we assume that the effective utility of the excitatory synapses of neuron j to the other neurons is proportional to the fraction of the available synaptic resources u. The decrease of neurotransmitters at the synapses and depression in release probability due to consecutive uses of neurotransmitters in previous spikes of the presynaptic neuron are the sources of STP. We assume by each spike of a presynaptic neuron, u is reduced by the factor qu and then recovers with the time constant τ_STP which is of order 100ms to a few seconds. Therefore, synaptic efficacy of the postsynaptic synapse of neuron j evolves as:

$\begin{array}{ll}\frac{d{u}_j}{dt}=\frac{1}{{\tau }_{STP}}(1-u_j)-q{u}_j{\displaystyle \sum _k}\delta (t-t_k^j)& (8)\end{array}$

3. Results

3.1. Long term synaptic plasticity by STDP tunes synaptic weights close to the balanced state

A typical neuron in the cortex has 10³ − 10⁴ synaptic connections with 80% of them of excitatory type and 20% of inhibitory type. On the other hand, even in the resting state, neurons on average have a non-zero firing rate with an average rate of 1Hz and their spike trains are very noisy with exponential inter-spike interval distribution indicating that the spiking of individual neurons is a Poisson point process. Yet another experimental fact about synaptic strength between neuron states is that usually, 10 − 20 presynaptic synchronous spikes suffice to bring a typical neuron to the firing threshold. If we take τ_m = 20ms as the membrane potential decay time constant, then during this time window a typical neuron receives 20 − 200 excitatory spikes, which are enough for the neuron to periodically spike at a very high rate. To avoid this, the inhibitory input in this time window should largely cancel the excitatory current. Therefore, for the average currents to maintain the average membrane potential below the threshold in order to avoid a high firing state and produce high variability in the spike trains, inhibitory and excitatory currents should be balanced. Dynamical balance of excitation and inhibition ensures a low level of activity, i.e., an asynchronous firing state. In the following, we present a synaptic plasticity rule which tunes the average synaptic weights to the balanced state.

We derive an equation for the evolution of the average and the variance of weights between excitatory and inhibitory neurons during the plasticity period.

Spike-timing dependent plasticity (Equation 7) changes synaptic weights on a very slow time scale compared to firing dynamics of the neurons, therefore, during a time period of [t, t + Δt] where Δt is long in comparison with the inter-spike time interval but small enough that the change in the weight w_ij of the synapse from neuron j to neuron i is infinitesimal, we can write:

$\begin{array}{l}\Delta w_{ij}={\displaystyle {\int }_t^{t+\Delta t}}{\displaystyle {\int }_0^{\infty }}{S}_j(s)S_i(s+\delta )f_{+}(w_{ij})K_{+}(\delta )d\delta ds\\ +{\displaystyle {\int }_t^{t+\Delta t}}{\displaystyle {\int }_0^{\infty }}{S}_i(s)S_j(s-\delta )f_{-}(w_{ij})K_{-}(\delta )d\delta ds& (9)\end{array}$

where S_i(t) and S_j(t) are the spike trains of the presynaptic and the postsynaptic neurons. Assuming that during this period, the firing rate of the output neuron is constant on average and there exist many pre- and post-synaptic spikes, we can write the mean change in the incoming synaptic weights to the neuron i as,

$\begin{array}{ll}\frac{{\langle \Delta w_{ij}\rangle }_j}{\Delta T}={\displaystyle {\int }_0^{\infty }}{\langle S_j(s)S_i(s+\delta )\rangle }_j{K}_{+}(\delta )d\delta -{\displaystyle {\int }_0^{\infty }}{\langle S_j(s)S_i(s-\delta )\rangle }_j{K}_{-}(\delta )d\delta & (10)\end{array}$

We want to investigate the evolution of the synaptic weights in the EI population in an asynchronous irregular state. Therefore, we assume that in the regime of spontaneous activity, neurons are firing as a Poisson process. Moreover, to estimate the cross-correlation of the pre- and the post-synaptic spike train we argue that the excitatory input to the cell has a positive correlation with preceding spikes in the target neuron. The magnitude of this excess correlation depends on the weight of the synapse and it is restricted to the time window before the firing of the postsynaptic neurons. With this in mind, we use the following approximation introduced in van Rossum et al. (2000) and Câteau and Fukai (2003) for the cross-correlations of spike trains to account for the causal contributions of presynaptic spikes to the postsynaptic ones of a synapse with the strength w_i:

$\begin{array}{l}\langle S_{pre}^E(s)S_{post}^E(s+\delta )\rangle ={\rho }_{pre}^E{\rho }_{post}+{\rho }_{pre}^E{w}_i(V_{Rexc}-V_i){\gamma }^E(\delta )\\ \langle S_{pre}^I(s)S_{post}^E(s+\delta )\rangle ={\rho }_{pre}^I{\rho }_{post}-{\rho }_{pre}^I{w}_i(V_i-V_{Rinh}){\gamma }^I(\delta )& (11)\end{array}$

Here, V_i is the voltage level of the postsynaptic neuron. As the second terms in both equations encode the excess correlation (anticorrelation) of the presynaptic excitatory(inhibitory) input preceding the firing at the postsynaptic neuron, we set γ^I(δ) = γ^E(δ) = 0 for δ < 0. For positive values of δ, this function which is independent of the rates and the weights of the synapses encodes the causal effect of the presynaptic spike which arrives δ units of time before firing of the postsynaptic neuron. Therefore, it is a decaying function of δ. Moreover, we have assumed the dependence on the weight of the synapse to be of a linear form, which is a good approximation in the regime of small synaptic weights. Inserting the above approximation and labeling the STDP kernels of Exc. to Exc. (EE) and Exc. to Inh. (IE) synapses as K^E and the STDP kernels of Inh. to Inh. (II) and Inh. to Exc. (EI) synapses as K^I, we can write the evolution of the average excitatory and inhibitory synaptic strength to the neuron i as

$\begin{array}{l}\frac{d\langle w_{ij}^E\rangle }{dt}=\langle {\rho }_j^E\rangle {\rho }_i({\overline{K}}_{+}^E-{\overline{K}}_{-}^E)+\langle {\rho }_j^E\rangle \langle w_{ij}^E\rangle (V_{Rexc}-\langle V_i\rangle )\overline{K_{+}^E{\gamma }^E}\\ \frac{d\langle w_{ik}^I\rangle }{dt}=\langle {\rho }_k^I\rangle {\rho }_i({\overline{K}}_{+}^I-{\overline{K}}_{-}^I)-\langle {\rho }_k^I\rangle \langle w_{ik}^I\rangle (\langle V_i\rangle -V_{Rinh})\overline{K_{+}^I{\gamma }^I}& (12)\end{array}$

Here, bars denote integrals of the kernels on the positive or negative real lines. In the population of sparsely connected and sufficiently homogeneous neurons, in terms of the number of connections of each neuron, and the regime of asynchronous homogeneous firing, i.e., when all the neurons fire with the same average rate but with a random phase of firing between them, the average weights evolve as

$\begin{array}{l}\frac{d{w}_{EE}}{dt}={\rho }_E^2\widehat{K^E}+{\rho }_E{w}_{EE}(V_{Rexc}-\langle V^E\rangle )\overline{K_{+}^E{\gamma }^E}\\ \frac{d{w}_{EI}}{dt}={\rho }_E{\rho }_I\widehat{K^I}-{\rho }_I{w}_{EI}(\langle V^E\rangle -V_{Rinh})\overline{K_{+}^I{\gamma }^I}\\ \frac{d{w}_{IE}}{dt}={\rho }_E{\rho }_I\widehat{K^E}+{\rho }_E{w}_{IE}(V_{Rexc}-\langle V^I\rangle )\overline{K_{+}^E{\gamma }^E}\\ \frac{d{w}_{II}}{dt}={\rho }_I^2\widehat{K^I}-{\rho }_I{w}_{II}(\langle V^I\rangle -V_{Rinh})\overline{K_{+}^I{\gamma }^I}& (13)\end{array}$

From the above equations, it is straightforward to see when $\widehat{K^E}:={\overline{K}}_{+}^E-{\overline{K}}_{-}^E<0$ and $\stackrel{-}{K^I}:={\overline{K}}_{+}^I-{\overline{K}}_{-}^I>0$, the stationary solutions satisfy:

$\begin{array}{ll}\frac{c_{EI}^{st}}{c_{EE}^{st}}=\frac{\widehat{K^I}\overline{K_{+}^E{\gamma }^E}}{\widehat{K^E}\overline{K_{+}^I{\gamma }^I}}=\frac{c_{II}^{st}}{c_{IE}^{st}}& (14)\end{array}$

We take the proportion of inhibitory synapses to excitatory synapses to be equal for both excitatory and inhibitory neurons, i.e., $\frac{k_{EI}}{k_{EE}}=\frac{k_{II}}{k_{IE}}$. The above condition brings the slopes of the excitatory and the inhibitory nullclines close to each other (Equation 6) leading to the intersection in the semi-linear regime and proportionality of excitation and inhibition:

$\begin{array}{ll}\frac{{\rho }_I^{st}}{{\rho }_E^{st}}\approx \frac{c_{EE}}{c_{EI}}& (15)\end{array}$

Figure 3 shows that STDP tunes the weights to this balanced inhibition dominated state in which ΔW = C_EEC_II − C_EIC_IE ≈ 0. As the system settles in this state, synaptic plasticity has a strong effect when neurons are in a higher (here the linear) firing regime. In this state, rates vary co-linearly according to the above equation. On the other hand, synaptic plasticity rules for w_II and w_EI, i.e., the second and fourth lines in equation (14), lead to a relation for stationary weights in the form of $\frac{c_{II}}{c_{EI}}=\frac{k_{II}{\rho }_I^{st}}{k_{EI}{\rho }_E^{st}}$. Comparing these last two equations, we arrive at $k_{II}{c}_{EE}^{st}=k_{EI}{c}_{II}^{st}$. Assuming k_II = k_EI, the mentioned relation adjusts the trace of the Jacobian at the fixed point in the linear section to be near zero. Therefore, the plasticity rule and the dynamics of the near-linear regime stabilize the system near the BT point in the long term. When the external input to Exc. and Inh. populations are finely tuned, STDP alone can lead to the emergence of avalanches in the system. In the following section, by introducing a short-term plasticity mechanism we aim to tune the system at the critical point in the presence of fluctuating input.

FIGURE 3

Figure 3. Effect of synaptic plasticity on the network with three different initial weight configurations when external excitatory input to both excitatory and inhibitory populations are of the same magnitude ρ_Ext = 150 Hz. (A–D) Evolution of average synaptic weights by STDP. (E) Change in the balance condition by STDP. Slopes of the Exc. and the Inh. nullclines approach each other under STDP in all three configurations. (F) The final state of the average neuron firing rates for these three networks lie below 1Hz. (G) Network activity for different clusters of neurons with different overall average inward synaptic weights. As STDP leads to an increase in the variance of the weight distribution, clusters with different overall connectivity strengths and correspondingly different average rates emerge.

3.2. Short-term plasticity tunes the network in a wide range of external input

In the following, first, we will discuss the adaptive role of short-term synaptic plasticity in bringing the network of the EI population to the avalanche regime. Afterward, we will discuss how internal or external noise close to the BT point can also cause the switch between the quiescent (Down) and the low firing (Up) states. We will discuss the Up-Down state transition by short-term depression can be achieved either through a switch between bi-stable states or by bringing the system close to the BT point by dampening the overall excitation. Here, we just consider the short-term plasticity of synapses between excitatory neurons. This type of plasticity might occur in other types of synapses as well, but we will not discuss this here. Because there exist numerous input synapses and we have assumed homogeneous connectivity, each neuron senses a large sample of the network activity and is connected with an overall average weight with a small variance to the excitatory neuron pool. With regard to these assumptions and structural homogeneity and based on Equation (8), we can write down the dynamic of the average synaptic weights to the neuron i in the state of the network with an excitatory population firing rate of magnitude ρ_E as:

$\begin{array}{ll}\frac{d{w}_{EE}}{dt}=\frac{w_{EE}^0-w_{EE}(t)}{{\tau }_{STP}}-w_{EE}(t)q{\rho }_E(t)& (16)\end{array}$

The rate equations for the EI population are of the form:

$\begin{array}{l}\frac{d{\rho }_E}{dt}=-\frac{1}{{\tau }_m}({\rho }_E(t)-f({\rho }_E(t),{\rho }_I(t),w_{EE}(t)))\\ \frac{d{\rho }_I}{dt}=-\frac{1}{{\tau }_m}({\rho }_I(t)-g({\rho }_E(t),{\rho }_I(t)))& (17)\end{array}$

Taking the time scale of short term plasticity to be much larger than the EI-network activity decay time constant, i.e., τ_STP >> τ_m, we can rewrite the dynamic in terms of fast, d/dt^f, and slow time, d/dt^s, evolution. Here, t_f = t/τ_m and t_s = t/τ_STP. Defining μ = t_s/t_f and ρ = $\left(\begin{array}{c}{\rho }^E\\ {\rho }^I\end{array}\right)$, we arrive at

$\begin{array}{l}\frac{d\rho }{d{t}^f}=-(\rho -f(\rho ,w_{EE}))\\ \frac{d{w}_{EE}}{d{t}^f}=\mu (w_{EE}^0-w_{EE})-q{\tau }_m{w}_{EE}{\rho }^E& (18)\end{array}$

This set of equations can have a stable fixed point or an oscillatory behavior.

The average synaptic efficacy in the stationary state with the average excitatory rate ${\rho }_E^{*}$ is:

$\begin{array}{ll}{\langle w_{EE}\rangle }_{St}=\frac{w_{EE}^0}{1+\tau q{\rho }_E^{*}}& (19)\end{array}$

In case there exists a fixed point or a stable limit cycle solution around this point in the (ρ_E, ρ_I, 〈W_EE〉_st) phase space, the system might settle down at this solution (Figure 4A). The dynamics of the EI population near this region (with the slow-fast assumption) can be written as

$\begin{array}{l}{\rho }_E^{st}=f(k_{EE}{\rho }_E^{st},k_{EI}{\rho }_I^{st},{\lambda }_E^{Ex},{\lambda }_I^{Ex},\frac{w_{EE}^0}{1+\tau q\langle {\rho }^{st}\rangle })\\ {\rho }_I^{st}=g(k_{EI}{\rho }_E^{st},k_{II}{\rho }_I^{st},{\lambda }_E^{Ex},{\lambda }_I^{Ex})& (20)\end{array}$

This mechanism is effective to bring the system close to the BT point. Short synaptic plasticity is a method of gain control that can bring the system from a wide range of input and initial states to the low activity background state. In Figure 5, we consider a system that is already tuned by STDP to the balanced state of weights (Equation 14) receiving various rates of external excitatory input to the excitatory population. In all cases, the system is initially away from the BT point. Without STP, the system shows a high firing rate oscillatory activity with an average rate of around 300Hz (with Nullcline configuration of Figure 7C). STP brings all of them closer to the avalanche regime. The average synaptic efficacy 〈u〉 in these cases does not oscillate significantly.

FIGURE 4

Figure 4. Output excitatory rates as a function of W_EE and the corresponding graphs for the average synaptic efficacy 〈W_EE〉_St at three values of q (dashed red curve belongs to the largest and the dashed green curve is for the lowest value). Based on the value of $W_{EE}^0$, two different scenarios can occur. In (A), by decreasing q, through a saddle-node bifurcation stable and unstable fixed points appear at low and high values of the rates. In (B), with higher $W_{EE}^0$, by decreasing q after the Hopf bifurcation of the low firing rate fixed point, an oscillatory solution for (u, ρ_out) emerges.

FIGURE 5

Figure 5. EI population with short term plasticity subjected to three different external excitatory rates ${\rho }_{Ext}^{exc}$ = [380, 280, 240]Hz in panels (A1–A3), (B1–B3), and (C1–C3) respectively. Other parameters of the model are W_EI = 1.5, W_II = 2, W_IE = 0.75, $W_{EE}^0=0.54$, ${\rho }_{Ext}^{inh}$= 150Hz, q = 0.3 and τ_STP = 10 * τ_syn. Left plots (A1, B1, and C1), show the excitatory population rates, middle plots (A2, B2 and C2), show the population average membrane potential and right plots (A3, B3 and C3), show the average excitatory synaptic efficacy 〈u〉.

When the external input rate is tuned very close to the BT point, where the quiescent fixed point and a low firing weakly (un)stable point lie close to each other, we see asynchronous avalanches of highly variable sizes (refer to Figures 6A, 7A). Without STP, we observe higher rates and less variable quiescent (Down) state time intervals (Figure 6B). Finite-size fluctuations kick the system out of the quiescent fixed point while STP plus fluctuations at the low rate fixed point drive the system back to the quiescent state. Average membrane potential and single neuron potentials in this case show a transition between two levels (refer to Figure 6B).

FIGURE 6

Figure 6. (A) System with the same parameters values as in Figure 5 except that here ${\rho }_{Ext}^{exc}=230Hz$. (B) EI population with the same parameters as in (A) but without STP. (C) Avalanche size distribution in a log-log plot for the network in (A) (blue) and two networks with smaller $W_{EE}^0$ which does not show scale free avalanches.

FIGURE 7

Figure 7. Excitatory (red) and inhibitory (blue) nullclines, i.e., solutions to Equation 4 for excitatory and inhibitory rates, respectively. Dashed lines are the linearized nullcline approximation with slopes of $k\frac{c_{EI}}{c_{EE}}$ (Exc.) and $k\frac{c_{II}}{c_{IE}}$ (Inh.) where k is a constant. STDP brings these slopes close to eachother (Equation 14). (A) Nullclines at the Avalanche regime: The volume of the basin of attraction of the stable quiescent state is small and a weakly stable or unstable fixed point at the intersection in the semi-linear regime of nullclines is close to the origin. (B) System with only a saddle fixed point with a limit cycle solution with medium firing rates. (C) System with only a high firing fixed point. (D) System with only a stable quiescent state fixed point.

Figure 6C shows the cumulative avalanche size distribution function in a log-log plot for avalanches as in Figure 6A. The slope of the linear regression line is close to 0.5 indicating that the avalanche size distribution function is a power law with the exponent τ = −1.5. In the Appendix 1, we have presented statistics of the avalanches including duration distribution function, shape collapse, and power law scaling of mean size vs. mean duration. The branching ratio for the final state of the system is shown in Figure 8B. For ρ_Ext = 230Hz, the branching ratio is slightly less than one which is in agreement with our prediction in the avalanche regime. Figure 8A shows excitatory and inhibitory stationary rates of the EI population subject to external rates in the range 200–500 Hz. STP leads to low firing rate states and prevents overactivation.

FIGURE 8

Figure 8. (A) The final excitatory (red) and inhibitory (blue) output rates for the system in Figures 5, 6A. STP works as a gain control mechanism. (B) The branching ratio in the network states shown in Figures 5, 6A is defined as the average number of post-synaptic neurons of a single presynaptic neuron which set to fire by receiving the presynaptic input spike. Higher external rates set more neurons close to the threshold and thus the branching factor increases. When the inflection point of the steady membrane potential distribution passes the membrane threshold in the steady-state firing rate regime this increment rate slows down leading to the concavity of the branching factor curve.

As the system resides closer to the low firing regime, the effective synaptic strength of Exc.-Exc. connections is lowered by STP (Figure 7B). On the longer time scale STDP transforms the combination of weights to move the system to the avalanche regime by aligning the slopes of linearized nullclines and setting the trace of Jacobian at the fixed point on a linear section close to zero. Figure 9, shows the cooperation of both STDP and STP brings the system initially away from the BT point to the critical avalanche regime. STP prevents high firing as well as being trapped in the quiescent state, therefore constraining the system to be in a low firing state which in a longer time scale leads to tunning in the vicinity of BT point by STDP.

FIGURE 9

Figure 9. (A) STDP and STP tune the weights to the avalanche regime by tuning the inhibitory feedback. (B) Evolution of average weights: W_II (black),W_EI (blue),W_IE (cyan), W_EE (red), $W_{EE}^{eff}=u{W}_{EE}$ (magnet), and $\Delta W=W_{EE}^{eff}{W}_{II}-W_{EI}WIE$ (green). (C) In the avalanche regime, both the determinant(green) and the trace(magnet) of the Jacobian of the fixed point at the linearized nullclines are nearly zero.

Another way that STP can cause a switch between two distinct firing states is in the EI population which possesses bi-stability. In this case, a change of u can make each of the bi-stable nodes unstable while the system resides near them. The decrease of u in the up-state makes the up-state fixed point unstable at some value of u(t) (and accordingly w_EE). Therefore, the system will jump to the remaining stable fixed point in a low or quiescent state. In the very low firing regime (the quiescent state), u will recover to its asymptotic value, and the average synaptic weight increases toward $w_{EE}^0$. If the quiescent state is unstable when u approaches its maximum value, we observe a transition to the high state. Moreover, if the volume of the basin of attraction of the quiescent fixed point is small, external and internal noise can also induce the transition to the high rate fixed point and the quiescent fixed point need not become unstable at $w_{EE}^0$. For high values of u, the up-state fixed point is the stable point of the fast system but an unstable point of the slow one. Therefore, following the slow path up-state loses stability and the fast system remains with only a stable low fixed point. The trajectory of the slow u is oscillatory in this case. Figure 10 shows both ways that STP can produce synchronous avalanche behavior in the system. When $W_{EE}=w^0$, the system is close to the constraints on the alignment of the semi-linear segments of the EI-nullclines which results in the presence of a high firing state as a unique fixed point of the system. In the high input rate case, ρ_Ext = 400Hz, corresponding to Figure 10A and the nullcline diagram of Figure 7C, due to STP, the system moves from a high state of activity to a limit cycle solution at lower firing rates. This final state is shown in Figure 10A and nullcline arrangements in this state are depicted in Figure 7B. Here, there is an unstable source in the linear branch sector which is surrounded by a limit cycle. Moreover, oscillations in 〈u〉 have a low amplitude because of the temporal averaging. On the other hand, the two bottom panels in Figure 10 are related to the situation of a switch between the high fixed point and the quiescent node. As shown in nullcline graphs in Figure 7C, at high synaptic efficacy the high firing state is the only stable fixed point, however, a high firing rate leads to a fast decline of the synaptic efficacy which brings the system to the state with the nullcline map of Figure 7D which has a stable quiescent fixed point. The final state activity, in this case, is composed of avalanches with a high rate of firing in a short time window. Decreasing the q factor can result in a longer up-state period. Also, 〈u〉 oscillates between two limits in these cases.

FIGURE 10

Figure 10. EI population with short term plasticity subjected to four different external excitatory rates ${\rho }_{Ext}^{exc}$ = [400, 330, 270, 225]Hz in panels (A1–A3), (B1–B3), (C1–C3), and (D1–D3) respectively. Other network parmeters are W_EI = 2, W_II = 2, W_IE = 0.75, $W_{EE}^0=0.74$, ${\rho }_{Ext}^{inh}$ = 150 Hz and STP parameters are q = 0.4 and τ_STP = 10 * τ_syn. Left plots (A1–D1): excitatory rates; middle plots (A2–D2): average membrane potential; right plots (A3–D3): average synaptic efficacy 〈u〉.

In this particular case, necessary conditions for up to down transitions are:

$\begin{array}{l}{k}_{EE}{W}_{EE}^0>k_{EI}{W}_{EI}\\ k_{EE}\frac{w_{EE}^0}{1+\tau q{\rho }^H}<k_{EI}{W}_{EI}& (21)\end{array}$

The first condition means that the slope of the excitatory nullcline is smaller than the inhibitory one, which indicates a stable high firing fixed point. The second condition states that at this high firing state the stationary weight is not accessible before a stability loss. The slope of nullclines increases by the decrease of effective W_EE which causes the high state to lose stability either through a Hopf or a saddle-node bifurcation.

Finally, the plots in Figure 10B depict the case where STP brings the system close to the BT point which shows low to medium size avalanches with higher variability. In summary, the transition from a quiescent state to a high firing state can be of two distinct types. One way is that by increasing W_EE, the quiescent fixed point and the unstable saddle move toward each other and in this way the basin of attraction of the quiescent fixed point shrinks and noise can initiate the escape from this fixed point to the higher firing state. This is the mechanism for alternation between quiescent periods and avalanche periods, corresponding to nullclines arrangements of Figures 7A,D with dynamics of Figure 6A. The other possibility is that a fixed point losses stability through a Hopf bifurcation either before or after the emergence of a saddle-node in the middle branch.

3.3. Continuum stochastic model of neuronal dynamics near BT point

We have shown that a single EI population in both regimes of sparse connectivity and all-to-all connectivity (Ehsani and Jost, 2022) exhibits scale-free avalanches near the critical point with Poisson firing of individual neurons. Next, we want to investigate if avalanche dynamics in weakly interconnected EI population persists. Our goal is to introduce a phenomenological stochastic field equation for the dynamics of population rates in the vicinity of the BT point.

We start with the model of Wilson and Cowan for the dynamics of the spatio-temporal mean fields of the excitatory and the inhibitory population rates, E(x, t) and I(x, t), in a 2D model of the cortex. Defining V(x, t)= $\left(\begin{array}{c}E(x,t)\\ I(x.t)\end{array}\right)$, D as the diffusion matrix and f : R² → R² as the gain function, one can write down the reaction-diffusion rates dynamics in the following form:

$\begin{array}{ll}\tau \frac{\partial V(x,t)}{\partial t}=D{\nabla }^{2}V(x,t)+f(V(x,t))& (22)\end{array}$

The ODE part of this equation is the dynamic of a single EI population. The corresponding low firing fixed point (E₀, I₀) is stable in a specific parameter regime. It usually loses stability either via a saddle-node or a Hopf bifurcation which leads to either a region of bi-stability of low and high firing states or the emergence of oscillations. However, still far away from the bifurcation point since the diffusion matrix is not a scalar multiple of the identity, Turing instabilities can occur in the system.

Besides Turing instabilities, it can also happen that the fixed point itself loses its stability at k = 0 through a Hopf bifurcation, where the real part of the eigenvalues of L becomes zero:

$\begin{array}{l}{\partial }_E{f}_E+{\partial }_I{f}_I=0\\ det(L)={\partial }_E{f}_E{\partial }_I{f}_I-{\partial }_I{f}_E{\partial }_E{f}_I>0& (23)\end{array}$

Furthermore, exactly at the BT point, we would also have det(L) = 0. Figure 11 shows the activity of 20 interconnected EI-populations each operating close to the BT point. Overall activity in this system is of synchronous avalanches type. Up-down state transitions also become synchronized. We can model weakly interconnected EI populations in the avalanche regime which shows oscillation of frequency ω_i as pulse-coupled oscillators and therefore investigate conditions on synchronization and traveling wave solutions. This analysis is out of the scope of the current work. Another approach consists in supplementing the macroscopic field equation with an appropriate noise term to derive the mesoscopic equation. As can be seen from Figure 11, the overall network activity is of avalanche type. This provides again evidence that avalanches are scale-free and occur in different temporal and spatial scales. However, our neural field model still lacks the internal finite-size fluctuation effects, inhomogeneities, and cross-correlation between individual neurons, and also inter-populations correlations.

FIGURE 11

Figure 11. Simulation result of 20 interconnected EI populations each of size N_E = 10, 000 arranged on a ring. Average synaptic weights between two different EI subnetworks decay with the phase difference of them. Parameters of synaptic weights in each EI population are W_EE = 0.6, W_EI = 2, W_II = 2, W_IE = 0.75, ${\rho }_{Ext}^E=230Hz$, and ${\rho }_{Ext}^I=150Hz$. (A) The average firing rate of the whole network. (B) Raster plot of the sub-population of 50 neurons in one single EI subnetwork. Although neurons fire at avalanche times, they do not fire in all of them, which leads to high variability in their interspike interval times. (C) The membrane potential of a single neuron is in a constant transition between a state close to the threshold and a state close to the resting potential. (D) The average membrane potential of an EI-population (red for Exc. and blue for Inh.). (E,F) The activity of two distinct EI subpopulations shows high variability in the sizes of avalanches in both of them.

To investigate the fluctuations around mean field rate dynamics, we study a Markov model inspired by the fact that close to the BT point, the firing of individual neurons is a Poisson point process. We consider a homogeneous network of size N in which temporal and spatial variances in the firing rates of neurons are minimal. In this network, fluctuations in the finite system firing rates in the steady-state will be proportional to $O(\frac{1}{\sqrt{N}})$. To model the finite-size stochastic effects, we need to write down dynamics of micro-state evolution that match the mean-field upon coarse-graining. As we have seen, the operating region of the EI population is around a low firing state where neurons fire with high variability of inter-spike intervals indicating that we can model their spiking as a Poisson process. In this regime of activity, we can write down the microscopic evolution of a model neuron with two active and inactive states. The transition rate α between the active and the inactive state should model the vanishing of the postsynaptic potentiation, and the rate of inactive to active transition depends on the input and is, therefore, denoted by f(i). We want to model the system in the statistical homogenous state, in which the probability that a neuron fires depends only on the number of active neurons and therefore is the same for every neuron in the population.

After using system size expansion method (refer to Appendix 2), we arrive at the following equations for variance of population rate at the stationary point of the macroscopic equation:

$\begin{array}{ll}\left(\begin{array}{ccc}{A}_{11}& 0& A_{12}\\ 0& A_{22}& A_{21}\\ A_{21}& A_{12}& A_{11}+A_{22}\end{array}\right)\left(\begin{array}{c}Var{(ϵ)}_{st}\\ Var{(i)}_{st}\\ Cov{(ϵ,i)}_{st}\end{array}\right)=-\frac{\alpha }{2}\left(\begin{array}{c}{\rho }_E^{st}\\ {\rho }_I^{st}\\ 0\end{array}\right)& (24)\end{array}$

which has the solution:

$\begin{array}{l}Var{(ϵ)}_{st}\approx c((A_{11}{A}_{22}-A_{21}{A}_{12}+A_{22}^2){\rho }_E+A_{12}^2{\rho }_I)\\ Var{(i)}_{st}\approx c((A_{11}{A}_{22}-A_{21}{A}_{12}+A_{11}^2){\rho }_I+A_{21}^2{\rho }_E)\\ Cov{(ϵ,i)}_{st}\approx -c(A_{11}{A}_{12}{\rho }_I+A_{21}{A}_{22}{\rho }_E)& (25)\end{array}$

with $c=\frac{-\alpha }{2(A_{11}+A_{22})(A_{11}{A}_{22}-A_{21}{A}_{12})}$. The average population rate and the fluctuation around the macroscopic state are:

$\begin{array}{l}\langle \frac{E}{N_E}\rangle ={\rho }_{E}Var(\frac{E}{N_E})=\frac{Var(ϵ)}{N_E}\\ \langle \frac{I}{N_I}\rangle ={\rho }_{I}Var(\frac{I}{N_I})=\frac{Var(i)}{N_I}& (26)\end{array}$

From equations (25) and (26), it can be seen that close to the bifurcation of the macroscopic equation, i.e., the BT point, where both trace and determinant of the Jacobian are close to zero, fluctuation magnitudes increase. The important point is that the population variance is linear in the rates.

Moreover, close to the BT region, the fixed point on the semi-linear regime has an eigenvalue close to zero. Dynamics along the slow manifold which provides the detailed temporal balance of inhibition and excitation enables us to write the instantaneous inhibition rate as the linear function of the excitatory rate. Therefore, we can write down the rate dynamics of the excitatory population in terms of a stochastic field equation when both inhibitory feedback and fluctuations are local. From equation (15), we know that near the BT point, there is a linear relation between rates, i.e., $I\approx \frac{c_{EE}}{c_{EI}}E$. Therefore, the average current to the excitatory population close to the BT point (Equation 6) can be written as γ = c_EEE − c_EII ≈ 0. The second derivative in the expansion of the gain function for the excitatory population from equation (22) in the region of low activity is

$\begin{array}{ll}\frac{1}{2}\frac{{\partial }^{2}f}{\partial E^2}{E}^2+\frac{{\partial }^{2}f}{\partial EI}EI+\frac{1}{2}\frac{{\partial }^{2}f}{\partial I^2}{I}^2& (27)\end{array}$

where we can use the following approximations for the gain function derivatives:

$\begin{array}{ll}\frac{{\partial }^{2}f}{\partial I^2}\propto W_{EI}^2,\frac{{\partial }^{2}f}{\partial E^2}\propto W_{EE}^2,\frac{\partial f}{\partial I\partial E}\propto -W_{EE}{W}_{EI}& (28)\end{array}$

By using the linear relation of the inhibitory and excitatory rates near the BT as the result of the projection of the dynamic to the slow manifold, we can replace inhibitory local field strength by a term linear in the local excitatory field. Besides, fluctuations in the average population activity, Equation (26), linearly depend on the rate. Therefore, we can write down the stochastic field equation for the excitatory rate in the region of small γ:

$\begin{array}{l}\frac{\partial E(x,t)}{\partial t}=\gamma E(x,t)+D\Delta E(x,t)-u{E}^2(x,t)+\psi (x,t)\\ \langle \psi (x,t)\psi (x^{\prime },t^{\prime })\rangle =2\sqrt{\frac{{\sigma }^2}{N}}E(x,t)\delta (x-x^{\prime })\delta (t-t^{\prime })& (29)\end{array}$

Here, u < 0 is the coefficient related to synaptic weights that can be explicitly derived by assuming a certain form of the gain function and proportionality of the rates. This stochastic partial differential equation after appropriate rescaling $E(x,t)=\frac{\sigma }{u\tau }S(x,t)$ agrees with the Langevin description of directed percolation which is of the following form with new transformed coefficients:

$\begin{array}{l}\frac{\partial S(x,t)}{\partial t}=({\gamma }^{\prime }+D^{\prime }\Delta )S(x,t)-u^{\prime }{S}^2(x,t)+\psi (x,t)\\ \langle \psi (x,t)\psi (x^{\prime },t^{\prime })\rangle =2u^{\prime }S(x,t)\delta (x-x^{\prime })\delta (t-t^{\prime })& (30)\end{array}$

At γ′ = 0, the above system shows an absorbing state phase transition. Thus, from any active state, the system relaxes by avalanches with a power-law size distribution to an inactive state.

In an isolated EI population, external drive to the inhibitory and the excitatory population should be present to counterbalance the dissipation by the leaking currents and thereby set the average membrane potential in these neurons at a state above the resting threshold. External excitatory input to the excitatory population is slightly higher than the external drive to the inhibitory population which leads to a slightly higher average membrane potential in the excitatory population. Furthermore, we can assume that the external spike train to each neuron is Poisson as well. The external drive by itself would not lead to significant firing in the individual neurons but the strengths of the internal connections between them are tuned in a way that bursts of activity occur in the excitatory population which is then followed by the inhibitory ones. The internal feedback inhibition is strong enough to kill the excitatory burst. In a slightly inhibition-dominated regime, we have sharp synchronous responses to the external input in a short time window. On the other hand, the network has a safe margin from an overly active state. In the absence of the input distinguished from random external noise, the system shows scale-free avalanches because of the maintenance of the inhibition-excitation balance. However, the external drive must compensate for the dissipation of the system to stay at or near the critical point. Without mechanisms like short-term plasticity, the external drive has to be fine-tuned for the system to show criticality. However, short-term plasticity in a network in which synaptic weights are already near a slightly inhibition dominated regime broadens the range of the external drive strength which leads to critical avalanches. Figure 12 shows final output rates for three different values of the external excitatory rates and $W_{EE}^0$ with STP. In all these cases, the operating point of the system is close to the bifurcation point of the quiescent state.

FIGURE 12

Figure 12. Solid curves are the output rates for three different external input strengths vs. W_EE: Blue (400Hz), red (310Hz), and black (220Hz). Dashed curves show average stationary synaptic weight, 〈W_EE〉, in the network with STP with different maximum synaptic efficacies: $W_{EE}^0=$1.3 (red), 0.9 (blue), 0.7 (green), and 0.65 (magenta). Intersections of the dashed and the solid curves are the fixed points of the EI network with STP for the corresponding control parameters. These fixed points are located in the low firing rate regime close to the avalanche region.

We can extend the short term synaptic depression of equation (8) to a continuum field equation by defining a field of excitatory synaptic efficacy Ω(x, t) ∝ 〈W_EE〉(x, t) with local dynamics of equation (16):

$\begin{array}{l}{\tau }_m\frac{\partial E(x,t)}{\partial t}=(-\alpha +\Omega (x,t))E(x,t)+D\Delta E(x,t)\\ -u{E}^2(x,t)+\psi (x,t)\\ \frac{\partial \Omega (x,t)}{dt}=\frac{1}{{\tau }_{STP}}({\Omega }_0-\Omega )-q\Omega E(x,t)\\ \langle \psi (x,t)\psi (x^{\prime },t^{\prime })\rangle =2\sqrt{\frac{{\sigma }^2}{N}}E(x,t)\delta (x-x^{\prime })\delta (t-t^{\prime })& (31)\end{array}$

Here, α represents both the decay of activity by leaky currents of the cells and the inhibition feedback which varies linearly with the excitatory rate. The dynamic excitatory synaptic strength brings the coefficient of the linear term to a value near zero (refer to Figure 13). This set of equations has a stationary synaptic efficacy solution of the value

$\begin{array}{ll}{\Omega }_{st}=\frac{{\Omega }_0}{1+q{\tau }_{STP}{E}_{st}}& (32)\end{array}$

On the other hand, in the active phase, the stationary homogeneous rate is :

$\begin{array}{ll}{E}_{st}=\frac{-\alpha +{\Omega }_{st}}{u}& (33)\end{array}$

Assuming |u| is a very small quantity and E_st is also small in the low firing rate regime then −α + Ω_st ≈ 0 and $E_{st}=\frac{1}{q{\tau }_{STP}}(\frac{{\Omega }_0}{\alpha }-1)$.

FIGURE 13

Figure 13. Intersection of curves given by equations (32) and (33). Setting Ω₀ to a value close to α by STDP and sufficient amount of synaptic depression leads to a stationary value of Ω_st very close to the critical value. The blue curve is associated with a higher value of q.

Long-term synaptic plasticity tunes Ω₀ so that the coefficient of the linear term is close to zero and a moderate level of short-term depression suffices to bring the system to the critical point. Altogether, equation (31) is the description of an EI interconnected spiking neuron network tuned to the critical point of balancing inhibition and excitation both by long-term synaptic plasticity and short-term synaptic depression. The system wanders around the phase transition point and shows avalanche dynamics with scale-free size and time distribution, the proportionality of inhibition and excitation, up and down state transitions of the membrane potential and population activity rates, and oscillations of order of the 10Hz resembling ubiquitous alpha-band oscillations.

4. Discussion

In this study, we have proposed a self-organizing model for the cortical dynamics which tunes the system to the regime of low firing avalanche dynamics corresponding to the ongoing intrinsic activity in the cortex. We showed that long-term synaptic plasticity by STDP tunes the synaptic weights to achieve the internal balance of inhibition and excitation. This effect does not depend on the exact shape of STDP kernels, however, the stability of results depends on the sign of the integral of the kernels. On the other hand, short-term depression of excitatory synapses can tune the system in response to the wide range of the strength of the external drive. We have not considered short-term plasticity for other types of synapses. However, we believe that under appropriate conditions on their respective strengths, we can observe the same qualitative regulation effects. Short-term homeostatic regulations acting at the time scale of neuronal dynamics are needed to self-tune the open non-conservative system subjected to varied external input at the critical point. On the other hand, long-term plasticity is responsible for bringing the system close to the bifurcation point. The interplay of these two types of synaptic plasticity is responsible for self-tuning to the critical state.

At the critical state, our system shows power law distribution functions for avalanches size and duration which matches the branching process and mean-field directed percolation process. In addition, scaling relation among the exponents and avalanche shape collapse is another proof that we are at the edge of a second-order phase transition. We have shown in Ehsani and Jost (2022) that close to the BT point, the branching factor is close to one because of the tight balance of inhibition and excitation and Poisson firing of neurons in the avalanche regime.

Feedback inhibition is proportional to the excitatory rate on the slow manifold i.e., I(t) = αE(t) where $\alpha \propto \frac{W_{EE}}{W_{EI}}$. By using this linear relation between rates, we rewrite the Wilson-Cowan field equation in the regime of low firing rate where the system is locally close to the BT point. Moreover, we showed that internal noise in local EI populations in this state is of Poisson type. Therefore, by knowing the noise variance and the fact that there is a tight temporal balance of inhibition and excitation, we showed that the excitatory field dynamics formulated as a stochastic field equation in the low firing rate regime matches the Langevin description of directed percolation at the absorbing state phase transition. However, there is not such an absolute absorbing state in our system. External drive to the neurons, which is also stochastic, would lead to the background level of firing and fluctuations in the network. In an open system, the external load has to be fine-tuned to compensate for the dissipation to remain at the critical point. Short-term depression of excitatory synapses allows this tuning for a wider range of external drives. Tuning the system at the critical point can be achieved by coupling fast population dynamics with slow adaptive gain and synaptic weight dynamics, which make the system wander around the phase transition point. Therefore, by introducing short-term and long-term synaptic plasticity, we have proposed a self-organized critical stochastic neural field model (Equation 31).

We have shown that in a weakly interconnected neuronal mass model the avalanche dynamics persist when local subpopulations are close to the BT critical point. However, a comprehensive study of the coupled field equation dynamics in a 2D or 3D network requires further study. Types of solution and dynamic repertoire of this coupled stochastic field equations have not been done in this study. Neural field models, without noise, can show vast dynamical solutions, including wave propagation in terms of a front solution in a bistable network (refer to Amari, 1977; Ermentrout, 1998), propagating pulses in an excitable medium (refer to Pinto and Ermentrout, 2001b), and spatially localized oscillations, spiral waves in the oscillatory regime of a local EI population (Troy and Shusterman, 2007), localized bump solutions (Pinto and Ermentrout, 2001a) and spatially periodic patterns called Turing patterns (Ermentrout and Cowan, 1979 also refer to Bressloff, 2011, for a comprehensive review). Analysis of the full stochastic version of neuronal field formulated in terms of coupled or non-coupled SPDE is a rather new area of study. The stochastic version of the continuum neural field for the excitatory population without STP has been discussed in Buice and Cowan (2007) and Bressloff (2019). Buice and Cowan used a coherent path formalism and the Doi-Peliti functional representation to translate the microscopic master equation to a path integral representation for activity fields (Buice and Cowan, 2007). One advantage of their method is that the study of scale invariance at criticality in the functional representation is possible. They proposed that the stochastic neural field equation for the excitatory system at a critical point can be written in the form of Langevin's description of directed percolation. Using coherent path formalism to investigate the critical behavior of the coupled system of Equation (31) can be illuminative for understanding the critical behavior of the system. Moreover, we have written down the dynamics of the excitatory field assuming a fast linear inhibitory feedback. Modeling the whole dynamic as a directed percolation of two dynamical variables is another possibility that we have not explored.

Data availability statement

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

Author contributions

ME and JJ designed the research. ME performed research and wrote the manuscript. JJ edited the manuscript. All authors reviewed the manuscript, contributed to the article, and approved the submitted version.

Funding

This study was funded by the International Max Planck Research School and the Max Planck Institute for Mathematics in Sciences in Leipzig.

Acknowledgments

Open access funding is provided by Max Planck Society. ME wants to thank International Max Planck Research School for funding his position during the research period.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

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

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2022.910735/full#supplementary-material

References

Amari, S. (1977). Dynamics of pattern formation in lateral-inhibition type neural fields. Biol. Cybern. 27, 77–87. doi: 10.1007/BF00337259

PubMed Abstract | CrossRef Full Text | Google Scholar

Beggs, J., and Plenz, D. (2003). Neuronal avalanches in neocortical circuits. J. Neurosci. 23, 11167–11177. doi: 10.1523/JNEUROSCI.23-35-11167.2003

CrossRef Full Text | Google Scholar

Benayoun, M., Cowan, J., van Drongelen, W., and Wallace, E. (2010). Avalanches in a stochastic model of spiking neurons. PLoS Comput. Biol 6, e1000846. doi: 10.1371/journal.pcbi.1000846

PubMed Abstract | CrossRef Full Text | Google Scholar

Bornholdt, S., and Roehl, T. (2003). Self-organized critical neural networks. Phys. Rev. E 67, 066118. doi: 10.1103/PhysRevE.67.066118

PubMed Abstract | CrossRef Full Text | Google Scholar

Bressloff, P. (2011). Spatiotemporal dynamics of continuum neural fields. J. Phys. A 45, 033001. doi: 10.1088/1751-8113/45/3/033001

PubMed Abstract | CrossRef Full Text | Google Scholar

Bressloff, P. (2019). Stochastic neural field model of stimulus-dependent variability in cortical neurons. PLoS Comput. Biol. 3, e1006755. doi: 10.1371/journal.pcbi.1006755

PubMed Abstract | CrossRef Full Text | Google Scholar

Brochini, L., de Andrade Costa, A., Abadi, M., Roque, A., Stolfi, J., and Kinouchi, O. (2016). Phase transitions and self-organized criticality in networks of stochastic spiking neurons. Sci. Rep. 6, 35831. doi: 10.1038/srep35831

PubMed Abstract | CrossRef Full Text | Google Scholar

Buice, M., and Cowan, J. (2007). Field-theoretic approach to fluctuation effects in neural networks. Phys. Rev. E 75, 051919. doi: 10.1103/PhysRevE.75.051919

PubMed Abstract | CrossRef Full Text | Google Scholar

Burkitt, A., Gilson, M., and van Hemmen, J. (2007). Spike-timing-dependent plasticity for neurons with recurrent connections. Biol. Cybern. 96, 533–546. doi: 10.1007/s00422-007-0148-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Câteau, H., and Fukai, T. (2003). A stochastic method to predict the consequence of arbitrary forms of spike-timing-dependent plasticity. Neural Comput. 15, 597–620. doi: 10.1162/089976603321192095

PubMed Abstract | CrossRef Full Text | Google Scholar

Chialvo, D. (2004). Critical brain networks. J. Phys. A Stat. Mech. Appl. 340, 756–765. doi: 10.1016/j.physa.2004.05.064

CrossRef Full Text | Google Scholar

Chialvo, D. (2010). Emergent complex neural dynamics. Nat. Phys. 6, 744–750. doi: 10.1038/nphys1803

CrossRef Full Text | Google Scholar

de Candia, A., Sarracino, A., Apicella, I., and de Arcangelis, L. (2021). Critical behaviour of the stochastic wilson-cowan model. PLoS Comp. Biol. 17, e1008884. doi: 10.1371/journal.pcbi.1008884

PubMed Abstract | CrossRef Full Text | Google Scholar

di Santo, S., Villegas, P., Burioni, R., and Muñoz, M. (2018). Landau-ginzburg theory of cortex dynamics: scale-free avalanches emerge at the edge of synchronization. Proc. Natl. Acad. Sci. U.S.A. 115, E1356-E1365. doi: 10.1073/pnas.1712989115

PubMed Abstract | CrossRef Full Text | Google Scholar

Ehsani, M., and Jost, J. (2022). Scale-free avalanches in excitatory-inhibitory populations of spiking neurons with conductance based synaptic currents [preprint]. Available online at: http://arxiv.org/abs/2203.07834

Google Scholar

Ermentrout, B. (1998). Neural networks as spatio-temporal pattern-forming systems. Rep. Prog. Phys. 61, 353. doi: 10.1088/0034-4885/61/4/002

CrossRef Full Text | Google Scholar

Ermentrout, G., and Cowan, D. (1979). A mathematical theory of visual hallucination patterns. Biol. Cybern. 34, 137–150. doi: 10.1007/BF00336965

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerstner, W., Kempter, R., van Hemmen, J., and Wagner, H. (1996). A neuronal learning rule for sub-millisecond temporal coding. Nature 383, 76–78. doi: 10.1038/383076a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Hahn, G., Petermann, D., Havenith, M., Yu, S., Singer, W., Plenz, D., et al. (2010). Neuronal avalanches in spontaneous activity in vivo. J. Neurophysiol. 104, 3312–3322. doi: 10.1152/jn.00953.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Hardstone, R., Poil, S., Schiavone, G., Jansen, R., Nikulin, V., Mansvelder, H. D., et al. (2012). Detrended fluctuation analysis: a scale-free view on neuronal oscillations. Front. Physiol. 3, 450. doi: 10.3389/fphys.2012.00450

PubMed Abstract | CrossRef Full Text | Google Scholar

Kempter, R., Gerstner, W., and van Hemmen, J. (2001). Intrinsic stabilization of output rates by spike-based hebbian learning. Neural Comput. 13, 2709–2741. doi: 10.1162/089976601317098501

PubMed Abstract | CrossRef Full Text | Google Scholar

Kistler, W., and van Hemmen, J. (1999). Short-term synaptic plasticity and network behavior. Neural Comput. 11, 1579–1594. doi: 10.1162/089976699300016151

PubMed Abstract | CrossRef Full Text | Google Scholar

Kistler, W., and van Hemmen, J. (2000). Modeling synaptic plasticity in conjuction with the timing of pre- and postsynaptic action potentials. Neural Comput. 12, 385–405. doi: 10.1162/089976600300015844

PubMed Abstract | CrossRef Full Text | Google Scholar

Levina, A., Herrmann, M., and Geisel, T. (2007). Dynamical synapses causing self-organized criticality in neural networks. Nat. Phys. 3, 857–860. doi: 10.1038/nphys758

CrossRef Full Text | Google Scholar

Levina, A., Herrmann, M., and Geisel, T. (2009). Phase transitions towards criticality in a neural system with adaptive interactions. Phys. Rev. Lett. 102, 118110. doi: 10.1103/PhysRevLett.102.118110

PubMed Abstract | CrossRef Full Text | Google Scholar

Linkenkaer-Hansen, K., Nikouline, V., Palva, V., and Ilmoniemi, R. (2001). Long-range temporal correlations and scaling behavior in human brain oscillations. J. Neurosci. 21, 1370–1377. doi: 10.1523/JNEUROSCI.21-04-01370.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Meisel, C., and Gross, T. (2009). Adaptive self-organization in a realistic neural network model. Phys. Rev. E80, 061917. doi: 10.1103/PhysRevE.80.061917

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, K. J., Sorensen, L., Ojemann, J. G., and den Nijs, M. (2009). Power-law scaling in the brain surface electric potential. PLoS Comput. Biol. 5, e1000609. doi: 10.1371/journal.pcbi.1000609

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, J., and Beggs, J. M. (2013). Attaining and maintaining criticality in a neuronal network model. Physica A 392, 1611–1620. doi: 10.1016/j.physa.2012.11.013

CrossRef Full Text | Google Scholar

Petermann, T., Thiagarajan, T., Lebedev, M., and Nicolelis, M. (2009). Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. Proc. Natl. Acad. Sci. U.S.A. 106, 15921–15926. doi: 10.1073/pnas.0904089106

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinto, D., and Ermentrout, G. (2001a). Spatially structured activity in synaptically coupled neuronal networks: lateral inhibition and standing pulses. SIAM J. Appl. Math. 62, 465. doi: 10.1137/S0036139900346465

CrossRef Full Text | Google Scholar

Pinto, D., and Ermentrout, G. (2001b). Spatially structured activity in synaptically coupled neuronal networks: traveling fronts and pulses. SIAM J. Appl. Math. 62, 206–225. doi: 10.1137/S0036139900346453

CrossRef Full Text | Google Scholar

Ribeiro, T., Copelli, M., Caixeta, F., Belchior, H., Chialvo, D., Nicolelis, M., et al. (2010). Spike avalanches exhibit universal dynamics across the sleep-wake cycle. PLoS ONE 5, e14129. doi: 10.1371/journal.pone.0014129

PubMed Abstract | CrossRef Full Text | Google Scholar

Rybarsch, M., and Bornholdt, S. (2014). Avalanches in self-organized critical neural networks: a minimal model for the neural soc universality class. PLoS ONE 9, e93090. doi: 10.1371/journal.pone.0093090

PubMed Abstract | CrossRef Full Text | Google Scholar

Scarpetta, S., and Candia, A. (2013). Neural avalanches at the critical point between replay and non-replay of spatiotemporal patterns. PLoS ONE 8, e64162. doi: 10.1371/journal.pone.0064162

PubMed Abstract | CrossRef Full Text | Google Scholar

Scarpetta, S., Giacco, F., Lombardi, F., and de Candia, A. (2013). Effects of poisson noise in a if model with stdp and spontaneous replay of periodic spatiotemporal patterns, in absence of cue stimulation. Biosystems 112, 258–264. doi: 10.1016/j.biosystems.2013.03.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Tagliazucchi, E., and Chialvo, D. (2013). Brain complexity born out of criticality. AIP Conf. Proc. 1510, 4. doi: 10.1063/1.4776495

CrossRef Full Text | Google Scholar

Troy, W., and Shusterman, V. (2007). Patterns and features of families of traveling waves in large-scale neuronal. SIAM J. Appl. Dyn. Syst. 6, 263–292 doi: 10.1137/06066638X

CrossRef Full Text | Google Scholar

Tsodyks, M., and Markram, H. (1997). The neural code between neocortical pyramidal neurons depends on neurotransmitter release probability. Proc. Natl. Acad. Sci. U.S.A. 94, 719–723. doi: 10.1073/pnas.94.2.719

PubMed Abstract | CrossRef Full Text | Google Scholar

van Rossum, M., Bi, G., and Turrigiano, G. (2000). Stable hebbian learning from spike timing-dependent plasticity. J. Neurosci. 20, 8812–8821. doi: 10.1523/JNEUROSCI.20-23-08812.2000

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeraati, R., Priesemann, V., and Levina, A. (2021). Self-organization toward criticality by synaptic plasticity. Front. Phys. 9, 619661. doi: 10.3389/fphy.2021.619661

PubMed Abstract | CrossRef Full Text | Google Scholar