Explosive, continuous and frustrated synchronization transition in spiking Hodgkin-Huxley neuronal networks: the role of topology and synaptic interaction
EExplosive, continuous and frustrated synchronizationtransition in spiking Hodgkin-Huxley neuronalnetworks: the role of topology and synaptic interaction
Mahsa Khoshkhou, Afshin Montakhab Department of Physics, College of Sciences, Shiraz Universiy, Shiraz 71946-84795, Iran
Abstract
Synchronization is an important collective phenomenon in interacting oscillatoryagents. Many functional features of the brain are related to synchronization ofneurons. The type of synchronization transition that may occur (explosive vs.continuous) has been the focus of intense attention in recent years, mostly in thecontext of phase oscillator models for which collective behavior is independentof the mean-value of natural frequency. However, synchronization properties ofbiologically-motivated neural models depend on the firing frequencies. In thisstudy we report a systematic study of gamma-band synchronization in spik-ing Hodgkin-Huxley neurons which interact via electrical or chemical synapses.We use various network models in order to define the connectivity matrix. Wefind that the underlying mechanisms and types of synchronization transitionsin gamma-band differs from beta-band. In gamma-band, network regularitysuppresses transition while randomness promotes a continuous transition. Het-erogeneity in the underlying topology does not lead to any change in the orderof transition, however, correlation between number of synapses and frequency ofa neuron will lead to explosive synchronization in heterogenous networks withelectrical synapses. Furthermore, small-world networks modeling a fine balancebetween clustering and randomness (as in the cortex), lead to explosive syn-chronization with electrical synapses, but a smooth transition in the case ofchemical synapses. We also find that hierarchical modular networks, such asthe connectome, lead to frustrated transitions. We explain our results based on
Preprint submitted to Journal of L A TEX Templates January 23, 2020 a r X i v : . [ c ond - m a t . d i s - nn ] J a n arious properties of the network, paying particular attention to the competitionbetween clustering and long-range synapses. Keywords:
Synchronization, Hodgkin-Huxley neuron, phasetransition, electrical and chemical synapses, complex networks
1. Introduction
The phenomenon of phase transition is an important part of modern statisti-cal physics with many applications in physical and biological systems [1, 2, 3, 4].It is generally believed that many naturally occurring systems self-organize tothe edge of a phase transition point which can lead to many functional advan-tages [5, 6, 7]. In the case of biological systems, such transitions are generallybelieved to be of the critical (continuous) nature associated with a critical point[8, 9] or an extended critical regime[10, 11]. Synchronization transition is aninteresting phase transition that might occur in some important systems suchas power grids [12], ecological systems [13], and seasonal epidemics spreading[14], with both a continuous as well as a more interesting explosive transition.It is also interesting to note that the critical brain hypothesis [9] had originallyassumed that the brain operates at the edge of an activity phase transition.However, recent theoretical [15, 16] as well as experimental [17] studies showthat the criticality may be associated with a synchronization transition. Thispossibility provides a stronger motivation to study various types of synchroniza-tion transition that may occur in network models of biological neurons.Phase synchronization also has a key role in large-scale integration, memory,vision and other cognitive tasks performed by the human brain [18, 19, 20, 21].In a healthy brain, synchrony must occur at a moderate level. Excessive synchro-nization leads to brain disorders like epilepsy or Parkinson, while schizophreniaand autism are related to deficit of synchronization among neurons [22]. Thus ahealthy brain is thought to be functioning at the edge of synchronization transi-tion between order and randomness [23, 16, 15]. From this perspective, a slightincrease in neural interactions might lead to a synchronization transition in lo-2al neural circuits. The type of resultant transition (continuous, explosive orfrustrated) is therefore important. For example, when the emerging transitionis a continuous one, then a small change in the interaction strength changesthe amount of synchronization slightly. But if the emerging transition is anexplosive one a small increase in the interaction strength may result in a sud-den emergence of global order in the neural circuit. Explosive synchronizationhas functional advantages if it occurs during a fast response, but it also hasdisadvantages if it occurs, for example, during an epileptic seizure.Brain oscillations are categorized in various frequency bands. For example,beta-band (13-30 Hz) are typically associated with cognitive task performance.Recently, we have provided a systematic study of beta-band synchronizationtransitions in network models of Izhikevich neurons [24] and showed that con-trary to the case of simple phase oscillators, biologically meaningful models ofneural dynamics exhibit synchronization transitions which depend on the aver-age firing frequency of neurons [24]. This difference is rooted in the fact thatphase oscillator dynamics has a single time-scale (the mean-value of naturalfrequencies) which can be re-scaled without having any significant influence onthe dynamics of the network [25], while biologically plausible neural dynamicstypically has more than one time-scale, e.g. refractory period. The frequency-dependent behavior can arise when one of these time-scales depends on a chang-ing parameter, while the other one does not, thus leading to changing ratio ofthe various time-scales [24]. In fact it was shown that the patterns of tran-sition changed significantly when one increased the average frequency to thegamma-band ( >
30 Hz).Gamma-band oscillations are also an important class of rhythms appearingduring a broad range of the brain activities [26], and have received a great dealof attention. Gamma-band oscillations have been observed in several corticalareas, as well as subcortical structures [27]. In sensory cortex, gamma powerincreases with sensory drive [28], cognitive tasks including feature binding[29],visual grouping [30], stimulus selection [31, 32] and attention [33]. In highercortex, gamma power is the dominant rhythm during working memory [34]3nd learning [35]. Also, it is reported that irregular gamma waves have beenobserved in pathologies such as Alzheimer [36].Our purpose here is to provide a systematic study of a biologically motivatedneuronal network. We therefore propose to study synchronization transition ingamma-band and seek the effect of synaptic interaction (chemical vs. electricalsynapses) as well as the topology of the network used on the ensuing transitiontype. However, Izhikevich neurons have a tendency to burst as opposed to spikewhen one increases the input in order to increase the frequency. Furthermore,increasing interaction strength in network of Izhikevich neurons also leads tobursting behavior. On the other hand Hodgkin-Huxley (HH) neurons have alarge stable spiking range in gamma frequencies [37]. We therefore use networkmodels of HH neurons in gamma band in order to study synchronization patternswhich emerge.Although synchronization of HH (or HH-type) neurons has been extensivelystudied before, e.g. in [38, 39, 40, 41, 42, 43, 44, 45, 46], a systematic studyof (the order of) synchronization transition has not been performed to the bestof our knowledge. In fact, much of such type of studies usually employ phaseoscillator models such as the Kuramoto model [47, 48]. Here, our emphasis isto ascertain the type of phase transition (e.g. continuous vs. explosive) thatmay occur in a collection of HH neurons and how that may depend on synapticinteraction and/or underlying structure (network) [49]. Surprisingly, we findthat one- and two-dimensional lattice networks of spiking HH neurons exhibitno transition. Instead they exhibit quasiperiodic partial synchronization as aresult of strong clustering which does not lead to global order due to lack oflong-range interactions. Random network structures like Erdos-Renyi (ER) andscale-free (SF) networks exhibit continuous transition with either electrical orchemical synapses, with no significant difference between SF and ER structures.However, small-world network with high clustering coefficient and long-rangeinteraction exhibits explosive (first-order) transition to synchronization whenneurons interact via electrical synapses, but exhibits continuous (second-order)transition when interacting via chemical synapses. Furthermore, we consider4he role of heterogeneity by introducing a correlation between frequency andthe degree of a given neuron. We find that while heterogeneity (in degree orfrequency) does not change the order of continuous transition, a correlationbetween the two can lead to explosive synchronization with electrical synapses,but not with chemical synapses. Finally, we show that hierarchical modular(HM) networks with both types of synapses exhibit frustrated synchronizationin an intermediate regime between disordered and ordered phases of the system.Some of the structures studied here have been studied in the beta-band andwill consequently be compared and contrasted. However, the case of correlatedheterogeneity as well as HM networks are just included in the current study andtheir counterparts in beta-band had not been studied in ref.[24]. Consequently,such results can be compared with those of phase oscillators independent offrequency.In the following section, we describe the model we use for our study. InSection (3), we describe our simulation details including the numerical methodsused. Extensive results of our numerical study are presented in Section (4), andwe close the paper with some concluding remarks in Section (5).
2. Model
Consider N Hodgkin-Huxley neurons on an arbitrary network. Electricalactivity of i th neuron of the network is described by a set of four nonlinearcoupled ordinary differential equations as follows[37]: C m dv i dt = I DCi + I syni − G Na m i h i ( v i − V Na ) − G K n i ( v i − V K ) − G L ( v i − V L ) (1) dm i dt = α m ( v i )(1 − m i ) − β m ( v i ) m i (2) dh i dt = α h ( v i )(1 − h i ) − β h ( v i ) h i (3) dn i dt = α n ( v i )(1 − n i ) − β n ( v i ) n i (4)5or i = 1 , , ..., N . Here v i is the membrane potential, m i and h i are variablesfor activation and inactivation of sodium current, and n i is the variable foractivation of potassium current [37]. α and β functions are the so-called ratevariables of HH neuron for each type of ionic currents and depend on the in-stantaneous membrane potential [37]. We use C m = 1 . G Na = 120, G K = 36, G L = 0 . V Na = 50, V K = −
77 and V L = − .
387 for the constant parameters[38]. I DCi is an external current which differs from a neuron to the other anddetermines dynamical properties of uncoupled HH neurons. It is shown that for I DC > . µA/cm a stable limit-cycle is the global attractor for a single HHneuron [50]. We choose values of I DCi randomly from a Poisson distributionwith mean value 10 . µA/cm . Therefore, intrinsic firing rates are non-identicaland most of the neurons spike regularly with gamma rhythms [26]. Here, we setthe mean intrinsic firing rate is f (cid:39)
75 Hz, unless otherwise stated.The term I syni in Eq.(1) represents synaptic current received by post-synapticneuron i . Functional form of this current depends on the synaptic type. For agap junction or an electrical synapse the synaptic current is [51]: I syni = 1 D i (cid:88) j g ji ( v j − v i ) (5)and if the synapse is chemical then [51]: I syni = 1 D i (cid:88) j g ji exp ( − t − t j τ s ) − exp ( − t − t j τ f ) τ s − τ f ( V − v i ) (6)where D i is in-degree of node i , g ji is the strength of synapse from pre-synapticneuron j to post-synaptic neuron i . Here we assumed that all existing synapseshave the same strength, viz g ji = ga ji , where g is the electrical conductance ofsynapse and a ji is the element of adjacency matrix of the underlying network.Also in Eq.(6) t j is the instance of last spike of pre-synaptic neuron j , τ s and τ f are the slow and fast synaptic decay constants and V is the reversal potentialof synapse which is equal to zero since we assumed that all synapses in our6ircuit are excitatory. In this study we take τ s = 1 . τ f = 0 . t j . Therefore, for example, one wouldexpect for a given value of coupling strength g , electrical synapses would providemore synchronization when compared to chemical synapses.In order to quantify the amount of phase synchronization in a neural popu-lation, we assign an instantaneous phase to each neuron as in [52]: φ i ( t ) = 2 π t − t mi t m +1 i − t mi (7)where t mi is the instant of m th spike of neuron i . Then we define a globalinstantaneous order parameter as: S ( t ) = 2 N ( N − (cid:88) i (cid:54) = j cos (cid:16) φ i ( t ) − φ j ( t )2 (cid:17) (8)The global order parameter S is the long-time-average of S ( t ) at the stationarystate and measures collective phase synchronization in oscillations of membranepotentials of all neurons, viz S = (cid:104) S ( t ) (cid:105) t . S is bounded between 0.5 and 1. Ifneurons spike out-of-phase, then S (cid:39) . S (cid:39)
1. For states with partial synchrony 0 . < S < S , we have also calculated the more com-monly used Kuramoto order parameter [47]: R ( t ) e iθ = 1 N (cid:88) j e iφ j ( t ) (9)with R = (cid:104) R ( t ) (cid:105) t where 0 ≤ R ≤ R = 0 indicates asynchronous, while R = 17ompletely synchronous, oscillations. Essentially, the same results are obtainedfor R as those obtained for S . However, from a statistical point of view R ( t )represents an average of N data points while S ( t ) represents an average of N ( N − / κ S = (cid:16) (cid:104) S ( t ) (cid:105) − (cid:104) S ( t ) (cid:105) (cid:104) S ( t ) (cid:105) (cid:17) / (10)or: κ R = (cid:16) (cid:104) R ( t ) (cid:105) − (cid:104) R ( t ) (cid:105) (cid:104) R ( t ) (cid:105) (cid:17) / (11)Such generalized susceptibilities are very useful tools in order to study phasetransitions in general, since critical systems are supposed to exhibit maximalfluctuations at the critical point, diverging in the thermodynamic N → ∞ limit.
3. Methods
We have scrutinized transition to phase synchronization in networks with N HH neurons interacting via two different synaptic types. We start by pro-viding a detailed description of our procedure. We first determine the networktopology by specifying elements of its adjacency matrix. These elements areeither zero or one depending on if the nodes are unconnected or connected,respectively. The links in our networks are symmetric. Synapses are also notplastic in this study. The strength of synapses is set with parameter g explainedin the previous section. After constructing each network, the synaptic type isdetermined. If synapses are supposed to be electrical, we use Eq.(5) to describesynaptic currents. While neurons are assumed to interact via chemical synapses,Eq.(6) is used. Next, we fix the values of I DCi and set the parameter g equalto zero. We then integrate Eqs.(1)-(4) using fourth order Runge-Kutta method8ith a fixed time step ∆ t = 10 − ms. Typically, much larger time steps areused in simulations of HH neurons, see for example [40, 41]. However, since longrelaxation times were required in our studies (particularly near the transitionpoints) we choose such a short time step in order to avoid the accumulationof errors. Using this small time step, we are able to specify t mi , the instant of m th spiking of each neuron with an accuracy of 10 − ms. Finally we obtain thephase of all neurons and calculate S ( t ) and R ( t ) at every time instant (Eqs.(8)and (9)). We allow the dynamics to progress for a long transient time (order of10 time steps) until the fluctuations in S ( t ) or R ( t ) reach a stationary state.After reaching stationary state, we run our simulation for another 2 × ms(2 × time steps) and evaluate the order parameter S and R by averaging S ( t ) and R ( t ) over this second interval. We next increase the value of g slightly(keeping all other conditions fixed) and repeat the whole process to evaluate S and R again. In this manner we obtain dependence of order parameter on cou-pling strength g , in each network topology and for each of the above mentionedsynaptic types. The initial condition of integration are random for g = 0 andthe system is evolved quasi-statically for larger g values. The synchronizationdiagrams that are reported here are results of averaging over five network re-alizations as well as other stochastic parameters. Our results are reported fortypically N ≈
4. Results
The first structure that we consider is regular network, a one dimensionalring of size N = 500 and z = 50 as well as a two-dimensional lattice of side L ( N = L × L ) with L = 22 and z = 16. z is the coordination number ofthe network. The results are shown in Figs.1(a-d). It is observed that increas-ing g does not lead to a transition in either case. It is somewhat surprising9 igure 1: Synchronization diagram of HH neurons on regular networks: (a) and(b) One-dimensional ring with electrical and chemical synapses. (c) and (d) Two-dimensional lattice with electrical and chemical synapses. (e)-(h) Raster plots of theone-dimensional ring with electrical synapses for four values of g . (i) Network activityfor the system in fully asynchronized (orange curve) and in quasiperiodic partiallysynchronized (green curve) states. (j) The Kuramoto order parameter R vs g for 1Drings with electrical synapses for different system sizes N . (k) and (l) The generalizedsusceptibilities κ R and κ S vs g for the same systems as in (j). z = 0 . N in each case. t = 0 indicates the beginning of stationary state. The synchronization diagrams andsusceptibility plots show the averaged results over five initial conditions.
10s one would expect a transition to synchrony for large g . We have thereforeinvestigated the raster plots for this system for different g values. Such rasterplots for one-dimensional ring with electrical synapses for four values of g areshown in Figs.1(e-h). Raster plots for rings with chemical synapses are qualita-tively the same as Figs.1(e-h). It is realized that imposing a small interactionamong neurons in a regular ring leads to formation of correlated regions. Thisis not unexpected since regular rings have high clustering coefficient [25]. In-creasing g slightly, regulates the phase of neurons on a local level. Since thereare no long-range synapses in the system, further increase of g could not van-ish phase lags among neurons belonging to far away areas of the network, butinstead results in the emergence of the so-called quasiperiodic partial synchro-nization. Quasiperiodic partial synchronization is denoted to the state of apopulation of interacting oscillators in which the system sets into a nontrivialdynamical regime where oscillators display quasiperiodic dynamics while col-lective observable of the system oscillate periodically [53, 54, 55]. A relevantcollective observable for a neural network is the network activity that is definedas A ( t ) = N (cid:80) Ni =1 v i ( t ). In Fig.1(i) we have plotted A(t) for a regular ring ofHH neurons for g = 0 .
00 when neurons are fully asynchronous (orange curve)and also for g = 0 .
80 which is where the network is in a quasiperiodic partialsynchronization state (green curve). It is seen that when neurons spike out oforder, A ( t ) fluctuates irregularly. But when the network is in state of quasiperi-odic partial synchronization neurons spike quasiperiodically and A ( t ) oscillatesperiodically. This state emerges in the rings from g (cid:39) .
20 and remains robustwhen g is increased further as the perspective of raster plots remain qualita-tively the same from g = 0 .
20 to g = 1 .
00 (or even for larger values of g whichare not shown here).For sake of comparison, in Fig.1(j) we show the R − g plot for 1D ringwith electrical synapses for three different system sizes N . The coordinationnumber in each system is set to be z = 0 . N . In light of these plots we findthat the synchronization diagram remains unaltered upon increasing the systemsize. Also, comparing Figs. 1(a) and 1(j) one verifies the equivalence of the11esults obtained based on the order parameters R and S , except for the morerefined statistics resulting from S . Moreover, the generalized susceptibilities κ R and κ S for the same systems as in Fig.1(j), are illustrated in Figs. 1(k)and 1(l), respectively. It is observed that increasing g does not lead to anydistinctive peak in κ R or κ s , confirming that no phase transition occurs in thissystems. Generalized susceptibilities for other regular networks studied here arequalitatively similar to Figs.1(j) and 1(l) (not shown).We note that one might suspect that the lack of transition observed in the1D lattice might be due to the low dimensional structure, similar to the lackof phase transition in, for example, 1D Ising model. That is why we have alsoperformed simulations for the 2D L × L lattice whose main results are shown inFig.1(c) and 1(d) which again show no transition. We note that raster plots ofthe 2D system are quantitatively the same as the 1D case with lesser coherence(due to smaller clustering) and that network oscillations, A ( t ), are also verysimilar to the 1D case (not shown). We next consider random networks with small-world effect but with muchsmaller clustering compared to regular networks. We consider ER network whichhas a homogeneous random structure as well as SF network which has a het-erogeneous random structure [25]. Such networks are constructed using a con-figurational model [25]. The networks size is N = 500. z = 50 for ER networkand z = 25 for SF network. Also the degree distribution function of SF net-work is P ( k ) ∼ k − γ with γ = 2 .
2. The results for synchronization transition ofHH neurons with electrical and chemical synaptic currents, on ER and SF areshown in Figs.2(a-d). It is observed that the system with both types of synapsesexhibits a continuous transition from asynchrony to synchrony. Raster plots ofspikes for the ER network with electrical synapses are illustrated in Figs.2(e-h). It is evident that since clustering coefficient is significantly reduced due torandomness (as compared to regular networks), neuronal clusters do not appearin the system. However, presence of a significant number of long-range con-12 igure 2: (a) and (b)Transition to phase synchronization in ER networks of spikingHH neurons with electrical and chemical synapses. (c) and (d) Transition to phasesynchronization in SF networks of spiking HH neurons with electrical and chemicalsynapses. (e)-(h) Raster plots for ER network with electrical synapses for variousvalues of g . (i) Kuramoto order parameter R vs g for ER networks of HH neuronswith electrical synapses for different system sizes N . (j) and (k) κ R and κ S vs g forthe same systems as in (i). z = 0 . N in each case. t = 0 indicate the beginningof stationary state. The synchronization diagrams and susceptibility plots show theaveraged results over five network realizations and initial conditions. g is increasedabove a certain threshold. Raster of spikes for other transitions are qualita-tively similar to those of Figs.2(e-h) (not shown). Looking at the value wherethe transition occurs, g t , for a given network, one concludes that synchroniza-tion is more conducive to electrical synapses than chemical synapses, i.e. g t isabout an order of magnitude smaller for electrical synapses. This makes senseas electrical synapses are known to be stronger than chemical synapses. On theother hand, the strong similarity between the results for ER and SF networks fora given synaptic type, including their corresponding value at transition, leadsone to conclude that the role of structural heterogeneity (SF network) is not animportant factor in influencing the type and shape of transition curves ( S − g plots).In Fig.2(i), we also show R − g plots of HH neurons with electrical synapseson ER networks with three system sizes N to be compared with with the S − g plot in Fig.2(a). Here, z = 0 . N in each system size. Furthermore, variationsof the generalized susceptibilities κ R and κ S upon increasing g for the systemsof Fig.2(i) are plotted in Figs.2(j) and 2(k), respectively. It is observed thatboth κ R and κ S show a specific peak at the transition point which grows withincreasing the system size. The behavior of the generalized susceptibilities fur-ther collaborates our order parameter results which indicate that our modeldoes not show synchronization transition for low dimensional systems (Fig.1),but exhibits definitive and continuous transition in a high dimensional structuresuch as complex networks (Fig.2).Regarding the results associated with figures 1 and 2, we can conclude twoimportant points: (i) the main results are unaltered upon increasing the systemsize N , and (ii) the synchronization diagrams exhibit qualitatively the same be-havior whether we employ R or S , except for the more refined statistics providedby S which enables us to determine the transition point clearly. Therefore, forthe rest of this paper we report the results only based on the order parameter S and for our largest available system size ( N (cid:39) .3. Small-world networks After considering regular networks with high clustering but large averagedistance on one hand, and highly random networks with strong small-worldeffect but negligible clustering on the other hand, we are interested in networksthat have high clustering coefficient, as well as small-world property. Therefore,we constructed Watts-Strogatz (WS) networks [25] with N = 500 and z =50 by random rewiring of two percent of links of a regular ring. This lowrewiring probability ( p = 0 .
02) allows the system to keep its large clusteringcoefficient while developing significantly low average distance (i.e. small-worldeffect). The resulting S − g curves for WS networks with electrical and chemicalsynapses are shown in Fig.3(a) and 3(b), respectively. Interestingly, we observea discontinuous (explosive) transition for the case of electrical synapses while acontinuous transition is observed for the chemical synapses. As we will discussshortly, this type of explosive synchronization is different in its mechanism thanthose seen for phase oscillators in heterogeneous networks such as [56, 57]. Thesynchronization transition is accompanied with a hysteresis loop if a backwardsweep in g is performed from the highly synchronized state. Therefore, as seenin Fig.3(a), not only the transition is explosive, the value of S is also history-dependent. This is to be contrasted with the case of chemical synapses in WSnetwork where increasing g leads to a continuous transition from asynchrony tosynchrony in neural spiking as is clear in Fig.3(b).Therefore, one- and two-dimensional regular networks produced no transi-tion, while random networks produced a continuous transition. However, small-world networks which lie somewhere between randomness and regularity exhibitexplosive synchronization (electrical) as well as continuous transition (chemical).The fact that transition type in WS network depends on the interaction typeis an interesting result and may be important from the point of view of neuro-science, since it has been reported that the brain networks at the microscopiclevel are similar to WS networks [58]. To elucidate the effect of topology andunderlying reason for different order of phase transitions, we display the rasterplots of HH neurons with electrical synapses on a WS network in Fig.3(c-f)15 igure 3: Transition to phase synchronization in WS networks of spiking HH neuronswith (a) electrical and (b) chemical synapses. Network size and coordination numberare N = 500 and z = 50, respectively and rewiring probability is p = 0 .
02. (c)-(f) Raster plots for the network with electrical synapses for various values of g ofthe system in forward direction. (g) and (h) Two raster plots for the same value of g inside the hysteresis loop ( g = 0 . t = 0 indicate the beginning of stationary state. The synchronization diagrams showthe averaged results over five network realizations and initial conditions. S = 0 . g will eventually leadto interactions among various parts of the network which eventually leads toglobal order in the system and therefore a phase transition, which was absentin models without long-range interaction. But, why do we observe an explosivesynchronization for electrical synapses but a continuous transition for chemicalsynapses? This has to do with the fact that long-range links provide stronginteractions for the case of electrical synapses as phase difference of farawayregions is considerable, while providing weak interaction for local interactionswhich are mostly synchronized. When nonlocal regions suddenly go in synchdue to strong electrical interactions a sudden jump in order parameter is ob-served. This is shown in Figs.3(d) and 3(e) for the two consecutive values of g ( g = 0 .
065 and g = 0 . g , one for the forwardbranch and one for the backward branch. Here, we note how small changes inglobal patterns of spikes can lead to significant change in the value of S . It might have been expected that heterogeneity in network structure as inthe SF network in Fig.2(b) would have led to a different transition pattern when17ompared to homogenous random network of ER. We note that the importanceof the role of heterogeneity in neural networks has attracted much attentionin recent literature. An important observation in regards to synchronizationin the Kuramoto model was that structural heterogeneity was not sufficient tolead to different transitional pattern, but a correlation between the frequencyand the degree of the node was the key element that would lead to explosivesynchronization in SF networks [56]. This means that the high frequency nodesin a network are also the highly connected nodes, while the low frequency nodesare sparsely connected. We note that the range of frequency in spiking HHneurons is relatively limited. However, one may attempt to make a heteroge-neous distribution even in this limited range. We have therefore studied a SFnetwork of size N = 500 with γ = 2 . k min = 7, k max = 47, z = 15. Wehave also produced the same distribution of frequencies with f min = 70 Hz and f max = 110 Hz and have studied the correlated ( f ∝ k ) and the uncorrelateddistribution of such frequencies. The results are shown in Fig.4. As is seen,the correlated case leads to explosive synchronization along with hysteresis inthe case of electrical synapses, but a smooth transition in the case of chemi-cal synapses. We also show the same results for the uncorrelated case whichindicates that the explosive synchronization is in fact due to the correlation be-tween the two heterogeneous distributions of degree and frequency, in the caseof electrical synapses.It is interesting to note that the mechanism for explosive synchronization isvery different here than that observed in WS network for electrical synapses.There, it was the combined effect of local order, which is achieved for low synap-tic weight g , and long-range order which sets in for large values of synapticweight, that leads to sudden order and explosive synchronization in the system,see raster plots in Fig.3. Here, in the case of correlated heterogeneity, the sys-tem is still essentially in a completely disordered phase just before the explosivesynchronization occurs. See g = 0 .
227 and g = 0 .
228 raster plots in Fig.4 whichare just before and after the explosive synchronization transition point. Thisindicates that the system truly goes through a sudden change from disordered18o ordered phase. The cause of such an explosive synchronization can easily beunderstood by looking at the ordered raster plots. One sees that in the syn-chronous phase the entire system is oscillating at the frequency of f ≈
110 Hzwhich is exactly the frequency of the only hub in the system. This indicates theessential role of the hub in this explosive synchronization. The entire networkmust adjust with the hub and once this happens an explosive synchronizationoccurs. This mechanism is very much similar to what happens in the case ofthe well-know explosive synchronization in the Kuramoto model [56]. However,we emphasize that the explosive synchronization we have observed only occursfor the stronger electrical synapses, and we did not observe any explosive syn-chronization for chemical synapses in the range of parameters studied here. Wefinally note that explosive synchronization occurs at much higher values of g when compared to the WS case and also exhibits a smaller hysteresis loop. Figure 4:
Synchronization diagram for the correlated scale-free network (a) electricalsynapses and (b) chemical synapses, and the corresponding uncorrelated case (c) and(d). Raster plots for various g are shown in parts (e)-(h) for the case of explosivesynchronization in panel (a). t = 0 indicate the beginning of stationary state. Thesynchronization diagrams show the averaged results over five network realizations andinitial conditions. .5. Hierarchical modular networks Figure 5: (a) Adjacency matrix of the HM network of size N = 512, coordinationnumber z = 15. There are 4 hierarchical levels. 16 modules in level 1, 8 modules inlevel 2, 4 modules in level 3 and 2 modules in level 4. (b) and (c) Synchronizationdiagram for electrical and chemical synapses, respectively. (d)-(i) Raster plots forneural network with electrical synapses for six different values of g . t = 0 indicatesthe beginning of stationary state. The synchronization diagrams show the averagedresults over five network realizations and initial conditions. So far we have investigated synchronization transition of spiking HH neu-rons for the most typical network topologies. However, it is believed that ahierarchical modular (HM) network is a more realistic representation of the ac-tual structural connectivity of the cortex on the large scale [58]. Therefore itis worthwhile to investigate synchronization of HH neurons in HM networks aswell. We construct networks with N = 512 nodes, z = 15 and 4 hierarchical20evels, see Fig.5(a). On the lowest level, the network is comprised of 16 modulesof 32 nodes where each node is connected to 10 randomly chosen other nodeswithin the same module. On the second level, there are 8 modules which areconstructed by connecting pairs of randomly chosen nodes belonging to twosmaller modules of the lower level. In this manner we construct each modulein a higher level by connecting members of two modules in the previous level.We should note that all links in the first level are inter-modular and all linksin higher levels are intra-modular. This network despite having considerableclustering also has small-world effect. S − g plots for HH neurons with electrical and chemical synapses for such aHM network is displayed in Figs.5(b) and 5(c), respectively. We observe that forboth types of synapses, there exists three regimes in the S − g plots. An asyn-chronous regime ( S = 0 .
5) for small values of g , a synchronous regime for largevalues of g and an intermediate regime between ordered and disordered phaseswhere S does not vary monotonically with increasing g , but reveals a fluctuat-ing behavior. Note that these fluctuations in order parameter are not due toinsufficient transient time because we have made sure that the system is in itsstationary state before taking measurement of S for each value of g . They alsodo not appear due to imprecise numerical integrations as we have taken muchcare in this regard. Emergence of this intermediate regime has been reported forthe Kuramoto model on human connectome network which has a hierarchicalmodular structure [59]. It is now believed that such intermediate regime is amanifestation of HM structure of the underlying network. As can be deducedfrom the raster plots in Figs.5(d-i), the HM structure leads to synchronizationwithin various modules which are themselves out of phase with various othermodules leading to relative synchrony (small S ) which nevertheless fluctuatesas various modules go in and out of phase with each other as we change thevalue of g . For example, for electrical synapses and g = 0 .
019 there is more syn-chronization than g = 0 .
020 as can clearly be seen why from the correspondingraster plots. Therefore, we observe the same type of frustrated synchronizationpatterns as in the Kuramoto model regardless of the synaptic interaction. We21lso note that, looking at the values of transition point g t , one sees strong simi-larity with fully random networks of SF and ER (see Fig.2). This indicates thatthe onset of synchronization here is also dictated by long-range links. But, oncesynchronization sets in, it is the strong clustering within various modules thatdictate the synchronization pattern for a range of g , before global order sets infor large enough g .
5. Concluding remarks
In a previous work, we studied beta-band synchronization in network modelsof spiking neurons. There, we showed that the type of synchronization transitionoccurring in a neural network depends on the firing rates of constituent neurons[24]. In this paper we have reported a systematic study of synchronization tran-sition in network models of spiking neurons in gamma-band. We employed HHneurons with electrical and chemical synapses. Our focus has been to charac-terize the combined effect of synaptic type and topological features on the typeof synchronization transitions that may occur. The mechanisms and patterns ofsynchronization transitions we obtained here for gamma-rhythms are distinctlydifferent from those we obtained for beta-rhythms in ref.[24]. For example, inthe beta-band in a one-dimensional lattice, we found a continuous transition forthe electrical synapses, while here for the gamma-band we observed no tran-sition for any synaptic type in one or two dimensions. Furthermore, here wefound smooth transitions for SF and ER networks in the gamma-band, whilepreviously we had observed explosive synchronization on such networks withelectrical synapses. On the other hand, here we observe explosive synchroniza-tion for the WS networks while in the beta-band we only saw a smooth transitionfor such networks. We also observed explosive synchronization in the case of SFnetworks with correlated heterogeneity which was not studied in the previousstudy for the beta-band.The underlying mechanisms leading to explosive synchronization in beta-band was rooted in the formation of anti-phase groups of neurons for intermedi-22te values of g and their sudden combination at a transition point [24]. This isdistinctly different from the mechanism that lead to explosive synchronizationin WS network of HH neurons or from the mechanism resulted in abrupt tran-sition in SF network of HH neurons with correlated heterogeneity. However,these three mechanisms of explosive synchronization have a common aspect.They all occur through electrical synapses. We have not observed explosivesynchronization in the case of chemical synapses. Our results highlight the factthat electrical synapses are more conducive to synchronization and can in factlead to entirely different transition patterns. This is in contrary to other studiesthat have concluded similar synchronization behavior for electrical and chemicalsynapses [38].We also note that it is interesting that our regular one and two dimensionallattice did not exhibit a transition which is what one would expect from thestudy of the Kuramoto model as it, too, does not exhibit a transition in low di-mensional systems [48]. However, the mechanism for such behavior are different,as we do observe considerable amount of order in our system with quasiperi-odic oscillations. We once again emphasize the key role of frequency as well assynaptic interactions in such studies, where in the beta-band, in one dimension,we had previously observed a continuous transition for the electrical synapseswhile no transition was seen for the chemical synapses.This brings us to emphasize the difference in the type of transitions we ob-served in WS networks. Electrical synapses, which are strong and fast, lead toexplosive synchronization while the slower and weaker chemical synapses lead toa smooth transition. This is particulary interesting as neuronal networks are ar-gued to by on the verge of a phase transition. This could, for example, be relatedto the fact that electrical synapses are useful in fast involuntary motor responsewhere a strong and fast collective action is desired, while a smooth transitionwith chemical synapses could be understood in terms of cortical neurons wheretoo much synchronization is deemed to be pathological [22].Furthermore, we investigated the role of correlated heterogeneity and foundthat in the case of electrical synapses one observes explosive synchronization23hile in chemical synapses a smooth transition occurs. This result could be in-teresting from two aspects. First, it shows that unlike what is generally believed,correlated heterogeneity does not always lead to explosive synchronization aschemical synapses showed a smooth transition. Secondly, it highlights the dis-tinctly different type of explosive synchronization that may occur in electricalsynapses. In WS network, explosive synchronization occurred after the systemgained a high degree of local order, but in the case of correlated heterogeneity,explosive synchronization was dictated by the role of the hub with no sign oforder in the system just before the transition occurred.We have also considered hierarchical modular networks which resulted inan intermediate regime between order and disorder. Such a behavior has beenpreviously shown to occur in the Kuramoto model [59] and our results indicatethat such frustrated transition is a more general property of neuronal systemsin HM networks, and it is furthermore independent of the type of synapticinteraction.We note that our choice of spiking HH neurons naturally limited our rangeof frequency to the gamma band. However, we observed the same type ofsynchronization patterns when we increased the natural spiking frequency ofthe HH neurons to the high gamma band up to f ≈
6. Acknowledgements
Support from Shiraz University research council is kindly acknowledged.This work has been supported in part by a grant from the Cognitive Sciencesand Technologies Council.
ReferencesReferences [1] H. B. Callen,
Thermodynamics and an Introduction to Thermostatistics ,Wiley, New York (1998).[2] S. K. Ma,
Modern theory of critical phenomena , Routledge (2018).[3] G. Weisbuch,
Complex systems dynamics , CRC Press, New York (2018).[4] A. Steyn-Ross and M. Steyn-Ross,
Modeling phase transitions in the brain ,Springer (2010).[5] D. Markovic and C Gros, Phys. Rep. , 41-74 (2014).[6] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. , 381 (1987).[7] P. Bak, K. Chen, and M. Creutz, Nature , 031001 (2018).[9] D.R. Chialvo, Nature physics (10), 744 (2010).2510] P. Moretti and M.A. Munoz, Nature Communications , 2521 (2013).[11] S. A. Moosavi, A. Montakhab, and A. Valizadeh, Sci. Rep. , 7107 (2017).[12] M. Rohden, A. Sorge, M. Timme, and D. Witthaut, Phys. Rev. Lett. (6), 064101 (2012).[13] B. Blasius, A. Huppert, and L. Stone, Nature , 6734, 354 (1999).[14] L. Stone, R. Olinky, and A. Huppert, Nature , 7135, 533 (2007).[15] S. di Santo, P. Villegas, R. Burioni, and M. A. Munoz, PNAS (7),1356-1365 (2018).[16] M. Khoshkhou and A. Montakhab, Front. Syst. Neurosci. (73) (2019).[17] A.J. Fontenele, et al. Phys. Rev. Lett. , 208101 (2019).[18] F. Varela, J. P. Lachaux, E. Rodriguez, and J. Martinerie, Nat. Rev. Neu-rosci. , 229-239 (2001).[19] A. Buehlmann and G. Deco, PLoS Comput. Biol. , e1000934 (2010).[20] P. Esir, A. Simonov, and M. Tsodyks, Front. Comput. Neurosci. , 21(2017).[21] J. Fell and N. Axmacher, Nat. Rev. Neurosci. , 105-118 (2001).[22] E. R. Kandel, J. H. Schwartz, and T. M. Jessell, Principles of Neural Sci-ence , McGraw-Hill, New York (2000).[23] D. Gireesh Elakkat and D. Plenz. PNAS , 7576 (2008).[24] M. Khoshkhou and A. Montakhab, Front. Comput. Neurosci. , 59 (2018).[25] C. Gros, Complex and Adaptive Dynamical Systems , 4th ed., Springer-Verlag Berlin Heidelberg (2015).[26] G. Buzsaki,
Rhythems of the Brain , Oxford University Press, New York(2006). 2627] X. Jia and A. Kohn, PLoS Biol. (4), e1001045 (2011).[28] J. A. Henrie and R. Shapley, J. Neurophysiol. , 479490 (2005).[29] P. Fries, D. Nikolic, and W. Signer, Trends Neurosci. (7), 309-316 (2007).[30] J. R. Vidal, M. Chaumon, J. Kevin, O. Regan, and C. Tallon-Baudy, J.Cog. Neurosci. (11), 1850-1862 (2006).[31] P. Berens, G. A. Keliris, A.S. Ecker, N. K. Logothetis, and A. S. Tolias,Front. Neurosci. , 199207 (2008).[32] J. Liu and W. T. Newsome, J Neurosci (1), 16-23 (2014).[35] E. P. Bauer, R. Paz and D. Pare, J. Neurosci. , 155168 (2006).[37] A. L. Hodgkin and A. F. Huxley, J. Physiol. , 500 (1952).[38] T. Perez, G. C. Garcia, V. M. Eguiluz, R. Vicente, G. Pipa, and C. Mirasso,PLoS ONE (5): e19900 (2011).[39] Y. Wang, D. T. W. Chik, and Z. D. Wang, Phys. Rev. E , 740 (2000).[40] O. Kwon and H. T. Moon, Phys. Lett. A , 319-324 (2002).[41] O. Kwon, K. Kim, S. Park, and H. T. Moon, Phys. Rev. E , 021911(2011).[42] Q. Wang, M. Perc, Z. Duan, and G. Chen, Phys. Lett. A , 5681-5687(2008).[43] C. Park, R. M. Worth, and L. Rubchinsky, Phy. Rev. E , 042901 (2011).2744] C. Zhou and J. Kurths, Chaos (1), 401-409 (2003).[45] T. de L. Prado, S. R. Lopes, C. A. S. Batista, J. Kurths, and R. L. Viana,Phys. Rev. E , 032818 (2014).[46] C. A. S. Batista, R. L. Viana, F. A. S. Ferrari, S. R. Lopes, A. M. Batista,and J. C. P. Coninck, Phys. Rev. E , 042713 (2013).[47] Y. Kuramoto, Lect. Notes Phys. , 420-422 (1975).[48] J. A. Acebron, L. L. Bonilla, C. J. Perez Vicente, F. Retort, and R. Spigler,Rev. Mod. Phys. , 1 (2005).[49] S. A. Moosavi and A. Montakhab, Phys. Rev. E , 3292 (1998).[51] A. Roth, A. and M. van Rossum, Computational Modeling Methods forNeuroscientists , Cambridge, MA: MIT Press (2009).[52] A. Pikovsky, M. Rosenblum, and G. Osipov, Physica D , 219 (1997).[53] C. van Vreeswijk, Phys. Rev. E , 5522 (1996).[54] M. Rosenblum and A. Pikovsky, Phys. Rev. Lett. , 064101 (2007).[55] R. Burioni, S. Santo, M. Volo, and A. Vezzani, Phys. Rev. E , 042918(2014).[56] J. Gomez-Gardenes, S. Gomez, A. Arenas, and Y. Moreno, Phys. Rev. Lett. , 128701 (2011).[57] I. Leyva, R. Sevilla-Escoboza, J. M. Buldu, I. Sendima-Nadal, J. Gomez-Gardenes, A. Arenas, Y. Moreno, S. Gomez, R. Jaimes-Reategui, and S.Boccaletti, Phys. Rev. Lett. , 168702 (2012).[58] O. Sporns, Networks of the Brain , The MIT Press, Cambridge (2011).[59] P. Villegas, P. Moretti, and M. A. Munoz, Sci. Rep. , 5990 (2014).2860] M.M. Asl, A. Valizadeh, P.A. Tass, Chaos , 106308 (2018).[61] S. A. Moosavi, A. Montakhab, and A. Valizadeh, Phy. Rev. E (2), 022304(2018).[62] S. A. Moosavi and A. Montakhab, Phys. Rev. E89