Critical Neural Networks Minimize Metabolic Cost

Abstract

Brain dynamics show a rich spatiotemporal behavior whose stability is neither ordered nor chaotic, indicating that neural networks operate at intermediate stability regimes including critical dynamics represented by a negative power-law distribution of avalanche sizes with exponent $\alpha =-1.5$. However, it is unknown which stability regimen allows global and local information transmission with reduced metabolic costs, which are measured in terms of synaptic potentials and action potentials. In this work, using a hierarchical neuron model with rich-club organization, we measure the average number of action potentials required to activate n different neurons (avalanche size). Besides, we develop a mathematical formula to represent the metabolic synaptic potential cost. We develop simulations variating the synaptic amplitude, synaptic time course (ms), and hub excitatory/inhibitory ratio. We compare different dynamic regimes in terms of avalanche sizes vs. metabolic cost. We also implement the dynamic model in a Drosophila and Erdos–Renyi networks to computer dynamics and metabolic costs. The results show that the synaptic amplitude and time course play a key role in information propagation. They can drive the system from subcritical to supercritical regimes. The later result promotes the coexistence of critical regimes with a wide range of excitation/inhibition hub ratios. Moreover, subcritical or silent regimes minimize metabolic cost for local avalanche sizes, whereas critical and intermediate stability regimes show the best compromise between information propagation and reduced metabolic consumption, also minimizing metabolic cost for a wide range of avalanche sizes.

1. Introduction

At present, there is not a consensus opinion about a general neural model that can explain the emergence of the dynamic richness of neural systems [1] in which a complex interaction occurs among inhibitory and excitatory neurons. The self-organized critical [2] (SOC) system hypothesis states that cell systems like the heart and the brain operate at intermediate stability regimens. Critical regimens are conformed by a high presence of local activity propagation, moderate probability of medium propagation, and low presence of global propagation. These properties allow the propagation of information, maximizing the dynamic pattern repertoire [3,4]. Ordered or subcritical systems are conformed by only local activation propagation in which the information cannot be transmitted to the entire system. On the contrary, supercritical or chaotic regimens are very sensitive to activity propagation promoting high global activation frequencies. Moreover, SOC systems are strongly related to long-range time correlations ($1/f$ noise [5]) and fractal fluctuations [6], indicating dynamic processes with memory, which is found in real neural networks [7]. SOC systems are on edge between order and chaotic regimes; this borderline represents a stability phase transition, which in statistic physics is represented by a power-law distribution of magnitude events [8]. Healthy neural networks display a negative power-law distribution of propagation activity or avalanche sizes: $p\left(\mathrm{size}\right)\approx {\left(\mathrm{size}\right)}^{-1.5}$ [9,10,11], where the avalanche size is the number of different neurons that fired during the spontaneous activation and propagation of neuron activity. It is called avalanche activity because of the reminiscent behavior of an initial snow mechanical failure that leads to subsequent failures that collect more snow. The power-law distribution of avalanche sizes $p\left(\mathrm{size}\right)\approx {\left(\mathrm{size}\right)}^{\alpha }$, $\alpha =-1.5$, indicates that as the size/propagation is larger, there is less probability of occurrence. If $\alpha$ exponent is more negative, the decay probability is more rapid, indicating a very low probability of global size avalanches and subcritical dynamics. In contrast, a more positive exponent indicates a probability function with a slow decay as the size increases indicating a high presence of global activity or supercritical regimens. A disruption in excitatory/inhibitory chemical concentrations leads to the loss of the power-law distribution [9]. Besides, brain activity and SOC models display a strong relation between $1/f$ noise and intermediate synchronization [12,13,14,15]. Finally, critical systems are robust against the interaction with external dynamics and environments [14,16].

On the other hand, the biological evolution provided the brain with a maximal computational power and minimum wiring cost [17], and an economical brain design must involve a compromise between optimal performance and reduced structural and functional cost [17,18]. Concerning neural structure, the organization of neural networks possess clustering and few long-distance links, properties called small-world [19,20], hierarchical structure involving a dominant small population of neurons that link the rest of neurons at different modules and scales [21], these highly connected neurons are called hubs, and they are connected to each other, forming a rich-club [22]. Rich-club property facilitates the information transmission between brain regions reducing the wiring cost [23,24]. Rich-club also promotes a wide range of dynamic patterns [25]. Moreover, a scale-free distribution of connectivity has been found in hippocampus’ cornu Ammonis 3 region [26]. A scale-free function is defined by the function $P\left(k\right)\approx k^{-\gamma }$, where k is the number of links, $p\left(k\right)$ is the probability of having k links, and $\gamma$ is the scaling exponent. A scale-free function indicates that the set of k values mathematically does not possess a finite mean when $\gamma \le 2$, or finite variance when $\gamma \le 3$, and their respective mean and variance values are not representative of the set because of the great influx of high k values. A scale-free distribution is conformed by most nodes having few links, whereas a tiny population of nodes has many links [27]. In this work, we use the Ravasz-Barabási network that posses hierarchical organization and scale-free connectivity distribution [28], involving clusters with dominant hubs.

In terms of functional cost, neurons used action potentials to transmit information within the neuron and send neurotransmitters to produce synaptic potentials and communicate with other neurons [29]. Action potentials and synaptic potentials represent $80\%$ of the energetic consumption [30].

1.1. Problem Statement

However, neural activity’s metabolic consumption has not been analyzed for the different propagation dynamic regimes and network structures. This work aims to measure the metabolic cost of the avalanche sizes, comparing avalanche sizes vs. spikes and avalanche sizes vs. synaptic cost. Moreover, it is necessary to propose a mathematical description of the metabolic cost of neural activity. Finally, the majority of neurons, around 85%, are excitatory, and the rest 15 % are inhibitory [31]. This suggests that the majority of low connectivity nodes are excitatory. However, the most connected nodes that represent a low percentage of neurons could be either a majority of excitatory or inhibitory neurons. We try to elucidate this problem by looking for different inhibition/excitation hub ratios that show intermediate and critical regimens.

1.2. Contributions

In this study, using a high-resolution neuron model, we variate the synaptic time course and amplitude to look for inhibitory/excitatory hub ratios that exhibit critical dynamics. We found that intermediate and critical regimens are present for a wide range of excitatory/inhibitory hub ratios from 25% to 75% of excitatory hub neurons. Next, we introduce a mathematical expression to represent the metabolic synaptic cost that depends on the synaptic amplitude, time course, and connectivity; then, we compare the different dynamic regimes in terms of avalanche size versus metabolic cost. We found that subcritical and critical dynamics show a power-law function between avalanche size and metabolic cost, minimizing its costs. However, subcritical dynamics show only local activity or small avalanche sizes, whereas critical dynamics show the best compromise between global communication and a low metabolic cost. In contrast, supercritical dynamics show the break of power-law function between the avalanche size and the metabolic cost, displaying the highest metabolic costs. Moreover, reduced connectivity promotes an increase in the energetic cost of avalanche sizes.

2. Methods

2.1. Izhikevich Neuron Model

The Izhikevich model [32] is used to reproduce the dynamics of neural networks. The model comprises a two-dimensional system of ordinary differential equations defined by the following equations:

$$\begin{array}{c}\hfill \frac{dv}{dt}=h(0.04v^2+5v+140-u+I+s),\end{array}$$

(1)

$$\begin{array}{c}\hfill \frac{du}{dt}=h\left(a(bv-u)\right),\end{array}$$

(2)

$$\begin{array}{c}\hfill \mathrm{if}v\ge 30\mathrm{mV},\mathrm{then}\left\{\begin{array}{c}v=c\hfill \\ u=u+d,\hfill \end{array}\right.\end{array}$$

(3)

where v represents the membrane potential of the neuron (${\mathrm{Na}}^{+}$ current), t is time, u represents a membrane recovery variable, which accounts for the repolarization and activation of ${\mathrm{K}}^{+}$ current. 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. Equation (3) states that if membrane potential v is equal to or greater than 30 mV, the neuron is spiking, and then the values v and u are reset to $v=c$ and $u=u+d$. The value I represents a noisy external input that receives each of the neurons, which represents a sensory stimulus, and the value s represents the sum of the incoming potentials of its nearest neighbors when they fired. Finally, h is the size of the time step in ms. It is considered the original diverse repertory of values [32] that for excitatory neurons: $a_i=0.02$, $b_i=0.2$, $I_i=5r_i$, $c_i=-65+15r_i$ and $d_i=8-6r_i$, where $r_i$ is a random variable uniformly distributed in the interval [0, 1], and i is the neuron index. If $r_i=0$, the neurons exhibit regular spiking, and if $r_i=1$, the neurons display chattering behavior. For inhibitory neurons: $a_i=0.02+0.08r_i$, $b_i=0.25-0.05r_i$, $I_i=2r_i$, $c_i=-65$ and $d_i=2$. Here, if $r_i=0$, the inhibitory neurons show a slow recovery with low-threshold spiking, however, if $r_i=1$, the neurons show fast-spiking. The 2nd-order Runge–Kutta method with time steps h of 0.1 ms is used to solve the system.

$$\begin{array}{c}\hfill A_1=h(0.04v^2+5v+140-u+I+s),\end{array}$$

(4)

$$\begin{array}{c}\hfill B_1=ha(bv-u),\end{array}$$

(5)

$$\begin{array}{c}\hfill A_2=h[0.04{(v+\frac{A_1}{2})}^2+5(v+\frac{A_1}{2})+140-(u+\frac{B_1}{2})+I+s],\end{array}$$

(6)

$$\begin{array}{c}\hfill B_2=ha[b(v+\frac{A_1}{2})-(u+\frac{B_1}{2})],\end{array}$$

(7)

$$\begin{array}{c}\hfill \frac{dv}{dt}=A_2,\end{array}$$

(8)

$$\begin{array}{c}\hfill \frac{du}{dt}=B_2,\end{array}$$

(9)

where $A_i$ and $B_i$, i = 1, 2, represent the increments in the Eurler’s method, the $A_2$ and $B_2$ increments evaluate the function in a midpoint between the present value and the next value, offering a more accurate approximation of the functions in Equations (1) and (2). This numerical solver has been previously used for the Izikevich model [33]. A high-resolution Izhikevich model with time step sizes of $h=0.1$ ms was implemented in previous studies [15,34]. However, the amplitude of synaptic potentials used is in the range of 50 mV or larger, which represents no biological plausibility because real synaptic amplitudes are in the range of $0.01<w<10$ mV [35]. The necessity to use these amplitudes is due to the synaptic transmission occurs only in one time step, namely synaptic time course, $\tau$, is 0.1 ms. Moreover, synaptic transmission is considered a slow-wave activity in comparison with the action potential [29]. Real synaptic time courses, $\tau$, are in the range of $0.2<\tau <1.5$ ms [36]. In this study, different $\tau$ are considered, from 0.1 ms to 3 ms.

2.2. Hierarchical Network and Rich-Club Organization

The neurons are located in a hierarchical scale-free network proposed by Ravasz and Barabási [28]. The connectivity degree distribution $P\left(k\right)$ of the model follows a negative power-law function $P\left(k\right)\sim k^{-\gamma }$, with $\gamma =2.1$, which is concordant with the degree exponent found for functional neural networks: $\gamma$ = 2.0–2.2 [21]. The first step consists in constructing a complete network of five linked nodes; next, four more replicas are created. Finally, we connect four nodes of each replica cluster to one node in the original cluster (hub node); this results in a network of 25 nodes, including a hub. The second step consists of 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 of a network with 125 nodes and five hubs (Figure 1). In this study, eight network replicas, similar to the one created in step 2, were used to form a network with N = 1000 nodes. We select a network size of 1000 because scale-free networks like Ravasz and Barabási network present autosimilitude, meaning that we expect the functional properties to be replicated at larger sizes. We also test a 5000 node network to look for critical dynamics and metabolic costs. The 1000 network size allows us to probe a wide variety of synaptic amplitudes and synaptic time courses using a C++ program in a conventional computer. From the hierarchical configuration, we consider that $25\%$ of all edges comprise reciprocal connections, and $50\%$ are unidirectional. The rest of the edges of the Ravasz and Barabási model ($25\%$) are no connected to avoid complete subgraphs and allow common neural motifs [35], which are different directional connectivity triangles. The selection of the different kinds of connections is taken at random. This partial connectivity allows the different clusters and hubs not to be fully connected, promoting independence in the individual activity. Besides, generally, complex dynamics is reached at intermediate connectivity [19]. Different levels of connectivity are suitable for simulation. It is also implemented simulations with different reduced connectivity to observe the consequences in the metabolic cost. Here, the computational cost depends on the structural cost (number of nodes and connections) and the simulated metabolic cost. Both kinds of costs involve more computations.

Next, we connect each pair of hubs with probability $\kappa$. This configuration allows hubs to communicate with each other. Two major types of hubs are defined: global hubs that are the most connected nodes (see Figure 1) and local hubs that are only connected to one cluster of 25 nodes and the rest of the hubs. The rich-club organization is a main structural property of neural network systems where nodes with high connectivity degree tend to connect each other [22]. In order to detect the rich-club phenomenon in our configuration, we use the normalized rich-club coefficient ${\mathsf{\Phi }}_{\mathrm{norm}}$ [37], which is defined as the number of edges between pairs of hubs normalized by the number of edges between pairs of hubs in a null model network. The null model comprises the same degree distribution with a random connectivity between pairs of nodes. Our previous work has shown that for $\kappa >0.5$ this network shows legitimate rich-club organization [15]. We included the concept of hub index to describe a structural detail. First, we assign an index i for all neurons: i = 1, …, 1000. Next, because each cluster at step 1 comprises 25 nodes with a hub; the hub indexes are spaced by 25, so each hub index hub${}_i$ is a multiple of 25: $i=25,50,75,\dots ,1000$, then additionally to $\kappa$ probability, all the hub interconnections need to satisfy |(hub)${}_i-$(hub)${}_j|<5^4$, which accounts biological plausibility promoting the presence of provincial hubs [22]. Finally, the parameter $\eta$ is the excitatory hub population. For simplicity, we consider that one global hub ($k=100$ in the Ravasz-Barabási model) is equivalent to 5 local hubs ($k=20$). For example, the proposed network comprises eight replicas (1000 nodes), which include eight global hubs and 32 local hubs. If $\eta =0.5$, there would be four excitatory global hubs and 16 excitatory local hubs, or five excitatory global hubs and 11 local hubs, etc. The selection of the excitatory behavior of each hub is taken at random, meeting the above condition. For the remaining low degree nodes, we set 85% as excitatory and 15% as inhibitory following natural neuron populations [31]. For the different combinations of $\kappa$, $\eta$, $\tau$, and w, we run 500 time evolutions with independent random noise seeds. Each time evolution comprises ${10}^5$ time steps. In total, this represents $5\times {10}^7$ time steps for each combination.

2.3. Goodness-of-Fit Test

Given the data sets of avalanche sizes, we want to know whether they are reliable power-law distributed, namely P(size) $\approx c\times$(size)${}^{\alpha }$, where size is the avalanche size, P(size) indicates the probability of the avalanche size occurrence and c is a constant. An exponent value $\alpha =-1.5$ is found in recordings of rat cortex and neuronal cultures [9,38]. This exponent indicates critical dynamics that also corresponds to the one found in the critical branching process [39]. For each combination of parameters of the network, we obtain the avalanche size distribution of the accumulated data from all independent evolutions. Next, using the Kolmogorov–Smirnov test [40], we compare segments of ${10}^4$ sizes with surrogate data that comprise the same length, scaling exponent $\alpha$. We used the Kolmogorov–Smirnov test because it compares two distributions with no limitation to Gaussian distributions, so this method is suitable for comparing power-law distributions. It is used ${10}^4$ segments because lower length segments may introduce a biased in the p-value [41]. Finally, we obtain the average p-value from all segments. The Kolmogorov–Smirnov test is very sensitive to the detection of power-law distributions [40], and a p-value > 0.1 can be used to detect reliable power-law distributions [41].

3. Results

3.1. Criticality in Rich-Club Neural Networks

First, we set all excitatory and inhibitory synaptic amplitudes w at 5 mV and −5 mV, respectively. We variate the synaptic time course $\tau$, rich-club connectivity $\kappa$, and excitatory hub population $\eta$. In Figure 2, we show the avalanche size distribution for different rich-club connectivity $\kappa$ and excitatory hub population $\eta$. Figure 2 shows three different synaptic time courses: $\tau =0.5$ ms (top), $\tau =1.1$ ms (middle), and $\tau =1.4$ ms (down). For $\tau =1.1$ ms and $\eta =0.3,0.5$, we observe distributions near to $\alpha =-1.5$ for sizes $>10$, and we used this exponent as a guide for the eye. As $\tau$ becomes larger, the frequency of large avalanches increases and the majority of avalanche distributions come from power-law exponents $\alpha <-1.5$ (top) to $\alpha >-1.5$ (down). It is worth noting that for $\eta >0.5$ and $\kappa \ge 0.5$, the distributions indicate a high frequency of global avalanches. Furthermore, rich-club connectivity ($\kappa =0.9$) promotes the occurrence of global avalanches for $\tau =0.5$ ms and $\tau =1.1$ ms.

In order to obtain more evidence about the occurrence of critical dynamics for the different parameters, we analyzed three cases of E/I hub ratios for different combinations of w and $\tau$. In Figure 3 is shown three heatmaps of avalanche $\alpha$ exponents for a rich-club configuration $\kappa =0.9$ and three excitation populations: $\eta =0.25$ (left), $\eta =0.5$ (middle), and $\eta =0.75$ (right). Black grids represent critical $\alpha$ exponents ($-1.65\le \alpha \le -1.35$) whose avalanche distribution satisfies p-value $>0.1$ (see Section 2), strongly suggesting reliable critical dynamics (see Section 2). Moreover, we also show the distributions with p-value $>0.005$ [42] to demonstrate that power-law functions are scarce. In Figure 3, we also included critical exponents with $0.005<$ p-value $<0.1$ (brown grids), indicating suggestive or near to critical regimes. In the three cases, subcritical regimes ($\alpha <-1.35$) are observed for synaptic time courses $\tau <1.0$ ms with low w. For these values, the synaptic resources are poor, and the information cannot be propagated to the entire system. On the contrary, if the synaptic time course is large enough $\tau >2.0$ ms, the system shows supercritical regimes ($\alpha >-1.65$), indicating a high occurrence of large avalanches.

For all the cases, we observe critical regimes (black grids) with different $\tau$ and w depending on the excitatory hub ratio $\eta$. It is important to note that for $\eta =0.25$ case (Figure 3 left), as the synaptic time increases, the system passes through a narrow phase transition. Conversely, when the majority of hubs are excitatory ($\eta =0.75$), the heatmap shows a larger region that satisfies $-1.65<\alpha <-1.35$ but not a p-value $>0.005$ (black and brown grids). Although the regions in the heatmaps for critical dynamics are small, the overall results indicate that the system can show reliable critical dynamics for a wide range of excitatory/inhibitory hub populations depending on the synaptic time course and the synaptic amplitude. For example, in $\eta =0.25$ case (Figure 3 left), the system shows critical dynamics at $0.8<\tau <1.1$ ms and $5<w<9$ mV. Moreover, for an intermediate excitation ratio (Figure 3 middle), the system shows criticality at $\tau =1.4$ ms and $3<w<4$ mV. Finally, at a high excitation hub ratio (Figure 3 right), it shows criticality at $\tau =1.5,1.6$ ms and $2<w<3$ mV. It is worth noting that the values of $\tau$ and w found for critical dynamics correspond to real values of synaptic time course $0.2<\tau <1.5$ ms [36] and synaptic amplitude $0.01<w<10$ mV [35]. We also observe that there is not a p-value $>0.1$ for $\alpha <-1.65$ and $\alpha >-1.35$, indicating that reliable power-law functions occur only in critical states; this result is concomitant with previous observations [43].

3.2. Metabolic Cost of Dynamic Regimes

In the previous section, we observe that both supercritical and critical dynamics can show a wide range of avalanche sizes. However, we want to know if both dynamics involve the same metabolic cost for specific avalanche sizes. We evaluate the two principal metabolic expenditures: action potentials and synaptic potentials [20,44]. The action potentials are represented by the number of times neurons are activated, where each neuron activation represents a spike. The synaptic cost of a single neuron activation is the sum of the outgoing synaptic potentials that depend on the number of outgoing links $k_{\mathrm{out}}$, synaptic time course $\tau$, and synaptic amplitude w (see Figure 4). As a result, the synaptic cost of a single neuron per spike is:

$$\begin{array}{c}\hfill S_i=\left|w\right|\frac{\tau }{h}{k}_{\mathrm{out}},\end{array}$$

(10)

where the absolute value of w is taken due to the negative synaptic potentials of inhibitory neurons. In this way, hub neurons (high $k_{\mathrm{out}}$) are high metabolic cost units. Next, we define the synaptic cost ($SC$) of a given avalanche as:

$$\begin{array}{c}\hfill SC=\sum _{t_i}^{t_f}(\sum _{i=1}^N{S}_i\left(t\right)),\end{array}$$

(11)

where the first sum represents all the time steps that comprise the avalanche and the second sum represents all the neurons in the network. Besides, since the sum of incoming links is equal to the sum of outgoing links, the $SC$ also can be represented as the sum of incoming synaptic potentials $s_i$. An ideal metabolic cost comprises no redundancy. As a result, the minimum cost configuration is an avalanche size equal to the number of spikes. Moreover, if the successive activation of neurons does not strongly depend on the $k_{\mathrm{out}}$, the minimum synaptic cost could be approximated as:

$$\begin{array}{c}\hfill S{C}_{\min }\left(\mathrm{size}\right)\approx {\overline{k}}_{\mathrm{out}}\left|w\right|\frac{\tau }{h}\times {\left(\mathrm{spikes}\right)}^{\beta },\end{array}$$

(12)

where ${\overline{k}}_{\mathrm{out}}$ is the mean value of $k_{\mathrm{out}}$ and $\beta =1$. The minimum expect cost is an avalanche size equal to the number of spikes and a linear relation between avalanche size and $SC$. A value $\beta >>1$ would indicate high redundancy and overactivation.

We compare the size of the avalanche versus the average amount of synaptic cost that the network needs. The top of Figure 5 shows the synaptic metabolic cost for subcritical, critical and supercritical configurations. We also consider three cases: left, $\eta =0.25$; middle $\eta =0.5$ and right $\eta =0.75$. Insets show a magnification of the largest avalanche sizes for critical dynamics in each case. For $\eta =0.25$: critical ($w=5.0$ mV, $\tau =1.1$ ms), subcritical ($w=4.0$ mV, $\tau =0.8$ ms), supercritical ($w=7.0$ mV, $\tau =1.6$ ms). For $\eta =0.5$: critical ($w=3.5$ mV, $\tau =1.3$ ms), subcritical ($w=2.5$ mV, $\tau =0.8$ ms), supercritical ($w=5.5$ mV, $\tau =1.8$ ms). For $\eta =0.75$: critical ($w=2.5$ mV, $\tau =1.5$ ms), subcritical ($w=2.0$ mV, $\tau =0.8$ ms), supercritical ($w=3.0$ mV, $\tau =2.4$ ms).

In the three $\eta$ cases, we observe that critical states show a power-law relation between the synaptic cost and the avalanche size: $SC\left(\mathrm{size}\right)\approx {\overline{k}}_{\mathrm{out}}\left|w\right|\frac{\tau }{h}\times {\left(\mathrm{spikes}\right)}^{\beta }$, where $\beta \approx 1.2$, indicating low activation redundancy (Equation (12)). This result is reminiscent of the real allometric metabolic consumption of the brain [20]. Moreover, subcritical regimes are unable to show large avalanches despite the excitatory hub population. On the contrary, supercritical regimes show large avalanche sizes in the three cases. This is attributed to the high amount of synaptic resources that lead to spontaneous activation and propagation. However, supercritical states show the highest metabolic costs: $\beta \approx 8.0$ for sizes $>700$, and they also lose the linear relation for large sizes. For critical regimes, we observed different maximum avalanche sizes depending on the excitatory hub population. For $\eta =0.25$, we observe that the maximum avalanche size is 300, this is due to 75% of hubs are inhibitory, so a majority of inhibitory hubs cannot propagate the information to the entire system. However, for $\eta =0.75$ (75% of excitatory hubs), the maximum size is around 600 because excitatory rich-club can propagate information to the majority of the neurons.

Next, we evaluate the number of spikes (representing action potentials) involved in the avalanche occurrence (see Figure 5 down). We observed that for subcritical and critical regimes, the spike cost also follows a power-law and linear relation with the avalanche size. Besides, supercritical regimes lose this linear relationship for large avalanches. As a guide for the eye, an ideal association: the number of spikes equal to avalanche size is represented by a solid line. This configuration involves no redundancy, namely, all the neurons spike only once during the avalanche occurrence. The results suggest that subcritical and critical regimes display a behavior near to a minimum redundancy. Moreover, critical regimes avoid high metabolic expenditures and exhibit the best compromise between network communication and low metabolic cost.

We also obtained the avalanche size, synaptic cost, and number of spikes for a system with N = 5000 to probe that critical dynamics remains for larger sizes. In Figure 6a we show avalanche size distribution for a system with $\eta =0.5$, $w=4.0$ mV and $\tau =1.3$ ms. The dash line represents a power-law function with $\alpha =-1.56$. This combination of $\eta$ and w values also corresponds to critical dynamics with N = 1000 (see Figure 3). We observe that the power-law distribution of avalanche sizes extends to size $=1250$. This value represents three times higher sizes in comparison with N = 1000 case (see Figure 5b,e). In Figure 6b,c we show the synaptic cost and number of spikes, respectively, of the critical configuration. We note that the synaptic cost and number of spikes follow power-law relations with the avalanche size. The dash line in Figure 6b represents a $\beta =1.2$ in correspondence with N = 1000 case. The dash line in Figure 6c represents an ideal relation with a number of spikes equal avalanche size. Synaptic cost and number of spikes indicate that critical dynamics in a larger system outperform supercritical dynamics (see Figure 5) using a much lower metabolic cost.

3.3. Drosophila Network

Next, in order to compare the behavior found in our model with a real neural structure, we implement the Izikevich neural model in a mesoscale connectome of the Drosophila medulla. Data was extracted from the Network Data Repository [45]. The Drosophila network shows hierarchical organization, small world, and rich-club properties [46]. These properties are also found in the hierarchical model. In Figure 7a we show the conectivity matrix of the Drosophila network. The network comprises 1472 nodes/neurons. Figure 7b shows the indegree and outdegree distribution of the network, we observed that both distributions follow a power-law function: $\gamma =1.4$ for indegree and $\gamma =1.7$ for outdegree distribution. In the simulation, we set a balanced excitation/inhibition ratio $\eta =0.5$ for the $4\%$ most connected population (rich-club). For the rest of the neurons, we set 85 % as excitatory and 15% as inhibitory. Figure 7c shows the heatmap of scaling exponents $\alpha$. We observe that reliable critical dynamics is reached at $0.9<\tau <1.1$ ms and $4.0<w<6.7$ mV. These values are similar to those found in critical dynamics in our model. Figure 7d shows avalanche size distributions for critical ($w=4.3$ mV and $\tau =1.1$ ms), subcritical ($w=3.6$ mV and $\tau =0.7$ ms), and supercritical dynamics ($w=7.8$ mV and $\tau =1.5$ ms), solid line represents a power-law function with $\alpha =1.56$. Figure 7e,f show the synaptic cost and number of spikes for three dynamics. We observe that in accordance with the Rabasz-Barabási model, critical dynamics keep power-law relations, $\beta \approx 1.2$, in synaptic cost vs. size (Figure 7e), and number of spikes similar to size (Figure 7f), indicating that the Drosophila network shows an economic metabolic design, together with a wide range of dynamic patterns.

3.4. Erdos–Renyi Network

To test that critical behavior is due to structural connectivity, i.e., hierarchical organization, scale-free, and rich-club connectivity, we probe the Erdos–Renyi network [47] in the Izhikevich model. For comparison, we set the same number of nodes and connections as in the Drosophila network. Figure 8a shows the Erdos–Renyi connectivity matrix. This structure, unlike that of Drosophila, shows no clustering. Figure 8b shows the indegree and outdegree connectivity. Figure 8c shows the heatmap of scaling exponents $\alpha$. We observe that there are not $\alpha$ exponents that meet p-value $>0.1$. Besides, the region of $\alpha >1.5$ is reached for lower values of $\tau$ and w in comparison with Drosophila and Hierarchical networks. This behavior may be attributed to the lack of clustering. The internal structure cannot stop the spread of activity; rather, the interconnected network promotes global activity. Figure 8d show the avalanche distributions for the w and $\tau$ pairs found in Drosophila network for critical ($w=4.3$ mV and $\tau =1.1$ ms), subcritical ($w=3.6$ mV and $\tau =0.7$ ms) and supercritical dynamics ($w=7.8$ mV and $\tau =1.5$ ms). We observe that the distributions do not follow power-law relations. For supercritical configuration ($w=7.8$ mV and $\tau =1.5$ ms combination) the avalanche sizes show either small size $<50$ or global avalanches size $>1000$. Figure 8e,f show the metabolic cost and number of spikes, respectively. An intermediate stability configuration also shows the best compromise between information propagation, dynamic repertoire, and economic metabolic cost.

3.5. Disconnected Network

We analyze the effects of high disconnection in the power-law scaling of metabolic costs. We disconnect the network for both low degree connectivity and rich-club connectivity $\kappa$. The synaptic cost and number of spikes for disconnected networks with 1/4 and 1/16 of the total links are shown in Figure 9. We set $\eta =0.5$ and $\tau =1.3$ ms. Since the disconnection configuration displays only small avalanches for $w<7$ mV (results not shown), we set $w=10$ mV for 1/4 case, and $w=50$ mV for 1/16 case to retrieve global avalanches. We observe that in both cases, the power-law relation is broken at a certain avalanche size depending on the disconnection degree. The larger the disconnection, the lower the cross-over size where the power-law relation brakes. These results are in accordance with previous experiments about the sensitivity of energetic cost in edge deletion [48]. We also observe that for small avalanche sizes, the synaptic cost and number of spikes are smaller than the critical case (represented by a solid line). However, global avalanches are much more expensive, with a higher exponent $\beta$ and a higher number of spikes. Although this configuration does not represent a minimum cost for global avalanches, and the disconnection promoted the break of the power-law, the metabolic costs follow Equation (12) where $k_{\mathrm{out}}$ is smaller, but the redundancy $\beta$ is higher.

3.6. Critical Neural Networks Minimize Metabolic Cost

Finally, to test if critical dynamics show the lowest synaptic cost and number of spikes for global avalanche sizes, sizes $>125$, we compare the synaptic cost (Figure 10a–c) and number of spikes (Figure 10d–f) of avalanche sizes among all pairs of w and $\tau$ with N = 1000. We plot the combination of w and $\tau$ that minimizes either the synaptic cost (Figure 10a–c) or number of spikes (Figure 10d–f) for a given avalanche size. We also consider three excitatory populations. As an additional criterion, we do not include avalanche size sets whose synaptic costs $>{10}^6$ mV or number of spikes $>1200$. These values correspond to cross-over points found for supercritical regimes at intermediate excitatory/inhibitory ratio: $\eta =0.5$ (Figure 5). Cross-over points were found using the method in [49]. White grids represent critical and near to critical exponents. For the three cases, we observe that the synaptic cost and number of spikes are minimized for intermediate values of $\tau$ and w, i.e., intermediate stability regimes, including critical and near to critical regimes. Similar avalanche sizes are minimized for different $\tau$ and w combinations, where the increment in the other replaces a decrease in one parameter. Both slightly subcritical (lower w and $\tau$) and slightly supercritical regimes (higher w and $\tau$) also minimize metabolic cost for different avalanche sizes. For a balanced hub excitation ratio, $\eta =0.5$, we observe that critical regimes represent the only region that minimizes both SC (Figure 10b) and number of spikes at the same time (Figure 10e).

4. Discussion

Neural networks need both excitatory and inhibitory hub neurons. Inhibitory hubs control the occurrence of overactivation that can promote neuropathologies like seizures or excitotoxicity. Moreover, inhibitory hubs have been found in the hippocampus of the rat [26] and in the connectome of the C. Elegans [50]. On the other hand, excitatory hubs not only facilitate the integration of local information but also send this information to other brain areas. In addition, the rich-club organization mainly supports the information processing in the brain [51]. In this work, we look for an excitatory/inhibitory hub ratio that can display critical dynamics variating the synaptic time course and synaptic amplitude. The neuronal model is used for three different network structures: a scale-free hierarchical network, Ravasz-Barabási; a random configuration, Erdos–Renyi; and a Drosophila network. It is important to mention the Izhikevich model’s limitations, such as no direct measure of conductance and no high-shaped fidelity of the action potential [52]. Besides, the model is susceptible to the integration time step. A reduction in the time step h requires adjusting the synaptic time course or amplitude [15]. However, the Izhikevich model allows for the simulation of a large number of time steps with reduced computational cost.

Rich-club configuration with clustering promotes a wide-range of avalanche sizes. In contrast, a random configuration exhibits either local or global propagation, representing a clear disadvantage in the dynamic repertoire. Moreover, connectivity is tightly related to metabolic cost: Equation (10), for example, the propagation activity of a neuron with high outdegree connectivity involves the release of neurotransmitters at multiple presynaptic terminals (Figure 4) but guarantees the spread of information. Previous studies showed that hub neurons/regions are related to high activity and metabolism, suggesting that connectivity k is an indicator of metabolism [53,54]. Besides, the quantity of neurotransmitters or the number of molecules released reflects a principal metabolic cost [44], which in our model is regulated by the amplitude of synaptic potentials w and synaptic time course $\tau$. In contrast, high levels of disconnection (Figure 9) promotes a supralinear expenditure of metabolic costs because the network cannot communicate and needs over-activity to reach partial global communication.

The results show that depending on the synaptic amplitude and synaptic time course, critical dynamics can coexist with a wide range of excitatory/inhibitory hub ratios. It is possible to find either a minority or majority of excitatory hubs in a rich-club configuration showing critical dynamics. Furthermore, this work found that critical dynamics and intermediate stability regimes are associated with a reduced and minimum metabolic cost, which suggests that criticality optimizes the metabolic sources to achieve global information transmission. In this way, critical dynamics could serve as a dynamic attractor because of its compromise between global communication and low metabolic cost. These results may explain previous results that neural networks are most of the time in a critical regime. In contrast, subcritical and supercritical regimes are present during a reduced time [11]. The results are consistent with the free energy brain hypothesis [55], which establishes that the brain operates under a low energy consumption, minimizes error prediction, and avoids disorder (supercritical states). Moreover, the results are also concordant with previous observations about the fact that fractal fluctuations related to critical states support information propagation in the body [56], and excitation/inhibition balance maximizes energy efficiency [57] and leads to criticality [58]. The results also indicate that critical networks arise from a balance between synaptic amplitude and synaptic time course. Finally, critical dynamics has been related to a maximum dynamic repertoire [3,4], high information storage [59,60], parametric dimension reduction [61], robustness [14,16,62], and computational efficiency [63]. These properties joined to a reduced and minimum metabolic cost suggest that criticality in biology systems might be shaped by the evolution to exhibit a common trade-off among all these advantages.

Code Availability

The source code of the neural network evolution (in C++), and the number of spikes vs. avalanche size and synaptic cost analysis (in python) is published at https://github.com/danielvelaguil/metabolic-consume-of-neural-networks (accessed on 5 February 2021).