1 Introduction
Self-organized critical states are found for many complex systems in nature,[1-2] from earthquakes to neuronal avalanches.[3-10] Several lines of evidence point to the existence of such critical states in brain activity.[11-22] The critical avalanches were observed by in vitro and in vivo experiments.[23-30] It was shown that neuronal avalanches may play an important role in cortical information processing and storage.[28,31-32] It was suggested that cortical networks that generate neuronal avalanches benefit from a maximized dynamic range.[33-37]
The excitable networks can be used as a simplified model to study the criticality in neural networks.[38-39] A model of an excitable network on Erdös-Rényi (ER) random graphs was proposed by Kinouchi and Copelli (KC model).[40] It is based on the branching process, which was used to understand the mean-field behavior of self-organized criticality.[41] It was shown that the criticality occurs at branching ratio is one and the dynamic range is maximized at critical point.
In realistic neural systems, the excitatory and inhibitory neurons coexist in networks.[42-44] The behavior of excitatory-inhibitory (EI) network is critical for understanding how neural circuits produce cognitive function,[45] and inhibition is always included in neural models.[19,39,46] The balance of excitation and inhibition has been shown that it plays an important role in the neuron discharge rate,[47] dynamic range of response,[48] and temporal structure of the spontaneous avalanche activity.[19,49] Prominent examples of balance include the intracellular recordings performed by Ferster[50] in simple and complex cells of the visual cortex in the cat. The inhibitory elements can be studied in the framework of KC model. In Ref. [51] network consists of the excitatory layer and the inhibitory layer which have the same coupling strength. The critical point is determined by the fraction of the excitatory elements.
The strength of inhibitory couplings is an attractive feature of neural networks. The inhibitory strength is several times larger than the excitatory strength. This is physiologically realistic, since inhibitory synapses are often closer to the cell body of postsynaptic neurons.[46] In Ref. [35] the criticality is changed dramatically by using antagonists to change the efficacy of synaptic couplings in experiment. It was shown that balanced excitation and inhibition establishes criticality, which maximizes the dynamic range. The change of the couplings was modeled by adjusting the branching ratio in an excitatory excitable network. However, little is known about the impact of changing the excitation and inhibition in EI network.
In this paper, we consider a model of EI network on the ER graphs using the KC model. We investigate the effect of inhibitory signal's strength on criticality and dynamic range. We show that the critical point is not affected by inhibitory signal's strength, and stronger inhibitory signal decreases the dynamic range. Moreover, to simulate the efficacy of antagonists, we change the network in three ways: removing excitatory/inhibitory nodes, deleting excitatory/inhibitory links, and weakening excitatory/inhibitory coupling strength. We show that the inhibition does not affect network criticality and sensitivity. We also conduct an analytic study of EI network to explain the effect of inhibition.
2 The Excitatory-Inhibitory Network Model
The EI excitable network model includes both inhibitory elements and excitatory elements.[48] The signal transferred from an excitatory node increases the probability that the neighbors of this node are excited, while the signal from an inhibitory node decreases this probability.
Each excitable element $i=1,\ldots, N$ has $n$ states: ${s}_i=0$ is the resting state, ${s}_i=1$ corresponds to excitation and the remaining ${ s}_i=2,\ldots,n-1$ are refractory states. To be precise, at discrete times $t=0,1,\ldots$ ($\Delta t=1$ ms) the states of the nodes ${s}_i$ are updated as follows:
(i) If node $i$ is in the resting state, ${ s}_i=0$, it can be respectively inhibited by inhibitory neighbor $j$ with probability $p^i_{ij}$, when node $j$ is in excited state $s_{j}=1$. In this case, the state of node $i$ will remain 0 in the next time step. Otherwise, it can be excited with probability $p^e_{il}$ by excitatory neighbor $l$ which is in excited states, or independently by an external stimulus with probability $\lambda$ ($\lambda=1-\exp(-\eta)$). The maximum transfer probabilities of excitatory and inhibitory nodes are denoted by $p^e_{\max}$ and $p^i_{\max}$. The probability $p^i_{ij}$ or $p^e_{ij}$ is a random variable with uniform distribution in the interval [0, $p^{i}_{\max}$] or [0, $p^{e}_{\max}$]. Each element receives external signals independently, that is, we have a Poisson process for each element.
(ii) After excitation the dynamics of nodes are deterministic: if ${s}_i=1$, then in the next time step its state changes to ${s}_i=2$, and so on until the state ${s}_i=n-1$ leads to the resting state ${s}_i=0$. So the node is a cyclic cellular automaton.
The system consists of $N$ elements, which are coupled by an ER random network. We randomly choose $N_e=f_eN$ nodes as excitatory elements, $0<f_e<1.0$. The rest $N_i=f_iN$ nodes as inhibitory ones, $f_i=1.0-f_e$. To build the random network with average degree $K$, $NK /2$ links are assigned to randomly chosen pairs of nodes.
The local branching ratio of excitatory nodes
$\sigma_j=\sum\limits^{K_j}_{i}p^e_{ij},$
corresponds to the average number of excitations created in the next time step by the excitatory node $j$, where each excitatory element $j$ is randomly connected to $K_j$ neighbors. The average branching ratio $\sigma=\langle\sigma_j\rangle$ is the relevant control parameter.
The network instantaneous activity is characterized by the density $\rho_t$ of active nodes $(s = 1)$ at a given time $t$. The average activity is defined as
$F=T^{-1}\sum\limits^{T}_{t=1}\rho_t, $
where $T$ is a large time window (of the order of $10^{3}$ step). As a function of the stimulus intensity $\eta$, networks have a minimum response $F_0$ and a maximum response $F_{\max}$. Variations in $\eta$ can be robustly coded by variations in $F$, discarding stimuli that are too weak to be distinguished from $F_0$ or too close to saturation. The range [$\eta_{0.1},\eta_{0.9}$] is found from its corresponding response interval [$F_{0.1},F_{0.9}$], where
$F_x=F_0+x(F_{\max}-F_0). $
The dynamic range $\Delta$ is defined as the range of stimuli (measured in dB)[40]
$\Delta=10\log(\eta_{0.9}/\eta_{0.1}).$
This choice of a $10%-90%$ interval is arbitrary, but is standard in the literature and does not affect our results.[40]
3 Simulations and Results
We first study how the strength of inhibitory couplings changes the transition of the network's average activity. The coupling strengths are determined by the values of $p^e_{\max}$ and $p^i_{\max}$. For excitatory couplings, we set $p^e_{\max}=2\sigma/K$. Then the mean value of $p^e$ is $\sigma/K$. To study the effects of the strength of inhibitory couplings in the model, we set the maximal inhibitory probability as $p^i_{\max}=2m\sigma/K$. To simulate that inhibitory coupling strength is larger than the excitatory coupling strength, we set $m>1$ in simulations. We input external stimulus for 1000 steps, then remove the stimulus and calculate the average activity $F_{\eta\rightarrow0}$.
Subcritical states have a negligible average activity, while supercritical states have a large average activity. As can be seen in the inset of Fig. 1, in subcritical state, the network becomes silent quickly. In supercritical state, network presents self-sustained activity. In the critical state, the duration of activity has a large variance. At this state, the average activity becomes into non-zero. We plot the response curves $F$ versus $\sigma$ for various inhibitory signal strength $m$. In Fig. 1, one can see that all the curves of $F$ show the transition at the same position. The critical value of the $\sigma$ is 1.25 which is agree with the value obtained in networks with $p^e_{\max} = p^i_{\max}$.[51] We obtained that the critical points are not changed as the strength of inhibitory couplings increases.
Fig. 1 (Color online) Response $F_{\eta\rightarrow0}$ versus $\sigma$ on EI networks for different $m$ with $N=10^5$, $K=10$, $n=5$, $f_e={4}/{5}$. The inset shows the instantaneous density of active sites for subcritical (square), critical (circle), and supercritical (triangle) states as function of time. Three different runs for each case are presented. |
We also investigate the effects of inhibitory strength on the response of EI network to stimulus. The average activity $F$ is used as the response of networks to stimulus. Figure 2 shows the relation between the response of the network and the stimulus strength for $m=5$. The results show that the EI networks can exhibit the same features of response curve as in excitatory networks. When $\sigma<\sigma_c$, the system is relatively insensitive. When $\sigma=0.85<\sigma_c$, the slope ($k$) of response curve is near to 1, when $\sigma=\sigma_c$, the slope comes to 0.5, sensitivity is enlarged because weak stimuli are amplified. As a result, the dynamic range increases monotonically with $\sigma$. When $\sigma>\sigma_c$, the spontaneous activity $F_0$ masks the presence of weak stimuli, the network maintains a large response in week stimuli.[40] Therefore dynamic range decreases. The response curves for different $m$ in critical states are shown in the inset of Fig. 2. The response curve is changed slightly. The effect of inhibitory strength $m$ on average activity $F$ may lead to slight change in the dynamic range.
Fig. 2 (Color online) Response curves (mean firing rate versus stimulus rate) from $\sigma=0.85$ to 1.65 (in intervals of 0.1) with $m=5$. Inset: response of the critical network with different $m$. $f_e=4/5$, $N=10^5$, $K=10$, $n=5$. |
We obtain the dynamic range versus the branching ratio for different inhibitory strength. In the simulation, we approximate $\eta_{0.9}$ to 1. The curves are shown in Fig. 3. As the same as excitatory networks, the dynamic range peaks at the critical point. One can see that, the dynamic range is decreased with the inhibitory strength $m$. As shown in inset, the error of dynamic does not overlap for different inhibitory coupling strength. We can conclude that increasing the inhibitory coupling strength do not affect the critical point, and causing slightly decrease of dynamic range. As the error of dynamic range is weak, it is not shown in later section.
Fig. 3 (Color online) Dynamic range versus branching ratio for different values of inhibitory coupling strength $m$. The points represent simulation results with $m=1$ (square), $m=5$ (circle), and $m=10$ (triangle) on EI networks. The lines correspond to theoretical results from Eqs. (20) and (21) with $m=1$ (black), $m=5$ (red), $m=10$ (blue). The inset presents the error of dynamic range for different inhibitory coupling strength. |
In order to modeling the experiments in Ref. [33] that using antagonists to change the excitation and inhibition, we change the network in three ways: removing a fraction ($r$) of excitatory or inhibitory nodes ($0<r<1$), deleting a fraction ($d$) of excitatory or inhibitory links ($0<d<1$), and weakening excitatory or inhibitory coupling strength by a fraction $w$ ($0<w<1$) on critical network.
Figure 4 shows that when reducing the excitation by each of the three methods, the dynamic range $\Delta$ decreases. Each method of reducing the excitation can recover the phenomenon in experiments.[33] As excitatory elements are removed, the change of the dynamic range is most dramatic. The second and third method of reducing the excitation has the same result.
Fig. 4 (Color online) The effect of decreasing the excitation (square) and inhibition (circle) on dynamic range. The triangle points represent the network's dynamic range in critical state. (a) Removing a fraction ($r$) of excitatory or inhibitory nodes. (b) Deleting a fraction ($d$) of excitatory or inhibitory links. (c) Weakening a fraction ($w$) of excitatory or inhibitory coupling strength in critical network. |
Most interestingly, none of the three methods of reducing inhibition can change the dynamic range. Therefore the model of EI excitable network cannot recover the experimental phenomenon of changing the inhibition. In the following we present an analytic treatment of the effect of inhibition on the criticality.
4 Analytical Results
4.1 Increasing the Inhibitory Coupling Strength
In the mean-filed approximation, the average degree $K$ is used to replace the degree $K_i$.[40,51-52] The excitatory probability $p^e_{ij}$ as the average value $\sigma/K$, and the inhibitory probability $p^i_{ij}$ as $m\sigma/K$. Therefore, the evolution of the response $F$ is described by the following mean-field map,
$F_{t+1}= P_{t}(0)\Big(1-\frac{m\sigma F_{t}}{K}\Big)^{f_{i}K} \\ \times\Big\{\eta +(1-\eta)\Big[1-\Big(1-\frac{\sigma F_{t}}{K}\Big)^{f_{e}K}\Big]\Big\}, $
where
$P_t(0)=1-(n-1)F_t, $
is the approximate probability that a randomly selected site is in the resting state ($s=0$). $(1-{m\sigma F_t}/{K})^{f_iK}$ is the probability that a resting node at time $t$ will not be inhibited at the next time step by any of its inhibitory neighbors, and $[1-(1-{\sigma F_t}/{K})^{f_eK}]$ is the probability of being excited by at least one of its excitable neighbors. In the stationary state, $F_{t+1}=F_{t}=F$. When external stimuli $\eta=0$, the response function $F$ is given by the solution of
$F=(1-(n-1)F)\Big(1-\frac{m\sigma F}{K}\Big)^{f_iK} \Big[1-\Big(1-\frac{\sigma F}{K}\Big)^{f_eK}\Big].$
Then we linearize the terms $(1-{m\sigma F}/{K})^{f_iK}$ and $(1-{\sigma F}/{K})^{f_eK}$ to the first order using binomial expansion,
$(1-\frac{m\sigma}{K}F\Big)^{f_iK}\simeq 1-f_i m \sigma F,$
$(1-\frac{\sigma}{K}F\Big)^{f_eK}\simeq 1-f_e \sigma F.$
The order parameter behavior is
$F=(1-(n-1)F)(1-f_im\sigma F)f_e\sigma F.$
In the right terms, the parameter associated with the inhibition ($f_i, m$) are multiplied with the second or third order of response $F$, while the parameter associated with the excitation ($fe$) can be multiplied with the first, second and third order of $F$. One solution of Eq. (10) is $F=0$, when $F>0$, Up to the first order, more terms that associated with the excitation are retained:
$f_e\sigma -f_if_em\sigma^2F-(n-1)f_e\sigma F =1.$
Then the response function $F$ is given by the solution of
$F=\frac{f_e\sigma-1}{f_if_em\sigma^2+(n-1)f_e\sigma}.$
When the control parameter $\sigma \rightarrow 1/f_e$, the activity $F \rightarrow 0$. Therefore, the critical point is
$\sigma_c = \frac{1}{f_e}.$
and $F=0$ when $\sigma<1/f_e$. The result does not depend on the average number of neighbours $K$. The parameter used in Fig. 1 is $f_e=4/5$, so the critical point is $\sigma_c=1.25$ for different values of $m$.
The analytic results show that the critical point is independent of the strength of the inhibitory signal. The theoretical result agrees with the simulations.
Next we analyze the effect of the inhibitory signal strength on the dynamic range. In the limit $F\rightarrow0$, the response function in Eq. (5) can be approximated by
$F=[1-(n-1)F]\exp{(-f_im\sigma F)}\{\eta+(1-\eta)[1-\exp{(-f_e\sigma F)}]\}. $
Given response $F=F_{0.1}$ we can obtain the stimulus $\eta_{0.1}$
$\eta_{0.1}=1-\exp{(f_e\sigma F_{0.1})}+\frac{F_{0.1}}{1-(n-1)F_{0.1}}\exp{[(mf_{i}+f_{e})\sigma F_{0.1}]}. $
For $\eta_{0.9}$, Eq. (14) is invalid because it applies only to $F\rightarrow0$. In order to simplify the analysis, we set $\eta_{0.9}$ to 1. It has been employed in Refs. [40, 51]. Using Eq. (4), we have
$\ \Delta=-10\log\Big[1-\exp{(f_e\sigma F_{0.1})}+\frac{F_{0.1}}{1-(n-1)F_{0.1}}\exp{[(mf_i+f_e)\sigma F_{0.1}\Big]}.$
We expand the terms $\exp{(f_e\sigma F_{0.1})}$ and $[(m f_i+f_e)\sigma F_{0.1}]$ in Eq. (16) to the second order
$exp(f_e \sigma F_{0.1}) \simeq1+f_e\sigma F_{0.1}+\frac{f_e^2\sigma^2F_{0.1}^2}{2},$
$exp{[(mf_i+f_e)\sigma F_{0.1}]}\simeq1+(mf_i+f_e)\sigma F_{0.1}+\frac{(mf_i+f_e)^2\sigma^2F_{0.1}^2}{2}.$
In the limit of $F\rightarrow0$, the dynamic range is
$\Delta=-10\log\Big[-f_e\sigma F_{0.1}+\frac{F_{0.1}}{1-(n-1)F_{0.1}}-\frac{f_e^2\sigma^2F_{0.1}^2}{2}+\frac{(mf_i+f_e)\sigma F_{0.1}^2}{1-(n-1)F_{0.1}}\Big].$
Up to the first order, $\Delta$ follows the equation,
$\Delta=-10\log\Big[-f_e\sigma F_{0.1}+\frac{F_{0.1}}{1-(n-1)F_{0.1}}+\frac{(mf_{i}+f_e)\sigma F_{0.1}^{2}}{1-(n-1)F_{0.1}}\Big].$
For a system with refractory time $n$, the maximal response $F_{\max}=1/n$. According to Eqs. (3) and (12), we have
$F_{0.1}=\frac{1}{10n}+\frac{9(f_e\sigma-1)}{10f_{i}f_{e}m\sigma^2+10(n-1)f_{e}\sigma}.$
When $\sigma\le1/f_e$,
$F_{0.1}=\frac{1}{10n}.$
We estimate the theoretical results of dynamic range using Eqs. (20) and (21). The inhibitory signal strength is in the term of $F_{0.1}^2$ which is relatively small in $\Delta$ as shown in Eq. (20). The theoretical results are shown as the curve in Fig. 3. The curves are in agreement with the simulations, that the dynamic range is taken to the maximum value in the critical state, and the increase of inhibitory strength can decrease the dynamic range.
Fig. 5 (Color online) The critical point and dynamic range of regular random network. (a) Response $F_{\eta\rightarrow0}$ versus $\sigma$ for different $m$. (b) Dynamic range versus branching ratio for different $m$. The point represents simulation results with $m=1$ (square), $m=5$ (circle), and $m=10$ (triangle). The lines correspond to the theoretical results with $m=1$ (black), $m=5$ (red), $m=10$ (blue). The $N=10000$, $n=2$, $f_e=4/5$, $K=30$. |
4.2 Reducing Excitation and Inhibition
For analysis of simulate the experiment in Ref. [33], we first remove a fraction ($r$) of excitatory nodes. The evolution of the average response $F$ is described by the following mean-field map,
$F_{t+1}=P_{t}(0)\Big(1-\frac{\sigma F_{t}}{K}\Big)^{f_{i}K} \\ \times \Big\{\eta +(1-\eta)\Big[1-\Big(1-\frac{\sigma F_{t}}{K}\Big)^{(1-r)f_{e}K}\Big]\Big\}.$
Removing the excitatory nodes decrease the probability of being excited by at least one of its remaining excitable neighbors. When external stimuli $\eta=0$, the response function $F$ is given by the following solution in the stationary state ($F_{t+1}=F_{t}=F$),
$F= (1-(n-1)F)\Big(1-\frac{\sigma F}{K}\Big)^{f_iK} \\ \times\Big[1-\Big(1-\frac{\sigma F}{K}\Big)^{(1-r)f_eK}\Big].$
Following the method used in Subsec. 4.1, the response function $F$ is given by the solution of
$F=\frac{(1-r)f_e\sigma-1}{(1-r)f_if_e\sigma^2+(n-1)(1-r)f_e\sigma}.$
When the control parameter $\sigma \rightarrow 1/(1-r)f_e$, the activity $F \rightarrow 0$. Therefore, the critical point is
$$ \sigma_c = {1}/{(1-r)f_e}. $$
$\sigma_c<1/f_e$, so network becomes into subcritical state.
Removing a fraction ($r$) of inhibitory nodes, the evolution of the average response $F$ is described as,
$F_{t+1}=P_{t}(0)\Big(1-\frac{\sigma F_{t}}{K}\Big)^{(1-r)f_{i}K} \\ \times\Big\{\eta +(1-\eta)\Big[1-\Big(1-\frac{\sigma F_{t}}{K}\Big)^{f_{e}K}\Big]\Big\}.$
A slight change in the probability of inhibition, changing the critical point to $\sigma_c = 1/f_e$, which equal to original network's branching ratio. The theoretical prediction fits the change of critical point. The response curves are presented in Fig. 6(a). As can be seen, removing the excitatory nodes increases the critical point. Therefore, the critical state becomes a subcritical state. However, removing the inhibitory nodes does not change the critical branching ratio.
Fig. 6 (Color online) Response curves (mean firing rate versus stimulus rate) for different fraction of decreasing the excitation (square) and inhibition (circle). The triangle points represent the network's mean firing rates in critical state. (a) In the case of removing a fraction ($r$) of excitatory nodes or inhibitory nodes, the $r=0.1,0.2$ are plotted separately. (b) In the case of deleting a fraction ($d$) of excitatory or inhibitory links, the $d=0.1,0.2$ are plotted separately. (c) In the case of weakening a fraction ($w$) of excitatory or inhibitory coupling strength, the $w=0.1,0.2$ are plotted separately. |
Deleting a fraction $d$ of excitatory/inhibitory links, leads to the following mean-filed map respectively:
$F_{t+1}=P_{t}(0)\Big(1-\frac{\sigma F_{t}}{K}\Big)^{f_{i}K}\Big\{\eta +(1-\eta)\Big[1-\Big(d+(1-d)\Big(1-\frac{\sigma F_{t}}{K}\Big)\Big)^{f_{e}K}\Big]\Big\},$
$F_{t+1}=P_{t}(0)\Big[d+(1-d)\Big(1-\frac{\sigma F_{t}}{K}\Big)\Big]^{f_{i}K}\Big\{\eta+(1-\eta)\Big[1-\Big(1-\frac{\sigma F_{t}}{K}\Big)^{f_{e}K}\Big]\Big\}.$
Their critical point is similar to removing nodes. When deleting the excitatory links, the critical point is $\sigma_c = {1}/{(1-r)f_e}$, when deleting the inhibitory links, the critical point is $\sigma_c = 1/f_e$. Figure 6(b) shows the network response curves as a function of branching ratio. The change of critical points is in agreement with predictions from our theoretical treatment.
As the excitatory/inhibitory coupling strength are weakened by a fraction $w$, the excitatory/inhibitory coupling strength, their average response $F$ is respectively given by the solution of:
$F_{t+1}=P_{t}(0)\Big(1-\frac{\sigma F_{t}}{K}\Big)^{f_{i}K}\Big\{\eta +(1-\eta)\Big[1-\Big(1-\frac{(1-w)\sigma F_{t}}{K}\Big)^{f_{e}K}\Big]\Big\},$
$F_{t+1}=P_{t}(0)\Big(1-\frac{(1-w)\sigma F_{t}}{K}\Big)^{f_{i}K}\Big\{\eta +(1-\eta) \Big[1-\Big(1-\frac{\sigma F_{t}}{K}\Big)^{f_{e}K}\Big]\Big\}.$
The corresponding critical points for them are $\sigma_c = 1/(1-r)f_e$, and $\sigma_c = 1/f_e$. This is consistent with the simulation results shown in Fig. 6(c).
Moreover, the dynamic range of EI network with reduced inhibition is given by the following solution
$\Delta=-10\log\Big[-f_e\sigma F_{0.1}+\frac{F_{0.1}}{1-(n-1)F_{0.1}}+\frac{((1-\alpha)f_{i}+f_e)\sigma F_{0.1}^{2}}{1-(n-1)F_{0.1}}\Big], $
where $\alpha = r, d, w$ for three methods of reducing inhibition respectively. The fraction of decreasing the inhibition is in the term of $F_{0.1}^2$, which is relatively small in $\Delta$.
Through the analytic treatment of three methods of reducing the efficacy of excitation and inhibition, we show that reducing the excitation of the network increases the critical value of $\sigma$, and changes the network in the critical state into the subcritical state. So the dynamic range is reduced. However, decreasing inhibition has no influence on critical point on EI network.
5 Discussions
In summary, we studied the role of inhibition in the criticality of EI network of excitable elements. We show that the critical point is not affected by inhibitory signal strength, and the increase of inhibitory signal strength will decrease the dynamic range slightly. Moreover, we use three different methods of reducing excitation and inhibition in order to simulating the experiment in Ref. [33]. Using a simulation and analysis, we show that reducing the excitation can make the critical state of the network into a subcritical state, thus dynamic range is decreased as the same as experiments. However reducing the inhibition cannot change network's criticality and dynamic range.
The results of the modeling study suggest that the experimental results in Ref. [33] cannot be understood by the change of the inhibition in excitable networks. We use the Poisson process to replace the dynamics of synapses, therefore time effect of spiking is ignored. Perhaps a better compromise between larger dynamic range and biological realism would be dynamics of synapses. Moreover, the process of stimulating nodes by neighbors causes the effect of inhibition is week. More sophisticated modeling studies are needed to enhance the effect of inhibition. The present model has the virtue of enabling analytical results that provide a benchmark for the performance of networks with other topologies. The scale-free small-world networks reflect important functional information about brain states.[54] The existence of functional hubs was observed by experiments, and perturbation of a single hub influenced the entire network dynamics.[55] In the future, verifying hypothesis that the inhibitory hubs may have an impact on criticality in a scale-free network would be a meaningful work.