Effects of neuronal variability on phase synchronization of neural networks
Kalel Luiz Rossi, Roberto Cesar Budzisnki, Joao Antonio Paludo Silveira, Bruno Rafael Reichert Boaretto, Thiago Lima Prado, Sergio Roberto Lopes, Ulrike Feudel
EE ff ects of neuronal variability on phase synchronization of neural networks K. L. Rossi a, ∗ , R. C. Budzinski a , J. A. P. Silveira a , B. R. R. Boaretto a , T. L. Prado a , S. R. Lopes a , U. Feudel b a Department of Physics, Universidade Federal do Paran´a, 81531-980 Curitiba, Brazil b Theoretical Physics / Complex Systems, ICBM, Carl von Ossietzky University Oldenburg, 26111 Oldenburg, Germany
Abstract
An important idea in neural information processing is the communication-through-coherence hypothesis, according to which com-munication between two brain regions is e ff ective only if they are phase-locked. Also of importance is neuronal variability, aphenomenon in which a single neuron’s inter-firing times may be highly variable. In this work, we aim to connect these two ideasby studying the e ff ects of that variability on the capability of neurons to reach phase synchronization. We simulate a network ofmodified-Hodgkin-Huxley-bursting neurons possessing a small-world topology. First, variability is shown to be correlated withthe average degree of phase synchronization of the network. Next, restricting to spatial variability - which measures the deviationof firing times between all neurons in the network - we show that it is positively correlated to a behavior we call promiscuity,which is the tendency of neurons to to have their relative phases change with time. This relation is observed in all cases we tested,regardless of the degree of synchronization or the strength of the inter-neuronal coupling: high variability implies high promiscuity(low duration of phase-locking), even if the network as a whole is synchronized and the coupling is strong. We argue that spatialvariability actually generates promiscuity. Therefore, we conclude that variability has a strong influence on both the degree and themanner in which neurons phase synchronize, which is another reason for its relevance in neural communication. Keywords:
Neural networks, Clustering, Variability, Phase synchronization, Hodgkin-Huxley
1. Introduction
The brain manages to robustly and e ffi ciently process hugeamounts of information [1, 2, 3]. Of the various di ff erent waysin which this information can be encoded, one then expects thatthey exhibit at least some degree of robustness. Two possiblecoding strategies are the rate and temporal codes, in which in-formation is encoded in the rate and in the timing of spikes,respectively [4]. Considering that neurons can have variationsin their inter-firing intervals [5, 6, 7, 8, 9, 10, 11], it is impor-tant to consider the interactions between this variability and thesupposed robustness of the two codes. The neuronal variabilitymay represent unwanted noise insofar as it causes fluctuationsin the firing rate [12]. In this case, for a constant input, onewould expect a constant firing rate coding it, which may in re-ality not happen because of the variability. On the other hand,variability may represent extra information since it a ff ects thetiming of spikes [12], in which case the di ff erences in the tim-ings, generated by the variability, can carry information.Besides these e ff ects on the codes, variability may also havean indirect e ff ect on information processing because it can influ-ence the phase synchronization of networks. This phenomenonof phase synchronization happens when neurons start their fir-ing at or su ffi ciently close to the same time. It is related to var-ious memory and conscious processes [13, 14, 15] and is alsoimportant in information processing, where it serves to increase ∗ Corresponding author the saliency of neural responses, thus facilitating their linking(binding) [13, 16, 17].Another important example relating oscillatory behavior toinformation coding is given by Fries [17]. It relies on the ob-servation that neuronal groups have a tendency to oscillate andthat these oscillations a ff ect the likelihood of spike output andsensitivity to synaptic input. The author proposes, that theseoscillations create windows of e ff ective communication, whenregions have their output and input sensitivity coinciding [17].Further, a phase can be defined for an oscillating time seriesto represent where in the oscillation cycle the series is. In thisway, coherently oscillating regions can be viewed as a phase-locked (their relative phase is kept constant). With this termi-nology, we can rephrase the communication-through-coherencehypothesis stated above to say that communication is e ff ectiveonly between phase-locked regions.The e ff ect that neural variability has on phase synchroniza-tion and phase-locking is still unclear. In this work, we in-vestigate this relation by simulating a network of temperature-dependent Hodgkin-Huxley-type neurons, proposed by Braun et al. [18], and henceforth called HB neurons, coupled in asmall-world and random topologies with excitatory chemicalsynapses.A network topology, also sometimes called connection scheme,refers to the structure in which neurons are connected. A small-world topology is characterized by having a short average pathlength (distance) between neurons, yet with high local cluster-ing [19]. From a neuronal information processing point of view, Preprint submitted to Neural Networks March 9, 2020 a r X i v : . [ q - b i o . N C ] F e b his is very appealing since it supports both segregated (local)and distributed (global) behaviors [20]. Small-world networksalso have been shown to have high e ffi cienciy, due to having ahigh complexity even with only few nonlocal connections [21],and to optimize wiring and energy costs [20, 22]. Supportingthis theoretical attractiveness, these topologies have been ob-served in a wide range of cases, such as in the nervous systemof C. elegans [19], the macaque monkey, the cat [23], and in thevertebrate brainstem reticular formation [24].Neurons in the HB model present bursts, which are rapidsequences of spikes followed by a long quiescent period. Burst-ing neurons are thought to be important in information process-ing [25]: compared to single spikes, bursts are more reliablytransmitted to postsynaptic neurons [26] and appear to trans-mit at least the same amount of information [27]. Bursts alsohave a greater ability to elicit responses on postsynaptic neu-rons [25, 28].Networks with the HB model have already had their syn-chronization characteristics studied in detail for di ff erent cou-pling strengths or temperatures in a small-world and scale-freetopologies [29, 30, 31, 32, 33, 34, 35, 36]. A transition fromdesynchronization, for very low coupling, to chaotic burst syn-chronization, for su ffi ciently strong coupling, was generally ob-served. The specific behavior between these two cases, how-ever, was seen to di ff er strongly for di ff erent temperatures [31,36].For some values of temperature, in the coupling region be-fore the transition to chaotic burst synchronization, the networkwas observed to be always desynchronized [31, 36]. Increas-ing the temperature, this always desynchronized region changesto one with a local maximum in the synchronization quanti-fier, meaning that with the increase in coupling the networkchanges from desynchronized to burst synchronized and then todesynchronized again before the chaotic burst synchronization[29, 33]. Further increase of the temperature makes the networkstart to synchronize even for very low coupling strengths. Forsome values of coupling, the synchronization is so strong thateven the spikes within bursts are synchronized [31]. That is,three temperature intervals were observed for which networkshave distinct synchronization behavior before an eventual tran-sition to chaotic burst synchronization.Studying the uncoupled neurons, it was observed [31] thatthere were also three di ff erent behaviors: for lower tempera-tures, neurons have various Inter-Burst Interval ( IBI ) values,indicative of chaotic behavior. With an increase of temperature,neurons then have two
IBI values. A further increase leads toonly one
IBI value. The three synchronization behaviors forweakly coupled networks are related to the three
IBI behaviorsfor uncoupled neurons: for weak coupling, networks with tem-perature values corresponding to chaotic behavior in uncoupledneurons did not synchronize. In this same regime, for tempera-tures in which uncoupled neurons had two
IBI s the coupled net-work had the local maximum in synchronization, as describedpreviously. Still for the same coupling strengths, for temper-atures in which uncoupled neurons had one
IBI the networksynchronized strongly.Therefore, it was concluded that, for these coupling strengths, the synchronization behavior of the coupled networks were linkedto the local dynamics of the uncoupled neurons [31].In the strong coupling regime, where networks achievedchaotic phase synchronization in a very similar way for all tem-peratures investigated, it was concluded that the forcing be-tween neurons is strong enough to synchronize the network re-gardless of what the uncoupled dynamics is. This also leadsto the remark that the individual dynamics only plays an im-portant role for weak coupling. Applying a pulsed current inthe network, the burst synchronization for weak coupling in thesecond temperature interval was successfully suppressed, whilestill maintaing the asymptotic synchronization [33].In this work, we extend some of these results and showthat the variability, measured by the standard deviation of theInter-Burst Intervals (
IBI ) divided by their mean value, cor-relates with the degree of synchronization, measured by theKuramoto Order Parameter [37], for all the temperature valuescorresponding to the di ff erent behaviors mentioned previously.We also show that there is a phenomenon we call promiscuity,which is the tendency of the relative phases between di ff erentneurons to change in time. The phases are calculated from thestart of each neuron’s firing (in the HB case, their bursts) andthe relative phase is just the di ff erence between the phases ofdi ff erent neurons. In this case, if promiscuity is zero, neuronsare phase-locked and the higher the promiscuity the less neu-rons stay in phase with each other.Promiscuity is shown to be di ff erent for network states withsimilar degrees of synchronization, but di ff erent coupling strengths- therefore allowing one to distinguish between them. We mea-sure promiscuity in two di ff erent ways. The first is by countinghow many neurons leave and enter the network cluster - definedhere as all neurons within one standard deviation of the meanvalue of each bursting event.The second way calculates how much, on average, the burststart times of di ff erent neurons drift in relation to each other.The higher this value, the higher the probability that the dif-ference between neurons’ bursts times changes over time. Thisaverage drift method has the advantage of being parameter-free,corroborates the results obtained by the clustering analysis andcan be easily applied to any time series of inter-burst or inter-spike intervals.The results obtained were for both small-world and randomnetworks and are very similar.We finally argue that variability generates promiscuity, dueto a simple statistical process, described in the conclusions.The paper is organized as follows: in section 2 we introducethe methods we use - including the neuronal model and the syn-chronization quantifiers -, we also define the two types of vari-ability, the clustering algorithm and the average drift measure.In section 3 we show the results and analysis, presenting, fi-nally, the conclusions in section 4.
2. Methods
To simulate the neurons, we used an adapted Hodgkin-Huxleymodel, proposed by Braun et al [18]. It consists of six di ff eren-2ial equations and introduces the influence of temperature. Themembrane potencial for the ith neuron at time t , V i ( t ), is givenby the master equation C M dV i ( t ) dt = − J i , Na − J i , K − J i , sd − J i , sr − J i , L − J i , coup , (1)in which C M is the specific membrane capacitance; J i , Na , J i , K ,and J i , L are the Sodium, Potassium, and Leakage current den-sities, respectively; J i , sd , J i , sr are current densities related to thesubthreshold depolarization (sd) and repolarization (sr); J i , coup is the current density due to the interneuronal coupling.The specific current densities J α , for α = { Na , K , sd , sr , L } ,are described by the Nernst potential E α of the respective ion orleak channels: J i , Na = ρ ¯ g Na a i , Na ( V i − E Na ) , (2) J i , K = ρ ¯ g K a i , K ( V i − E K ) , (3) J i , sd = ρ ¯ g sd a i , sd ( V i − E sd ) , (4) J i , sa = ρ ¯ g sa a i , sa ( V i − E sa ) , (5) J i , L = ¯ g L ( V i − E L ) . (6)The g α are maximum specific conductances, ρ is the first scalefactor, which serves to introduce the temperature dependence,and is given by ρ = ρ ( T − T ) / T . (7)The a α are variables responsible for the channel activations,whose temporal evolutions are da i , Na dt = φτ Na ( a i , Na , ∞ − a i , Na ) , (8) da i , K dt = φτ K ( a i , K , ∞ − a i , K ) , (9) da i , sd dt = φτ sd ( a i , sd , ∞ − a i , sd ) , (10) da i , sa dt = φτ sa ( − η J i , sd − γ a i , sa ) (11)where the parameters η and γ describe the increase and de-crease, respectively, of the concentration of Ca + ions. Thesecond scale factor φ is φ = φ ( T − T ) / T . (12)The a i ,α, ∞ ( α = { Na , K , sd } ) are activation functions which de-pend on the membrane potential through the equations a Na , ∞ = + exp[ − s Na ( V i − V )] , (13) a K , ∞ = + exp[ − s K ( V i − V )] , (14) a sd , ∞ = + exp[ − s sd ( V i − V )] . (15)All previously unmentioned terms are constants, given intable 1. Finally, the coupling current J i , coup is J i , coup = g c εν ( V i − E syn ) (cid:88) j ∈ Γ i r j ( t ) , (16) where ε is the coupling parameter, which controls the strengthof the coupling; ν is the normalization factor, defined as themaximum degree of connectivity of the network; g c is a unitaryparameter with conductance units; E syn is the synaptic reversalpotential, whose value is taken to ensure that the synapse isalways excitatory; Γ i is the set of neurons connected to the i -thneuron; r j is the coupling variable and has a temporal evolutiongiven by [38] dr j dt = (cid:32) τ r − τ d (cid:33) − r j + exp[ − s ( V j − V )] − r j τ d . (17)We note that the coupling current can be written as J coup = ¯ gP ( V − E syn ) [4], so that we can identify ¯ g = g c ε/ν as themaximum conductivity of the postsynaptic membrane and thesummation as the fraction of bound postsynaptic receptors P .Thus, r j can be seen as the fraction of bound receptors of the i -th neuron that can be activated by the j -th neuron.The integration is done using the CVODE solver [39], witha 12th order Adams-Moulton predictor-corrector method, ab-solute and relative tolerances of 10 − , and maximum time step h = .
01 ms. The initial conditions are given randomly in theinterval [ − ,
0] for V and [0 ,
1] for the other variables. Sim-ulations are run for 25 s, 15 s of which are considered to betransient dynamics and disregarded.The parameter values were obtained from [18], with two al-terations done for convenience, which do not alter the dynam-ical behavior: (i) temperatures were rescaled, by changing T ,so that the relevant behavior lies in a range plausible for mam-mals (around 37 ◦ C); (ii) characteristic times were rescaled sothat the frequency of bursts was around 10 Hz. Since the systemis invariant under time translations, this change only a ff ects thevalues of the burst times: except for these values, the dynamicbehavior is exactly the same.In Fig. 1 we show a typical membrane potential (blue line)for an uncoupled neuron and the corresponding burst start times(orange circles). t (ms) − − − V ( m V ) Figure 1: Membrane potential V of an uncoupled neuron, represented as theblue line, for T = ◦ C and the burst start times, represented as the orangecircles. able 1: Parameter values of the constants for the Huber-Braun neuron model [18]. Membrane capacitance C m = . µ F / cm Maximum conductances (mS / cm ) g Na = . g K = . g sd = . g sa = . g l = . g c ≡ . τ Na = . τ K = . τ sd = τ sa = τ r = . τ d = . E Na = E K = − E sd = E sa = − E l = − V = − V = − V = − E syn = ρ = . φ = . T = ◦ C T = ◦ C s Na = .
25 mV − η = .
012 cm / µ A s sd = .
09 mV − γ = . s K = .
25 mV − s = . − Network parameters: N = N = ν = We used the small-world topology proposed by Newmanand Watts [40], in which a regular network of size N and k = N − kN of con-nections is added, thus totalling N connections. This numberof added connections is chosen in such a way that the aver-age path length of the network is low (comparable to a randomtopology) and the clustering coe ffi cient is high (much higherthan random), so that this network is locally clusterized, butstill neurons are not so far away (in a number of connectionssense) from each other [19].The network generated with this algorithm is directed, mean-ing that if neuron i is connected j , then j is not necessarilyconnected to i . This is chosen to be so because the synapsesare chemical and therefore not symmetrical. The network alsohas no self-loops, since we consider that neurons don’t connectback to themselves.Another algorithm, proposed by Watts and Strogatz [19]also has been used to generate the small-world networks, withsimilar results being obtained. Phase synchronization is a state in which the oscillators ofthe system have similar phases, but may have di ff erent ampli-tudes [41]. To quantify the degree of synchronization of a neu-ral network, we first introduce the phase θ i ( t ) for each neuron.We define it as starting in 0, increasing by 2 π every time a burststarts and being a linear interpolation in between [42] θ i ( t ) = π k i + π t − t k , i t k + , i − t k , i , ( t k < t < t k + ) (18)where t k , i is the time at which the k -th burst of the i -th neuronocurred, called the burst start time.We remark that this definition works for all parameter val-ues studied in this work. For coupling strengths higher than theones used, single spikes start appearing a significant distancebetween adjacent bursts and so a problem arises when defininga phase. In the cases we studied, however, this is not significantbecause these events are rare and the distances are small, so theisolated spike are considered to belong to the previous burst. The measurement of the phase synchronization is done viathe Kuramoto order parameter [37] R ( t ) = N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N (cid:88) j = e i θ j ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (19)where i is the imaginary unit here. The quantifier R is 0 ifthe network is completely desynchronized (each neuron has an-other neuron with opposite phase) and is 1 if the network iscompletely synchronized in phase (every neuron has the samephase). We may take the time average, (cid:104) R (cid:105) = n t f (cid:88) t = t R ( t ) , (20)where t is the transient time and t f is the total simulation timeand n = ( t f − t ) / h is the number of steps, with h being the timestep. We define the k -th Inter-Burst Interval (IBI) of the i -th neu-ron as the di ff erence between its k th and ( k + IBI k , i = t k + , i − t k , i . (21)The variability is measured by the coeficient of variability CV , defined as [43, 44] CV = σ ( IBI ) (cid:104) IBI (cid:105) , (22)where (cid:104) IBI (cid:105) is the time average taken over the IBIs of all neu-rons and σ ( IBI ) is the standard deviation of the IBIs, whichmay be calculated in 2 di ff erent ways: (1) by taking the spatialaverage of the deviation of a neurons’ IBI over time, denoted CV t ; (2) by taking the time average of the deviation over space,denoted CV s : σ t = N N (cid:88) i = σ (cid:18) k max (cid:88) k = IBI k , i (cid:19) (23) σ s = k max k max (cid:88) k = σ (cid:18) N (cid:88) i = IBI k , i (cid:19) , (24)in which k max denotes the total number of IBI s analyzed.4 .5. Clustering Analysis
In a phase synchronized network, neurons’ burst times tendto be distributed closely around a mean value for every burstingevent. In panel (b) of Fig. 2 are depicted the burst start timesof all neurons in a network for one single event for T = ◦ Cand ε = . . . . . C o un t (a)
700 800 900 t (ms) N e u r o n (b) Figure 2: Panel (b) depicts the network raster plot for T = ◦ C and ε = . To quantify this observation, we first define a cluster to bethe set of neurons that are within one standard deviation fromthe mean value of the bursting event. In Fig. 2, the mean valuecorresponds to the blue line, and the cluster to the green area.We are interested in studying how the relative bursting timesof neurons change with time. In this way, it is appropriate tostudy the changes in the cluster composition throughout time.With this in mind, we generalize the cluster definition to con-sider the behavior for multiple bursting events in the followingway: identify the clusters C k for every bursting event k for T events and consider only neurons who are in all clusters. Thisintersection C ∩ C ∩ · · · ∩ C T will be called a T -cluster andthe collection of T bursting events is called a cluster event.This can be summarized in the following algorithm1. Identify the approximate time of the bursting event througha local maximum on the histogram of the burst times.2. Choose the burst times that are su ffi ciently close (withinone half of a bursting period) for each neuron.3. Take the mean value of the burst times obtained above,resulting in the mean burst time.4. Collect all neurons whose burst time is within one stan-dard deviation from the mean burst time in a cluster.5. Repeat the above for T events.6. Take the intersection of all the clusters obtained above toget the T -cluster. The number of neurons inside each cluster C k is naturallydefined as the cluster size CS k . Once all clusters are identi-fied, we select two adjacent clusters, described by the sets C k and C k + , respectively, and count the number L k of neurons thatleave the cluster (i.e. the size of the set C k \ C k + ). Finally, takingthe average over various bursting events, we obtain the meancluster size CS and the mean number of neurons that leave theclusters L . A downside of this clustering analysis is its neces-sity to identify what we call a bursting event and identifyingwhich of a neuron’s burst start times correspond to this event.This relies on a certain minimum degree of synchronization ofthe network: if the network is unsu ffi ciently synchronized, theanalysis cannot be applied. In that case, another method needsto be considered. We introduce such a method in the next sub-section (2.6). In this section, we define a quantity to measure if the tempo-ral distances between neurons’ bursts remain locked or changein time. To do this, we start with the time t k , i of the start of k -th burst of the i -th neuron, and then calculate the temporaldistance of each pair of neurons ( i , j ) for each event k δ ki j = | t k , i − t k , j | . (25)Then we calculate the absolute di ff erence of this distancebetween sucessive events for each pair ( i , j ), ∆ li j = | δ ki j − δ k − i j | . (26)The temporal average of ∆ i j , (cid:104) ∆ i j (cid:105) , measures the tendencyof the pair ( i , j ) to drift away from each other across time. Weaverage the result over all pairs of neurons, resulting in ∆ ≡(cid:104) ∆ i j (cid:105) , which is termed the average drift of the network. Thisaverage drift tendency ∆ measures how much, on average, thetemporal distances between neurons’ firings change. If it is low,neurons are locked together, the di ff erence between their burststart times remaining fixed. If it is high, neurons are not phase-locked.
3. Results and discussions
The main scenario of synchronization and the neurons’ vari-ability as a function of the coupling strength ε is depicted inFig. 3. The degree of synchronization, measured by the time-averaged Kuramoto order parameter (Eq. 20), is shown in thefirst row of Fig. 3 for di ff erent values of temperature T .The coupling strength values are in the interval [0 , . T = ◦ C, in which amonotonic transition to phase synchronization is observed, sim-ilar to that noticed in Kuramoto Oscillators [37]. In panel (b),for T = ◦ C, there is a local maximum of the synchronizationquantifier in the weak coupling regime, characterizing a non-monotonic transition [29]. Panel (c) corresponds to T = ◦ C.In this case, the local maximum of (cid:104) R (cid:105) is replaced by a global5 . . . . . h R i (a) (b) (c) . . . . . . C V (d) (e) (f) CV t CV s ε I B I ( m s ) (g) ε (h) ε (i) − − − − − − − − − l n ( λ ) Figure 3: The scenario of synchronization and its relation to the neurons’ variability as a function of coupling strength for di ff erent values of temperature. The firstrow corresponds to the mean Kuramoto order parameter ( (cid:104) R (cid:105) ); the second to the coe ffi cients of variability ( CV t and CV s ); the third to the IBI bifurcation diagram.The columns correspond, from left to right, to temperatures T = ◦ C, T = ◦ C, and 40 ◦ C. Di ff erent behaviors are observed: in the first column ( T = ◦ C), atraditional transition from non-synchronized to a phase-synchronized state is observed, with variabilities being high. Increasing the temperature to T = ◦ C and T = ◦ C, in the second and third columns, leads the system to a di ff erent scenario, with a non-monotonic synchronization level phenomena. In these cases, thelocal maxima of (cid:104) R (cid:105) happens just where variabilities ( CV t and CV s ) are low. Results are given by an average over 10 initial conditions. maximum, in which spike synchronization can be observed forweak coupling, as reported in [31].We notice that the network has very di ff erent behaviors forweak coupling, ranging from desynchronization to burst syn-chronization to even complete synchronization by a change oftemperature T . For high coupling strengths, roughly in [0 . , . T , due to the networkhaving chaotic phase synchronization [31].The second row of Fig. 3 shows the variabilities, as definedin section 2.4: temporal variability ( CV t ) is in blue and spatialvariability ( CV s ) is in orange. Panel (d) ( T = ◦ C) depictshigh values of CV for small values of ε , which decreases inthe high coupling regime. Panel (e) ( T = ◦ C) depicts a CV minimum for the weak coupling regime, which is followed bya maximum. Increasing the coupling further, CV decreases, butnot below the previous minimum. In panel (f) ( T = ◦ C),we see a significant decrease in the two variabilities. Followingthe periodic behavior in the uncoupled case, the variabilties arevery close to zero for very weak coupling ( ε (cid:46) . CV s start to increase and a corresponding decrease in the degree of phase synchronization (cid:104) R (cid:105) is observed in panel(c). This culminates in a maximum of the CV s, which happenstogether with a minimum in (cid:104) R (cid:105) .These results show a clear correlation for the weak couplingregime ( ε (cid:46) . CV t or CV s ) is associatedwith low degree of synchronization ( (cid:104) R (cid:105) ), as in T = ◦ C andlow variability is associated with high synchronization (for both T = ◦ C and T = ◦ C). This relation is especially clear forthe last two temperatures, in which regions with almost zerovariability are the ones with almost complete synchronization.The third row of Fig. 3 shows the bifurcation diagrams forthe Inter-Burst Intervals of neurons in the network. The colorof each point is determined by the logarithmic density ln( λ ) ofthe corresponding IBI . The density λ is equal to the number oftimes the IBI was observed divided by the total number of
IBI sobserved. The colorbar, shown on the right, has values in theinterval [ − , − IBI s, par-6icularly in the width of the probability density function, lead-ing to larger standard deviations or in the location of its max-imum value leading to a shift in the mean value. The changesin variability are related to changes in the width and each ofthe transitions described above in the slope of variability can beassociated to such changes in the width of the
IBI
PDF.An important detail consists in the behavior for tempera-tures T = ◦ C , ◦ C, in which the network desynchronizesas coupling strength is increased. The analysis suggests that, inthese cases, the coupling strength begins to be strong enough toinduce a more chaotic behavior (reflected in the increase of thevariability), which desynchronizes the network.A further increase in the coupling strength makes the net-work assume a chaotic burst synchronization. This result is inline with previous observations [31, 36] that the local, individ-ual behavior of the neurons is very relevant to the global, cou-pled behavior for weak coupling. For higher coupling strengthsthat is no longer the case, since forcing is strong enough to bedominant. Here, we arrive at this conclusion from the variabil-ity statistical analysis, a di ff erent approach than the one usedpreviously: instead of using the individual membrane potentialsor the mean field, we analyse just the burst times. To get more information about the synchronization charac-teristics of the network, Fig. 4 depicts the raster plot of theburst starting times for the synchronization states of the net-work in the weak and strong coupling regimes. In these cases,after the transient time ( t ), we identify the cluster, using theprotocol defined in section 2.5, for the first bursting event. Allneurons inside the first cluster (green) are painted blue and alloutside are painted red. For subsequent events, this coloringscheme is mantained, which means that a blue (red) dot in Fig.4 represents a neuron that was (not) inside the cluster in the firstevent (green area). In this context, panel (a) depicts the networkfor T = ◦ C and ε = . ε = .
09 (a chaotic phase-synchronized state). In this case, there is very strong minglingbetween neurons (blue and red dots are mixed).A similar behavior is observed for the network with T = ◦ C. Panel (c) of Fig. 4 depicts the state of synchronizationin the weak coupling regime ( ε = . ε = . T = ◦ C where there is a mixture betweenthe blue and red dots in a very similar way to the behavior of anetwork with the same coupling and T = ◦ C. The di ff erencein temperature does not seem to a ff ect the behavior for the highcoupling regime.Given the definition of promiscuity, introduced in section3.2, the cases stated previously in Fig. 4, ε = . (a) (b) (c) t (ms) (d) N e u r o n Figure 4: Raster plots of the burst start times for di ff erent values of temperature T and coupling strength ε . Blue points correspond to neurons that were in thefirst cluster (defined with T =
1) of the graph and red squares to neurons thatwere not. The region painted green displays the interval within which neuronsare considered to be in the cluster, which is defined by the standard deviationof the bursting start times. Here, panels (a) and (b) correspond to the case of T = ◦ C where ε = . ε = .
09, respectively. Panels (c) and (d)represent the case of T = ◦ C for ε = . ε = .
09, respectively. Forboth cases, in the region of smaller coupling strength neurons tend to stay inthe cluster, while for the higher coupling neurons tend to mingle more. Time t shown is the time after 15 s of transient dynamics. els (a) and (c)) have low promiscuity (neurons are more phase-locked) and ε = . The previous observations, depicted in Fig. 4, can be mademore quantitative by identifying the T -clusters (as defined insection 2.5) for all clustering events and calculating CS / N (theaverage cluster size CS normalized by the network size N ) asa function of the number T of bursting events. This is donein Fig. 5, in which panel (a) represents the behavior for T = ◦ C and panel (b) for T = ◦ C. In both panels, blue linescorrespond to ε = . ε = .
09, representative of the strong coupling synchronization (seeFig. 3).The results of Fig. 5 show a similar behavior of CS / N for both temperatures considered: in the weak coupling regime(blue lines), the normalized size of the cluster tends to remainhigh even for high values of T , contrary to the strong couplingregime (orange lines), in which CS / N tends to vanishing val-ues as the number of bursting events is increased. These resultsindicate that neurons tend to stay in the cluster for longer times(i.e. promiscuity is lower) for the weak coupling cases. This7 . . . . . . C S / N (a) T . . . . . . C S / N (b) ε = ε = Figure 5: Cluster size CS normalized by the network size N as a function ofthe number of bursting events T for T = ◦ C (panel a) and T = ◦ C (panelb). The cluster is defined following section 2.5 and its size’s dependence on T refletcs how many neurons remain together for T bursting events. Blues lines( ε = . ε = . phenomenon is especially clear for T = ◦ C, in which theblue line has zero decline, indicating that all neurons in the T -cluster stay together through at least 15 subsequent events.It is worth remembering that networks in the strong cou-pling regime have the same or even higher values of the av-erage Kuramoto order parameter (cid:104) R (cid:105) , in comparison to oneswith weak coupling. Even still, they have higher promiscuity,showing that higher degree of synchronization does not implyin longer duration of phase-locking (i.e. lower promiscuity).The curves for strong coupling ( ε = .
09) are very similiar forthe two temperatures, reflecting the fact that dynamics in thiscase is really due to the coupling, not to the local dynamics(which is a ff ected by the temperatures) [31].We extend the previous analysis for all coupling strengthsin Fig. 6, which depicts the normalized cluster size CS / N (firstcolumn) and the average normalized number of neurons thatleave subsequent clusters L / CS (second column) as a functionof the coupling strength ε . Each row corresponds to di ff erentvalues of temperature: T = ◦ C , ◦ C , ◦ C from top to bot-tom, respectively. Furthermore, di ff erent values of the numberof bursting events T are considered, represented in di ff erentline colors: blue for T =
1, orange for T =
2, green for T = T =
10. The light gray areas are regions where it isnot possible to define cluster by using the approach defined inthe section 2.5.The di ff erence in behavior is much more stark for weak cou-pling strengths ( ε < . T = ◦ C (panels(a) and (b)) there is no defined cluster since the network is notsu ffi ciently synchronized. For T = ◦ C (panels (c)), there is ahigh number of neurons in the cluster, with a peak at the localmaximum of the Kuramoto Order Parameter, shown in Fig. 3. Increasing the number T of events, the normalized cluster size CS / N decreases, but not by much, indicating that most neuronsstay in the cluster for at least T =
10 events, which is reflectedin the low value of L / CS (panel (d)). Still for weak coupling,for T = ◦ C (panel (e)) almost the whole network is clustered( CS / N ≈ ε (cid:46) . CS / N as T is increased, which is corroborated by thevalue of L / CS , in panel (f) and is a quantitative representationof what was observed in the raster plot of Fig. 4. This region isprecisely the one with almost zero spatial variability ( CV s ) (cf.Fig. 3).For higher values of coupling, ( ε (cid:38) . T = ◦ C,the decrease with T is much more significant, reflecting theincrease of CV s . In this regime the behavior is very similarfor all temperatures: neurons group together around a cluster,but do not remain inside for more than, at most, a few burstingevents: that is, the duration of phase locking is very short.In short, comparing regions of weak coupling with ones ofhigh coupling, we note that the degree of synchronization (mea-sured by (cid:104) R (cid:105) ) may be similar, but the promiscuity is very dif-ferent: networks with weaker coupling are less promiscuous inthese cases.Furthermore, we note that, once again, promiscuity and vari-ability follow the same tendencies, as can be seen by comparing CV s in Fig. 3 with L / CS in Fig. 6. At last, a parameter-free approach is used to verify the re-sults about the promiscuity of a neural network. As definedin the section 2.6 the average temporal drift ( ∆ ) measures theaverage changes in the distance between burst times.Fig. 7 shows ∆ as a function of the coupling strength ε fordi ff erent temperature, as considered in main scenario (Fig. 3).For T = ◦ C (blue line), the drift starts high and decreaseswith an increase of coupling, indicating higher promiscuity forweaker coupling strengths. This was not observed with the clus-tering method, as depicted in Fig. 6 since synchronization inthis region is so weak that a cluster could not be defined.For T = ◦ C (orange line), the average temporal drift ∆ has a minimum in the region of weak coupling, coinciding withthe local maximum of (cid:104) R (cid:105) (Fig. 3). A coupling strength in-crease makes ∆ go up until a maximum is reached, correspond-ing to the local minimum of (cid:104) R (cid:105) . Further increase of couplingdecreases ∆ , but the final value is still higher than in the lowcoupling regime.Similar phenomena are observed for T = ◦ C (green line):there is lower promiscuity (lower ∆ ) for weak coupling ( ε (cid:46) . ∆ ) for high coupling ( ε (cid:38) . ∆ ) has a non-monotonicdependence on the coupling strength ( ε ) and this dependenceprofile changes with the temperature T . However, the relationbetween promiscuity ∆ and the spatial variability CV s , observedin Fig. 3 is maintained in all cases: the two quantities have thesame tendencies. This coincidence of tendencies between the8 . . . . . . C S / N (a) L / C S (b) . . . . . . C S / N (c) L / C S (d) .
00 0 .
02 0 .
04 0 .
06 0 . ε . . . . . . C S / N (e) .
00 0 .
02 0 .
04 0 .
06 0 . ε L / C S (f) T = T = T = T = Figure 6: Normalized cluster size ( CS / N ) (first column) and normalized number of neurons that leave the clusters ( L / CS ) (second column) as a function of couplingstrength ε for temperatures T = ◦ C, T = ◦ C, T = ◦ C (first, second and third rows, respectively). The line colors are varied according to the number T of events needed for neurons to be grouped in the cluster: blue lines are for T =
1, orange for T =
2, green for T =
5, and red for T =
10. Regions in whichno cluster could be defined are depicted in light gray. CS / N has a higher value, and smaller decrease as T increases, in the weak coupling regions, where spatialvariability is smaller. The lower decrease in these regions indicates that neurons remain in the cluster for longer times. A similar information can be observed for L / CS , in which the quantifier is lower for weak coupling, meaning that fewer neurons leave the cluster in this case. Also, the correlation between L / CS and thespatial variability is notable. All values are averages over 10 di ff erent initial conditions, with the errorbars corresponding to the standard deviation over them. drift ∆ and the spatial variability CV s calls into question if theyare equivalent by construction. This is in fact not the case: it ispossible for ∆ to change, but for CV s to remain constant.Very similar results were also observed in random networkswith the same number N of connections, created through thealgorithm proposed in [45].
4. Conclusions
The Inter-Burst Interval (
IBI ) variabilities were measuredin a network of modified-Hodgkin-Huxley neurons and shownto a ff ect the phase synchronization of the network. Firstly, forsmall coupling strengths, the degree of phase synchronizationwas observed to be inversely proportional to both temporal andspatial variabilities. This correlation is observed even for dif-ferent behaviors: from no synchronization to almost completesynchronization.Secondly, spatial variability was shown to be positively cor-related with neuronal promiscuity, the tendency of neurons tochange their relative phases in time. That is, for low variability,the neuronal promiscuity is low: neurons tend to stay togetherthroughout time; for high variability, promiscuity is high: neu-rons tend to mingle.It is striking that, in general, promiscuity is not directly re-lated to the degree of synchronization (it was observed only for small coupling strengths). As such, there can be highly syn-chronized networks with low promiscuity, but also ones withhigh. In fact, this high promiscuity is present in all investigatedtemperatures for high coupling strengths, which is counterintu-itive: in those cases, neurons that are strongly coupled are alsoless locked together.Since, according to the communication-through-coherencehypothesis [17], absence of phase coherence may prevent com-munication between neuronal groups, then we note that thisrelation between variability and phase-locking is important tokeep in mind in the study of neural communication.We argue that the relation between spatial variability andpromiscuity is actually causal, in that the former actually gen-erates the latter. This occurs simply because the higher thevariability, the higher the pool of available IBI s, so the higherthe chance neurons will choose di ff erent IBI and therefore theharder it is for neurons to stay locked together (ie the higherthe promiscuity). This e ff ect is statistical in nature, so it shouldbe present to some extent in networks with spatial variabilityin the inter-firing (spikes or bursts) intervals. In fact, it was al-ways observed for this neuronal model, with small-world andrandom topologies, as shown in the paper and in the appendix,respectively.It is worth noting that some regions of di ff erent couplingstrengths, but same degree of phase synchronization, exhibited9 .
00 0 .
02 0 .
04 0 .
06 0 . ε ∆ ( m s ) T = ◦ C T = ◦ C T = ◦ C Figure 7: Average drift ∆ as a function of the coupling strength ε . The blue, or-ange and green lines refer to T = ◦ C, T = ◦ C, and T = ◦ C, respectively.The average drift is a second measure of promiscuity and, as such, corroboratesthe results obtained via the clustering method in Fig. 6. It has the advantagesof being parameter-free and working for all degrees of synchronization in thenetwork. Results are averaged over 10 initial conditions, with errorbars equalto the standard deviation over them. di ff erent dynamical behaviors (reflected in their promiscuity),illustrating that the synchronization behavior of a network goesbeyond just the degree of synchronization.We conclude that spatial variability a ff ects not just the de-gree of phase synchronization, but also its form: it generatespromiscuity, influencing the duration of phase-locking. We sug-gest, then, that spatial variability is a useful quantity to considerin the study of neuronal communication. Acknowledgements
This study was financed in part by the Coordenac¸ ˜ao de Aper-feioamento de Pessoal de N´ıvel Superior - Brasil (CAPES) -Finance Code 001. The authors also acknowledge the supportof Conselho Nacional de Desenvolvimento Cient´ıfico e Tec-nol´ogico, CNPq - Brazil, grant number 302785 / ReferencesReferences [1] S. Den`eve, A. Alemi, R. Bourdoukan, The brain as an e ffi cient and robustadaptive learner, Neuron 94 (2017) 969–977.[2] T. B. DeMarse, D. A. Wagenaar, A. W. Blau, S. M. Potter, The neu-rally controlled animat: Biological brains acting with simulated bodies,Autonomous Robots 11 (2001) 305–310.[3] J. v. Neumann, The Computer and the Brain, Yale University Press, NewHaven, CT, USA, 1958.[4] P. Dayan, L. F. Abbott, Theoretical Neuroscience: Computational andMathematical Modeling of Neural Systems, The MIT Press, 2005.[5] M. N. Shadlen, W. T. Newsome, The variable discharge of cortical neu-rons: Implications for connectivity, computation, and information coding,Journal of Neuroscience 18 (1998) 3870–3896.[6] R. Nogueira, S. Lawrie, R. Moreno-Bote, Neuronal variability as a proxyfor network state, Trends in Neurosciences 41 (2018) 170 – 173. [7] M. P. Nawrot, C. Boucsein, V. R. Molina, A. Riehle, A. Aertsen, S. Rotter,Measurement of variability dynamics in cortical spike trains, Journal ofNeuroscience Methods 169 (2008) 374 – 390.[8] R. Vogels, W. Spileers, G. A. Orban, The response variability of striatecortical neurons in the behaving monkey, Experimental Brain Research77 (1989) 432–436.[9] M. M. Churchland, B. M. Yu, J. P. Cunningham, L. P. Sugrue, M. R.Cohen, G. S. Corrado, W. T. Newsome, A. M. Clark, P. Hosseini, B. B.Scott, D. C. Bradley, M. A. Smith, A. Kohn, J. A. Movshon, K. M. Arm-strong, T. Moore, S. W. Chang, L. H. Snyder, S. G. Lisberger, N. J. Priebe,I. M. Finn, D. Ferster, S. I. Ryu, G. Santhanam, M. Sahani, K. V. Shenoy,Stimulus onset quenches neural variability: a widespread cortical phe-nomenon, Nature Neuroscience 13 (2010) 369 EP –.[10] D. Tolhurst, J. Movshon, A. Dean, The statistical reliability of signalsin single neurons in cat and monkey visual cortex, Vision Research 23(1983) 775 – 785.[11] P. Kara, P. Reinagel, R. Reid, Low response variability in simultaneouslyrecorded retinal, thalamic, and cortical neurons, Neuron 27 (2000) 635 –646.[12] R. B. Stein, E. R. Gossen, K. E. Jones, Neuronal variability: noise or partof the signal?, Nature Reviews Neuroscience 6 (2005) 389 EP –.[13] J. Fell, N. Axmacher, The role of phase synchronization in memory pro-cesses, Nature Reviews Neuroscience 12 (2011) 105 EP –.[14] R. Gaillard, S. Dehaene, C. Adam, S. Clmenceau, D. Hasboun,M. Baulac, L. Cohen, L. Naccache, Converging intracranial markers ofconscious access, PLOS Biology 7 (2009) 1–21.[15] S. Dehaene, L. Charles, J.-R. King, S. Marti, Toward a computational the-ory of conscious processing, Current Opinion in Neurobiology 25 (2014)76 – 84.[16] W. Singer, Neuronal synchrony: A versatile code for the definition ofrelations?, Neuron 24 (1999) 49–65.[17] P. Fries, A mechanism for cognitive dynamics: neuronal communicationthrough neuronal coherence, Trends in Cognitive Sciences 9 (2005) 474– 480.[18] H. A. Braun, M. T. Huber, M. Dewald, K. Schafer, K. Voigt, Computersimulations of neuronal signal transduction: the role of nonlinear dynam-ics and noise, International Journal of Bifurcation and Chaos 8 (1998)881–889.[19] D. J. Watts, S. H. Strogatz, Collective dynamics of small-worldnetworks,Nature 393 (1998) 440–442.[20] D. S. Bassett, E. Bullmore, Small-world brain networks, The Neurosci-entist 12 (2006) 512–523.[21] V. Latora, M. Marchiori, E ffi cient behavior of small-world networks,Physical Review Letters 87 (2001) 198701.[22] J. Karbowski, Optimal wiring principle and plateaus in the degree ofseparation for cortical neurons, Physical Review Letters 86 (2001) 3674–3677.[23] C. C. Hilgetag, G. A. Burns, M. A. O’Neill, J. W. Scannell, M. P. Young,Anatomical connectivity defines the organization of clusters of corticalareas in the macaque monkey and the cat, Philosophical transactions ofthe Royal Society of London. Series B, Biological sciences 355 (2000)91–110.[24] M. D. Humphries, K. Gurney, T. J. Prescott, The brainstem reticular for-mation is a small-world, not scale-free, network, Proceedings. Biologicalsciences 273 (2006) 503–511.[25] H. A. Swadlow, A. G. Gusev, The impact of ’bursting’ thalamic impulsesat a neocortical synapse, Nature Neuroscience 4 (2001) 402–408.[26] J. E. Lisman, Bursts as a unit of neural information: making unreliablesynapses reliable, Trends in Neurosciences 20 (1997) 38 – 43.[27] P. Reinagel, D. Godwin, S. M. Sherman, C. Koch, Encoding of visualinformation by lgn bursts, Journal of Neurophysiology 81 (1999) 2558–2569.[28] J. Csicsvari, H. Hirase, A. Czurko, G. Buzski, Reliability and state de-pendence of pyramidal cellinterneuron synapses in the hippocampus: anensemble approach in the behaving rat, Neuron 21 (1998) 179 – 189.[29] R. C. Budzinski, B. R. R. Boaretto, T. L. Prado, S. R. Lopes, Detection ofnonstationary transition to synchronized states of a neural network usingrecurrence analyses, Physical Review E 96 (2017) 012320.[30] T. L. Prado, S. R. Lopes, C. A. S. Batista, J. Kurths, R. L. Viana, Synchro-nization of bursting hodgkin-huxley-type neurons in clustered networks,Physical Review E 90 (2014) 032818.
31] R. C. Budzinski, B. R. R. Boaretto, T. L. Prado, S. R. Lopes, Tempera-ture dependence of phase and spike synchronization of neural networks,Chaos, Solitons & Fractals 123 (2019) 35–42.[32] C. A. S. Batista, R. L. Viana, F. A. S. Ferrari, S. R. Lopes, A. M. Batista,J. C. P. Coninck, Control of bursting synchronization in networks ofhodgkin-huxley-type neurons with chemical synapses, Physical ReviewE 87 (2013) 042713.[33] B. R. R. Boaretto, R. C. Budzinski, T. L. Prado, J. Kurths, S. R. Lopes,Protocol for suppression of phase synchronization in hodgkin–huxley-type networks, Physica A: Statistical Mechanics and its Applications 528(2019) 121388.[34] K. Xu, J. P. Maidana, S. Castro, P. Orio, Synchronization transition inneuronal networks composed of chaotic or non-chaotic oscillators, Sci-entific reports 8 (2018) 8370.[35] B. R. R. Boaretto, R. C. Budzinski, T. L. Prado, S. R. Lopes, Suppres-sion of phase synchronization in scale-free neural networks using externalpulsed current protocols, Mathematical and Computational Applications24 (2019) 46.[36] B. R. R. Boaretto, R. C. Budzinski, T. L. Prado, J. Kurths, S. R. Lopes,Neuron dynamics variability and anomalous phase synchronization ofneural networks, Chaos: An Interdisciplinary Journal of Nonlinear Sci-ence 28 (2018) 106304.[37] Y. Kuramoto, Chemical oscillations, waves, and turbulence, volume 19,Springer Science & Business Media, 2012.[38] A. Destexhe, Z. F. Mainen, T. J. Sejnowski, An e ffi cient method for com-puting synaptic conductances based on a kinetic model of receptor bind-ing, Neural Computation 6 (1994) 14–18.[39] A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E.Shumaker, C. S. Woodward, SUNDIALS: Suite of nonlinear and dif-ferential / algebraic equation solvers, ACM Transactions on MathematicalSoftware (TOMS) 31 (2005) 363–396.[40] M. E. J. Newman, D. J. Watts, Scaling and percolation in the small-worldnetwork model, Physical Review E 60 (1999) 7332–7342.[41] A. Pikovsky, M. G. Rosenblum, J. Kurths, Synchronization, A UniversalConcept in Nonlinear Sciences, Cambridge University Press, Cambridge,2001.[42] M. V. Ivanchenko, G. V. Osipov, V. D. Shalfeev, J. Kurths, Phase synchro-nization in ensembles of bursting oscillators, Physical Review Letters 93(2004) 134101.[43] W. Softky, C. Koch, The highly irregular firing of cortical cells is in-consistent with temporal integration of random epsps, Journal of Neuro-science 13 (1993) 334–350.[44] Z. A. M. Stevens, Charles F., Input synchrony and the irregular firing ofcortical neurons, Nature Neuroscience 1 (1998).[45] P. Erd¨os, A. R´enyi, On the evolution of random graphs, Publications ofMathematics Institute Hungarian Academy Science 5 (1960) 17–61.algebraic equation solvers, ACM Transactions on MathematicalSoftware (TOMS) 31 (2005) 363–396.[40] M. E. J. Newman, D. J. Watts, Scaling and percolation in the small-worldnetwork model, Physical Review E 60 (1999) 7332–7342.[41] A. Pikovsky, M. G. Rosenblum, J. Kurths, Synchronization, A UniversalConcept in Nonlinear Sciences, Cambridge University Press, Cambridge,2001.[42] M. V. Ivanchenko, G. V. Osipov, V. D. Shalfeev, J. Kurths, Phase synchro-nization in ensembles of bursting oscillators, Physical Review Letters 93(2004) 134101.[43] W. Softky, C. Koch, The highly irregular firing of cortical cells is in-consistent with temporal integration of random epsps, Journal of Neuro-science 13 (1993) 334–350.[44] Z. A. M. Stevens, Charles F., Input synchrony and the irregular firing ofcortical neurons, Nature Neuroscience 1 (1998).[45] P. Erd¨os, A. R´enyi, On the evolution of random graphs, Publications ofMathematics Institute Hungarian Academy Science 5 (1960) 17–61.