Critical synchronization and 1/f noise in inhibitory/excitatory rich-club neural networks

In recent years, diverse studies have reported that different brain regions, which are internally densely connected, are also highly connected to each other. This configuration seems to play a key role in integrating and interchanging information between brain areas. Also, changes in the rich-club connectivity and the shift from inhibitory to excitatory behavior of hub neurons have been associated with several diseases. However, there is not a clear understanding about the role of the proportion of inhibitory/excitatory hub neurons, the dynamic consequences of rich-club disconnection, and hub inhibitory/excitatory shifts. Here, we study the synchronization and temporal correlations in the neural Izhikevich model, which comprises excitatory and inhibitory neurons located in a scale-free hierarchical network with rich-club connectivity. We evaluated the temporal autocorrelations and global synchronization dynamics displayed by the system in terms of rich-club connectivity and hub inhibitory/excitatory population. We evaluated the synchrony between pairs of sets of neurons by means of the global lability synchronization, based on the rate of change in the total number of synchronized signals. The results show that for a wide range of excitatory/inhibitory hub ratios the network displays 1/f dynamics with critical synchronization that is concordant with numerous health brain registers, while a network configuration with a vast majority of excitatory hubs mostly exhibits short-term autocorrelations with numerous large avalanches. Furthermore, rich-club connectivity promotes the increase of the global lability of synchrony and the temporal persistence of the system.


Methods
Hierarchical network and rich-club organization. The neurons are located in a hierarchical scale-free network proposed by Ravasz and Barabási 81 . The connectivity degree distribution P(k) of the model follows a negative power-law function ∼ γ − P k k ( ) , with γ = . 2 1. The first step consists in constructing a cluster of five linked nodes (a complete network); creating four replicas, and finally connecting four nodes of each replica cluster to one node in the first cluster (hub node); this results in a network of 25 nodes including a hub. The second step consists in replicating the first step four more times, and connecting the resulting 16 peripheral nodes to the hub node proposed in step one; the output consists in a network with 125 nodes and 5 hubs (Fig. 1a). The algorithm presented can be repeated indefinitely, where each step increases the number of nodes by a factor of 5. In this study, 5 network replicas, similar to the one created in step 2, were used to form a network with 625 nodes (see Fig. 1b). In our case, we consider that all edges are bidirectional. Specifically, we assumed that there are two weights for each bidirectional edge, one incoming and one outgoing link, and that each neuron sends an independent synaptic potential. All excitatory outgoing links activate and all inhibitory outgoing links send negative potential, which is in accordance with Dale's law 82 . We also performed simulations establishing only unidirectional links to detect changes in the dynamics (see Supplementary Material). Next, each pair of hubs is connected with probability κ; such a deliberate setup allows hubs to communicate with each other. Two major types of hubs are defined: global hubs (with connections of ≥ k 100 nodes) and local hubs (only connected to one group of 25 nodes and the rest of the hubs ≤ ≤ k 20 4 4). A previous report has identified a lack of links between hubs in this network and its importance in modeling neural networks 83 . The rich-club organization is a main structural property of neural network systems where nodes with high connectivity degree tend to connect each other 7 . In order to detect the rich-club phenomenon in our configuration, we use the normalized rich-club coefficient Φ norm 84 , which is defined as the number of edges between pairs of hubs normalized by the number of edges in a null model network. The null model comprises the same degree distribution with random connectivity between pairs of nodes. For Φ > 1 norm , the configurations show legitimate rich-club organization. Figure 1c shows rich-club coefficients of the model for different values of κ and connectivity degree k.
Izhikevich neuron model. We use the integrate-and-fire neuron model presented by Izhikevich 85 . The principal advantage of the model over others is that it reproduces many dynamic features of real neurons with low computational cost 86 . The model is a two-dimensional system of ordinary differential equations defined by the following equations: where h is the size of the time step, v represents the membrane potential of the neuron, the notation =  v dv dt / , u represents a membrane recovery variable, which emulates the sodium -potassium pump. The parameters a and b represent the time scale and the sensitivity of the recovery variable u, respectively. The parameters c and d, represent the after-spike reset of v and u respectively. The value I represents a noisy thalamic input that is received for each of the neurons, and the value s represents the sum of the incoming potentials of its nearest neighbors when they fired. In our simulation, we consider a diverse repertoire of values for the parameters. For excitatory neurons: , the inhibitory neurons show a slow recovery with low-threshold spiking; however if = r 1 i , the neurons show fast spiking. In the hierarchical model, the hubs represent 4% of all nodes; we variate the proportion of excitatory-inhibitory hub neurons to observe its dynamics, so we set the 5 most connected neurons or global hubs as either inhibitory (case 1) or excitatory (case 2). Next, we variate the proportion of excitatory-inhibitory local hubs and define the value η as the probability of a local hub becoming inhibitory. For the rest of the nodes (96% of the total), we set 80% of them as excitatory and 20% as inhibitory.
In order to generate a time series of the evolution, we define the state of the system as, where N is the total number of units. We consider time evolutions of the system comprising 625 neurons. After a transient period (8000 time steps) the state of the system is monitored for 10 4 additional time steps. time resolution and synaptic weight. It is well known that the best individual accuracy of the Izhikevich model is reached for time steps, h, less than 0.5 ms [87][88][89] . However, the original Izhikevich model 85 uses the forward Euler integration with time steps of 0.5 and 1 ms, and other authors have used similar time steps for achieving collective and critical dynamics 77,80,90,91 (see Table 1). A very recent study by Pauli et al. 92 , focused specifically on the reproducibility of the Izhikevich model, demonstrates that the model is highly sensitive to the integration time step. Namely, they observed dynamic changes in network dynamics while they were decreasing the time resolution from 1.0 ms to 0.1 ms. The authors also offer a guide to reproduce global dynamics at time steps of 0.1 ms, which is based on the increment of the synaptic weight s in relation to the corresponding average firing rate. Here we adopt Pauli et al. 's approach by looking for an average excitatory firing rate of 5 spikes per second (spks/s), which is an intermediate firing rate close to the ones used in the original Izhikevich model (2-7 spks/s) 93 . In addition, we use the 2nd-order Runge-Kutta (RK-2) method to implement the numerical integration of the differential equations with a time resolution of 0.1 ms (see Results Section for additional discussion about other time resolutions). Using the RK-2 with a time resolution of 0.1 ms can be considered as a high-quality simulation 89 .
We also provide an online source of the model's implementation, where the time evolution, a simulation in mp4 format, and the detrended fluctuation analysis can be obtained from a standalone script 94 .
Detrended Fluctuation Analysis. Detrended fluctuation analysis (DFA) is a reliable method to detect long-range time correlations in nonstationary time series 95 . DFA considers the following steps. First, the signal (in this case the state of the system, S(t)) is integrated; the resulting series (y i ) is divided into windows of size n and for each window, a straight line is fitted to the points (y n ). Next, the root-mean-square fluctuation is computed of the detrended sequence within each window: is present, then the correlation exponent α characterizes the original signal. It is known that α = . 0 5 corresponds to white noise (uncorrelated signal), α = .
1 5 corresponds to Brownian noise, and α = 1 corresponds to a long-range correlated process or 1/f noise. In this work, we obtain F(n) and α for the state of the system S(t) under different configurations.
Global Lability of synchronization. The aim of this work is to measure the avalanche activity among sets in the system in terms of excitatory/inhibitory hub proportion and rich-club connectivity. Here, we use the global lability of synchronization method 53 , which is suitable for the detection of local and global events of synchrony. An advantage of this measure over the standard avalanche analysis is that the definition of the avalanche size is based on a quiescent initial condition in the network, so some kinds of sustained firing or supercritical behaviors are difficult to measure using the avalanche analysis. Moreover, the lability of global synchronization allows us to detect a wide range of stability regimes. This measure is based on the rate of change in the number of phase-locked signals pairs between successive time steps. We generate 125 signals from the average local state (of membrane potentials v i ) of clusters of 5 nodes. As in the original lability work, we assume that two signals are synchronized at certain time step if they satisfy two conditions: the absolute value of the phase difference, , must be smaller than π/4; an index synchronization parameter, γ t ( ) ). The phase difference between two signals s i (t) and s j (t) is given by 96 represent the Hilbert transform of s i (t) and s j (t), respectively. The Hilbert transform is defined as 96 : where P.V. indicates the Cauchy principal value. The phase synchronization index is used to meet the second condition 96,97 : where v is the size of the time window. We set ν = 50 time steps. For γ = 0 i j , , there is a large frequency difference during the time window, while for γ = 1 i j , , the phase difference is constant during the whole interval, indicating full synchronization. Next, we count the number of pairs of signals that are synchronized at time t: Finally, from the statistics represented by M(t), we calculate the square of the difference in the number of synchronized pairs between two successive time steps, representing the lability of global synchrony 53 : 2 A large value of  indicates a global synchronization/desynchronization process, whereas a small value indicates a local change in synchronization. Previous studies reported that systems operating at critical regimens, such as Ising and Kuramoto models, show negative power-law distributions of  values 53,98 , and similar distributions were also found for human brain functional signals 53 .

Results
Integration time step, firing rate and synaptic weight. To determine how robust our results are, first we address the question about the effects of changes of integration time steps on firing rate and on synaptic weight. Figure 2 shows the average firing rate of the network (Fig. 1b) vs. the synaptic weight for case 1 (κ = . 0 5 and η = 0), for different time resolutions. The value of the parameter η was selected to have a balanced inhibition/ excitation hub activity, where all global hubs are inhibitory and all local hubs are excitatory. We observe that as the size of the time step increases, the model needs less synaptic weight in order to display the same firing rate. For instance, for time steps of 0.1 ms, a firing rate of 5 spks/s is reached approximately at 40 mV of synaptic weight. It is important to mention that the injected current is applied for a single time step. For this reason, the different curves are scaled versions of each other, with a scaling factor given by the ratio of the time resolutions. Thus, for h = 0.1 ms, the synaptic weight needs to be large to recover the dynamics of the original model.
Furthermore, other studies have found that synaptic strength and connectivity (degree) play a key role in driving the system toward criticality. As the synaptic strength increases, the system shows more synchronization 90 , and more connectivity is associated with more activity in the system 90 . Lastly, in order to achieve critical dynamics, synaptic weight and connectivity have to be large enough, otherwise, the system shows subcritical behavior 90 . In Table 1 we present the main characteristics of representative spiking neuron models, which are able to reproduce either global or critical dynamics.
For the first two models presented by Izhikevich, we notice that the decrease of incoming connectivity in the 2006 model was replaced by an increase in synaptic strength. A similar strategy is used in Massobrio's study. In general, it is observed that for smaller time steps, the synaptic weight is bigger. A high time resolution reduces global activity, however, more incoming connectivity and more synaptic strength have been used to increase global activity. time evolution and temporal autocorrelation. In order to characterize the signature of dynamical evolutions of specific configurations cases, where the hub activity plays an important role, first we investigate the time evolution of the hierarchical network model at step 1 governed by the Izhikevich model.
In Fig. 3 we show representative time evolutions of the simplest configuration in our model. The network comprises 1 cluster with 25 units (one cluster in Fig. 1a). Each network posses 5 randomly selected inhibitory neurons, 19 excitatory neurons and 1 hub. For the case when the hub is inhibitory (Fig. 3a,b), we observe that the activation of neurons is transmitted to the hub, and the hub sends a negative feedback and rapidly stops the propagation of the firings. In contrast, for the case when the hub is excitatory (Fig. 3c,d), the activation of the excitatory hub promotes the activation of all units, giving rise to fully synchronized states. Next, we investigate the time evolution of the hierarchical network model at step 2 governed by the Izhikevich model. Figure 4a depicts the time evolution of the neural model where the most-connected node (node 125) and two local hubs (node 50 and 100) are inhibitory. In the figure, dash lines surround the temporal evolution of hubs, and also divides the 5 clusters of 25 neurons. For example, neuron 25 is a local excitatory hub which is connected with 20 neurons with the neuron indexes in the range 1-24.
We observe that the avalanche activity is present in terms of temporal activity of clusters, which are not stable and this fact contributes to the emergence of 1/f fluctuations in the S(t) dynamics (Fig. 4b,c). Previous findings indicate that there are several real and simulated systems operating at critical regimes which exhibit 1/f dynamics. The majority of these systems do not posses a long-tail connectivity distribution to show 1/f dynamics. However, the evidence of functional and structural connectivity suggest that the real topology of brain networks is characterized by the presence of hub regions, thus it is important to study the emergence of 1/f dynamics in the context of hub interconnectivity and functionality.
The activation of clusters (comprising 25 nodes) is mainly initiated by small clusters (comprising 5 nodes) with no hub nodes, then the activation is transmitted and amplified to other small clusters by excitatory hubs. This qualitative behavior is quite similar to the findings reported by Luccioli et al. 68 , where the excitatory hubs are   8 . It is noteworthy that during the time evolution, the firing activity between the different clusters is correlated, suggesting that excitatory hubs play a central role in information transmission between clusters. Besides, we observe diverse activity patterns. In the literature, critical dynamics are related to high information transmission and the existence of many attractor states 8,55,99 . On the other hand, the global inhibitory hub controls the occurrence of large avalanches (Fig. 4a), suggesting the important role that inhibitory hubs (GABA neurons) play in giving rise to critical dynamics. To our knowledge, there is some evidence that integrate-and-fire neurons can display 1/f fluctuations [100][101][102] . Here Fig. 4 shows that a mixed population of inhibitory and excitatory hubs and an intermediate rich-club connectivity promote the deployment of 1/f fluctuations. In order to test the robustness of the dynamics with respect to the number of neurons in the network, we repeat the temporal analysis with a set of 5 clusters with 125 units each. Now 5 replicas of the hierarchical network at step 2 with 625 nodes are considered, where 5 of them are global hubs and 20 are local hubs (see Fig. 1b). Two different cases are considered: global hubs are inhibitory (case 1) and global hubs are excitatory (case 2). The time evolution of case 1, where global hubs are inhibitory (κ = . 0 75 and the fraction of local inhibitory hubs η = . 0 75) is shown in Fig. 5a. The time evolution shows that the activation of big clusters (comprising 125 nodes) is correlated with other big clusters within and outside the 5 replicas (see Supplementary Material for the behavior of the correlation exponents in terms of the number of clusters). Interestingly, the state of the system exhibits 1/f fluctuations (Fig. 5b,c). We observe that the dynamics corresponding to a mixed population of inhibitory and excitatory local hubs with inhibitory global hubs is reminiscent of the case 1 with 1 cluster.
The time evolution of case 2, where global hubs are excitatory (κ = . 0 15 and η = . 0 9) is shown in Fig. 6a. We observe a global correlated behavior, where the vast majority of neurons in the big clusters are synchronized. We also observe that 1/f fluctuations get destroyed and the state of the system S(t) displays Brownian-like noise with α ≈ .
1 38 (Fig. 6b,c). To strengthen our results, we systematically calculate the scaling temporal exponent α in terms of the rich-club connectivity κ and the excitatory/inhibitory hub proportion η. The phase space of the scaling temporal exponent α is shown in Fig. 7. For the hierarchical network in case 1 (Fig. 7a), the phase space mainly exhibits 1/f noise (α ≈ . 1 0). However, for η < . 0 3 and κ > . 0 4, we observe a region where the scaling exponent is α > . 1 3, indicating Brownian-like fluctuations. When the number of excitatory hubs is dominant (low value of η) and the rich-club connectivity increases, the system becomes more activated, which is reflected in the increment of the α exponent. For intermediate values of both κ and η, global inhibitory hubs control the production of global avalanches and this fact is reflected in the deployment of 1/f noise, also in agreement with the results described in Fig. 4. In contrast, for case 2 (Fig. 7b), the phase space mostly exhibits exponent values close to Brownian noise (α ≈ . 1 5). Notably, when the hierarchical network structure is destroyed but the same incoming degree distribution is preserved, the model still displays a variety of dynamical behaviors, but with a more widespread presence of Brownian-like noise, indicating a more activated behavior with short-range correlations (see Supplementary Material). In addition, considering only unidirectional links promotes the decrement of α exponents values and more presence of 1/f noise (see Supplementary Material). The fact that the phase space gets slightly altered by Global lability of synchrony. Figure 8 shows the probability density of lability values for case 1 and different values of η and κ. We observe that the probabilities follow a power-law , with exponent values within the range δ . ≤ ≤ . 0 5 0 9. For the estimation of the power-law exponent, we used the method proposed by Hanel et al. 103 , which is suitable for the estimation of power-law exponents that are located within the range δ < ≤ . 1 0 and η = . 0 25, the corresponding exponent is δ = . ± . 0 83 0 08. Interestingly, for these intermediate values of κ and η we also observe long-range correlated temporal dynamics in Fig. 7a. Moreover, as the fraction of inhibitory neurons increases (a larger η), the frequency of lability events decreases, indicating a reduction in the occurrence of local and global synchronization changes. For η = 0 (mostly excitatory hubs) and high values of κ, the δ exponents lie within the range from 0.66 to 0.58, indicating a legitimate power-law behavior. These values of κ and η correspond to Brownian fluctua-   Fig. 9, we show the probability density of lability values for case 2. We observe that for low hub interconnectivity, κ = .
0 25, the majority of distributions follow a power-law function with exponent δ ≈ . 0 5. These distributions are similar to the ones found in case 1 (η = .
0 0 and κ ≥ . 0 75), suggesting that for these configurations the synchronization exponent δ ≈ . 0 5 is associated with Brownian fluctuations in the state of the system (see Fig. 7a). Moreover, when the majority of hubs are excitatory and the rich-club connectivity is high (κ = . 1 0 and η = . 0 0), the lability distribution displays an almost flat behavior with exponent δ ≈ . ± . 0 24 0 06, where local and global phase changes take place with similar probabilities.

Discussion
Our current understanding of brain function is directly linked to network topology and dynamics occurring within the brain, which have been explored in detail by network neuroscience. However, limitations imposed by recording tools -which provide low quality readings of the brain -continue to generate gaps of our understanding of how brain networks operate. Our simple neural model with scale-free topology is able to display a variety of dynamical patterns, levels of synchronization and temporal autocorrelation in terms of a hub's inter-connectivity and excitatory/inhibitory proportions. Our results show that 1/f fluctuations are closely related to a power-law distribution of global phase changes, which emerge when there is a mixed population of inhibitory and excitatory hubs for broad interval values of rich-club connectivity. However, when the majority of hubs are excitatory and the rich-club connectivity is high enough, the system displays Brownian fluctuations, indicating a high frequency of global phase changes. These results are consistent with previous findings about the fact that identified hub neurons in Human and C. elegans neural networks are principally inhibitory (GABA interneurons) 24,104 . It is worth mentioning that similar critical behavior displaying 1/f noise, critical avalanche distribution and critical synchronization have been also detected in regular and small-world configurations 78,105,106 , and that in many cases the existence of 1/f noise depends on the balance between excitatory and inhibitory interaction 77,80 . However, much evidence suggests that real neural networks posses units that exert more influence than others, and that this influence is reflected in the connectivity of the units and synaptic strength. Our study based on the high resolution Izhikevich model comprises large synaptic weights that are higher than those reported in real electrophysiological recordings. Further works should incorporate the problem of the duration of synaptic release in order to reduce the weights and generate a model with more realistic parameter values. Furthermore, our results are concordant with previous findings about the fact that a scale-free distribution of degree connectivity promotes critical distribution of avalanches; configurations with only excitatory hubs are unable to show criticality 90 and 1/f noise appears when a mixed population of inhibitory and excitatory hubs are present 77 . Moreover, the disruption of rich-club connectivity promotes the desynchronization of the system and the reduction of the temporal autocorrelation exponent, while the total shift from inhibitory to excitatory behavior of hubs destroys 1/f fluctuations and power-law synchronization. Our results also show that the disruption of the hierarchical connectivity promotes the increase of global activity and the presence of Brownian fluctuations, on the other hand, the unidirectional link configuration facilitates the decrement of α exponents and the existence of 1/f noise.
Although different mechanisms and models have been proposed to generate 1/f noise [72][73][74][75] , it is important to remark that few models are able to generate long-range correlations with critical power-law synchronization based on large number of communicating units or neurons 77,78 . In these models, the topology and a balanced excitation/inhibition population of units seem to play a important role for an efficient functional communication between different parts of the network 80,90 . In the context of neurophysiological systems, mental illnesses like epilepsy, schizophrenia and Alzheimer have been associated to alterations in the rich-club connectivity or the failure of hub responses. In addition, we demonstrate that, in the context of our model, ingredients like rich-club connectivity, inhibitory global hubs and excitatory local hubs are relevant to generate a vast repertoire of spatio-temporal patterns, including 1/f fluctuations and critical synchronization. In the context of brain studies, which focus on network topology and dynamics, our results indicate that the rich-club connectivity and the hub inhibitory/excitatory population are important properties, which may help to understand the variety of temporal correlations and synchronization levels reported in brain dynamics.