Beta-Rhythm Oscillations and Synchronization Transition in Network Models of Izhikevich Neurons: Effect of Topology and Synaptic Type

Despite their significant functional roles, beta-band oscillations are least understood. Synchronization in neuronal networks have attracted much attention in recent years with the main focus on transition type. Whether one obtains explosive transition or a continuous transition is an important feature of the neuronal network which can depend on network structure as well as synaptic types. In this study we consider the effect of synaptic interaction (electrical and chemical) as well as structural connectivity on synchronization transition in network models of Izhikevich neurons which spike regularly with beta rhythms. We find a wide range of behavior including continuous transition, explosive transition, as well as lack of global order. The stronger electrical synapses are more conducive to synchronization and can even lead to explosive synchronization. The key network element which determines the order of transition is found to be the clustering coefficient and not the small world effect, or the existence of hubs in a network. These results are in contrast to previous results which use phase oscillator models such as the Kuramoto model. Furthermore, we show that the patterns of synchronization changes when one goes to the gamma band. We attribute such a change to the change in the refractory period of Izhikevich neurons which changes significantly with frequency.

1. Introduction

Synchronization is an important collective phenomenon that may emerge in locally interacting physical and biological oscillatory systems (Pikovsky et al., 2001; Barahona and Pecora, 2002; Motter et al., 2005; Arenas et al., 2008). Neural tissue of central nervous system can generate oscillatory activity in various scales from individual neuron firing to macroscopic oscillations in large neural ensembles (Engel et al., 2001; Varela et al., 2001; Buzsaki and Draguhn, 2004; Jensen and Lisman, 2005; Buzsaki, 2006). Macroscopic rhythmic activity which is observed in electroencephalography (EEG) recordings, is believed to occur due to emergence of synchronization in oscillations of constituent neurons. Synchronization of neural activity has a fundamental role in brain functions such as vision, memory, action, perception, information transfer, thought, and so on (Kandel et al., 2000; Steinmetz et al., 2000; Fell and Axmacher, 2001; Fries et al., 2001; Varela et al., 2001; Uhlhass et al., 2008; Jiang et al., 2012).

Neural oscillations have been documented to cover a broad spectrum of frequencies. These oscillations are observed widely in every level of central nervous system and are usually categorized into five frequency bands: delta 0.5−3.5 Hz, theta 4−7 Hz, alpha 8−12 Hz, beta 13−30 Hz, and gamma >30 Hz (Buzsaki, 2006; Rosanova et al., 2009). Beta rhythms are associated with normal wakeful consciousness states and appear when one is alert, attentive or when a person is engaged in problem solving or decision making. Beta waves are also associated with the activities of motor cortex (Baker, 2007; Engel and Fries, 2010).

Synchronization in neural population has been in the focus of intense experimental and theoretical research recently. See (Salinas and Sejnowski, 2001; Bennett and Zukin, 2004; Zhou et al., 2006; Mikkelsen et al., 2013; Sadeghi and Valizadeh, 2014; Matias et al., 2015; Esfahani et al., 2016; Pazó and Monthrió, 2016) for a few examples. Although beta-band activities have a significant role in brain functions, they have attracted less attention than other frequency bands (Engel and Fries, 2010). This is all the more important as many fundamental functions of the brain are associated with such oscillations. For example, synchronization transition is an important issue. From a theoretical point of view, synchronization in a neuronal network occurs as one increases synaptic strength. How this transition occurs is of fundamental importance. Generally, the transition can occur either as a continuous transition or a discontinuous (explosive) manner. If continuous, a small change can lead to small changes in systems response; however, if explosive, a small change can lead to dramatic changes in system's response. In addition to the type of synaptic interaction, the role of network topology is of key issue in determining the order of synchronization transition. In this paper, we intend to investigate the effect of network topology and synaptic type on synchronization phase transition in populations of spiking neurons with none-identical intrinsic frequencies in beta band. Specifically, we will focus on the order of the emerging phase transition for various network structures and different synaptic interactions. It is believed that normal brain activity requires it to be close to a phase boundary (a critical point) which consequently provide access to both synchronous and asynchronous oscillations with small change in the input (Beggs and Timme, 2012; Hesse and Gross, 2014; Moosavi et al., 2017). Hence, it is important to know whether the emerging synchronization transition is continuous or abrupt.

It is usual to evoke phase oscillators to characterize transition properties of neural oscillations. See (Shimokawa and Shinomoto, 2006; Cumin and Unsworth, 2007; Maistrenko et al., 2007; Kitzbichler et al., 2009; Gómez-Gardeñes et al., 2010; Botcharova et al., 2014; Timms and English, 2014) for some examples. While this choice offers many computational and analytic advantages, it suffers from some drawbacks. For example, it is not possible to consider a biologically realistic dynamical model as a phase oscillator, since many important features such as realistic synaptic interaction, are not easily implemented in phase oscillator models such as the Kuramoto model. Also, the spiking patterns of real neurons with wide range of frequencies are washed out in phase oscillator models. We therefore propose to study neuronal dynamics according to the Izhikevich model (Izhikevich, 2003) which is obtained by reducing some biological aspects of Hodgkin-Huxley (HH) neuron using bifurcation methods (Izhikevich, 2007). This model is computationally simpler than HH neuron, but is still biologically plausible.

To describe the functional form of synaptic interactions, we use two experimentally documented synaptic types: electrical synapses or gap junctions and chemical synapses (Kopell and Ermentrout, 2004; Simon et al., 2005; Sohal and Huguenard, 2005; Roth and van Rossum, 2009; Pérez et al., 2011; Kuo et al., 2016). These two types of interactions will appear as distinct expressions for synaptic currents to be added to the Izhikevich neurons. To describe the structure of synaptic interaction, we couple neurons via a network. It is well-known that network connectivity can have strong effects on patterns of collective behavior such as synchronization (Watts and Strogatz, 1998; Strogatz, 2001). It is believed that key elements such as small-world effect, clustering, and heterogeneity are of fundamental importance effecting the general collective behavior of a network. We therefore propose to study various network structures starting with a regular ring with high clustering and no randomness. We next consider small-world networks which provide a balance between high clustering and small-world effect. We also consider the more random structures such as Erdos-Renyi (homogeneous) and scale-free (heterogeneous) networks with low clustering but dense long-range synapses.

Our main results are as follows: (i) we find that electrical synapses are more conducive to synchronization than chemical synapses, leading to explosive synchronization in beta band in random networks. (ii) we find that the effect of clustering is far more important than small-world effect in determining the order of transition. (iii) we find that patterns of synchronization are distinctly different in beta band from the corresponding transitions in the high frequency gamma band.

2. Methods

To construct a neural circuit, we consider N Izhikevich neurons on an arbitrary network with a specific (symmetric) adjacency matrix A. The electrical activity of each neuron of this ensemble is described by a set of two ordinary nonlinear coupled differential equations (Izhikevich, 2003):

$\begin{array}{ll}\frac{d{v}_i}{dt}=0.04v_i^2+5v_i+140-u_i+I_i^{DC}+I_i^{syn}& (1)\end{array}$

$\begin{array}{ll}\frac{d{u}_i}{dt}=a(b{v}_i-u_i)& (2)\end{array}$

with the auxiliary after-spike reset:

$\begin{array}{ll}if{v}_i\ge 30,then{v}_i\to cand{u}_i\to u_i+d& (3)\end{array}$

for i = 1, 2, …, N. Here v_i is the membrane potential and u_i is the membrane recovery variable. When v_i reaches its apex (v_max = 30 mV), voltage and recovery variable are reset according to Equation (3). The term $(0.04v_i^2+5v_i+140)$ has been chosen so that v has mV unit and t has ms units (Izhikevich, 2003). In addition a, b, c, and d are four adjustable parameters in this model. Tuning these parameters, Izhikevich neuron is capable of reproducing about twenty different intrinsic firing patterns observed in real neurons (Izhikevich, 2006; Izhikevich and Edelman, 2008). In this paper we set a = 0.02, b = 0.2, c = −65, and d = 8, which corresponds to regular spiking pattern (Izhikevich, 2003, 2006).

The term $I_i^{DC}$ is an external current which determines intrinsic firing rate of uncoupled Izhikevich neurons. Regularly spiking Izhikevich neurons exhibits a Hopf bifurcation at I^DC = 3.78 (Kim and Lim, 2013). We choose values of $I_i^{DC}$ randomly from a Poisson distribution with mean value 10. Thus the intrinsic firing rates f_i lay in beta band and are different from one neuron to the other. The term $I_i^{syn}$ in Equation (1) denotes synaptic current received by post-synaptic neuron i. If the synapse is electrical, the synaptic current is (Roth and van Rossum, 2009; Pérez et al., 2011):

$\begin{array}{ll}{I}_i^{syn}=\frac{1}{D_i}{\displaystyle \sum _j}{g}_{ji}(v_j-v_i)& (4)\end{array}$

and if the synapse is chemical then (Roth and van Rossum, 2009; Pérez et al., 2011):

$\begin{array}{ll}{I}_i^{syn}=\frac{1}{D_i}{\displaystyle \sum _j}{g}_{ji}\frac{exp(-\frac{t-t_j}{{\tau }_s})-exp(-\frac{t-t_j}{{\tau }_f})}{{\tau }_s-{\tau }_f}(V_0-v_i)& (5)\end{array}$

where D_i is in-degree of node i, g_ji is the strength of synapse from pre-synaptic neuron j to post-synaptic neuron i. g_ji = gA_ji, where g is the electrical conductance of synapse and A_ji is the element of adjacency matrix of the underlying network (Gros, 2015). A_ji = 1 if nodes j and i are connected and A_ji = 0, otherwise. Also in Equation (5) τ_s = 1.7 and τ_f = 0.2 are the slow and fast synaptic decay constants (Roth and van Rossum, 2009), t_j is the instance of last spike of pre-synaptic neuron j and V₀ is the reversal potential of synapse which is equal to zero since we assumed that all synapses in our circuit are excitatory. We only consider networks which are composed of one giant cluster, and thus no isolated nodes or clusters exist.

Our main goal is to study the role of various network properties on synchronization patterns that may emerge. The leading network properties we consider are clustering coefficient (C) and average path length (L). Also, the existence of hubs in heterogeneous networks are thought to play an important role in synchronization. We consider a regular ring with high clustering and large average path length (large-world-effect), a slightly random network which preserves clustering but has small-world effect, as well as two random networks which exhibit small clustering but strong small-world effect: a homogeneous Poissonian network with no hub as well as a heterogeneous scale-free network with hubs. Details of the networks used and comparison with corresponding theoretical values are summarized in Table 1.

TABLE 1

Table 1. Theoretical values of clustering coefficient C and its corresponding measured values for the four network structures we have used in this study, first and second columns.

We integrate the dynamical equations using fourth-order Runge-Kutta method with a time step of 0.01 ms in order to obtain v_i(t). We typically evolve the entire system for a long time and make sure that the system has reached a stationary state. We then perform our measurements and calculations. We obtain the instants of firings of all neurons and then assign an instantaneous phase to each neuron between each pairs of successive spikes, as in Pikovsky et al. (1997):

$\begin{array}{ll}{\varphi }_i(t)=2\pi \frac{t-t_i^m}{t_i^{m+1}-t_i^m}& (6)\end{array}$

where $t_i^m$ is the instant of m^th spike of neuron i. We define a global instantaneous order parameter:

$\begin{array}{ll}S(t)=\frac{2}{N(N-1)}{\displaystyle \sum _{i\ne j}c}o{s}^2\left(\frac{{\varphi }_i(t)-{\varphi }_j(t)}{2}\right)& (7)\end{array}$

where the sum is over all pairs of neurons in the system whether they are connected or not. The global order parameter S is the long-time-average of S(t) in the stationary state of the system (S = 〈S(t)〉_t) and measures the collective phase synchronization in neuronal oscillations. S is bounded between 0.5 and 1. If neurons spike out-of-phase, then S ≃ 0.5, where they spike completely in-phase S ≃ 1 and for states with partial synchrony 0.5 < S < 1. Synchronization transition is displayed in S−g plots where transition is expected to occur at a given value of g. We note that we have also calculated the Kuramoto order parameter and have found identical results as the ones calculated using Equation (7). The relevant codes have been shared in public domain at figshare.com.

3. Results

In Figure 1A, we have plotted the gain function of regularly spiking Izhikevich neuron. Gain function of a neuron shows the dependence of firing rate on the external stimulating current (Gerstner et al., 2014). It is seen that the neuron shows type II excitability in this parameter regime and is capable of generating regular spikes with a broad range of intrinsic frequencies from theta-band to gamma-band. This range is more diverse than the possible range of firing rates of HH neuron (Lee et al., 1998). Also for an illustration, we have plotted the time dependence of the electrical and chemical synaptic currents which an exemplary neuron in a network receives from its neighbors at the beginning of stationary state in Figures 1B,C, respectively. Here g = 0.15 and neurons of the circuit are unsynchronized. We note that the pattern of electrical synaptic current is very different from that of the chemical synaptic current. For one thing, it is an order of magnitude stronger (15 vs. 0.5). Secondly, they are dispersed and act more as a pulse as opposed to fluctuating current due to chemical synapses. Thus it is expected that electrical synapses have more impact on emergence of synchronization in the system.

FIGURE 1

Figure 1. (A) Gain function of Izhikevich neuron which shows the dependence of firing rate of an uncoupled neuron to the external current. (B) Electrical [and (C) Chemical] synaptic current which an exemplary neuron in a network receives for g = 0.15. In this case the neurons of the network are unsynchronized. t = 0 indicates the beginning of stationary state.

Next we focus on synchronization transition in network models. We construct networks of size N = 1, 000 and coordination number (average connectivity) z = 50, unless otherwise stated. Also we set the values of $I_i^{DC}$ so that the intrinsic firing rates f_i are in beta-band and have mean-value f = 〈f_i〉 ≃ 22 Hz, unless otherwise stated. Synchronization diagrams for regular ring of Izhikevich neurons with electrical and chemical synaptic interactions are illustrated in Figures 2A,B, respectively. It is observed that the network with electrical synapses exhibits a continuous transition to phase synchronization, while no transition occurs in the network with chemical synapses. Investigation of raster plots of the system with electrical synapses (Figure 2C) reveals that when synaptic interaction is weak, neurons spike out-of-phase. Note that the mean firing rate is f ≃ 22 Hz, and therefore each neuron should fire about seven times in the 300 ms window that is illustrated here. Increasing g slightly, leads to two neural groups each of which contains neurons that spike partially coherently but the members of two groups spike anti-phase with respect to each other. See g = 0.15 in Figure 2C. When we increase g, the phases of a number of members in one group gradually match the phases of the members of the other group. Hence the order parameter of synchronization increases continuously from S = 0.5 to higher values and the neural network exhibits a continuous transition to phase synchronization. In case of ring with chemical synapses, no anti-phase groups form. Since neuronal interactions are local, increasing synaptic strength, leads to the formation of wave-like pattern in order of neuronal spikes (Figure 2D). Although increasing g enhances local coherence in neuronal oscillations and each neuron has a small phase lag with adjacent neurons in the circuit, there exists no global order in the network and S ≃ 0.5, see Figure 2D. We cannot increase g to arbitrary large values, as after a certain value (near g ≃ 0.6) neurons start bursting instead of regular spiking (Kim and Lim, 2013). We therefore conclude that spiking Izhikevich neurons with chemical synapses with purely local interactions lead to local order without any long-range order necessary for a phase transition. However, the effect of long-range interaction may change this picture.

FIGURE 2

Figure 2. Synchronization diagram of Izhikevich neurons on a regular ring for: (A) circuit with electrical synapses and (B) circuit with chemical synapses. (C) Raster plots of the system of panel (A) for four different values of g. (D) Raster plots of the system in (B) for four different values of g. The mean firing rate is f ≃ 22 Hz here, and t = 0 indicates the beginning of stationary state. t is measured in units of ms.

To examine this we consider transition to phase synchronization in Watts-Strogatz (WS) small-world networks (Watts and Strogatz, 1998). In Figure 3A, we have plotted the variation of clustering coefficient C and average distance L when we rewire the previous ring with different probabilities p and found that for p = 0.01 the resulting network has significant small-world effect and clustering coefficient, simultaneously. Figures 3B,C show the synchronization diagrams of Izhikevich neurons with electrical and chemical synapses in WS networks with p = 0.01. It is seen that the resultant synchronization diagrams are similar to S−g plots of regular ring except for a different (larger) transition point in circuits with electrical synapses. Investigation of the raster plots of WS neural networks (not shown) reveals that the underlying reason for the observed synchronization transitions is exactly the same as the reason explained for regular ring above. Since our regular ring and WS network have approximately the same value of C but distinctly different values of L, similarity of the transitions which they produce indicates that the clustering coefficient (and not the small-world effect) is the main topological factor that plays an important role in the resulting transitions.

FIGURE 3

Figure 3. (A) Dependence of clustering coefficient C and average path length L on rewiring probability p in WS networks with N = 1, 000 and z = 50. C and L are normalized with C₀ and L₀ which are the clustering coefficient and average path length of a regular ring (p = 0), respectively. (B,C) Synchronization diagram of Izhikevich neurons on WS networks with p = 0.01, for electrical and chemical synaptic interactions. The mean firing rate is f ≃ 22 Hz here.

In order to examine the role of clustering coefficient further, we investigate more random topologies, viz, Erdos-Renyi (ER) network (Erdos and Renyi, 1959) and scale-free (SF) network with small average path length and negligible clustering coefficient (Gros, 2015). Both these networks have random structures with ER being homogeneous and SF being heterogenous exhibiting hubs which are thought to be to play an important role in synchronization phenomena (Gomez-Gardenes et al., 2011). Also, both these networks exhibit small-world effect, while the existence of hubs in SF networks leads to a relatively smaller L for a fixed z and N. Dependence of order parameter S on coupling strength g for an ER network of Izhikevich neurons with electrical synapses is illustrated in Figure 4A. The network exhibits a first-order or explosive transition to phase synchronization, with a large hysteresis loop, as neurons spike with beta rhythms. Note that the transition is truly explosive as S jumps directly to its maximum value immediately at the transition point, indication full synchrony in the network. Explosive synchronization is a novel phenomenon and has attracted much attention recently. Different mechanisms have been reported for generation of such type of synchronization transition so far (Gomez-Gardenes et al., 2011; Leyva et al., 2012; Ji et al., 2013; Skardal and Arenas, 2014; Zhang et al., 2015), and the key role played by heterogeneity has been in focus in this regard.

FIGURE 4

Figure 4. Synchronization diagram of Izhikevich neurons on ER network: (A) S−g plot in forward and backward evolution of the system with electrical synapses. (B) S−g plot for the system with chemical synapses. (C) Raster plots for the system with electrical synapses in forward direction. (D) Raster plots of the system with chemical synapses. The mean firing rate is f ≃ 22 Hz and t = 0 indicates the beginning of stationary state. t is measured in units of ms.

Seeking the underlying reason of this explosive transition, we investigate raster plots of this neural circuit for different values of g. Four such raster plots for forward evolution of the system are shown in Figure 4C. We find that neurons spike out-of-order initially. As g is increased slightly, the neurons in the system are organized into two distinct groups in which members of each group spike almost coherently, as they oscillate anti-phase with the other group. Further increase of g regulates neuronal phases in each group but the phase lag between two groups remain robust. Therefore there exists no global phase coherence in the system and S = 0.5, see g = 0.33 in Figure 4C. There exists a transition point for which these two anti-phase groups abruptly join together leading to complete phase coherence. Hence the order parameter suddenly jumps from S = 0.5 to S = 1, see Figure 4A and g = 0.34 in Figure 4C.

When interaction among neurons is mediated via the softer chemical synapses (Figure 4B), anti-phases groups do not form in the neural network. Since the clustering coefficient of ER network is negligible (C = 0.054 here) and long-range interaction is significant, wave-like patterns in neuron spikes do not appear. See raster plots in Figure 4D. The gradual increase of g subsequently leads to gradual increase in global order in the system leading to a continuous transition at which global order appear in the system, see Figure 4D. Further increase of g leads to increase of S as more and more neurons align their phases. Therefore Izhikevich neurons with mean firing rate f ≃ 22 Hz produce continuous transition to phase synchronization when they interact via chemical synapses on an ER network.

Next, we ask whether heterogeneity in SF networks can change the picture obtained from ER networks above. Figure 5 displays synchronization diagrams of Izhikevich neurons on SF networks. Here we have generated uncorrelated SF networks (Catanzaro et al., 2005) with coordination number z = 20 and degree distribution function P(k)~k^−γ with γ = 3. Smaller z is necessary here in order to give real meaning to heterogeneity needed in our study for SF networks. Note, that despite using smaller z the network still displays significant small-world effect and small clustering coefficient, see Table 1. S−g plots of Izhikevich neurons with mean firing rate f ≃ 22 Hz on SF networks with electrical and chemical synapses are illustrated in Figures 5A,B, respectively. Interestingly, it is observed that the resulting synchronization diagrams are essentially exactly the same as the results obtained for ER network. We therefore conclude that the existence of hubs does not have a significant effect in the synchronization pattern in the parameter regime we have focused for Izhikevich neurons. Furthermore, the dramatic change in clustering coefficient of random networks (ER or SF) lead to decidedly different synchronization pattern when compared to clustered networks such as regular ring or WS network.

FIGURE 5

Figure 5. Synchronization diagram of Izhikevich neurons on SF networks with N = 1, 000, z = 20 and degree distribution function P(k)~k^−γ with γ = 3. (A) Phase transition diagram for the system with electrical synapses in forward and backward evolution of the system. (B) Phase transition diagram for the system with chemical synapses. The mean firing rate is f ≃ 22 Hz.

When synchronization transition is studied in a population of phase oscillators such as Kuramoto model, results are independent of the mean value of frequency distribution. Therefore we can switch to a rotating frame of reference where the mean value of frequencies is zero (Gros, 2015). In contrast we found that the resulting synchronization transitions which we obtain in neural circuits depend on the mean frequency of firing. For example, in Figure 6 we have illustrated synchronization diagrams and raster plots of Izhikevich neurons with firing frequencies in high gamma band (f = 〈f_i〉 ≃ 70 Hz) on ER network. Comparing these results with plots of Figure 4, it is observed that while neurons oscillate with high gamma rhythms, electrical synapses lead to a continuous transition to phase synchronization (rather than the explosive transition in beta frequencies) and chemical synapses do not lead to any synchronization in the system (as opposed to a continuous transition in beta frequencies). Also investigation of raster plots shows that while neurons fire with high gamma frequencies, interactions via electrical synapses do not results in anti-phase synchronization. (Compare Figure 4C and Figure 6C). This is a curious result that needs further investigation.

FIGURE 6

Figure 6. Synchronization diagram of Izhikevich neurons with spiking frequencies in the high gamma band on ER networks: Phase transition diagram for the system with (A) electrical synapses and (B) chemical synapses. (C) Raster plots for the system with electrical synapses. The mean firing rate is f ≃ 70 Hz and t = 0 indicates the beginning of stationary state.

In order to justify the frequency-dependent behavior of our neural networks, we illustrate spike trains of an individual Izhikevich neuron for two different values of firing frequencies f = 22 Hz and f = 70 Hz in Figure 7A. The horizontal solid green line at −55 mV indicates the threshold for firing, while the solid red line at −60 mV indicates the resting potential. One can see that the firing pattern of the two neurons are exactly the same, i.e., the dynamics above threshold are identical. However, the dynamics below the resting potential is decidedly different, as the lower frequency beta oscillation takes much longer to reach resting potential. Note that the hyperpolarization is stronger in the beta regime and the relative refractory period (time during which the system remains below resting potential) is clearly longer. This time scale τ_ref which renders the neuron to be relatively unexcitable is an important factor. In other words, while changing the frequency of Izhikevich neurons does not change the time scale of firing, it has a strong influence on the refractory period. This relative change of the time scales (firing vs. refractory) can provide an explanation for why anti-phase oscillations (and consequently explosive synchronization) occur in low frequency regime but not in the high frequency regime. In fact, existence of anti-phase oscillations have been attributed to separation of time scales in models of epidemic spreading (He and Stone, 2003). In Figure 7B, we plot the refractory period of Izhikevich neurons as function of frequency in the beta and gamma regime. One sees that in the Izhikevich neuron, the refractory period can become considerably long as one lowers the frequency of oscillations. While such a behavior may be an artifact of the model, one can see that many other neuronal dynamics models also exhibit similar behavior, i.e., a long time associated with slow increase in potential at low frequencies. Therefore, one might suspect that anti-phase oscillations and explosive synchronization might be associated with other generic neuronal models as well.

FIGURE 7

Figure 7. (A) Spike train of an Izhikevich neuron with spiking frequency 22 Hz (solid black curve) and with frequency 70 Hz (dashed blue curve). The solid green line indicate the threshold of firing and solid red line is the rest potential. (B) Dependence of refractory period τ_ref of Izhikevich neuron on spiking frequency. t = 0 indicates the beginning of stationary state.

4. Discussion

Synchronization transition in a network of oscillators has attracted much attention in recent years. The Kuramoto model has been used extensively in this regard with important implications for neural networks. However, it is a very crude approximation to consider neurons as phase oscillators. In this work we have studied synchronization transition in network models of biologically plausible neurons. We used Izhikevich neurons in beta frequency range coupled with electrical and chemical synapses on various network structures. We found that stronger electrical synapses are more conducive to synchronize than chemical synapses, regardless of network structure. We also found that electrical synapses can lead to anti-phase synchronization while no such behavior was seen for chemical synapses. As far as network structure was concerned, we found that the clustering coefficient, and not the small-world effect, is the key topological factor that determines the order of synchronization transition. When we introduced short-cuts into the regular ring, no significant change in the pattern of transition was observed. However, when random networks with small clustering coefficient were considered, synchronization patterns were significantly different. Additionally, heterogeneity in network structure did not play an important role as ER and SF results were identical. We note that anti-phase synchronization leading to explosive synchronization transition is a new and interesting mechanism which has not been reported in the existing literature to the best of our knowledge. The standard mechanisms reported in the literature for explosive synchronization are associated with heterogeneity and disorder. The reported results in this work becomes more interesting when we note that explosive synchronization occurred in the beta band and was not observed in the gamma band, i.e., it is frequency-dependent and therefore of dynamical origin, as opposed to the more widely-studied structural underpinnings. We note the fact that beta and gamma rhythms have different synchronization patterns has been reported before (Kopell et al., 2000). Such a frequency dependent behavior in synchronization patterns seem important and deserves further investigation (Khoshkhou and Montakhab, submitted).

Author Contributions

AM and MK conceived the project. MK carried out the computations. AM and MK analyzed the data. AM and MK wrote the paper.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Support from Shiraz University research council is kindly acknowledged. This work has been supported in part by a grant from the Cognitive Sciences and Technologies Council. We would also like to thank Alireza Valizadeh for many useful conversations throughout this work.

References

Arenas, A., Diaz-Guilera, A., Kurths, J., Morena, Y., and Zhou, C. (2008). Synchronization in complex networks. Phys. Rep. 469, 93–133. doi: 10.1016/j.physrep.2008.09.002

CrossRef Full Text | Google Scholar

Baker, S. N. (2007). Oscillatory interactions between sensorimotor cortex and the periphery. Curr. Opin. Neurobiol. 17, 649–655. doi: 10.1016/j.conb.2008.01.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Barahona, M., and Pecora, L. M. (2002). Synchronization in small-world systems. Phys. Rev. Lett. 89:054101. doi: 10.1103/PhysRevLett.89.054101

PubMed Abstract | CrossRef Full Text | Google Scholar

Beggs, J. M., and Timme, N. (2012). Being critical of criticality in the brain. Front. Physio. 3:163. doi: 10.3389/fphys.2012.00163

PubMed Abstract | CrossRef Full Text | Google Scholar

Bennett, M. V., and Zukin, R. S. (2004). Electrical coupling and neuronal synchronization in the mammalian brain. Neuron 41, 495–511. doi: 10.1016/S0896-6273(04)00043-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Botcharova, M., Farmer, S. F., and Berthouze, L. (2014). Markers of criticality in phase synchronization. Front. Syst. Neurosci. 8:00176. doi: 10.3389/fnsys.2014.00176

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsaki, G. (2006). Rhythems of the Brain. New York, NY: Oxford University Press.

Buzsaki, G., and Draguhn, A. (2004). Neuronal oscillations in cortical networks. Science, 304, 1926–1929. doi: 10.1126/science.1099745

PubMed Abstract | CrossRef Full Text | Google Scholar

Catanzaro, M., Boguñá, M., and Pastor-Satorras, R. (2005). Generation of uncorrelated random scale-free networks. Phys. Rev. E 71, 027103. doi: 10.1103/PhysRevE.71.027103

PubMed Abstract | CrossRef Full Text | Google Scholar

Cumin, D., and Unsworth, C. P. (2007). Generalising the Kuramoto model to the study of neuronal systems. Physica D 226, 181–196. doi: 10.1016/j.physd.2006.12.004

CrossRef Full Text | Google Scholar

Engel, A. K., and Fries, P. (2010). Beta-band oscillations signalling the status quo? Curr. Opin. Neurobiol. 20, 156–165. doi: 10.1016/j.conb.2010.02.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Engel, A. K., Fries, P., and Singer, W. (2001). Dynamic predictions: oscillations and synchrony in top-down processing. Nat. Rev. Neurosci. 2, 704–716. doi: 10.1038/35094565

PubMed Abstract | CrossRef Full Text | Google Scholar

Erdos, P., and Renyi, A. (1959). On random graphs. Publ. Math. 6, 290–297.

Google Scholar

Esfahani, Z. G., Gollo, L. L., and Valizadeh, A. (2016). Stimulus-dependent synchronization in delayed-coupled neuronal networks. Sci. Rep. 6:23471. doi: 10.1038/srep23471

PubMed Abstract | CrossRef Full Text | Google Scholar

Fell, J., and Axmacher, N. (2001). The role of phase synchronization in memory processes. Nat. Rev. Neurosci. 2, 105–118. doi: 10.1038/nrn2979

CrossRef Full Text | Google Scholar

Fries, P., Reynolds, J. H., Rorie, A. E., and Desimone, R. (2001). Modulation of oscillatory neuronal synchronization by selective visual attention. Science 291, 1560–1563. doi: 10.1126/science.1055465

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerstner, W., Kistler, W. M., Naud, R., and Paniuski, L. (2014). Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition. Cambridge: Cambridge University Press.

Google Scholar

Gomez-Gardenes, J., Gomez, S., Arenas, A., and Moreno, Y. (2011). Explosive synchronization transition in scale-free networks. Phys. Rev. Lett. 106:128701. doi: 10.1103/PhysRevLett.106.128701

CrossRef Full Text | Google Scholar

Gómez-Gardeñes, J., Zamora-Lopez, G., Moreno, Y., and Arenas, A. (2010). From modular to centralized organization of synchronization in functional areas of the cerebral cortex. PLoS ONE 5:e12313. doi: 10.1371/journal.pone.0012313

PubMed Abstract | CrossRef Full Text | Google Scholar

Gros, C. (2015). Complex and Adaptive Dynamical Systems. Berlin; Heidelberg: Springer-Verlag.

He, D., and Stone, L. (2003). Spatio-temporal syncronization of recurrent epidemics. Proc. R. Soc. Lond. B 270, 1519–1526. doi: 10.1098/rspb.2003.2366

PubMed Abstract | CrossRef Full Text | Google Scholar

Hesse, J., and Gross, T. (2014). Self-organized criticality as a foundamental property of nerural systems. Front. Sys. Neurosci. 8:166. doi: 10.3389/fnsys.2014.00166

PubMed Abstract | CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2003). Simple model of spiking neurons. IEEE Trans. Neural Netw. 14, 1569–1572. doi: 10.1109/TNN.2003.820440

PubMed Abstract | CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2006). Polychronization: computation with spikes. Neural Comput. 18, 245–282.

PubMed Abstract | Google Scholar

Izhikevich, E. M. (2007). Dynamical Systems in Neuroscience. Cambridge, MA: MIT Press.

Google Scholar

Izhikevich, E. M., and Edelman, G. M. (2008). Large-scale model of mammalian thalamocortical systems. Proc. Natl. Acad. Sci. U.S.A. 105, 3593–3598. doi: 10.1073/pnas.0712231105

PubMed Abstract | CrossRef Full Text | Google Scholar

Jensen, O., and Lisman, J. E. (2005). Hippocampal sequence-encoding driven by a cortical multi-item working memery buffer. Trends Neurosci. 28, 67–72. doi: 10.1016/j.tins.2004.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, P., Peron, T. K. D. M., Menck, P. J., Rodregues, F. A., and Kurths, J. (2013). Cluster explosive synchronization in complex networks. Phys. Rev. Lett. 110:218701. doi: 10.1103/PhysRevLett.110.218701

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, J., Dal, B., Peng, D., Zhu, C., Liu, L., and Lu, C. (2012). Neural synchronization durong face-to-face communication. J. Neurosci. 32, 16064–16069. doi: 10.1523/JNEUROSCI.2926-12.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Kandel, R. E., Schwartz, J. H., and Jessell, T. M. (2000). Principles of Neural Science. New York, NY: McGraw-Hill.

Kim, S. Y., and Lim, W. (2013). Coupling-induced population synchronization in an excitatory population of subthreshold Izhikevich neuron. Cogn. Neurodyn. 7, 495–503. doi: 10.1007/s11571-013-9256-y

CrossRef Full Text | Google Scholar

Kitzbichler, M. G., Smith, M. L., Christensen, S. R., and Bullmore, E. (2009). Broadband criticality of human brain network. PLoS Comput. Biol. 5:e1000314. doi: 10.1371/journal.pcbi.1000314

PubMed Abstract | CrossRef Full Text | Google Scholar

Kopell, N., and Ermentrout, B. (2004). Chemical and electrical synapses perform complementary roles in the synchronization of interneuronal networks. Proc. Natl. Acad. Sci. U.S.A. 101, 15482–15487. doi: 10.1073/pnas.0406343101

PubMed Abstract | CrossRef Full Text | Google Scholar

Kopell, N., Ermentrout, G. B., Whittington, M. A., and Traub, R. D. (2000). Gamma-rhythms and beta-rhythms have different synchronization properties. Proc. Natl. Acad. Sci. U.S.A. 97 1867–1872. doi: 10.1073/pnas.97.4.1867

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuo, S. P., Schwartz, G. W., and Rieke, F. (2016). Nonlinear spatiotemporal integration by electrical and chemical synapses in the retina. Neuron 90, 2. doi: 10.1016/j.neuron.2016.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S. G., Neiman, A., and Kim, S. (1998). Coherence resonance in a Hodgkin-Huxley neuron. Phys. Rev. E 57:3292. doi: 10.1103/PhysRevE.57.3292

CrossRef Full Text | Google Scholar

Leyva, I., Sevilla-Escoboza, R., Buldu, J. M., Sendima-Nadal, I., Gomez-Gardenes, J., Arenas, A., et al. (2012). Explosive first-order transition to synchronization in networked chaotic oscillators. Phys. Rev. Lett. 108:168702. doi: 10.1103/PhysRevLett.108.168702

CrossRef Full Text | Google Scholar

Maistrenko, Y. L., Lysyansky, B., Hauptmann, C., Burylko, O., and Tass, P. A. (2007). Multistability in the Kuramoto model with synaptic plasticity. Phys. Rev. E 75:066207. doi: 10.1103/PhysRevE.75.066207

PubMed Abstract | CrossRef Full Text | Google Scholar

Matias, F. S., Carelli, P. V., Mirasso, C. R., and Copelli, M. (2015). Self-organized Near-zero-lag synchronization induced by spike-timing dependent plasticity in cortical populations. PLoS ONE 10:e0140504. doi: 10.1371/journal.pone.0140504

PubMed Abstract | CrossRef Full Text | Google Scholar

Mikkelsen, K., Imparato, A., and Torcini, A. (2013). Emergence of slow collective oscillations in neural networks with spike-timing dependent plasticity. Phys. Rev. Lett. 110:208101. doi: 10.1103/PhysRevLett.110.208101

PubMed Abstract | CrossRef Full Text | Google Scholar

Moosavi, S. A., Montakhab, A., and Valizadeh, A. (2017). Refractory period in network models of excitable nodes: self-sustaining stable dynamics, extended scaling region and oscillatory behavior. Sci. Rep. 7: 7107. doi: 10.1038/s41598-017-07135-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Motter, A. E., Zhou, C. S., and Kurths, J. (2005). Enhancing complex-network synchronization. Europhys. Lett. 69, 334. doi: 10.1209/epl/i2004-10365-4

CrossRef Full Text | Google Scholar

Pazó, D., and Monthrió, E. (2016). From quasiperiodic patial synchronization to collective chaos in populations of inhibitory neurons with delay. Phys. Rev. Lett. 11006:238101. doi: 10.1103/PhysRevLett.116.238101

CrossRef Full Text

Pérez, T., Garcia, G. C., Eguiluz, V. M., Vicente, R., Pipa, G., and Mirasso, C. (2011). Effect of the topology and delayed interactions in neural networks synchronization. PLoS ONE 6:e19900. doi: 10.1371/journal.pone.0019900

CrossRef Full Text | Google Scholar

Pikovsky, A., Rosenblum, M., and Osipov, G. (1997). Phase synchronization of chaotic oscillators by external driving. Phys. D 104:219. doi: 10.1016/S0167-2789(96)00301-6

CrossRef Full Text | Google Scholar

Pikovsky, A., Rosenbum, M., and Kurths, J. (2001). Synchronization: A Universal Concept in Non-linear Sciences. New York, NY: Cambridge University Press.

Google Scholar

Rosanova, M., Casali, A., Bellina, V., Resta, F., Mariotti, M., and Massimini, M. (2009). Natural frequencies of human cortocothalamic circuits. J. Neurosci. 29, 7679–7685. doi: 10.1523/JNEUROSCI.0445-09.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Roth, A., and van Rossum, M. (2009). Computational Modeling Methods for Neuroscientists. Cambridge, MA: MIT Press.

Google Scholar

Sadeghi, S., and Valizadeh, A. (2014). Synchronization of delayed coupled neurons in presence of inhomogeneity. J. Comput. Neurosci. 36, 55–66. doi: 10.1007/s10827-013-0461-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Salinas, E., and Sejnowski, T. J. (2001). Correlated neuronal activity and the flow of neural information. Nat. Rev. Neurosci. 2:539. doi: 10.1038/35086012

PubMed Abstract | CrossRef Full Text | Google Scholar

Shimokawa, T., and Shinomoto, S. (2006). Inhibitory neurons can facilitate rhythmic activity in a neural network. Phys. Rev. E 73:066221. doi: 10.1103/PhysRevE.73.066221

PubMed Abstract | CrossRef Full Text | Google Scholar

Simon, A., Oláh, S., Molnár, G., Szabadics, J., and Tamás, G. (2005). Gap-junctional coupling between neurogliaform cells and various interneuron types in the neocortex. J. Neurosci. 25, 6278–6285. doi: 10.1523/JNEUROSCI.1431-05.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Skardal, P. S., and Arenas, A. (2014). Disorder induces explosive synchronization. Phys. Rev. E 89:062811. doi: 10.1103/PhysRevE.89.062811

PubMed Abstract | CrossRef Full Text | Google Scholar

Sohal, V. S., and Huguenard, J. R. (2005). Inhibitory coupling specifically generates emergent gamma oscillations in diverse cell types. roc. Natl. Acad. Sci. U.S.A. 102, 18638–18643. doi: 10.1073/pnas.0509291102

PubMed Abstract | CrossRef Full Text | Google Scholar

Steinmetz, P. N., Roy, A., Fitzgerald, P. J., Hsiao, S. S., Johnson, K. O., and Niebur, E. (2000). Attention modulates synchronized neuronal firing in primate somatosensory cortex. Nature 404, 187–190. doi: 10.1038/35004588

PubMed Abstract | CrossRef Full Text | Google Scholar

Strogatz, S. H. (2001). Exploring complex network. Nature 410, 268. doi: 10.1038/35065725

CrossRef Full Text

Timms, L., and English, L. Q. (2014). Synchronization in phase-coupled Kuramoto oscillator networks with time delay and synaptic plasticity. Phys. Rev. E 89:032906. doi: 10.1103/PhysRevE.89.032906

CrossRef Full Text | Google Scholar

Uhlhass, P. J., Pipa, G., Lima, B., Melloni, L., Neuenschwander, S., Nikolic, D., and Singer, W. (2008). Neural synchronization in cortical networks: history, concept and current status. Front. Integ. Neurosci. 3:17. doi: 10.3389/neuro.07.017.2009

CrossRef Full Text

Varela, F., Lachaux, J. P., Rodriguez, E., and Martinerie, J. (2001). The brainweb: phase synchronization and large-scale integration. Nat. Rev. Neurosci. 2, 229–239. doi: 10.1038/35067550

PubMed Abstract | CrossRef Full Text | Google Scholar

Watts, D. J., and Strogatz, S. H. (1998). Collective dynamics of small-world networks. Nature 393:440. doi: 10.1038/30918

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Boccaletti, S., Guan, S., and Liu, Z. (2015). Explosive synchronization in adaptive and multilayer networks. Phys. Rev. Lett. 114:038701. doi: 10.1103/PhysRevLett.114.038701

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, C., Zemanova, L., Hilgetag, C. C., and Kurths, J. (2006). Hierarchical organization unveiled by functional connectivity in complex brain networks. Phys. Rev. Lett. 97:238103. doi: 10.1103/PhysRevLett.97.238103

PubMed Abstract | CrossRef Full Text | Google Scholar