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.


. 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 scalefree 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(Levina et al., , 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 −1 ) 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. (2007Levina et al. ( , 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 andBeggs (2013) anddi 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.  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 .
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 shortterm 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 . /fncom. .
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.
. Materials and methods . . 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 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 (2) 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 0 will change the inhibitory or excitatory ion channel conductance of the postsynaptic neuron for t > t 0 according to: 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.

. . 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 conductancebased currents introduced in this section. Each excitatory neuron in the population is randomly connected to k EE = N Exc 100 = 200 excitatory and k EI = k EE 4 inhibitory neurons and each inhibitory neuron is connected to k IE = k EE and k II = 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( k EE ) spikes are adequate for firing. Assuming homogeneity in the population as we have discussed in the introduction we can build a meanfield equation for the excitatory and inhibitory population in this sparse network, assuming each neuron receives input with the same statistics.
. . 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 .
/fncom. . At equilibrium, population rates satisfy a system of fixed point equations of the form: where c xy = ck xy w xy (V R y − 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 0 is a constant that depends on V rest , V th , the maximal rates and the SD of the input. V R y 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: τ 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 ρ E Ext 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: where d is a constant equal to . 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 .
/fncom. . 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.

. . 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 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 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 E max and w I max . We use the model of Tsodyks and Markram (1997) for shortterm 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: .

Results
. . Long term synaptic plasticity by STDP tunes synaptic weights close to the balanced state A typical neuron in the cortex has 10 3 − 10 4 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 .
/fncom. . 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: (9) 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, (10) 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 : 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 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 From the above equations, it is straightforward to see when We take the proportion of inhibitory synapses to excitatory synapses to be equal for both excitatory and inhibitory neurons, i.e., k EI k EE = 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: Figure 3 shows that STDP tunes the weights to this balanced inhibition dominated state in which W = C EE C II − C EI C 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 c II c EI = k II ρ st I k EI ρ st E . Comparing these last two equations, we arrive at k II c st EE = k EI c st II . 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.
. . 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: The rate equations for the EI population are of the form: 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 .

FIGURE
Output excitatory rates as a function of W EE and the corresponding graphs for the average synaptic e cacy 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 , two di erent 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 , by decreasing q after the Hopf bifurcation of the low firing rate fixed point, an oscillatory solution for (u, ρ out ) emerges.
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 ρ * E is: 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 slowfast assumption) can be written as 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.
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 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.  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.
. /fncom. .  Another way that STP can cause a switch between two distinct firing states is in the EI population which possesses bistability. In this case, a change of u can make each of the bistable 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 0 EE . 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 0 EE . 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 .
/fncom. .  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. In this particular case, necessary conditions for up to down transitions are: 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.

. . 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) = E(x, t) I(x.t) , D as the diffusion matrix and f : R 2 → R 2 as the gain function, one can write down the reaction-diffusion rates dynamics in the following form: The ODE part of this equation is the dynamic of a single EI population. The corresponding low firing fixed point (E 0 , I 0 ) 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: 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 pulsecoupled 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.
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 To model the finite-size stochastic effects, we need . /fncom. . 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: which has the solution: with c = −α 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: 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 semilinear 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 ) can be written as γ = c EE E − c EI I ≈ 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 1 2 where we can use the following approximations for the gain function derivatives: 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 γ : 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) = σ uτ S(x, t) agrees with the Langevin description of directed percolation which is of the following form with new transformed coefficients: 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 0 EE with STP. In all these cases, the operating point of the system is close to the bifurcation point of the quiescent state.
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 .

FIGURE
Intersection of curves given by equations ( ) and ( ). Setting to a value close to α by STDP and su cient 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.
dynamics of equation (16): 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 On the other hand, in the active phase, the stationary homogeneous rate is : 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 = 1 qτ STP ( 0 α − 1).
Long-term synaptic plasticity tunes 0 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 longterm 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.

. 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 shortterm 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. Shortterm 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 α ∝ 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 . /fncom. . 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 andCowan, 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.