Self-organized criticality (SOC) and stochastic oscillations (SOs) are two theoretically contradictory phenomena that are suggested to coexist in the brain. Recently it has been shown that an accumulation-release process like sandpile dynamics can generate SOC and SOs simultaneously. We considered the effect of the network structure on this coexistence and showed that the sandpile dynamics on a small-world network can produce two power law regimes along with two groups of SOs—two peaks in the power spectrum of the generated signal simultaneously. We also showed that external stimuli in the sandpile dynamics do not affect the coexistence of SOC and SOs but increase the frequency of SOs, which is consistent with our knowledge of the brain.
The brain, one of the most interesting and complex dynamical systems, is hypothesized to be a self-organized critical system due to the observed apparent power law in neurophysiological data (Bressler & Kelso, 2001; Beggs & Plenz, 2003; Rubinov & Sporns, 2010; Schuster, 2014). In a dynamical system with many interacting components, complexity arises in a state between ordered and disordered states (Beggs & Timme, 2012). By defining an order parameter in a dynamical system, there would be a phase transition between ordered and disordered states. The point at which phase transition occurs is called the critical point (Goldenfeld, 1992; Yeomans, 1992). Self-organized criticality (SOC) is a property of natural complex dynamical systems that work near the critical point (Bak, Tang, Wiesenfeld, 1987; Beggs & Plenz, 2003). Many variables in systems with the SOC property can be described by a power law function (Goldenfeld, 1992; Yeomans, 1992; Stanley, 1999; Nishimori & Ortiz, 2010). Power law functions show no characteristic scale, and therefore the data have a fractal structure.
Some studies report that networks of neurons can produce neuronal avalanches with power law distributions (Beggs & Plenz, 2003; Friedman et al., 2012). A neuronal avalanche is a sequence of neuronal activities in which the time interval between those activities is less than a specific constant (Benayoun, Cowan, van Drongelen, & Wallace, 2010).
Brain activities are characterized by stochastic oscillations (SOs; Buzsáki, 2006). These oscillations belong to repetitive activities with an inherently typical scale believed to be one of the essential features of the brain (Buzsáki, 2006).
The thought-provoking question is, How does the brain produce these two types of inconsistent properties, SOC and SO, while one indicates a scale-free behavior and the other demonstrates a typical scale? Computational modeling approach may help in this regard. As an example, Wang et al. (2016) have tried to answer this question by drawing on a new viewpoint on the prototype sandpile model—the Bak-Tang-Wiesenfeld (BTW) model.
The sandpile model can easily be interpreted as a simplified model of neuronal dynamics (Abbott, 1999; Priesemann et al., 2014; Zarepour, Niry, & Valizadeh, 2015). Neurons act like leaky integrators with nonlinear interactions (Dayan & Abbott, 2003). Each individual neuron integrates the received synaptic potentials from other neurons. If the potential reaches a threshold, it may excite another neuron and may trigger an avalanche. This process can be seen in many other complex systems with an accumulation-release process (Malamud, Morein, & Turcotte, 1998; Paczuski, Maslov, & Bak, 1996; Harris, 2002). The evolution rule of each cell in sandpile model is similar to the evolution rule of a nonleaky integrate-and-fire neuron (Abbott, 1999; Priesemann et al., 2014). This model acts like a complex nonequilibrium system with a simple interaction of agents. In the sandpile model, the individual components with specific thresholds are interacting to produce so-called avalanches and progress toward a critical state (Bak, 1996). Wang et al. showed that in a regular, small, square lattice, sandpile dynamics with weak external stimuli can produce SO and SOC simultaneously.
However, the brain network structure is not regular, and some studies report that structural and functional brain networks have a small-world property (Sporns, Chialvo, Kaiser, & Hilgetag, 2004; Bassett & Bullmore, 2006; Reijneveld, Ponten, Berendse, & Stam, 2007; Stam & Reijneveld, 2007). In a small-world network, the characteristic path length is small in comparison with a regular network, and the clustering coefficient is large in comparison to a random one. It may be concluded that a regular topology like a square lattice is not appropriate for modeling brain networks. Studies on neurodynamical models have shown that different topologies of neuronal connections can generate different dynamical modes (Jirsa & Kelso, 2000; Sporns & Tononi, 2001; Sporns, 2004). Also, it has been shown that small-world topology is correlated with cognitive performance (Stam, Jones, Nolte, Breakspear, & Scheltens, 2007). Moreover, there are reports that the brain has a small modular structural network (a few hundred neurons) and densely connected spatial clusters that imply a high clustering coefficient (Mountcastle, 1997; Buxhoeveden & Casanova, 2002; Hilgetag & Kaiser, 2004). It seems that the study of sandpile dynamics on a network with both the small-world property and a high clustering coefficient, both brain network features, will be useful to understand normal brain behavior.
In recent research, some modifications have been imposed on the BTW model to study the small-world property in self-organized critical systems with square lattice structures (De Arcangelis & Herrmann, 2002; Chen, Chiao, Lee, Cheng, & Wu, 2008; Hoore & Moghimi-Araghi, 2013), but the clustering coefficient is almost zero in square lattices.
In this letter, we study the effect of network structure on the coexistence of SOs and SOC in a model of a brain network with sandpile dynamics. We take a ring lattice network and rewire it with the standard Watts-Strogatz protocol (Watts & Strogatz, 1998) to achieve networks with small-worldness and high clustering. We propose a range of network parameters to have both the SOC and SOs in our suggested brain dynamical model.
2.1 Network Topology
We create a small-world network using the Watts-Strogatz algorithm (Watts & Strogatz, 1998) with nodes, each of which interacts with its neighbors. A binary adjacency matrix defines the neighborhood of each node. In the first step, the network is a ring lattice with the same degree for all nodes, which is similar to a one-dimensional lattice with periodic boundary for a small average degree (. Next, links were rewired by a fixed rewiring probability (. If is too small (around 0.001), the network is regular, with a high characteristic path length and a high clustering coefficient. When is large (around 1), the network will be a random one, with a low characteristic path length and low clustering coefficient. In a specific range of (0.01 0.2), the network has a small-world property with a low characteristic path length and high clustering coefficient.
In order to study the effect of the network structure on the coexistence of SO and SOCs, we change the rewiring probability (, average degree (, and network size ( and evaluate the existence of SOC and SOs for different sets of parameters.
2.2 Network Dynamics
We apply a sandpile model on different types of networks. At first, the height of the sand in all nodes is zero; then a constant value ( is added to the height of randomly selected node. If the height of sand in that node ( exceeds its threshold (, the node will be unstable and collapse. Therefore, its height is subtracted by , and the height of its neighbors is increased by one. If there is no site to collapse, another node will be selected randomly, and its height will be increased by .
The height of each node in a sandpile model can be increased by two procedures: collapsing its neighbors as an internal factor or adding sand ( as an external stimulus.
By collapsing one node, the neighbors also may exceed their threshold and collapse to their neighbors and so on. We define an avalanche as a sequence of collapses in different nodes. The size of an avalanche is the number of collapsed nodes. In a critical system, the distribution of avalanche size shows power law behavior. In each time step, we can record the number of the superthresholded sites (sites that exceed their threshold) in a time series and extract frequency properties of this signal by its power spectrum (PS). Conceptually, this time series is similar to neural population recording.
2.3 Distribution Fitting
For this study, we focused on the effect of network structure and external stimuli on the coexistence of SOs and SOC. The network parameters in this research are the rewiring probability (, average degree (, and network size (. Average degree ranges from 4 to 32 and rewiring probability is chosen as 0.001, 0.01, 0.05, 0.1 and 0.5. Network size ranges from 400 to 2500. The external stimuli parameter (i.e., is chosen as 0.1, 0.25, 0.5, 0.75, 1, 2, 3, and 4.
3.1 Coexistence of SOC and SOs
In order to study the effects of small-worldness, the distribution functions of avalanche size and power spectrum of superthresholded sites time series (PSS) are depicted in Figure 1. In this figure, the parameters are plotted for three different types of networks with different average degrees and = 2500: (1) a regular network ( = 0.001), (2) a small-world network ( = 0.1), and (3) a relatively random network ( = 0.5). Previous studies showed that BTW dynamics cannot generate SOC in one-dimensional systems (Bak et al., 1987; Dhar, 1999; Ruelle & Sen, 1992), but exhibited criticality in high-dimension lattices (Bak et al., 1987; Bak, Tang, & Wiesenfeld, 1988), although there are extended sandpile dynamics with critical behavior in one dimension (Pinho, Prado, & Salinas, 1997; Bhaumik & Santra, 2016). However, Figures 1a and 1d show that in a regular ring lattice network with the low mean degree ( = 2 or 4), there is neither SOC nor SOs. It can be explained as follows. In our model, leaky nodes are distributed throughout the network, which is more physiologically plausible. Thus, there is more leakage than in other one-dimensional models, and dissipation of sand is more likely. Therefore, the probability of large avalanches will drop quickly (see Figure 1a). A regular ring lattice with a large mean degree is not one-dimensional (Dehmer, Mowshowitz, & Emmert-Streib, 2013); it will have mean-field properties, and the threshold of sites will be larger and the amount of moving sand in each toppling is large. Thus, more sites can contribute to making large avalanches. In this case, avalanche distributions are bimodal and demonstrate the existence of two types of avalanches: regular avalanches with partial power law distribution and large avalanches with exponential decrease.
As illustrated in Figures 1b and 1e, in a small-world network, both power law behavior and sharp peaks in PSS emerge in denser networks (i.e., networks with larger . In the case of random networks (see Figures 1c and 1f), power law behavior and sharp peaks in the PSS exist for any average degree of the networks. Also, there are two peaks in the PSS (see Figures 1e and 1f). For larger values of two peaks can be generated in a smaller .
3.2 Two Regimes in Avalanche Dynamics and Oscillations
As shown in Figures 1b and 1e, in small-world networks, the avalanche size distribution has two power law regimes and PSS has two peaks in some cases. To elucidate these phenomena, two distinct power laws with exponents and are illustrated in Figure 2a. The first power law has an ordinary exponent in the Abelian sandpile model ( = 1) and the second one has a mean-field exponent ( = 1.5) affected by long-range connection. The probability of contribution of these long-range bonds in the smaller avalanches is lower, and small avalanches remain local; thus their statistics are similar to the BTW model. However, the grains of sand in large avalanches can spread throughout the network by long-range connections (their existence is a characteristic of small-world structure) and reach leaky nodes more easily; hence, dissipation would be larger. Therefore, the large avalanches will occur less frequently (i.e., the power law should have a larger exponent). The inset plot in Figure 2a shows that the two exponents remain unchanged in small-world networks (0.05 0.3).
The two regimes reach together at the crossover point, and its corresponding size ( is calculated by estimating the slope of two regimes with the MLE method. Since the crossover point belongs to both regimes and the second regime is affected by long-range bonds, would be a function of . Previous studies proposed that changes as a power of — (Lübeck & Usadel, 1997; Bhaumik & Santra, 2013, 2016; Hoore & Moghimi-Araghi, 2013). The crossover size for a small-world regime is sketched in a bilogarithmic coordinate in the inset of Figure 2b and shows a straight line with a slope = 0.85, which is consistent with other studies (Hoore & Moghimi-Araghi, 2013; Bhaumik & Santra, 2016). By increasing , the two regimes converge and finally become one power law in the case of random networks.
The PSS of small-world networks shows two peaks, indicating two modes of network activity: a sharp peak in low frequencies and a smaller one in higher frequencies (see Figure 2c). To demonstrate the relation between large avalanches and sharp peaks in PSS, the correlation between the sizes of successive avalanches is shown in the inset of Figure 2c. The first values are negative since after the large avalanches, the smaller ones occur. The value of correlation has a significant peak for a specific lag (, which is quite distinct from the correlation of the shuffled version. This means that the large avalanches probably occur after the avalanche because the large values in correlation are associated with larger values in the given data. Corresponding frequency to (i.e., ) can be calculated using the average time between avalanches. The calculated frequency ( is consonant with the frequency of the sharp peak in PSS (see the lower inset in Figure 2c). For example, in a small-world network with size 2500 and = 16 and = 4, the specific lag is = 49, the corresponding frequency is = (1.84 0.23) 10, and the frequency of the sharp peak in PSS is = (1.49 0.12) 10. Therefore, the origin of the sharp peak is the correlation of large avalanches.
3.3 External Stimuli Affect SOs through Waiting Times
In this study, the external stimuli, are applied only in the quiescent state to the system; thus, there would be a separation of timescales (STS). STS is a necessary condition for a model to be self-organized critical (Bak et al., 1987; Drossel & Schwabl, 1992; Clar, Drossel, & Schwabl, 1996; Dickman, Muñoz, Vespignani, & Zapperi, 2000; Pruessner, 2012; Priesemann et al., 2014). Therefore, we study the effects of external stimuli and STS on the coexistence of SOC and SOs. The distributions of avalanche size, avalanche duration, waiting time (WT), and PSS for different in a small-world network ( = 0.1) are illustrated in Figures 3a, 3c, and 3e, respectively. In this case, for each , the avalanche distribution of size and duration is the power law. Changing external stimuli does not affect avalanche distributions in the investigated range—that is, z equal to 0.1 up to 4 (see Figures 3a and 3c); thus, power law behavior is invariant of this value.
Figure 3c shows that WT is exponentially distributed, as seen previously (Boffetta, Carbone, Giuliani, Veltri, & Vulpiani, 1999). This indicates that there is no correlation between the occurrence time of successive avalanches. It is natural that growth in leads to a decrease in WT and STS, and we show that average of WTs is inversely proportional to the external stimulus strength (see the inset of Figure 3c). Furthermore, an increase in can displace the peak in PSS to larger values. Hence, network activities will occur more frequently.
To show that the emergence of oscillations is not a consequence of WT, zero values in the superthresholded time series, which represent WT, are excluded from the time series, and the PS of the reduced data is sketched in Figures 3c and 3d for different and as inset plots. The existence of peaks in the PS of reduced data shows that the emergence of SOs does not result from WT. Since the location of the peak is independent of in the reduced time series (see the inset of Figure 3e), one can conclude that changes the frequency of oscillations through changing WT.
Figures 3b and 3d show that power law behavior is not affected by network size; instead, the network size imposes a cutoff on the distribution of size and duration, which reflects the system's finite size (Kadanoff, Nagel, Wu, & Zhou,1989; Ktitarev, Lübeck, Grassberger, & Priezzhev, 2000).
However, Figure 3f shows that network size can change the location of the peak in PSS. The impact of external stimuli and network size on the frequency of the PSS peak is illustrated in Figure 4. It shows that peak frequency will increase by increasing the . For smaller network sizes, the frequency shift is faster, and the peak frequency has larger values. This figure demonstrates how peak frequency depends on network size and external stimuli.
Excitability is a key feature of neurons by which they interact with other neurons through synaptic connections. Membrane potential has a rest value, usually 50 and 90 mV. The membrane has some voltage-gated ion channels, which are responsible for producing the cell impulse (i.e., action potential). When the cell is stimulated up to a threshold, voltage-gated sodium channels open, and the potential reaches up to 30 mV. Then voltage-gated potassium channels open, and the membrane potential returns to its rest value. This accumulation-release process is a primary procedure in neuronal activities.
The neuronal accumulation-release process could be considered the same as the sandpile model proposed to explain SOC. In a sandpile model with a small system, there is a temporal correlation between large avalanches (avalanches comparable with system size; Wang et al., 2016). This temporal correlation could be explained as follows. After a global release in the system, it takes time to accumulate for another large avalanche. By considering the number of superthresholded sites in each time step, there is a time series that its PS will have a peak at low frequencies (Wang et al., 2016). This peak in PS shows the existence of SOs in a sandpile model. In this study, we investigated the effects of network structure and external stimuli level on the coexistence of SOC and SOs by simulating a sandpile dynamic.
The small-sized systems in our study may be interpreted as small subregions of the brain because their nonleaky neurons are usually similar. Some of these neurons have connections with neurons outside the subregion, which represent boundary neurons. Because we consider only the intraconnections, the boundary neurons should be low-degree nodes in comparison to other neurons in the subregion. If we want dissipation of sand in the boundaries of the system (like the BTW model), it seems that the leak term should be assigned to low-degree nodes. Also, the nodes with more connections to the outside are potentially more capable of moving sand to the outside and causing more dissipation. Therefore, we assigned an inversely proportional leakiness to some low-degree nodes.
4.1 Effects of Structure
Previously it has been shown that SOC exists in a sandpile model on an ordinary square lattice that is a type of regular network (Wang et al., 2016). However, our results showed that in a ring lattice, SOC does not exist for any value of the average degree (see Figures 1a and 1d). It may be explained as follows. In a regular ring lattice, there are no long-range interactions, and dissipations are more effective. Therefore, activities have no chance of propagating throughout the network; the probability of large avalanches will decrease, and the distribution deviates from the power law. Some reports show that during an epileptic seizure, functional and structural brain networks change to a regular structure with a larger characteristic path length and larger clustering coefficient (Koepp & Woermann, 2005; Ponten, Bartolomei, & Stam, 2007; Schindler, Bialonski, Horstmann, Elger, & Lehnertz, 2008; Kramer et al., 2010). On the other hand, Meisel, Storch, Hallmeyer-Elgner, Bullmore, and Gross (2012) showed that the avalanche size distribution deviates from the power law during an epileptic seizure. Our model suggests that deviation of avalanche size distribution from the power law during an epileptic seizure may be related to the regular structure of brain networks during this period. Furthermore, Stam et al. (2007) showed that Alzheimer's disease (AD) is correlated with a loss of the small word property in brain networks. Therefore, we may expect that the dynamics of neuronal avalanches in Alzheimer's disease deviates from the power law.
In small-world networks, the coexistence of SOs and SOC emerges in networks with a high average degree, 8 (see Figures 1b and 1e). This is consistent with some studies suggesting that brain networks have densely modular structures (Mountcastle, 1997; Hilgetag & Kaiser, 2004) and small-world property (Sporns et al., 2004; Bassett & Bullmore, 2006; Reijneveld et al., 2007; Stam & Reijneveld, 2007). In the case of small-world networks, avalanche size distribution has two distinct power law regimes: one with the BTW (sandpile model on a regular lattice) slope ( 1) and the other with another slope related to the small-world property and mean-field solution ( 1.5) (see Figure 2a). The second regime is in consonance with experimental data showing that the exponent of neuronal avalanche size distribution is 1.5 (Mountcastle, 1997; Beggs & Plenz, 2003; Gireesh & Plenz, 2008). Furthermore, we showed that crossover avalanche size is a power law function of rewiring probability with slope 0.85.
In a small-world network with 16 there are two peaks in the power spectrum (see Figure 2c), related to two frequency bands. In the brain, SOs in different frequency bands are related to a range of cognitive and perceptual processes (Varela, Lachaux, Rodriguez, & Martinerie, 2001; Steriade, 2006). Different frequency bands can coexist in the brain and interact with each other (Steriade, 2006; Csicsvari, Jamieson, Wise, & Buzsáki, 2003). The interaction of different frequency bands is an interesting topic and a broad area of research that is known as cross-frequency coupling (CFC). CFC demonstrates that different groups of neurons are coordinated by distinct frequencies (Canolty, Cadieu, Koepsell, Knight, & Carmena, 2012). CFC is a general mechanism for all brain rhythms (Schroeder & Lakatos, 2009; Buzsáki & Watson, 2012) and is related to behavioral performance during learning (Tort, Komorowski, Manns, Kopell, & Eichenbaum, 2009), neuronal communication, and information encoding (Lisman, 2005; Jensen & Colgin, 2007). Our model has more frequency content by the coexistence of two rhythms and their possible interactions. This result was not shown in Wang et al. (2016) in which the model exhibited just one broad peak in PS.
4.2 Effects of External Stimuli
Slow oscillations are obviously related to brain resting state with low external stimuli (i.e., closed eyes) (Kirschfeld, 2005). These slow oscillations have been known as alpha oscillations (8–13 Hz). The amplitude of alpha oscillations in the power spectrum of electroencephalography (EEG) will decrease and vanish when the eyes open (Kirschfeld, 2005). Also, large-scale synchronization in high-frequency bands [(beta (12–30 Hz) and gamma (30–70 Hz))] is associated with large external stimuli and cognitive tasks (Varela et al., 2001). Proper brain dynamical models that generate brain rhythms should be affected by changing external stimuli. The external stimuli can indicate the effect of sensory inputs (out of the neural system) or spontaneous internal activities in the brain (Priesemann et al., 2014).
In this work, we studied the effect of external stimuli on the rhythm variations of the proposed neuronal model. We showed that an external stimulus ( does not affect criticality (see Figure 3a). However, changing leads to changes in the peak frequency in the PSS (see Figure 4). Increasing causes nodes to reach the threshold faster. In this case, after a global release in the system due to a large avalanche, it takes less time for another large avalanche to happen. In fact, large avalanches occur over shorter time intervals, and therefore the specific frequency of the system will increase. For larger networks, this growth is slower (see Figure 4). This frequency is independent of and and depends highly on and (see Figure 4). This result is biologically plausible since in cognitive tasks, the brain has more inputs and the external stimuli increase. Also, cognitive tasks are related to high-frequency activities (gamma rhythms) in the brain (Varela et al., 2001). In our suggested sandpile-based model of neuronal networks and the brain, increasing external stimuli change the frequency of activities and increases it.
The specific frequencies and the increment of frequencies have been shown previously for 1 in an ordinary square lattice (Wang et al., 2016). Here, we showed that it also happens for 1 in a more complex network. Moreover, our proposed model has more frequency content, which makes it more suitable to describe biological and cognitive facts.
Although we assume that each node is representative of a neuron, it is possible that each node represents an area in the brain. Actually, this study can be applied to different scales of the brain. It is worth mentioning that the sandpile model is a very simple model with equal amounts of for all nodes; in a real neuronal network, this value is dynamic and changes for each node over time. As future work, we suggest modeling the effect of heterogeneity of in the model.
We are grateful to Saman Moghimi-Araghi for sharing his valuable comments on this research.
A.S. and M.J. contributed equally to this work.