Network Events on Multiple Space and Time Scales in Cultured Neural Networks and in a Stochastic Rate Model

Cortical networks, in-vitro as well as in-vivo, can spontaneously generate a variety of collective dynamical events such as network spikes, UP and DOWN states, global oscillations, and avalanches. Though each of them has been variously recognized in previous works as expression of the excitability of the cortical tissue and the associated nonlinear dynamics, a unified picture of the determinant factors (dynamical and architectural) is desirable and not yet available. Progress has also been partially hindered by the use of a variety of statistical measures to define the network events of interest. We propose here a common probabilistic definition of network events that, applied to the firing activity of cultured neural networks, highlights the co-occurrence of network spikes, power-law distributed avalanches, and exponentially distributed ‘quasi-orbits’, which offer a third type of collective behavior. A rate model, including synaptic excitation and inhibition with no imposed topology, synaptic short-term depression, and finite-size noise, accounts for all these different, coexisting phenomena. We find that their emergence is largely regulated by the proximity to an oscillatory instability of the dynamics, where the non-linear excitable behavior leads to a self-amplification of activity fluctuations over a wide range of scales in space and time. In this sense, the cultured network dynamics is compatible with an excitation-inhibition balance corresponding to a slightly sub-critical regime. Finally, we propose and test a method to infer the characteristic time of the fatigue process, from the observed time course of the network’s firing rate. Unlike the model, possessing a single fatigue mechanism, the cultured network appears to show multiple time scales, signalling the possible coexistence of different fatigue mechanisms.


Introduction
The spontaneous activity of excitable neuronal networks exhibits a spectrum of dynamic regimes ranging from quasi-linear, small fluctuations close to stationary activity, to dramatic events such as abrupt and transient synchronization.Understanding the underpinnings of such dynamic versatility is important, as different spontaneous modes may imply in general different state-dependent response properties to incoming stimuli and different information processing routes.
Cultured neuronal networks offer a controllable experimental setting to open a window into the network excitability and its dynamics, and have been used intensively for the purpose.
In addition, recent observations in-vitro (and later even in-vivo) revealed a rich structure of network events ('avalanches') that attracted much attention because their spatial and temporal structure exhibited features (power-law distributions) reminiscent of what is observed in a 'critical state' of a physical system (see e.g.[7,8], and [9,10] and references therein).Generically, an avalanche is a cascade of neural activities clustered in time; while there persist ongoing debate on the relation between observed avalanches and whatever 'criticality' may mean for brain dynamics [11], understanding their dynamical origin remains on the agenda.
Quasi-synchronous NS, avalanches and small activity fluctuations are frequently coexisting elements of the network dynamics.Besides, as we will describe in the following, in certain conditions one can recognize network events which are clearly distinct from the mentioned network events, which we name here as 'quasi-orbits'.
The abundant modeling literature on the above dynamical phenomena has been frequently focused on specific aspects of one of them; on the other hand, getting a unified picture is made often difficult by different assumptions on the network's structure and constitutive elements and, importantly, by different methods used to detect the dynamic events of interest.
In the present paper we define a common probabilistic criterium to detect various coexisting dynamic events (NS, avalanches and quasi-orbits) and adopt it to analyze the spontaneous activity recorded from both cultured networks, and a computational rate model.
Most theoretical models accounting for NS are based on an interplay between network self-excitation on one side, and on the other side some fatigue mechanism provoking the extinction of the network spike.For such a mechanism two main options, up to details, have been considered: neural 'spike-frequency adaptation' [3,12] and synaptic 'short-term depression' (STD) [4,5,[13][14][15][16].Despite their differences, both mechanisms share the basic property of generating an activity-dependent self-inhibition in response to the upsurge of activity upon the generation of a NS, the more vigorously, the stronger the NS (i.e. the higher the average firing rate).In this paper, we will mainly focus on STD, stressing the similarities of the two mechanisms, yet not denying their possibly different dynamic implications.
While STD acts as an activity-dependent self-inhibition, the self-excitability of the network depends on the balance between synaptic excitation and inhibition; investigating how such balance, experimentally modifiable through pharmacology, influences the dynamics of spontaneous NSs is interesting and relevant as a step towards the identification of the 'excitability working point' in the experimental preparation.
To study the factors governing the co-occurrence of different network events and their properties we adopt a rate model for the dynamics of the global network activity, that takes into accounts finite-size fluctuations and the synaptic interplay between one excitatory and one inhibitory population, with excitatory synapses being subject to STD.
On purpose we implicitly exclude any spatial topology in the model, which is meant to describe the dynamics of a randomly connected, sparse network, since we intend to expose the exquisite implications of the balance between synaptic excitation and inhibition, and the activity-dependent self-inhibition due to STD.In doing this, we purposely leave out not only the known relevance of a topological organization [9,17,18], but also the role of cliques of neurons which have been proposed to play a pivotal role in the the generation of NS as functional hubs [19], as well as the putative role of 'leader neurons'.
We perform a systematic numerical and analytical study of NSs for varying excitation/inhibition balance.The distance from an oscillatory instability of the mean-field dynamics (in terms of the dominant eigenvalue of the linearized dynamics) largely appears to be the sole variable governing the statistics of the inter-NS intervals, ranging from a very sparse, irregular bursting (coefficient of variation CV ∼ 1) to a sustained, periodic one (CV ∼ 0).The intermediate, weakly synchronized regime (CV ∼ 0.5), in which the experimental cultures are often observed to operate, is found in a neighborhood of the instability that shrinks as the endogenous fluctuations in the network activity become smaller with increasing network size.
Moreover, the model robustly shows the co-presence of avalanches with NS and quasi-orbits.The avalanche sizes are distributed according to a power-law over a wide region of the excitation-inhibition plane, although the crossing of the instability line is signaled by a bump in the large-size tail of the distribution; we compare such distributions and their modulation (as well as the distributions of NS) across the instability line with the experimental results from cortical neuronal cultures; again the results appear to confirm that neuronal cultures operate in close proximity of an instability line.
Taking advantage of the fact that the sizes of both NS and quasi-orbits are found to be significantly correlated with the dynamic variable associated with STD (available synaptic resources) just before the onset of the event, we developed a simple optimization method to infer, from the recorded activity, the characteristic time-scales of the putative fatigue mechanism at work.We first tested the method on the model, and then applied it to in-vitro recordings; we could identify in several cases one or two long time-scales, ranging from few hundreds milliseconds to few seconds.
Weak or no correlations were found instead between avalanche sizes and the STD dynamics; this suggests that avalanches originate from synaptic interaction which amplifies a wide spectrum of small fluctuations, and are mostly ineffective in eliciting a strong self-inhibition.

Experimental data
Cortical neurons were obtained from newborn rats within 24 hours after birth, following standard procedures [2].Briefly, the neurons were plated directly onto a substrate-integrated multielectrode array (MEA).The cells were bathed in MEM supplemented with heat-inactivated horse serum (5%), glutamine (0.5 mM), glucose (20 mM), and gentamycin (10 µg/ml) and were maintained in an atmosphere of 37 • C, 5% CO2/95% air in a tissue culture incubator as well as during the recording phases.The data analyzed here was collected during the third week after plating, thus allowing functional and structural maturation of the neurons.MEAs of 60 Ti/Au/TiN electrodes, 30 µm in diameter, and spaced 200 µm from each other (Multi Channel Systems, Reutlingen, Germany) were used.The insulation layer (silicon nitride) was pretreated with poly-D-lysine.All experiments were conducted under a slow perfusion system with perfusion rates of ∼100 µl/h.A commercial 60-channel amplifier (B-MEA-1060; Multi Channel Systems) with frequency limits of 1-5000 Hz and a gain of 1024× was used.The B-MEA-1060 was connected to MCPPlus variable gain filter amplifiers (Alpha Omega, Nazareth, Israel) for additional amplification.Data was digitized using two parallel 5200a/526 analog-to-digital boards (Microstar Laboratories, Bellevue, WA.)Each channel was sampled at a frequency of 24000 Hz and prepared for analysis using the AlphaMap interface (Alpha Omega.)Thresholds (8× root mean square units; typically in the range of 10-20 µV ) were defined separately for each of the recording channels before the beginning of the experiment.The electrophysiological data is freely accessible for download at https://technion.academia.edu/ShimonMarom/Neural-Network-Data-(Marom's-Lab).

Network rate dynamics
A set of Wilson-Cowan-like equations [20] for the spike-rate of the excitatory (ν E ) and the inhibitory (ν I ) neuronal populations lies at the core of our dynamic mean-field model: where τ E and τ I represent two characteristic times (of the order of few to few tens of ms), and Φ is the gain function of the input currents, I E and I I , that in turn depend on ν E and ν I .We chose Φ to be the transfer function of the leaky integrate-and-fire neuron under the assumptions of Gaussian, uncorrelated input of mean µ and infinitesimal variance σ 2 [21]: where τ V is the membrane time constant, τ ARP is a refractory period, and V rest , V reset , and V thresh are respectively the rest, the post-firing reset, and the firing-threshold membrane potential of the neuron (we assume the membrane resistance R = 1).The values of the parameters are shown in Table 1.The model incorporates the non-instantaneous nature of synaptic transmission in its simplest form, by letting the νs being low-pass filtered by a single synaptic time-scale τ : One can regard the variables νs as the instantaneous firing rates as seen by post-synaptic neurons, after synaptic filtering.The form of Eq. 3 and our choice of τ values (see Table 1) implicitly neglects slow NMDA contributions and is restricted to AMPA and GABA synaptic currents.Thus, the input currents I E and I I in Eq. 1 will be functions of the rates νs through these filtered rates; with reference to Eq. 2, the model assumes the following form for the mean and the variance of the current I E (the expressions for I I are similarly defined): where the n E and n I are the number of neurons in the excitatory and inhibitory population respectively; c is the probability of two neurons being synaptically connected; J EE (J EI ) is the average synaptic efficacy from an excitatory (inhibitory) pre-synaptic neuron to an excitatory one, σ 2 J is the variance of the J-distribution; finally, an external current is assumed in the form of a Poisson train of spikes of rate ν ext driving the neurons in the network with average synaptic efficacy J ext .In Eq. 4 r E (t) (0 < r E < 1) is the fraction of synaptic resources available at time t for the response of an excitatory synapse to a pre-synaptic spike; the evolution of r E evolves according to the following dynamics, which implements the effects of short-term depression (STD) [22,23] into the network dynamics: where 0 < u STD < 1 represents the (constant) fraction of the available synaptic resources consumed by an excitatory postsynaptic potential, and τ STD is the recovery time for the synaptic resources.Finally, for a network of n neurons, we introduce finite-size noise by assuming that the signal the synapses integrate in Eq. 3 is a random process ν n of mean ν; in a time-bin dt, we expect the number of action potentials fired to be a Poisson variable of mean n ν(t) dt; Eq. 3 will thus become: Putting all together, the noisy dynamic mean-field model is described by the following set of (stochastic) differential equations: complemented by Eqs. 2, 4, and 6.
Spike-frequency adaptation (SFA) (only included in simulations of Fig. 9) is introduced by subtracting a term to the instantaneous mean value of the proportional to the instantaneous value of the variable c E , that simply integrates ν n E : with a characteristic time τ SFA .This additional term aims to model an afterhyperpolarization, Ca 2+ -dependent K + current [24,25].In this sense, c E can be interpreted as the cytoplasmic calcium concentration [Ca 2+ ]), whose effects on the network dynamics are controlled by the value of the "conductance" g SFA .Simulations are performed by integrating the stochastic dynamics with a fixed time step dt = 0.25 ms.
In the following, by "spike count" we will mean the quantity ν(t) n dt.

Network events detection
For the detection of network events (NSs, quasi-orbits, and avalanches) we developed a unified approach based on Hidden Markov Models (HMM) [26].
Despite HMM have been widely used for temporal pattern recognition in many different fields, to our knowledge few attempts have been made to use them in the context of interest here [27,28].For the purpose of the present description, we just remind that a HMM is a stochastic system that evolves according to Markov transitions between "hidden", i.e. unobservable, states; at each step of the dynamics the visible output depends probabilistically on the current hidden state.Such models can be naturally adapted to the detection of network events, the observations being the number of detected spikes per time bin, and the underlying hidden states, between which the system spontaneously alternates, being associated with high or low network activity ('network event -no network event').A standard optimization procedure adapts then the HMM to the recorded activity sample by determining the most probable sequence of hidden states given the observations.The two-step method we propose is based on HMM, has no user-defined parameters, and automatically adapts to different conditions.
In the first step, the algorithm finds the parameters of the two-state HMM (one low-activity state, representing the quasi-quiescent periods, and one high-activity state, associated with network events) that best accounts for a given sequence of spike counts -the visible states in the HMM; such fitting is performed through the Baum-Welch algorithm [26].In the second step, the most probable sequence for the two alternating hidden levels, given the sequence of spike counts and the fitted parameters, is found through the Viterbi algorithm.Network events are identified as the periods of dominance of the high activity hidden state.
In order to retain only the most significant events a minimum event duration is imposed; such threshold is self-consistently determined as follows.The second step of the algorithm is applied to a "surrogate" time-series obtained by randomly shuffling the original one; calling x the vector listing the durations of the found events, and x 75 its (estimated) 75th percentile, the threshold is determined as log 0.25 × 10 3 ( x x>x 75 − x 75 ) + x 75 .The motivation behind this procedure lies in a distribution of event durations for surrogate data consistently found to have a roughly exponential large-value tail; the average of such tail-distribution is estimated, and the threshold is then set to a value such as to have a probability P = 10 −3 of finding a surrogate event of duration larger than this value.The procedure was found to give results that are more consistent for different realizations of the surrogate data-set than just taking the 99.9-th percentile of the detected event durations.
For the results presented in Figs. 2 and 3, a further analysis is carried out.The distribution of network event sizes is often found (both in simulations and in experimental data) to be well approximated by a sum of an exponential and a normal distribution (see Results for further discussion): The parameters of the two distributions and their relative weight 0 ≤ p 0 ≤ 1 are estimated by maximizing the log-likelihood on the data.A threshold for the event size is determined as the value having equal probability of being generated by either the exponential or the normal distribution.In Figs. 2 and 3 NSs are defined as events having size larger than this threshold.In those cases in which a threshold smaller than the peak of the normal distribution could not be determined, no threshold was set.
As already remarked, we used essentially the same algorithm for detecting NS/quasi-orbits and avalanches.The only significant difference is that, in the case of avalanches, the emission probability of the low-activity hidden state is kept fixed during the Baum-Welch algorithm to p(n) δ n0 (δ ij is the Kronecker delta; p(n) is the probability of emitting n spikes in a time-bin).Thus the lower state is constrained to a negligible probability of outputting non-zero spike-counts, conforming to the intuition that in between avalanches the network is (almost) completely silent.More precisely, we set p(1) = 10 −6 n , where n is the average number of spikes that the network emits during a time-bin dt.After the modified Baum-Welch first step, avalanches are determined, as above, by applying the Viterbi algorithm; no threshold is applied in this case, neither to the avalanche duration nor to its size.
Simulations and data analysis have been performed using custom-written mixed C++/ MATLAB (version R2013a, Mathworks, Natick, MA) functions and scripts.

Size distribution for quasi-orbits
Consider a generic planar linear dynamics with noise: where A is 2×2 real matrix, and ξ = (ξ(t), 0) is a white noise with ξ(t)ξ(t ) = δ(t−t ).We here assume that the system is close to a Hopf bifurcation; in other words that the matrix A has complex-conjugated eigenvalues λ ± = λ + i λ, with λ < 0 and | λ| λ.By means of a linear transformation, the system can be rewritten as: with σ x and σ y constants determined by the coordinate transformation.Making use of Itō's lemma to write: and summing the previous two equations, we find for the square radius

As long as λ
| λ|, it is physically sound to make the approximation: for 0 ≤ t ≤ T = 2 π/ λ and then to average the variance of the noise over such period to get: in order to rewrite Eq. ( 13) as: Such stochastic differential equation is associated with the Fokker-Planck equation: that admits an exponential distribution as stationary solution: that is, a Rayleigh distribution for w:

Results
In the following, we will study a stochastic firing-rate model and make extensive comparison of its dynamical behavior with the activity of ex-vivo networks of cortical neurons recorded through a 60-channel multielectrode array.
The stochastic firing-rate model consists of two populations of neurons, one excitatory and one inhibitory, interacting through effective synaptic couplings; excitatory synaptic couplings follow a dynamics mimicking shortterm depression (described after the Tsodyks-Markram model, [22]).We adopted the transfer function of the leaky integrate-and-fire neuron subject to white-noise current with drift [21] as the single population input-output function; moreover the activity of each population is made stochastic by adding realistic finite-size noise.Working with a noisy mean field model allows in principle to easily sweep through widely different network sizes and, more importantly, allows us to perform the stability analysis.
To start the exploration that follows, we chose a reference working point where the model's dynamics has a low-rate fixed point (2 − 4 Hz) just on the brink of an oscillatory instability or, in other words, where the dominant eigenvalue λ of the dynamics, linearized around the fixed point, is complex with null real part.The model network (Fig. 1, panel A) shows in proximity of this point a dynamical behavior qualitatively similar, in terms of population spikes, to what is observed in ex-vivo neuronal networks (Fig. 1, panel B).deterministically when the network is again in the condition of generating a new NS.Weak excitability, on the other hand, leads to rare NSs, approaching a Poisson statistics (CV INSI 1), since excitability is so low that fluctuations are essential for recruiting enough activation to elicit a NS, with STD playing little or no role at the ignition time; below an "excitation threshold", NSs disappear.

Excitation-inhibition balance and NS statistics
The solid lines in Fig. 2 are derived from the linearization of the 5dimensional dynamical system (see Eq. 7), and are curves of iso-λ, where λ is the dominant eigenvalue of the Jacobian: λ = 0 Hz (black line, signaling a Hopf bifurcation in the corresponding deterministic system), λ = 3.5 Hz (red line), and λ = −3.5 Hz (blue line).Values of CV found in typical cultured networks are close to model results near the bifurcation line λ = 0 Hz.We observe, furthermore, that such lines roughly follows iso-INSI and iso-CV INSI curves, suggesting that a quasi one-dimensional representation might be extracted.
We show in Fig. 3 INSI (panel A) and CV INSI (panel B) against λ for the same networks (circles) of Fig. 2, and for a set of larger networks (N = 8000 neurons, squares) that are otherwise identical to the first ones, pointwise in the excitation-inhibition plane.The difference in size amounts, for the new, larger networks, to weaker endogenous noise entering the stochastic dynamics of the populations' firing rates.The points are seen to approximately collapse onto lines for both sets of networks, thus confirming λ as the relevant control quantity for INSI and CV INSI .It is seen that, for the smaller networks, INSI and CV INSI depend smoothly on λ, due to finite-size effects smearing the bifurcation.Also note the branch of points (marked by superimposed crosses) for which λ = 0 and then no oscillatory component is present, corresponding to points in the extreme top-left region of the planes in Fig. 2. For the set of larger networks, the dependence of INSI and CV INSI on the λ is much sharper, as expected given the much smaller finite-size effects; this shrinks the available region, around the instability line, allowing for intermediate, more biologically plausible values of CV INSI .
We remark that NSs are highly non-linear and relatively stereotyped events, typical of an excitable non-linear system.The good predictive power of the linear analysis for the statistics of INSI signals that relatively small fluctuations around the system's fixed point, described well by a linear analysis, can ignite a NS.

A spectrum of network events
Our mean-field, finite-size network is a non-linear excitable system which, to the left of the Hopf bifurcation line, and close to it, can express different types of excursions from the otherwise stable fixed point.Large (almost stereotyped for high excitation) NSs are exquisite manifestations of the nonlinear excitable nature of the system, ignited by noise; the distribution of NS size (number of spikes generated during the event) is relatively narrow and symmetric.
Noise can also induce smaller, transient excursions from the fixed point which can be adequately described as quasi-"orbits" in a linear approximation.In fact, noise induces a probability distribution in the maximum distance w from the fixed point reached by the system in each of such events.Such distribution can be computed (see Methods, Eq. 18), and it is a Rayleigh distribution (p(w) ∝ w exp(−α w 2 )).On the other hand, we found a high correlation between w and the duration of the event.Therefore the event size (the 'area' below the firing rate time profile during the excursion from the fixed point) is expected to scale as w 2 , so that it should be exponentially distributed.Fig. 4, panel A, shows the activity of a simulated network (blue line) alongside with detected network events.At the single event level the distinction between a NS and a quasi-orbit cannot be other than statistical.From the best fit for the expected size distribution (an exponential plus a Gaussian distribution) a threshold for the event size can be determined to separate events that are (a-posteriori ) more probably quasi-orbits from the ones that are more probably NSs (for details, see Methods).Following such classification, the green line in Fig. 4, panel A, marks the detection of two NSs (first and third event) and two quasi-orbits (second and fourth event).
As one moves around the excitation-inhibition plane, to the left of the  4. Algorithms for network events detection.Panel A: total network activity from simulation (blue line) with detected NS/quasi-orbits (green line) and avalanches (red line).Four large events (green line) are visible; the first and third are statistically classified as network spikes; the other smaller two as quasi-orbits.Note how network spikes and quasi-orbits are typically included inside a single avalanche.Panel B: a zoom over 0.5 seconds of activity, with discretization time-step 0.25 ms, illustrates avalanches structure (red line).bifurcation line, the two types of events contribute differently to the overall distribution of network event sizes.Qualitatively, the farther from the bifurcation line, the higher the contribution of the small, "quasi-linear" events.This fact can be understood by noting that the average size of such events is expected to scale as 1/| λ|, where λ is the real part of the dominant eigenvalue of the (stable) linearized dynamics (see Methods, Eq. 17).The av-erage size is furthermore expected to scale with the amount of noise affecting the dynamics, thus the contribution of quasi-linear events is also expected to vanish for larger networks.
To take into account the coexistence of the two types of events, we fitted the distribution of burst sizes with the sum of an exponential and a Gaussian distributions (see Methods), which turns out to be acceptably good over wide regions of the excitation and inhibition plane (see Fig. 5, panels A and B).The same fit seems to adequately describe the analogous distributions from experimental data, which show a balanced contribution of the two types of events (see Fig. 5, panels C and D).
As mentioned in the introduction, avalanches are cascades of neural activities clustered in time (see Methods for our operational definition; examples of different methods used in the literature to detect avalanches can be found in [7,[29][30][31]).Fig. 4, panel A and panel B, shows an example of the structure of the detected avalanches (red lines) in the model network.
We extracted avalanches from simulated data, as well as from experimental data.For simulations, we choose data corresponding to three points in the (w exc , w inh ) plane of Fig. 2, with constant w inh = 1 and increasing w exc , with the rightmost falling exactly over the instability line (black solid line in Fig. 2).Three experimental data sets were extracted from different periods of a very long recording of spontaneous activity from a neural culture; each data set is a 40-minute recording.
In Fig. 6 we show (in log-log scale) the distribution of avalanche sizes for the three simulated networks (top row) and the three experimental (bottom row) data sets (blue dots); red lines are power-law fits [32].
From the panels in the top row we see that the distributions are well fitted, over a range of two orders of magnitude, by power-laws with exponents ranging from about 1.5 to about 2.2, consistent with the results found in [7].Note that in the cited paper the algorithm used for avalanche detection is quite different from ours, and the wide range of power-law exponents is related to their dependence on the time-window used to discretize data.In [33] (adopting yet another algorithm for avalanche detection), both the shape of the avalanche distribution and the exponent vary depending on using pharmacology to manipulate synaptic transmission, over a range compatible with our model findings; notably, they find the slope of the power-law to be increasing with the excitability of the network, which is consistent with our modeling results.
Panels B and C of Fig. 6  A broad spectrum of synchronous network events: simulations vs ex-vivo data.Distribution of event sizes for simulations (panel A, corresponding to the point (w exc , w inh ) = (0.86, 0.9) in Fig. 2; panel B, (w exc , w inh ) = (0.86, 0.75)) and two ex-vivo datasets (panel C and D).Note the logarithmic scale on the y-axis.The network of panel A has higher total inhibition and lies farther from the oscillatory instability line than the network of panel B. Theoretical analysis leads to recognize in the broad spectrum of events two families, network spikes and quasi-orbits, that differ in behavior both with respect to noise and in terms of how close the network is to an oscillatory instability; the continuous lines are fits of the theoretical distribution of event sizes, a sum of an exponential (for quasi-orbits) and a Gaussian (for NS) distribution (see Methods).The vertical lines mark the probabilistic threshold separating NS and quasi-orbits.
Avalanche Size p(AS) Panels A-C: mean-field simulations, with fixed inhibition w inh = 1.and increasing excitation (w exc = 0.9, 0.94, 1).The distributions are well fitted by power-laws; panel B and C clearly show the buildup of 'bumps' in the high-size tails, reflecting the increasing contribution from network spikes and quasi-orbits in that region of the distribution.Panels D-F from ex-vivo data, different ∼ 40-minute segments from one long recording; power-laws are again observed, although fitted exponents cover a smaller range; in panels E and F , bumps are visible, similar to model findings.The similarity between the theoretical and experimental distributions could reflect changes of excitatory/inhibitory balance in time in the experimental preparation.Since all the three simulations lay on the left of or just on the bifurcation line (black line in Fig. 2), the shown results are compatible with the experimental network operating in a slightly sub-critical regime.
high-size tails, increasing with the self-excitation of the network; this is understood as reflecting the predominance of a contribution from NS and possibly quasi-orbits in that region of the distribution, on top of a persisting wide spectrum of avalanches.Also this feature is consistent with the findings of [33].
Turning to the plots in the bottom row of Fig. 6, we observe the following features: power-laws are again observed over two decades and more; in panels E and F , bumps are visible, similar to model findings; power-law exponents cover a smaller range just above 2.
While the sequence of plots in two rows (modeling and experiment) clearly shows similar features, we emphasize that experimental data were extracted from a unique long recording, with no intervening pharmacological manipulations affecting synaptic transmission; on the other hand, it has been suggested [34] that a dynamic modulation of the excitatory/inhibitory balance can indeed be observed in long recordings; although our model would be inherently unable to capture such effects, it is tempting to interpret the suggestive similarity between the theoretical and experimental distributions in Fig. 6 as a manifestation of such changes of excitatory/inhibitory balance in time, of which the theoretical distributions would be a 'static' analog.If this interpretation is correct, this means that the experimental preparation operates below, and close, to an oscillatory instability; on the other hand, contrary to NS, the appearance of avalanches does not seem to be exquisitely related to a Hopf bifurcation, rather they seem to generically reflect the nonlinear amplification of spontaneous fluctuations around an almost unstable fixed point -a related point will be mentioned in the next section.We also remark that we obtain power-law distributed avalanches in a (noisy) mean-field rate model, by definition lacking any spatial structure; while the latter could well determine specific (possibly repeating) patterns of activations (as observed in [17]), it is here suggested to be not necessary for power-law distributed avalanches.

Inferring the time-scales
The fatigue mechanism at work (STD in our case) is a key element of the transient network events, in its interplay with the excitability of the system.While the latter can be manipulated through pharmacology, STD itself (or spike frequency adaptation, another -neural -fatigue mechanism) cannot be directly modulated.It is therefore interesting to explore ways to infer relevant properties of such fatigue mechanisms from the experimentally accessible information, i.e. the firing activity of the network.We focus in the following on deriving the effective (activity-dependent) time scale of STD from the sampled firing history.
The starting point is the expectation that the fatigue level just before a NS should affect the strength of the subsequent NS.We therefore measured the correlation between r (fraction of available synaptic resources) and the total number of spikes emitted during the NS (NS size) from simulations.We found that the average value of r just before a NS is an effective predictor of the NS size, the more so as the excitability of the network grows.
Based on the r-NS size correlation, we took the above "experimental" point of view, that only the firing activity ν is directly observable, while r is not experimentally accessible.Furthermore, the success of the linear analysis for the inter-NS interval statistics (due to the NS being a low-threshold very non-linear phenomenon), suggests that without assuming a specific form for the dynamics of the fatigue variable f , we may tentatively adopt for it a generic linear integrator form, of which we want infer the characteristic time-scale τ To do this, first we reconstruct f (t) from ν(t) for a given τ * ; then we set up an optimization procedure to estimate τ * optim , based on the maximization of the (negative) f -NS size correlation (a strategy inspired by a similar principle was adopted in [11]).Fig. 7, panel A, shows an illustrative example of how the correlation peaks around the optimal value.As a reference, the dotted line marks the value below which 95% of the correlations computed from surrogate data fall; surrogate data are obtained by shuffling the values of f at the beginning of each NS.
We remark that in this analysis we use both NS and quasi-orbit events (which are both related to the proximity to a Hopf bifurcation.)This is reasonable since we expect to gain more information about the anti-correlation between f and NS size by including both types of large network events.
For each point of the excitation-inhibition plane, we performed such an optimization procedure.We expect the optimal values τ * optim to depend on the average activity ν of the network for the different points.Indeed, when the fatigue variable follows the Tsodyks-Markram model of STD (which of course was actually the case in the simulations), linearizing the dynamics of r around a fixed point r ( r = 1/(1 + u STD ν τ STD )), r behaves as a simple linear integrator with a time-constant:  19) and the size of the immediately subsequent network spike plotted against the time-scale τ * of the low-pass integrator (continuous line).The correlation presents a clear (negative) peak for an 'optimal' value τ * optim = 0.58 s of the low-pass integrator; such value is interpreted as the effective time-scale of the putative slow self-inhibitory mechanism underlying the statistics of network events -in this case, short-term synaptic depression (STD); as a reference, the dotted line marks the value computed for surrogate data (see text).Panel B: for each point in the (w exc , w inh )-plane (see Fig. 2), τ * optim vs average network activity; the continuous line is the best fit of the theoretical expectation for STD's effective time-scale (Eq.20); the fitted values for the STD parameters τ STD and u STD are consistent with the actual values used in simulation (τ STD = 0.8 s, u STD = 0.2). that depends both on u STD and ν .The optimal τ * values across the excitation-inhibition plane against ν are plotted in Fig. 7, panel B (dots).The solid line is the best fit of τ STD and u STD from Eq. 20, which are consistent with the actual values used in simulations.
This result is suggestive of the possibility of estimating from experiments the time-scale of an otherwise inaccessible fatigue variable, by modeling it as a generic linear integrator, with a "state dependent" time-constant.
Fig. 8 shows the outcome of the same inference procedure for two segments of experimental recordings.The plot in panel A is qualitatively similar to panel A in Fig. 7: although the peak is broader and the maximum correlation (in absolute value) is smaller, the τ * peak is clearly identified and statistically significant (with respect to surrogates, dotted line), thus suggesting a dominant time scale for the putative underlying, unobserved fatigue process.However, Fig. 8, panel B, clearly shows two significant peaks in the correlation plot; it would be natural to interpret this as two fatigue processes, with time scales differing by an order of magnitude, simultaneously active in the considered recording segment.Correlation between low-pass filtered network activity f (see Eq. 19) and the size of the immediately subsequent network spike plotted against the timescale τ * of the low-pass integrator for two experimental datasets (different periods -about 40 minutes each -in a long recording).The plot in panel A is qualitatively similar to the simulation result shown in panel A of Fig. 7: a peak, although broader and of smaller maximum (absolute) value, is clearly identified and statistically significant (with respect to surrogate data, dotted line).Panel B shows two significant peaks in the correlation plot, a possible signature of two concurrently active fatigue processes, with time scales differing by roughly an order of magnitude.
To test the plausibility of this interpretation, we simulated networks with simultaneously active STD and spike-frequency adaptation (SFA, see Methods).Fig. 9 shows the results of time scale inference for two cases sharing the same time scale for STD (800 ms) and time scale of SFA differing by a factor of 2 (τ SFA = 15 and 30 s respectively).In both cases the negative correlation peaks at around τ * 500 ms; this peak is plausibly related to the characteristic time of STD, consistently with Fig. 7.The peaks at higher τ * s, found respectively at 11 and 18 s, roughly preserve the ratio of the corresponding τ SFA values.
This analysis provides preliminary support to the above interpretation τ SFA = 30 s.In both cases the correlation presents a STD-related peak at around τ * 500 ms (τ STD = 800 ms), consistently with Fig. 7.The peaks at higher τ * s, found respectively at 11 and 18 s, roughly preserve the ratio of the corresponding τ SFA values.
of the double peak in Fig. 8, right panel, in terms of two coexisting fatigue processes with different time scales.We also checked to what extent the avalanche sizes were influenced by the immediately preceding amount of available synaptic resources r, and we found weak or no correlations; this further supports the interpretation offered at the end on the previous section, that avalanches are a genuine manifestation of the network excitability which amplifies a wide spectrum of small fluctuations.

Discussion
Several works recently advocated a key role of specific network connectivity topologies in generating 'critical' neural dynamics as manifested in powerlaw distributions of avalanches size and duration (see [18,35]).Also, it has been suggested that 'leader neurons', or selected coalitions of neurons, play a pivotal role in the onset of network events (see e.g.[19,36,37]).While a role of network topology, or heterogeneity in neurons' excitability, is all to be expected, we set out to investigate what repertoire of network events is accessible to a network with the simplest, randomly sparse, connectivity, over a wide range of excitation-inhibition balance, in the presence of STD as an activity-dependent self-inhibition.In the present work we showed that network spikes, avalanches and also large fluctuations we termed 'quasi-orbits' coexist in such networks, with various relative weights and statistical features depending on the excitation-inhibition balance, which we explored extensively, including the role of finite-size noise (irregular synchronous regimes in balanced excitatory-inhibitory networks has been studied in [29]).We remark in passing that the occurrence of quasi-orbits is exquisitely related to the proximity to a Hopf bifurcation for the firing rate dynamics; on the other hand, the occurrence of NS and, presumably, avalanches, does not necessarily require this condition: for instance, NS can occur in the proximity of a saddlenode bifurcation, where the low-high-low activity transitions derive from the existence of two fixed points, the upper one getting destabilized by the fatigue mechanism (see e.g.[38,39]).We also remark that, with respect to the power-law distribution of avalanches, it is now widely recognized that while criticality implies power-law distributions, the converse is not true, which leaves open the problem of understanding what is actually in operation in the neural systems observed experimentally (for a general discussion on the issues involved, see [40]).In the present work, we do not commit ourselves to the issue of whether avalanches could be considered as evidence of Self-Organized Criticality.
On the methodological side, in the present work, we contribute a unified probabilistic model for detection of NS and avalanches, which we believe makes the study of their co-existence easier to develop and more consistent.The more so, given that the very definition of avalanches and the consequent detection algorithms vary significantly across published papers on the subject, which makes it difficult to compose the reported results of different studies in a global picture.
Besides, under an assumption which, in the simulation tests, turned out to be a posteriori a fairly good one, we assumed the (experimentally not accessible) fatigue mechanism to be described by a generic linear integrator, but with a characteristic time scale dependent on the average network activity, and in this way, based only on the observed firing rate, through a simple optimization process we could provide an estimate of the characteristic time of the fatigue mechanism.We demonstrated a preliminary application on experimental data from ex-vivo cultured networks, which gave encouraging results.

Figure 1 .Figure 2 .
Figure 1.Time course of the network firing rate.Panel A: noisy mean-field simulations; panel B: ex-vivo data.Random large excursions of the firing rate (network spikes and quasi-orbits) are clearly visible in both cases.

Figure 3 .
Figure 3. Stability analysis of the linearized dynamics captures most of the variability in the inter-network-spike interval (INSI) statistics.INSI (panel A) and CV INSI (panel B) vs the real part λ of the dominant eigenvalue of the Jacobian of the linearized dynamics, for two networks that are pointwise identical in the excitation-inhibition plane, except for their size (circles: 200 neurons, as in Fig. 2; squares: 8000 neurons).The data points almost collapse on 1-D curves when plotted as functions of λ, leading effectively to a "quasi one-dimensional" representation of the INSI statistics in the (w exc , w inh )-plane.The region in which the INSIs are neither regular (CV INSI ∼ 0) nor completely random (CV INSI 1), as typically observed in experimental data, shrinks for larger networks.The crosses superimposed on the data points mark a null imaginary part λ.
Figure 5.A broad spectrum of synchronous network events: simulations vs ex-vivo data.Distribution of event sizes for simulations (panel A, corresponding to the point (w exc , w inh ) = (0.86, 0.9) in Fig.2; panel B, (w exc , w inh ) = (0.86, 0.75)) and two ex-vivo datasets (panel C and D).Note the logarithmic scale on the y-axis.The network of panel A has higher total inhibition and lies farther from the oscillatory instability line than the network of panel B. Theoretical analysis leads to recognize in the broad spectrum of events two families, network spikes and quasi-orbits, that differ in behavior both with respect to noise and in terms of how close the network is to an oscillatory instability; the continuous lines are fits of the theoretical distribution of event sizes, a sum of an exponential (for quasi-orbits) and a Gaussian (for NS) distribution (see Methods).The vertical lines mark the probabilistic threshold separating NS and quasi-orbits.

Figure 6 .
Figure 6.Avalanche size distribution: simulations vs ex-vivo data.Panels A-C: mean-field simulations, with fixed inhibition w inh = 1.and increasing excitation (w exc = 0.9, 0.94, 1).The distributions are well fitted by power-laws; panel B and C clearly show the buildup of 'bumps' in the high-size tails, reflecting the increasing contribution from network spikes and quasi-orbits in that region of the distribution.Panels D-F from ex-vivo data, different ∼ 40-minute segments from one long recording; power-laws are again observed, although fitted exponents cover a smaller range; in panels E and F , bumps are visible, similar to model findings.The similarity between the theoretical and experimental distributions could reflect changes of excitatory/inhibitory balance in time in the experimental preparation.Since all the three simulations lay on the left of or just on the bifurcation line (black line in Fig.2), the shown results are compatible with the experimental network operating in a slightly sub-critical regime.

Figure 7 .
Figure 7. Slow time-scales inference procedure: test on simulation data.Panel A: correlation between low-pass filtered network activity f (see Eq.19) and the size of the immediately subsequent network spike plotted against the time-scale τ * of the low-pass integrator (continuous line).The correlation presents a clear (negative) peak for an 'optimal' value τ * optim = 0.58 s of the low-pass integrator; such value is interpreted as the effective time-scale of the putative slow self-inhibitory mechanism underlying the statistics of network events -in this case, short-term synaptic depression (STD); as a reference, the dotted line marks the value computed for surrogate data (see text).Panel B: for each point in the (w exc , w inh )-plane (see Fig.2), τ * optim vs average network activity; the continuous line is the best fit of the theoretical expectation for STD's effective time-scale (Eq.20); the fitted values for the STD parameters τ STD and u STD are consistent with the actual values used in simulation (τ STD = 0.8 s, u STD = 0.2).

Figure 8 .
Figure 8. Slow time-scales inference procedure on ex-vivo data.Correlation between low-pass filtered network activity f (see Eq.19) and the size of the immediately subsequent network spike plotted against the timescale τ * of the low-pass integrator for two experimental datasets (different periods -about 40 minutes each -in a long recording).The plot in panel A is qualitatively similar to the simulation result shown in panel A of Fig.7: a peak, although broader and of smaller maximum (absolute) value, is clearly identified and statistically significant (with respect to surrogate data, dotted line).Panel B shows two significant peaks in the correlation plot, a possible signature of two concurrently active fatigue processes, with time scales differing by roughly an order of magnitude.

Figure 9 .
Figure 9. Slow time-scales inference procedure on simulation data with STD and spike-frequency adaptation.Correlation between lowpass filtered network activity f (see Eq.19) and the size of the immediately subsequent network spike plotted against the time-scale τ * of the low-pass integrator.In this case, the mean-field model includes, besides short-term depression (STD), a mechanism mimicking spike-frequency adaptation.Panel A: spike-frequency adaptation with characteristic time τ SFA = 15 s.Panel B: τ SFA = 30 s.In both cases the correlation presents a STD-related peak at around τ * 500 ms (τ STD = 800 ms), consistently with Fig.7.The peaks at higher τ * s, found respectively at 11 and 18 s, roughly preserve the ratio of the corresponding τ SFA values.