{'fraction_non_alphanumeric': 0.03409806567701305, 'fraction_numerical': 0.009671614934772829, 'mean_word_length': 4.909356725146199, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
New hallmarks of criticality in recurrent neural networks

A rigorous understanding of brain dynamics and function requires a conceptual bridge between multiple levels of organization, including neural spiking and network-level population activity. Mounting evidence suggests that neural networks of cerebral cortex operate at criticality. How operating near this network state impacts the variability of neuronal spiking is largely unknown. Here we show in a computational model that two prevalent features of cortical single-neuron activity, irregular spiking and the decline of response variability at stimulus onset, are both emergent properties of a recurrent network operating near criticality. Importantly, our work reveals that the relation between the irregularity of spiking and the number of input connections to a neuron, i.e., the in-degree, is maximized at criticality. Our findings establish criticality as a unifying principle for the variability of single-neuron spiking and the collective behavior of recurrent circuits in cerebral cortex.


I. Introduction
Linking the evolutionary-derived multi-scale organization of the brain to neural dynamics and computation represents a major challenge in systems neuroscience [1]. Among the vast spectrum of spatial and temporal scales of brain activity, two experimentally accessible levels of brain organization are (i) the single-neuron spiking and (ii) the population activity of the network in which the neurons are embedded to various degrees. Single-neuron spiking in cerebral cortex is characterized by statistical properties, such as irregular spiking [2] and reduced variability during sensory stimulation [3]. Population activity is characterized by complex spatiotemporal activity, including scale-free activity [4], which is predicted to occur for a network state near criticality [5]. These observations at two adjacent levels of brain organization raise the question, to what extent the network state controls the variability of single-neuron spiking? Irregular spiking, defined as the mean coefficient of variation (CV) being larger than one, is known as one of the most widespread features of cortical activity in vivo [2,6]. Theoretical studies have linked the prevalence of irregular spiking to the fluctuations of synaptic inputs at the sub-threshold regime [7,8]. Despite the popularity of this hypothesis, here we provide an alternative scenario, by which the onset of irregular spiking can simply emerge at the transition between two phases of order and chaos. Further, using a computational model we show that even at the presence of other mechanisms for irregular spiking (CV > 1), neuronal activity shows maximum irregularity when the network resides at criticality. Moreover, in addition to irregular spiking, we show that criticality also gives rise to pronounced decline in neural variability after the onset of stimulus, as well as maximum correlation between a number of single-neuron properties such as their CV's and firing rates. Our findings provide us with robust measures of critical dynamics that at the same time could further our understanding about the implications of criticality for brain dynamics and function.

II. Results
To explore the impact of critical dynamics on the statistics of single-neuron spiking, we used a model network consisting of excitatory binary probabilistic neurons with sparse connectivity and external inputs (see fig.1a and Methods). This model permits the use of the maximum eigenvalue λ of the transition probability matrix P ij as a control parameter to tune the network near the critical point [9,10]. Such tuning results in characteristic avalanche size distributions for the subcritical (λ < 1), critical (λ = 1), and supercritical (λ > 1) regime ( fig.1  b-d).
To investigate the hypothesized impact of the network state on the statistics of neuronal spiking, we simulated the network activity for different values of λ and quantified the single-neuron spiking statistics using the CV of the inter-spike-interval distributions ( fig.1 e). The CV is defined as the ratio of the standard deviation and the mean of the inter-spike-interval distribution for a given neuron. The irregular spiking is basically characterized by CV > 1, whereas CV = 1 is considered as Poisson spiking. It turned out that at the presence of a constant slow drive (see Methods) the CV values distributed around a mean greater than one, thus indicating irregular spiking. In contrast, small deviations of the network state towards either the subcritical or the supercritical regime resulted in CV values distributed around a mean of 1 or less. In summary, when tuning the network from the subcritical to the supercritical state, the mean CV peaked near λ = 1 ( fig.1f); with increasing network size, the peak moved closer towards λ = 1 and becomes narrower as well. This observation suggests that, at the large-size limit, the irregular spiking is an emergent property of recurrent networks operating at the critical point (λ = 1). In other words, the scale-free fluctuations in network activity at criticality translate into irregular single-neuron spiking.
The above results could be confounded by the finite network size. Importantly, the width of the mean CV as a function of λ decreased with increasing network size . (e) Inter-spike-interval CV distributions of simulated spike trains from a network of N = 5000 neurons with 3% connectivity and external inputs η = 1/(5N ) to all neurons. Simulations were conducted for the subcritical (λ = 0.9, blue), critical (λ = 1.02, red), and supercritical (λ = 1.06, green) network states. (f ) The average CV as a function of the maximum eigenvalue λ of the transition probability matrix P ij for three network sizes. Connectivity was 3% and external input was 1/(5N ) to all neurons. The curves were based on 13 values of λ within the range shown and were smoothed using Matlab spline. (g) The average population coupling and the change in mean response as a function of λ for a network of N = 5000 neurons with 3% connectivity and external inputs of strength η = 1/(5N ) applied to all neurons. To compute the change in mean response, we increased the external input strength by a factor of 10 half-way through the simulation, i.e., from η = 1/(5N ) to η = 2/N . and the peak location moved closer to λ = 1 ( fig.1f). To further test whether the deviation of the peak from λ = 1 is indeed due to the finite size effect, we looked at two major characteristics of criticality, namely maximum correlations and mean response. In order to see how the average correlation among neurons behaves in terms of the control parameter λ, we computed a commonly-used measures of coordinated network activity that is known to be maximized at criticality. One such quantity is the average population coupling, which represents a measure of the overall level of correlated fluctuations within the network [11] (see Methods).
Comparing the CV and the average population coupling reveals that they both peak at the same point (close to λ = 1), also similar to the change in mean response to an increase in external input. This demonstrates that the peak for CV coincides with the critical point, which is characterized by maximum average population coupling and change in mean response. Therefore, maximum CV could be regarded as a hallmark of critical dynamics for recurrent neural networks.
The observed onset of irregular spiking near criticality ( fig.1f) has an intuitive explanation in the extreme limits of connection strength. In the limit of weak connections (λ < 1), spiking is largely driven by the Poisson external input alone, thus resulting in Poisson like spiking. On the other hand, neurons become mostly active leading to more regular spiking in the limit of very strong connections (λ > 1). However, it is at criticality (λ ≈ 1) that the scale-free fluctuations in network activity translate to single-neuron irregular spiking ( CV > 1). Furthermore, although there can be alternative mechanisms for irregular spiking, such as synchronous inputs (see discussion), we show that the impact of the network state on the statistics of single-neuron spiking is largely robust with respect to structured external inputs; the irregularity of spiking is always maximized near network criticality even at the presence of other sources of irregular spiking ( fig.S1,fig.S2,fig.S3). The observed impact of the network state on the statistics of single-neuron spiking ( fig.1e,f) raised the question to what extent this networkto-neuron impact is regulated by a neuron's in-degree, i.e., the number of input connections to a neuron. To address this question, we took advantage of the distribution of in-degrees provided by a model network with sparse and random connectivity.
We found that near criticality (λ ≈ 1), a neuron's CV tended to increase with increasing in-degree ( fig.2a). Importantly, however, this correlation between a neuron's activity statistics (CV) and its connectivity (in-degree) changed drastically when tuning the network state away from criticality. In the subcritical regime only weak correlation was found and in the supercritical regime the relation reversed, i.e., a neuron's CV tended to decrease with increasing in-degree ( fig.2b). In general, when tuning the network state through the critical regime, the relation between a neuron's CV and its in-degree transforms from small correlation in the subcritical regime (λ < 1), to maximized correlation at the critical point (λ ≈ 1), and to anti-correlation in the supercritical The relationship that solely near criticality is the irregular spiking maximally reflective of the neuron's indegree, suggests a novel measure to test the criticality hypothesis in electrophysiological experiments. It is of practical importance that in a network dominated by excitatory neurons, a neuron's firing rate scales with its in-degree and that this relation is independent of the network state. Thus, spike train recordings from a population of neurons can be informative about the network state: a maximum correlation between a neuron's CV and its firing rate is indicative of a critical network state, whereas weaker correlation or anti-correlation is indicative of the subcritical or supercritical network state, respectively ( fig.2d-f).
In comparison, the relation between a neuron's "population coupling" (see Methods) and its in-degree is less decisive about the network state, as population coupling increases with a neuron's in degree for all three network states ( fig.2g, h). However, this relation is also dominant for the critical network state ( fig.2i). Consistent with these model results, recent experimental data from in vivo recordings showed increased population coupling with increasing synaptic inputs [11]. However, unless the network state is manipulated and a maximum in the correlation between the population coupling and the in-degree is determined, that data set is not informative about the network state.
In addition to the irregular spiking evaluated above, a decline of response variability at stimulus onset is a prominent cortical phenomenon [3]. The important question that how an external input suppresses the variability of ongoing activity has remained intriguing [12,13]. To address this question in the context of the network state, we simulated spiking activity of the model network near criticality (λ ≈ 1) in response to repeated increases in the external input. We quantified the across-trial firing rate variability in terms of the Fano factor, which is a measure of the across-trial variability in the number of spikes in relation to the trial-averaged mean number of spikes during a given window of time (see Methods). The Fano factor has a value of approximately 1 for repeated Poisson spike trains. Recordings from cortical neurons show a Fano factor above 1 for ongoing activity and a significant drop during external stimulation [3]. Our simulations of a model network in the critical state reproduced this drop in Fano factor from a high value during ongoing activity to a lower value (still above 1) during external stimulation ( fig.3a). In contrast, in the subcritical and the supercrit-ical network state, the Fano factor remained unchanged when switching from ongoing to evoked activity. This result demonstrates that the experimentally observed drop in Fano factor could naturally emerge as a characteristic feature of the critical network state. The result is qualitatively robust with respect to the details of the external inputs ( fig.S4).
Since the Fano factor is defined as a ratio of the variance divided by the mean, the declining Fano factor at criticality could result trivially from an increased mean spike count accompanied by a weak dependence of the variance on the network state during external stimulation. Evaluating the changes in the spike count and the variance separately ruled out this possibility ( fig.3b). The "mean spike count" and the variance increased with external stimulation for all three network states. However, only near criticality, when the network is maximally sensitive to external stimuli, did the mean spike count increase more drastically, thus yielding a reduced Fano factor. In conclusion, our model simulation suggests that the experimentally observed drop in Fano factor after stimulus onset [3] can be the result of the cortical circuit operating at criticality.

III. Discussion
A wealth of evidence indicates that neural networks of cerebral cortex operate at criticality [4,14,15]. Here, we showed that (i) irregular spiking [2,6,16], (ii) its relation to the neuron's in-degree, and (iii) the decline of response variability at stimulus onset [3], are all emergent properties of a recurrent network operating at criticality. A number of separate dynamical, biophysical, and structural mechanisms have been proposed to generate the observed irregular spiking [12,[17][18][19][20][21]. The significance of our work resides in part in establishing criticality as one unifying principle for both the collective behavior of cortical circuits and the statistics of singleneuron spiking. Indeed, experimental evidence for the predicted coexistence of irregular spiking and criticality has recently been provided in recordings of ongoing cortical activity in vivo [6]. While the coexistence of irregular spiking and power-law avalanche size distributions has been demonstrated in more complex model networks [22][23][24][25][26], our work extends qualitatively beyond these important earlier studies in three fundamental dimensions. First, the choice of a network of excitatory probabilistic integrate-and-fire model neurons allows the precise analytic evaluation of the network state via a single control parameter λ, i.e., the maximum eigenvalue of the transfer probability matrix.
This approach avoids the need to rely on the cumbersome and less precise avalanche analysis to evaluate the network state. Second, our discovery that the relation between the irregularity of spiking and the neuron's in-degree/firing rate is maximized at criticality provides an important new measure of criticality. Furthermore, this relation establishes a valuable conceptual link between criticality and the important field of network theory, where a node's in-degree is a basic system parameter. Third, as the observed decline in response variability is regarded as an essential mechanism to enhance response fidelity to stimuli [3], our discovery of its relation to network criticality offers a starting point toward unraveling the possible roles of critical dynamics in neural coding. In conclusion, it will be interesting to see to what extent the presented findings will generate a paradigm shift in the study of criticality of neural systems: our results build a much-needed bridge between critical dynamics and neural coding, and provide novel and robust measures to test the criticality hypothesis itself.

IV. Methods
We simulated a model network consisting of excitatory binary probabilistic model neurons with sparse connectivity and external inputs. Network size ranged from 5000 to 20000 model neurons. The strength of the connection from neuron j to neuron i is quantified in terms of the transition probability Pij , which is the probability that a spike in neuron j causes a spike in neuron i in the next simulation time step. For a network of N neurons and an average connectivity K, each neuron is connected to N − 1 other neurons with probability K/N . For each (on average) K(N − 1) connections among neurons a Pij is assigned by drawing a random number from uniform distribution in the interval [0 2 K ]. With sufficiently large this yields a network with a normally distributed connectivity with average K and a transition matrix Pij with maximum (absolute value) eigenvalue of 1. In order to deviate the network from the critical point we can simply multiply Pij by a factor smaller(greater) than one. The binary state Xi(t) of neuron i denotes whether the model neuron spikes (Xi(t) = 1 or does not spike (Xi(t) = 0 at time t. At each time step, the state of all neurons were updated synchronously according to the following update rule: where ηi(t) of neuron i quantifies the probability of that neuron spiking only due to external input, ξ(t) is a random number in [01] drawn from a uniform distribution, and Θ is the step function. In addition to the update rule, a refractory period of two time-steps was implemented. The external input ηi(t) was chosen to be smaller than the transition probability Pij , which itself was small for large networks, Pij ∼ 1/N . Because of the weak external inputs, we employed the approximation 1 − η ≈ 1 in the above update rule. The maximum eigenvalue λ of the transition probability matrix Pij describes the network state at the thermodynamic limit: λ < 1 denotes subcritical regime, λ ≈ 1 denotes the near critical regime and λ > 1 denotes the supercritical regime. However, for finitesized networks we evaluate the exact critical point based on the peak of average population coupling (see below). The Pij values were drawn from a uniform distribution and then scaled by constant to reach the desired maximum eigenvalue λ. The constant external input ηi(t) was either constant.
The external input ηi(t) was either constant or was modeled as a binary Poisson process followed by smoothing with a Gaussian filter with a width of 20 time steps and multiplied by an amplitude factor η0 between 0 and 1. The synchronous input (fig.S1) was simulated with replicating a single binary Poisson process smoothed with a Gaussian filter with a width of 100 time-steps. In order to apply some variability among the stimuli received by different neurons each smoothed Poisson process was multiplied by a factor of η0 + 0.2ǫ where ǫ was drawn from a normal distribution.
We analyzed the simulated spike trains with respect to five complementary statistical measures: neuronal avalanches, coefficient of variation, population coupling, change in mean response, and Fano factor. Following commonly-used procedures [27], a neuronal avalanche was defined as an episode of continuous (each time step) network spiking, framed by time steps of no spikes. Avalanche size was taken as the number of spiking neurons. The single-neuron spike train variability was quantified using the coefficient of variation CV ≡ σisi/µisi, defined as the ratio of the standard deviation σisi and the mean µisi of the inter-spike-interval (ISI) distribution for a given neuron. We managed simulation times to be sufficiently long to ensure stable CV values. The coordination within the network was quantified using the population coupling , which is defined as the zero-lag cross-correlation between the spike train Xi(t) of neuron i and the remaining network activity Ni(t) = j =i Xj (t) from all other spike trains, i.e., where the angled brackets indicate a time average [11]. Averaging the population coupling over many neurons within a large network yields the average population coupling , which represents a measure of the overall level of coordination within the network. The change in mean response was computed as the difference in the mean spike counts for external inputs of η = 1/(5N ) and η = 2/N . When quantifying the spike train variability in the context of repeated experimental situations, it is convenient to use the Fano factor, which, for a given time window and repeated trials, is defined as the ratio of the variance and mean in the number of spikes. We chose a window of 200 time-steps (sliding in increments of 20 time steps) and repeated trials 2000 times. We calculated the average Fano factor from 60 randomly sampled model neurons. Fig. S1: Irregular spiking at criticality in a recurrent network with synchronous external inputs. (a) The temporal structure and strength of the external input η(t) to 10% of the neurons in a recurrent model network of 5000 neurons and 1% connectivity. The external input η(t) was generated from Poisson pulses of rate 10/N , smoothed by a Gaussian filter of width 100 time-steps and amplitude of 0.2(1 + ǫ), where ǫ is drawn from a normal distribution (see Methods). This synchronous input was added to a background constant external input of 1/ (10N ). (b) Inter-spike-interval CV distributions of simulated spike trains for the subcritical (λ = 0.95, blue), critical (λ = 1.09, red), and supercritical (λ = 1.09, green) network state. At the critical regime the spike trains show highest irregularity, which is indicated by the peak of the CV distribution located near 1.3. (c, d) The inter-spike-interval CVs from simulated spike trains versus the neuron's in degree (c) and its normalized rate (d) for the three network states. (e, f ) The population coupling from simulated spike trains versus the neuron's in degree (e) and its normalized rate (f) for the three network states. (g) The population coupling vs a neuron's CV for the three network states. (h) The Spearman correlation coefficients between CV and in-degree (rate), population coupling and in-degree (rate), population coupling and CV are all maximized near criticality.