Abstract

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.

1  Introduction

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  Methods

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 .

Therefore, the height of the randomly selected node changes as
formula
2.1
When node exceeds its threshold, we will have
formula
2.2
formula
2.3

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.

In our model, the threshold ( for each node is equal to the degree of the node ( plus a leak term (:
formula
2.4
It is common to choose the boundary sites of the regular lattice as the leaky nodes (Wang et al., 2016). Therefore, defining the leak term in the small-world network can be done as follows. First, in an ordinary square lattice of size , the ratio of boundary nodes ( to total nodes is
formula
2.5
Therefore, in all the networks used in this study, nodes were selected as the boundary nodes. Second, the degree of boundary nodes in the square lattice is less than the other nodes. Thus, in this study, low-degree nodes should be chosen as leaky boundary nodes. Third, among boundary nodes, the corner nodes that have the smallest degree have the maximum value of leak term. Hence, the amount of leak term is assigned inversely proportional to the degree of each selected node. With this procedure, a node with the lowest degree will have the maximum leak. We set the maximum level of the leak to a fixed value. With a fixed maximum leak, the observed properties in different topologies become more comparable to other topologies in our work and also with previous studies in which the effect of small-worldness is studied (De Arcangelis & Herrmann, 2002; Chen et al., 2008; Hoore & Moghimi-Araghi, 2013).
In order to keep the maximum leak larger than (which is up to 4), the maximum leak is set to 5. Consequently, the leak term for every selected boundary node is
formula
2.6

2.3  Distribution Fitting

For evaluating the exponent of power law distributions, we use maximum likelihood estimator (MLE), an appropriate method for fitting distributions. The likelihood for a typical distribution is defined as
formula
2.7
where is the distribution and is the value of the th data. The , , … are the parameters of the distribution. The size and duration of avalanches are discrete values. Therefore, in a discrete power law distribution, the logarithm of the likelihood is
formula
2.8

The proper choice of leads to the maximum likelihood (Clauset, Shalizi, & Newman, 2009; Marshall et al., 2016).

3  Results

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.

Figure 1:

Distribution of avalanche size and duration (inset); with = 1600 and = 1 for (a), = 0.001, (b) = 0.1, and (c) = 0.5. PSS for (d) = 0.001, (e) = 0.1, and (f) = 0.5. The critical behavior and sharp peak could not be generated in a regular network (panels a and d) but emerged in small-world networks (panels b and e) for higher mean degrees. These properties exist for any value in the random networks (panels c and f).

Figure 1:

Distribution of avalanche size and duration (inset); with = 1600 and = 1 for (a), = 0.001, (b) = 0.1, and (c) = 0.5. PSS for (d) = 0.001, (e) = 0.1, and (f) = 0.5. The critical behavior and sharp peak could not be generated in a regular network (panels a and d) but emerged in small-world networks (panels b and e) for higher mean degrees. These properties exist for any value in the random networks (panels c and f).

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).

Figure 2:

Two power law regimes and two peaks in small-world networks. (a) Avalanche size distribution with two distinct linear parts with different slopes. The inset shows the exponents of two parts for different in the range of small-world networks (0.05 P 0.3). (b) The effect of rewiring probability on crossover of two regimes. The inset shows the power law behavior of in the range of small-world networks. (c) Two peaks in PSS. The upper inset shows a correlation between the sizes of successive avalanches and their random shuffled version. The lower inset shows the frequency of peak ( and corresponding frequency to correlation of large avalanches (.

Figure 2:

Two power law regimes and two peaks in small-world networks. (a) Avalanche size distribution with two distinct linear parts with different slopes. The inset shows the exponents of two parts for different in the range of small-world networks (0.05 P 0.3). (b) The effect of rewiring probability on crossover of two regimes. The inset shows the power law behavior of in the range of small-world networks. (c) Two peaks in PSS. The upper inset shows a correlation between the sizes of successive avalanches and their random shuffled version. The lower inset shows the frequency of peak ( and corresponding frequency to correlation of large avalanches (.

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 3:

(a) Distribution of avalanche size for = 1600, = 32, = 0.1 for different . The power law behavior is independent of . (b) Distribution of avalanche size with = 24, = 0.1, and = 1. The cutoff increases by network size. (c) Distribution of avalanche duration and WTs for parameters in panel a. Durations are independent of external stimuli, while WTs depend highly on Inset: Average of duration and WTs. (d) Distribution of duration and WTs for parameters in panel b. Distribution of WTs and power law behavior of durations are independent of . Inset: Average of duration and WTs. Average of durations increases by network size. (e) PSS for parameters in panel a. Inset: PS of reduced time series by removing WTs. Emergence of SOs in PSS does not result from WTs. (f) PSS for parameters in panel b. Inset: PS of reduced time series for different sizes. The larger networks have smaller peak frequencies.

Figure 3:

(a) Distribution of avalanche size for = 1600, = 32, = 0.1 for different . The power law behavior is independent of . (b) Distribution of avalanche size with = 24, = 0.1, and = 1. The cutoff increases by network size. (c) Distribution of avalanche duration and WTs for parameters in panel a. Durations are independent of external stimuli, while WTs depend highly on Inset: Average of duration and WTs. (d) Distribution of duration and WTs for parameters in panel b. Distribution of WTs and power law behavior of durations are independent of . Inset: Average of duration and WTs. Average of durations increases by network size. (e) PSS for parameters in panel a. Inset: PS of reduced time series by removing WTs. Emergence of SOs in PSS does not result from WTs. (f) PSS for parameters in panel b. Inset: PS of reduced time series for different sizes. The larger networks have smaller peak frequencies.

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.

Figure 4:

Peak frequency in PSS for = 16 and = 0.1, with different size and . The peak frequency grows by increasing external stimuli. For smaller networks, this growth is faster, but the peak location is independent of network structure and depends on only network size and external stimuli.

Figure 4:

Peak frequency in PSS for = 16 and = 0.1, with different size and . The peak frequency grows by increasing external stimuli. For smaller networks, this growth is faster, but the peak location is independent of network structure and depends on only network size and external stimuli.

4  Discussion

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.

Acknowledgments

We are grateful to Saman Moghimi-Araghi for sharing his valuable comments on this research.

References

Abbott
,
L. F.
(
1999
).
Lapicque's introduction of the integrate-and-fire model neuron (1907)
.
Brain Research Bulletin
,
50
(
5
),
303
304
.
Bak
,
P.
(
1996
).
How nature works: The science of self-organized criticality
.
Nature
,
383
(
6603
),
772
773
.
Bak
,
P.
,
Tang
,
C.
, &
Wiesenfeld
,
K.
(
1987
).
Self-organized criticality: An explanation of the 1/f noise
.
Physical Review Letters
,
59
(
4
),
381
.
Bak
,
P.
,
Tang
,
C.
, &
Wiesenfeld
,
K.
(
1988
).
Self-organized criticality
.
Physical Review A
,
38
(
1
),
364
.
Bassett
,
D. S.
, &
Bullmore
,
E. D.
(
2006
).
Small-world brain networks
.
Neuroscientist
,
12
(
6
),
512
523
.
Beggs
,
J. M.
, &
Plenz
,
D.
(
2003
).
Neuronal avalanches in neocortical circuits
.
Journal of Neuroscience
,
23
(
35
),
11167
11177
.
Beggs
,
J. M.
, &
Timme
,
N.
(
2012
).
Being critical of criticality in the brain
.
Frontiers in Physiology
,
3
,
163
.
Benayoun
,
M.
,
Cowan
,
J. D.
,
van Drongelen
,
W.
, &
Wallace
,
E.
(
2010
).
Avalanches in a stochastic model of spiking neurons
.
PLoS Comput. Biol.
,
6
(
7
),
e1000846
.
Bhaumik
,
H.
, &
Santra
,
S. B.
(
2013
).
Critical properties of a dissipative sandpile model on small-world networks
.
Physical Review E
,
88
(
6
),
062817
.
Bhaumik
,
H.
, &
Santra
,
S. B.
(
2016
).
Dissipative stochastic sandpile model on small-world networks: Properties of nondissipative and dissipative avalanches
.
Physical Review E
,
94
(
6
),
062138
.
Boffetta
,
G.
,
Carbone
,
V.
,
Giuliani
,
P.
,
Veltri
,
P.
, &
Vulpiani
,
A.
(
1999
).
Power laws in solar flares: Self-organized criticality or turbulence?
Physical Review Letters
,
83
(
22
),
4662
.
Bressler
,
S. L.
, &
Kelso
,
J. S.
(
2001
).
Cortical coordination dynamics and cognition
.
Trends in Cognitive Sciences
,
5
(
1
),
26
36
.
Buxhoeveden
,
D. P.
, &
Casanova
,
M. F.
(
2002
).
The minicolumn hypothesis in neuroscience
.
Brain
,
125
(
5
),
935
951
.
Buzsáki
,
G.
(
2006
).
Rhythms of the brain.
New York
:
Oxford University Press
.
Buzsáki
,
G.
, &
Watson
,
B. O.
(
2012
).
Brain rhythms and neural syntax: Implications for efficient coding of cognitive content and neuropsychiatric disease
.
Dialogues Clin. Neurosci.
,
14
(
4
),
345
367
.
Canolty
,
R. T.
,
Cadieu
,
C. F.
,
Koepsell
,
K.
,
Knight
,
R. T.
, &
Carmena
,
J. M.
(
2012
).
Multivariate phase–amplitude cross-frequency coupling in neurophysiological signals
.
IEEE Transactions on Biomedical Engineering
,
59
(
1
),
8
11
.
Chen
,
C. C.
,
Chiao
,
L. Y.
,
Lee
,
Y. T.
,
Cheng
,
H. W.
, &
Wu
,
Y. M.
(
2008
).
Long-range connective sandpile models and its implication to seismicity evolution
.
Tectonophysics
,
454
(
1
),
104
107
.
Clar
,
S.
,
Drossel
,
B.
, &
Schwabl
,
F.
(
1996
).
Forest fires and other examples of self-organized criticality
.
Journal of Physics: Condensed Matter
,
8
(
37
), 6803.
Clauset
,
A.
,
Shalizi
,
C. R.
, &
Newman
,
M. E.
(
2009
).
Power-law distributions in empirical data
.
SIAM Review
,
51
(
4
),
661
703
.
Csicsvari
,
J.
,
Jamieson
,
B.
,
Wise
,
K. D.
, &
Buzsáki
,
G.
(
2003
).
Mechanisms of gamma oscillations in the hippocampus of the behaving rat
.
Neuron
,
37
(
2
),
311
322
.
Dayan
,
P.
, &
Abbott
,
L. F.
(
2003
).
Theoretical neuroscience: Computational and mathematical modeling of neural systems
.
Journal of Cognitive Neuroscience
,
15
(
1
),
154
155
.
De Arcangelis
,
L.
, &
Herrmann
,
H. J.
(
2002
).
Self-organized criticality on small world networks
.
Physica A
: Statistical Mechanics and Its Applications,
308
(
1
),
545
549
.
Dehmer
,
M.
,
Mowshowitz
,
A.
, &
Emmert-Streib
,
F.
(Eds.). (
2013
).
Advances in network complexity
.
Hoboken, NJ
:
Wiley
.
Dhar
,
D.
(
1999
).
Studying self-organized criticality with exactly solved models
.
arXiv:cond-mat/9909009
.
Dickman
,
R.
,
Muñoz
,
M. A.
,
Vespignani
,
A.
, &
Zapperi
,
S.
(
2000
).
Paths to self-organized criticality
.
Brazilian Journal of Physics
,
30
(
1
),
27
41
.
Drossel
,
B.
, &
Schwabl
,
F.
(
1992
).
Self-organized criticality in a forest-fire model
.
Physica A
: Statistical Mechanics and Its Applications,
191
(
1–4
),
47
50
.
Friedman
,
N.
,
Ito
,
S.
,
Brinkman
,
B. A.
,
Shimono
,
M.
,
DeVille
,
R. L.
,
Dahmen
,
K. A.
, …
Butler
,
T. C.
(
2012
).
Universal critical dynamics in high resolution neuronal avalanche data
.
Physical Review Letters
,
108
(
20
),
208102
.
Gireesh
,
E. D.
, &
Plenz
,
D.
(
2008
).
Neuronal avalanches organize as nested theta-and beta/gamma-oscillations during development of cortical layer 2/3
.
Proceedings of the National Academy of Sciences
,
105
(
21
),
7576
7581
.
Goldenfeld
,
N.
(
1992
).
Lectures on phase transitions and the renormalization group
.
Boulder, CO
:
Westview Press
.
Harris
,
T. E.
(
2002
).
The theory of branching processes
.
Courier Corporation
.
Hilgetag
,
C. C.
, &
Kaiser
,
M.
(
2004
).
Clustered organization of cortical connectivity
.
Neuroinformatics
,
2
(
3
),
353
360
.
Hoore
,
M.
, &
Moghimi-Araghi
,
S.
(
2013
).
Critical behavior of a small-world sandpile model
.
Journal of Physics A
: Mathematical and Theoretical,
46
(
19
),
195001
.
Jensen
,
O.
, &
Colgin
,
L. L.
(
2007
).
Cross-frequency coupling between neuronal oscillations
.
Trends in Cognitive Sciences
,
11
(
7
),
267
269
.
Jirsa
,
V. K.
, &
Kelso
,
J. S.
(
2000
).
Spatiotemporal pattern formation in neural systems with heterogeneous connection topologies
.
Physical Review E
,
62
(
6
),
8462
.
Kadanoff
,
L. P.
,
Nagel
,
S. R.
,
Wu
,
L.
, &
Zhou
,
S. M.
(
1989
).
Scaling and universality in avalanches
.
Physical Review A
,
39
(
12
),
6524
.
Kirschfeld
,
K.
(
2005
).
The physical basis of alpha waves in the electroencephalogram and the origin of the “Berger effect.”
Biological Cybernetics
,
92
(
3
),
177
185
.
Koepp
,
M. J.
, &
Woermann
,
F. G.
(
2005
).
Imaging structure and function in refractory focal epilepsy
.
Lancet Neurology
,
4
(
1
),
42
53
.
Kramer
,
M. A.
,
Eden
,
U. T.
,
Kolaczyk
,
E. D.
,
Zepeda
,
R.
,
Eskandar
,
E. N.
, &
Cash
,
S. S.
(
2010
).
Coalescence and fragmentation of cortical networks during focal seizures
.
Journal of Neuroscience
,
30
(
30
),
10076
10085
.
Ktitarev
,
D. V.
,
Lübeck
,
S.
,
Grassberger
,
P.
, &
Priezzhev
,
V. B.
(
2000
).
Scaling of waves in the Bak-Tang-Wiesenfeld sandpile model
.
Physical Review E
,
61
(
1
),
81
.
Lisman
,
J.
(
2005
).
The theta/gamma discrete phase code occurring during the hippocampal phase precession may be a more general brain coding scheme
.
Hippocampus
,
15
(
7
),
913
922
.
Lübeck
,
S.
, &
Usadel
,
K. D.
(
1997
).
Numerical determination of the avalanche exponents of the Bak-Tang-Wiesenfeld model
.
Physical Review E
,
55
(
4
),
4095
.
Malamud
,
B. D.
,
Morein
,
G.
, &
Turcotte
,
D. L.
(
1998
).
Forest fires: An example of self-organized critical behavior
.
Science
,
281
(
5384
),
1840
1842
.
Marshall
,
N.
,
Timme
,
N. M.
,
Bennett
,
N.
,
Ripp
,
M.
,
Lautzenhiser
,
E.
, &
Beggs
,
J. M.
(
2016
).
Analysis of power laws, shape collapses, and neural complexity: New techniques and Matlab support via the NCC toolbox
.
Frontiers in Physiology
,
7
.
Meisel
,
C.
,
Storch
,
A.
,
Hallmeyer-Elgner
,
S.
,
Bullmore
,
E.
, &
Gross
,
T.
(
2012
).
Failure of adaptive self-organized criticality during epileptic seizure attacks
.
PLoS Comput. Biol.
,
8
(
1
), e1002312.
Mountcastle
,
V. B.
(
1997
).
The columnar organization of the neocortex
.
Brain
,
120
(
4
),
701
722
.
Nishimori
,
H.
, &
Ortiz
,
G.
(
2010
).
Elements of phase transitions and critical phenomena.
Oxford
:
Oxford University Press
.
Paczuski
,
M.
,
Maslov
,
S.
, &
Bak
,
P.
(
1996
).
Avalanche dynamics in evolution, growth, and depinning models
.
Physical Review E
,
53
(
1
),
414
.
Pinho
,
S. T. R.
,
Prado
,
C. P. C.
, &
Salinas
,
S. R.
(
1997
).
Complex behavior in one-dimensional sandpile models
.
Physical Review E
,
55
(
3
),
2159
.
Ponten
,
S. C.
,
Bartolomei
,
F.
, &
Stam
,
C. J.
(
2007
).
Small-world networks and epilepsy: graph theoretical analysis of intracerebrally recorded mesial temporal lobe seizures
.
Clinical Neurophysiology
,
118
(
4
),
918
927
.
Priesemann
,
V.
,
Wibral
,
M.
,
Valderrama
,
M.
,
Pröpper
,
R.
,
Le Van Quyen
,
M.
,
Geisel
,
T.
, …
Munk
,
M. H.
(
2014
).
Spike avalanches in vivo suggest a driven, slightly subcritical brain state
.
Frontiers in Systems Neuroscience
,
8
.
Pruessner
,
G.
(
2012
).
Self-organised criticality: Theory, models and characterisation
.
Cambridge
:
Cambridge University Press
.
Reijneveld
,
J. C.
,
Ponten
,
S. C.
,
Berendse
,
H. W.
, &
Stam
,
C. J.
(
2007
).
The application of graph theoretical analysis to complex networks in the brain
.
Clinical Neurophysiology
,
118
(
11
),
2317
2331
.
Rubinov
,
M.
, &
Sporns
,
O.
(
2010
).
Complex network measures of brain connectivity: Uses and interpretations
.
Neuroimage
,
52
(
3
),
1059
1069
.
Ruelle
,
P.
, &
Sen
,
S.
(
1992
).
Toppling distributions in one-dimensional Abelian sandpiles
.
Journal of Physics A: Mathematical and General
,
25
(
22
), L1257.
Schindler
,
K. A.
,
Bialonski
,
S.
,
Horstmann
,
M. T.
,
Elger
,
C. E.
, &
Lehnertz
,
K.
(
2008
).
Evolving functional network properties and synchronizability during human epileptic seizures
.
Chaos: An Interdisciplinary Journal of Nonlinear Science
,
18
(
3
),
033119
.
Schroeder
,
C. E.
, &
Lakatos
,
P.
(
2009
).
Low-frequency neuronal oscillations as instruments of sensory selection
.
Trends in Neurosciences
,
32
(
1
),
9
18
.
Schuster
,
H. G.
(
2014
).
Criticality in neural systems
.
D.
Plenz
, &
E.
Niebur
(Eds.).
Hoboken, NJ
:
Wiley
.
Sporns
,
O.
(
2004
). Complex neural dynamics. In
V. K.
Jirsa
&
J. A. S.
Kelso
(Eds.),
Coordination Dynamics: Issues and Trends
(pp.
197
215
).
Berlin
:
Springer
.
Sporns
,
O.
,
Chialvo
,
D. R.
,
Kaiser
,
M.
, &
Hilgetag
,
C. C.
(
2004
).
Organization, development and function of complex brain networks
.
Trends in Cognitive Sciences
,
8
(
9
),
418
425
.
Sporns
,
O.
, &
Tononi
,
G.
(
2001
).
Classes of network connectivity and dynamics
.
Complexity
,
7
(
1
),
28
38
.
Stam
,
C. J.
, &
Reijneveld
,
J. C.
(
2007
).
Graph theoretical analysis of complex networks in the brain
.
Nonlinear Biomedical Physics
,
1
(
1
),
3
.
Stam
,
C. J.
,
Jones
,
B. F.
,
Nolte
,
G.
,
Breakspear
,
M.
, &
Scheltens
,
P.
(
2007
).
Small-world networks and functional connectivity in Alzheimer's disease
.
Cerebral Cortex
,
17
(
1
),
92
99
.
Stanley
,
H. E.
(
1999
).
Scaling, universality, and renormalization: Three pillars of modern critical phenomena
.
Reviews of Modern Physics
,
71
(
2
),
S358
.
Steriade
,
M.
(
2001
).
Impact of network activities on neuronal properties in corticothalamic systems
.
Journal of Neurophysiology
,
86
(
1
),
1
39
.
Steriade
,
M.
(
2006
).
Grouping of brain rhythms in corticothalamic systems
.
Neuroscience
,
137
(
4
),
1087
1106
.
Tort
,
A. B.
,
Komorowski
,
R. W.
,
Manns
,
J. R.
,
Kopell
,
N. J.
, &
Eichenbaum
,
H.
(
2009
).
Theta–gamma coupling increases during the learning of item–context associations
.
Proceedings of the National Academy of Sciences
,
106
(
49
),
20942
20947
.
Varela
,
F.
,
Lachaux
,
J. P.
,
Rodriguez
,
E.
, &
Martinerie
,
J.
(
2001
).
The brainweb: Phase synchronization and large-scale integration
.
Nature Reviews Neuroscience
,
2
(
4
),
229
239
.
Wang
,
S. J.
,
Ouyang
,
G.
,
Guang
,
J.
,
Zhang
,
M.
,
Wong
,
K. M.
, &
Zhou
,
C.
(
2016
).
Stochastic oscillation in self-organized critical states of small systems: Sensitive resting state in neural systems
.
Physical Review Letters
,
116
(
1
),
018101
.
Watts
,
D. J.
, &
Strogatz
,
S. H.
(
1998
).
Collective dynamics of “small-world” networks
.
Nature
,
393
(
6684
),
440
442
.
Yeomans
,
J. M.
(
1992
).
Statistical mechanics of phase transitions
.
Gloucestershire
:
Clarendon Press
.
Zarepour
,
M.
,
Niry
,
M. D.
, &
Valizadeh
,
A.
(
2015
).
Functional scale-free networks in the two-dimensional Abelian sandpile model
.
Physical Review E
,
92
(
1
),
012822
.

Author notes

A.S. and M.J. contributed equally to this work.