Neurobiologically realistic determinants of self-organized criticality in networks of spiking neurons

Abstract

Self-organized criticality refers to the spontaneous emergence of self-similar dynamics in complex systems poised between order and randomness. The presence of self-organized critical dynamics in the brain is theoretically appealing and is supported by recent neurophysiological studies. Despite this, the neurobiological determinants of these dynamics have not been previously sought. Here, we systematically examined the influence of such determinants in hierarchically modular networks of leaky integrate-and-fire neurons with spike-timing-dependent synaptic plasticity and axonal conduction delays. We characterized emergent dynamics in our networks by distributions of active neuronal ensemble modules (neuronal avalanches) and rigorously assessed these distributions for power-law scaling. We found that spike-timing-dependent synaptic plasticity enabled a rapid phase transition from random subcritical dynamics to ordered supercritical dynamics. Importantly, modular connectivity and low wiring cost broadened this transition, and enabled a regime indicative of self-organized criticality. The regime only occurred when modular connectivity, low wiring cost and synaptic plasticity were simultaneously present, and the regime was most evident when between-module connection density scaled as a power-law. The regime was robust to variations in other neurobiologically relevant parameters and favored systems with low external drive and strong internal interactions. Increases in system size and connectivity facilitated internal interactions, permitting reductions in external drive and facilitating convergence of postsynaptic-response magnitude and synaptic-plasticity learning rate parameter values towards neurobiologically realistic levels. We hence infer a novel association between self-organized critical neuronal dynamics and several neurobiologically realistic features of structural connectivity. The central role of these features in our model may reflect their importance for neuronal information processing.

Figure 1. Hierarchically modular connectivity and spike-timing-dependent plasticity.

(a) An illustrative connectivity matrix of a hierarchical modular network. This network consists of sixteen -neuron modules, organized into four hierarchical levels. Squares in the connectivity matrix outline the nesting of hierarchical level (small orange squares) inside hierarchical level (large yellow squares). In the present study we considered networks of hierarchical levels and neurons. (b) An illustration of the synaptic plasticity rule used in the study. (c) Weight frequency distributions for the STDP rule with soft bounds (used in most simulations). (d) Weight frequency distributions for the STDP rule with soft bounds and reduced learning rate (used in some simulations). (e) Weight frequency distributions for the STDP rule with hard bounds (used in some simulations). Error bars represent the standard deviation.

Figure 2. Properties of hierarchically modular connectivity.

(a) Power-law, exponential and linear density scaling functions of the networks used in this study. (b) Dimensionless measures of wiring cost for each density scaling function in (a). The wiring cost was computed by equating synaptic cost with hierarchical-level number, and averaging the cost over all synapses. Hence, synapses in hierarchical level were assigned a cost of , synapses in hierarchical level were assigned a cost of , etc. Approximate values of the network (c) clustering coefficient and (d) characteristic path length across a range of randomizations of each hierarchical topology. Color-coding is the same as in (a).

Figure 3. Relationship between neuron spikes and module spikes.

(a) An illustrative scatter plot of the numbers of neuron spikes and module spikes at sampled intervals. Integer spike numbers were jittered by the addition of uniformly distributed random numbers between and . (b) Correlation coefficients between numbers of neuron spikes and numbers of module spikes as a function of network randomization. Correlations were computed from five-minute spike rasters. Error bars represent the standard error of the mean from simulations.

Figure 4. Phase transition from subcritical to supercritical network dynamics.

(a) Illustrative ordered (left), intermediate (center) and random (right) nonhierarchical connectivity matrices. Nonhierarchical networks are characterized by a homogeneous density of between-module excitatory connections (Figure 2a, blue lines) (b) Cumulative probability distributions of avalanche sizes, avalanche durations and inter-avalanche intervals emergent on nonhierarchical networks. Subcritical dynamics (concave distributions) correspond to less randomized networks, supercritical dynamics (convex distributions) correspond to more randomized networks, while critical dynamics (linear-like distributions, in bold) occur between these two extremes. Gray, pink and light blue distributions correspond to random networks (concave distributions) and ordered networks (convex distributions). (c) An illustrative module spike raster of critical dynamics.

Figure 5. Relationship between spike-timing-dependent synaptic plasticity and network dynamics.

(a) Fluctuations of within-module synaptic weights over a 20 second period. Synaptic weights were rank-ordered and assigned a rank-specific color at the first sampled time step. At subsequent steps, weights were re-ranked and therefore reordered, but the color-coding remained fixed. The mixing of colors hence represents fluctuations in rank positions. Stable synaptic weight distributions allowed the inference of weight fluctuations from these rank fluctuations. Weights were sampled at intervals. (b) Illustrative fluctuations in the number of module spikes (top) and in the mean within-module excitatory synaptic weights (bottom), recorded over a minute period from a single module. Module spikes were binned at second intervals, and synaptic weights were sampled at second intervals. (c) Cumulative probability distributions of avalanche sizes, avalanche durations and inter-avalanche intervals and (d) an illustrative module spike raster of dynamics emergent on nonhierarchical network topologies in Figure 4a, with frozen synaptic weights.

Figure 6. Relationship between hierarchical modularity, wiring cost and network dynamics.

Statistical significance of power-law distributions of avalanche sizes (black) and durations (red) as a function of network randomization, for (a) nonhierarchical and hierarchical (b) power-law, (c) exponential and (d) linear density scaling functions. Gray lines show the threshold for power-law scaling. Error bars represent the standard error of the mean from simulations.

Figure 7. Illustrative power-law distributions of avalanche sizes and durations.

Cumulative probability distributions of avalanche (a) sizes and (b) durations for the optimal power-law (), exponential () and linear () density scaling functions. Gray and pink distributions correspond to random networks (concave distributions) and ordered networks (convex distributions). (c) Mean exponents of statistically significant power-law distributions of avalanche sizes (black) and durations (red), as a function of network randomization. Error bars represent the standard error of the mean from simulations.

Figure 8. Role of modularity and low wiring cost in emergence of self-organized critical network dynamics.

(a) An illustrative lattice connectivity matrix with the optimal power-law () density scaling function. Presumed network modules are shown in yellow, while inhibitory synapses are shown in red. Note the presence of significant numbers of intermodule inhibitory synapses. (b) A variant of the connectivity matrix in (a), with modularity of inhibitory neurons. (c) Statistical significance of power-law distributions of avalanche sizes (black) and durations (red) for the lattice (top), the lattice with restored module spike rate (center), and the modified lattice, in which inhibitory neurons were arranged into explicit modules. Gray lines show the threshold for power-law scaling. Error bars represent the standard error of the mean from simulations.

Figure 9. Robustness of self-organized critical network dynamics.

Statistical significance of power-law distributions of avalanche sizes (black) and durations (red) as a function of network randomization for the optimal power-law () density scaling function associated with (a) changes in external current (default ), (b) changes in conduction delays (default delays are uniformly distributed between and ), (c) weaker postsynaptic response (default ) and slower STDP learning rate (default ) (d) doubling of module size to neurons and network size to neurons and reductions in postsynaptic response, STDP learning rate and external current (e) changes from soft to hard STDP weight bounds, and (f) removal of inhibitory synapses. Gray lines show the threshold for power-law scaling. Error bars represent the standard error of the mean from simulations.