Loss of synchronization in complex neuronal networks with delay
Judith Lehnert, Thomas Dahms, Philipp Hövel, Eckehard Schöll
aa r X i v : . [ c ond - m a t . d i s - nn ] N ov epl draft Loss of synchronization in complex neuronal networks withdelay
Judith Lehnert , Thomas Dahms , Philipp Hövel , and Eckehard Schöll (a) Institut für Theoretische Physik, Technische Universität Berlin, Hardenbergstraße 36, 10623 Berlin, Germany Bernstein Center for Computational Neuroscience Berlin, Philippstraße 13, Haus 2, 10115 Berlin, Germany
PACS – Synchronization, nonlinear dynamics
PACS – Neural network
PACS – Complex system
Abstract – We investigate the stability of synchronization in networks of delay-coupled excitableneural oscillators. On the basis of the master stability function formalism, we demonstrate thatsynchronization is always stable for excitatory coupling independently of the delay and couplingstrength. Superimposing inhibitory links randomly on top of a regular ring of excitatory coupling,which yields a small-world-like network topology, we find a phase transition to desynchronizationas the probability of inhibitory links exceeds a critical value. We explore the scaling of the criticalvalue in dependence on network properties. Compared to random networks, we find that small-world topologies are more susceptible to desynchronization via inhibition.
Introduction. –
Studies of complex networks havesparked tremendous scientific activities in many researchareas and the analysis of network topologies in real-worldsystems has become a field of large interest. For instance,there is evidence that neuronal networks on the level of sin-gle neurons coupled through synapses or gap junctions, aswell as on the level of cortex areas and their pathways ex-hibit the small-world (SW) properties [1,2]. The high clus-tering coefficient of the SW networks enhances local com-munication efficiency, while the small shortest path lengthenables efficient global communication [3]. Thus, the SWarchitecture is optimal for processing and transmission ofsignals within and between brain areas. However, the syn-chronizability of small-world networks depends in a deli-cate way upon the network topology [4]. Next to thisstructural aspect, inhibition plays a prominent role inmany neural processes [5]. Without an inhibitory mecha-nism, excitation in a compound system would not decay,but spread through the whole network, finally leading topersistent spiking of all neurons. Thus, encoding and pro-cessing of information would be impossible.In this Letter, we combine both fundamental aspects– inhibition and SW property – in order to emphasizethe important interplay of excitation and inhibition incomplex networks. We start with a regular ring network (a)
E-mail: [email protected] that consists of purely excitatory links with delay. Thus,it exhibits strong and stable synchronization. Depend-ing on the initial conditions, both isochronous and clus-ter synchronization are possible implying multistability.Additional inhibitory connections, which we include in aSW-like manner [1, 6], result in a loss of synchronization.A similar transition was reported for phase oscillators inRef. [7] for unidirectional rings, but the effect of inhibi-tion upon excitable systems could not be treated by thatmodel. For the node dynamics, we consider a genericmodel to demonstrate the fundamental relevance and im-portance of our findings in the field of neuroscience. In thisarea, synchronization can be related to cognitive capaci-ties [8] as well as to pathological conditions, e.g., epilepsy[9]. A better understanding of the loss of synchronizationwill eventually lead to future therapeutic treatments [10].Throughout this Letter, we consider a network of N delay-coupled FitzHugh-Nagumo (FHN) oscillators. TheFHN system describes neuronal dynamics by a two-variable model [11]. Because of its simplicity it can beconsidered as a paradigmatic model of excitable systems,which also occur in several other natural contexts rangingfrom cardiovascular tissues to the climate system [12, 13].p-1. Lehnert et al. Fig. 1: (Color online) Dynamics in the synchronization man-ifold in dependence on the coupling strength C and delay τ .The gray scale (color code) indicates the period of spiking os-cillations T , the black region corresponds to fixed-point dy-namics. Left and right insets show time series of the activator u s for ( C = 0 . , τ = 1 ) and ( C = 1 . , τ = 1 . ), respectively.Parameters: a = 1 . , ǫ = 0 . . Here, the network dynamics is described by ǫ ˙ u i = u i − u i − v i + C N X j =1 G ij [ u j ( t − τ ) − u i ] (1a) ˙ v i = u i + a, (1b)where u i and v i denote the activator and inhibitor vari-ables of the nodes i = 1 , . . . , N , respectively. The param-eter a determines the threshold of excitability. A singleFHN oscillator is excitable for a > and exhibits self-sustained periodic firing beyond the Hopf bifurcation at a = 1 . Here, we will focus on the excitable regime with a = 1 . . The time-scale parameter ǫ is chosen as ǫ = 0 . . C is the coupling strength. G = { G ij } , i, j = 1 , . . . , N ,denotes the coupling matrix that determines the topol-ogy of the network. In the following we will assume unityrow sum of G . This ensures that each neuron receives thesame input if the network is synchronized. The delay time τ takes into account the finite propagation speed of an ac-tion potential. We investigate complete synchronizationwith ( u i ( t ) , v i ( t )) = ( u s ( t ) , v s ( t )) ≡ x s ( t ) for i = 1 , . . . , N ,which is also known as zero-lag or isochronous synchro-nization. This state is a solution of qs. (1) and reducesthe system’s dynamics to ˙ x s = F ( x s ) + C H [ x s ( t − τ ) − x s ( t )] (2)with F ( x ) = (cid:16) [ u − u / − v ] /ǫu + a (cid:17) and the matrix H = (cid:0) /ǫ
00 0 (cid:1) .The N − constraints of complete synchronization de-fine a two-dimensional synchronization manifold (SM) inthe N -dimensional phase space.As we operate in the excitable regime, the dynamics onthis SM, in particular the period, will depend on the choiceof the coupling parameters C and τ as depicted in fig. 1. The grayscale (color code) corresponds to the period T ofthe oscillations on the SM, which we find to follow T = τ + δ with δ ≪ τ accounting for a short activation time [14].For small coupling strength C the incoming signal is notsufficient to trigger oscillations (black region). For smalldelay times τ consecutive spikes run into the refractoryphase of the previous one, which prevents oscillations aswell. From here on we consider C and τ sufficiently largesuch that the coupling induces oscillations. Stability analysis. –
In the following we address thequestion whether the oscillatory solution on this mani-fold is transversely stable. The master stability function(MSF) [15] allows us to quantify this transversal stability.It can be calculated as largest Lyapunov exponent Λ fromeq. (1) linearized around eq. (2): ˙ ζ ( t ) = [ D F ( x s ( t )) − C H ] ζ ( t ) + ( α + iβ ) H ζ ( t − τ ) . (3)Here, D F denotes the Jacobian of F . The idea of the MSFis to calculate the stability of a synchronized solution foran arbitrary topology matrix G . For this purpose, theparameter α + iβ represents a continuous parametrizationof { Cν i } , where ν i , i = 1 , . . . , N , are the eigenvalues of G . In the same sense, the vector ζ is a generalization ofthe variational vectors transformed to the correspondingeigensystem. In the ( α , β )-plane the MSF typically givesrise to regions with negative Λ . If all rescaled transversaleigenvalues { Cν i } of a given network are located withinthis stable region, perturbations from the SM will decayexponentially and the synchronized dynamics will be sta-ble. Due to the unity row sum condition, G will alwayshave one eigenvalue ν = 1 . This longitudinal eigenvalueis associated with perturbations within the SM and isnot relevant for the stability of synchronization. Λ( Cν ) determines the type of dynamics on the SM. For peri-odic dynamics in the SM, as in the present case, we have Λ( Cν ) = 0 .Figure 2 depicts the MSF for the network of FHN oscil-lators given by eqs. (1). Dark (blue) colors mark the stableregion. As an illustration the rescaled eigenvalues { Cν i } of a bidirectionally coupled ring ( N = 8 ) are shown as redsymbols. The corresponding coupling matrix is given by G i,i +1 mod N = G i,i − N = 1 ( i = 1 , . . . , N ) and zerootherwise. The rescaled longitudinal eigenvalue Cν = C is depicted by a black (red) square. All rescaled transver-sal eigenvalues (black (red) circles) lie inside the stableregion indicating that the synchronization of the bidirec-tionally coupled ring is stable. Shape of stability region. –
The MSF must be cal-culated for each combination of C and τ . Although dif-ferent C and τ lead to quantitatively different Lyapunovexponents Λ , the shape of the stable regions remains qual-itatively very similar. In particular, it is in very good ap-proximation given by the circle S ((0 , , C ) with center atthe origin and radius C (dotted circle in fig. 2) indepen-dent of the specific values of C and τ . The rotational sym-metry has recently been proved generally for large τ [16].p-2oss of synchronization in complex neuronal networks with delay -0.3 0 0.3 α -0.3 0 0.3 β -0.5 0 0.5 Λ -0.3 0 0.3 α -0.3 0 0.3 β -0.5 0 0.5 Λ -0.3 0 0.3 α -0.3 0 0.3 β -0.5 0 0.5 Λ -0.3 0 0.3 α -0.3 0 0.3 β -0.3 0 0.3 α -0.3 0 0.3 β -0.3 0 0.3 α -0.3 0 0.3 β -0.3 0 0.3 α -0.3 0 0.3 β -0.3 0 0.3 α -0.3 0 0.3 β -0.3 0 0.3 α -0.3 0 0.3 β -0.3 0 0.3 α -0.3 0 0.3 β -0.3 0 0.3 α -0.3 0 0.3 β (b)(c)(a) Fig. 2: (Color online) (a) Master stability function for anetwork of FHN systems given by eqs. (1). Dotted curve: S ((0 , , C ) . Red circles (square): Rescaled transversal (longi-tudinal) eigenvalues Cν i of a bidirectionally coupled ring with N = 8 nodes. Parameters: a = 1 . , ǫ = 0 . , C = 0 . , τ = 1 .(b) Scheme of a bidirectional regular network ( N = 20 , k = 2 ),and (c) a random network ( N = 20 , fixed number of links kN ) with excitatory coupling (gray (green) arrows) on whichinhibitory links (black (red) arrows) are superimposed. Only for small τ and C the stable region is slightly largerthan the circle and shows a bulge around α = − C , β = 0 [17]. The positive α -axis is always intersected at α = C ,which corresponds to Λ( Cν ) = 0 as discussed above. Forany choice of τ and C that leads to periodic dynamics onthe SM, the circle S ((0 , , C ) serves as a lower bound forthe stability boundary. See the appendix for an analyticderivation of this circle S ((0 , , C ) in the limit of largecoupling strength and as a lower bound for all couplingstrengths. We conclude that the stability of the synchro-nized periodic dynamics, if such a solution exists, dependsonly on the topology and neither on the coupling strengthnor on the delay time. Excitatory coupling. –
For excitatory coupling, i.e., G ij ≥ , all eigenvalues of G are located inside the sta-ble region. Using Gershgorin’s circle theorem [18], whichgives an upper bound of the eigenvalues, and the constantrow sum assumption, all Gershgorin circles ( i = 1 , . . . , N ),centered at G ii with radius P j = i G ij = 1 − G ii becauseof the unity row sum, lie inside the unit circle. Thus allrescaled eigenvalues { Cν i } , i = 1 , . . . , N , are located in-side S ((0 , , C ) , i.e., inside the stable region. Networkswith purely excitatory coupling will always exhibit stablesynchronization. Inhibitory coupling. –
As a consequence of thisresult, desynchronization can only be achieved by intro-ducing negative entries in the coupling matrix G , i.e.,inhibitory coupling between neurons. This inhibition isa crucial feature in neural processes, e.g., to overcomeunwanted synchronization associated with pathologicalstates.Particularly, we consider the following variation of theWatts-Strogatz SW network [1, 6]: (i) Start with a one-dimensional ring of N nodes, where every node is con- f p k k=6 k=12 k=18 k=24 k=30 Fig. 3: (Color online) Fraction of desynchronized networks f vs the probability of additional inhibitory links p for N = 100 . k varies from 6 to 30. Thin black curve: Example fit to f ( p ) ( p c = 0 . , b = 186 ) for k = 24 . Number of realizations:500 for each value of k . Parameters as in fig. 1. nected by excitatory links to its k neighbors on eitherside. (ii) For each of the kN links of the network addan inhibitory link with probability p connecting two ran-domly chosen nodes. (iii) Do not allow self-coupling ormore than one link between any pair of nodes. (iv) Nor-malize the entries of the coupling matrix G by dividingeach row by the absolute value of its row sum. In the casethat the row sum of the i th row is negative we set G ii = 2 to ensure unity row sum. Figure 2(b) illustrates such a SWnetwork for N = 20 and k = 2 , where gray (green) andblack (red) arrows indicate excitatory and inhibitory cou-pling, respectively. For each realization of such a network,we determine the stability of synchronization by checkingwhether the full eigenspectrum Cν i of the coupling ma-trix is contained in the stable region S ((0 , , C ) . Herebywe compute the fraction f of desynchronized networks.Figure 3 shows f as a function of p for different couplingranges k . This Figure is virtually identical for all delaytimes, that is, for all parameters within the color shadedarea of Fig.1. To obtain this Figure we made use of thecircular shape of the stable region of the master stabilityfunction. Only for very small delays or coupling strength,the stable region is slightly larger and thus the shape ofthe curves shown in Figure 3 might be shifted slightly tolarger values of p .For fixed k a steep transition between synchronizationand desynchronization takes place as p approaches a crit-ical value p c . This critical value p c and the steepness b of the transition can be fitted with a sigmoidal func-tion f ( p ) = 1 / [ e − b ( p − p c ) + 1] . Figure 4(a) depicts thecritical probability p c for f ( p c ) = 0 . in dependence on k/N for different network sizes. It can be seen thatfor SW networks p c follows a linear relation p c ( k/N ) =1 . k/N − . independently of the network size N . Fig-ure 4(c) shows the steepness b as a function of the networksize demonstrating that the transition becomes increas-ingly sharp as N increases. This indicates a first-ordernonequililibrium phase transition in the thermodynamiclimit [23].To verify whether this phase transition and especially itsindependence of the network size is common in networkswith inhibitory links or unique to the SW structure, wep-3. Lehnert et al. p c k/N(a) 00.250.65 10 20 40 50 60 p c k(b)05001500 0 250 750 1000 b N(c) b N p c 〈 Ι 〉 N(d)
Fig. 4: (Color online) Critical value p c for different networksizes, (a) in dependence on k/N , (b) in dependence on k : SWnetworks N = 60 (dark (blue) circles), N = 100 (dark (purple)squares), and N = 500 (turquoise triangles). Random networks N = 60 (black (red) crosses), N = 100 (light (orange) circles),and N = 500 (lightgray (yellow) squares). (c) Steepness b forSW (circles) and random (squares) networks vs. N for constant k/N = 0 . . Inset: Blow-up for random networks. (d) p c vs. N for k = 10 for a SW (black (red) filled circles) and a random(black (red) empty squares) network. Number of inhibitorylinks h I i vs. N for constant k for a SW (gray (green) emptycircles) and a random (gray (green) filled squares) network.Number of realizations: 500. construct a different network for comparison: The regu-lar excitatory network is replaced by a random networkwith fixed number of excitatory links kN equivalently tothe regular network used before, as shown in fig. 2(c).We only consider realizations where this underlying ex-citatory network is fully connected. The construction ofthe inhibitory links then follows steps (ii-iv) as above. Wefind that a phase transition to desynchronization still oc-curs with critical probabilities of inhibitory links p c as de-picted in fig. 4(a) by black (red) crosses, gray (orange)circles, and lightgray (yellow) squares for N = 60 , ,and , respectively. We observe, however, that the val-ues of p c are higher, i.e., the random network can tol-erate more inhibitory links than the SW network beforedesynchronizing. Furthermore, the function p c ( k/N ) is nolonger independent of the network size N . Instead, p c isa function of log( k ) as can be seen in fig. 4(b), where p c is plotted in dependence on k for the different net-work sizes N . Figure 4(d) depicts p c in dependence on N for constant k = 10 for a random (black (red) emptysquares) and for a SW network (red (gray) circles). Forrandom networks, p c is independent of N for sufficientlylarge N , while for SW networks it approaches zero. Re-call that p c is the mean value of the ratio of inhibitory toexcitatory links. Thus, we conclude that in SW networkswith increasing network size N but same local structure(constant k ) an infinitesimally small ratio of inhibition toexcitation is needed to prevent synchronization, while in a random network even for very large networks only anon-vanishing ratio impairs synchronization. We find in aSW network with constant k that the mean value of thenumber of inhibitory links h I i := p c kN causing desyn-chronization scales as h I i = 1 . k − . kN for small N and approaches zero for large N (see fig. 4(d) green(gray) empty circles). In contrast, in a random network h I i is proportional to N (see fig. 4(d) green (gray) filledsquares), i.e., an increasing number of inhibitory links isneeded. This difference to SW networks can be under-stood in an intuitive way: In a SW, any added inhibitorylink is part of a shortest path for many pairs of nodes,as it shortens the mean path length considerably with re-spect to the underlying regular ring. In the random net-work, however, where the mean path length is relativelylow even without added shortcuts, only few node pairs willgain shorter paths by adding inhibitory links. Consider-ing the dynamics on a network, perturbations from thesynchronized state spread along the shortest paths first,changing the response of the receiving node, and informa-tion flow along longer paths will reach the receiving nodeonly at a later time and will not influence the change ofthe initial response. In conclusion, if a large fraction ofthe inhibitory links is part of the shortest paths - like inthe small-world topology superimposed to a regular ring -these inhibitory shortcuts become dominant. Conclusions. –
We have shown how the interplay ofexcitatory and inhibitory couplings leads to desynchro-nization in networks of neural oscillators. The desynchro-nization is achieved via a phase transition from a com-pletely synchronized state. This can be seen as a first steptowards an understanding of the robustness of differentstates of synchrony, e.g., cluster synchronization, in arbi-trary networks with weighted links or distributed delays.Note that for appropriate network topologies the frame-work of the MSF presented above can indeed be extendedto cluster synchronization where the oscillators synchro-nize in M clusters with a constant phase lag π/M be-tween subsequent clusters [19]. The corresponding SM is M dimensional. Hence, M longitudinal eigenvalues ex-ist. The MSF, however, is again very well approximatedby the circle S ((0 , , C ) and thus, we observe multistabil-ity between zero-lag and cluster synchronization.Excitable systems can be classified into type-I and type-II excitability [13, 20]. In addition to the generic type-IIFitzHugh-Nagumo model used in this paper, we have con-sidered the normal form of a saddle-node bifurcation onan invariant circle (SNIC) as a generic model of type-Iexcitability [21]. For sufficiently large delay times andcoupling strength the MSF is again given by the cir-cle S ((0 , , C ) implying that the previously obtained re-sults persist. In particular, the same phase transition oc-curs. This indicates that the phenomena observed hereare generic for any excitable system. Appendix: Analytic approximation of the stabil-ity region. –
The numerical calculation of the masterp-4oss of synchronization in complex neuronal networks with delaystability function has shown that S ((0 , , C ) is a lowerbound and a very good approximation of the stable regionfor all τ and C . As τ and C increase the approxima-tion becomes even better. A Taylor expansion as donein Ref. [22] for the investigation of time-delayed feedbackcontrol of an unstable periodic orbit gives analytic insightin the problem. This analysis is very general and doesnot use the specific form of the local dynamics in termsof the FHN model. It only assumes that the synchro-nized dynamics is oscillatory with period T . Using a Flo-quet ansatz ζ = e (Λ+ i Ω) t Q ( t ) with the periodic function Q ( t ) = Q ( t + T ) in Eq. (3) yields (Λ + i Ω) Q ( t ) + ˙ Q ( t ) (4) = ( D F − C H ) Q ( t ) + ( α + iβ ) e − (Λ+ i Ω) τ HQ ( t − τ ) . Λ + i Ω is the Floquet exponent, whose real part coincideswith the Lyapunov exponent in the case of a periodic orbit.Assume T = τ . In the case of the FHN system this is anapproximation since the period of the oscillations differsby a small activation time δ ≪ τ from the delay time τ following T = τ + δ . Then Q ( t − τ ) can be substituted by Q ( t ) : (Λ + i Ω) Q ( t ) + ˙ Q ( t ) (5) = ( D F ) Q ( t ) + [ − C + ( α + iβ ) e − (Λ+ i Ω) τ ] | {z } κ HQ ( t ) . We expand the solution Γ( κ ) = Λ + i Ω of the eigenvalueproblem defined by Eq. (5) in a Taylor approximation: Γ( κ ) = Γ(0) + Γ ′ (0) κ + O ( κ ) . (6)Using Γ(0) ≡ λ + iω and Γ ′ (0) ≡ χ ′ + iχ ′′ we obtain Λ + i Ω = λ + iω + ( χ ′ + iχ ′′ )[ − C + ( α + iβ ) e − (Λ+ i Ω) τ ] . (7)Note that κ = 0 if ( α, β ) = ( C, corresponding to thedynamics within the synchronization manifold. Thus thefirst term in the Taylor approximation corresponds to theGoldstone mode, i.e., λ + iω = 0 for ( α, β ) = ( C, . Equa-tion (7) then becomes Λ + i Ω = χ ′ [ − C + ( α + iβ ) e − (Λ+ i Ω) τ ] . (8)Here, we assume χ ′′ = 0 . Separating Eq. (8) into real andimaginary part leaves us with Λ = χ ′ {− C + e − Λ τ [ α cos(Ω τ ) + β sin(Ω τ )] } , Ω = χ ′ e − Λ τ [ − α sin(Ω τ ) + β cos(Ω τ )] . (9)Equation (9) can be solved numerically yielding the cir-cular stability region. On the border of the stability thereal part of the Floquet exponent vanishes. Using Λ = 0 in Eq. (9) yields after algebraic manipulations: α b = − Ω sin(Ω τ ) χ ′ + C cos(Ω τ ) ,β b = Ω cos(Ω τ ) χ ′ + C sin(Ω τ ) , (10) where α b and β b denote the values of α and β , respectively,on the bounder of stability. Finally we obtain α b + β b = C + Ω χ ′ . (11)Obviously α b + β b > C holds, demonstrating that S ((0 , , C ) is a lower bound for the stable region. Forlarge C the term C on the right hand side dominates.Thus, the boundary of stability is very well approximatedby S ((0 , , C ) for large coupling strength. ∗ ∗ ∗ This work was supported by DFG in the framework ofSFB 910. PH acknowledges support by the BMBF underthe grant no. 01GQ1001B (
Förderkennzeichen ). REFERENCES[1] D. J. Watts and S. H. Strogatz, Nature , 440 (1998).[2] O. Sporns, G. Tononi, and G. M. Edelman, Cereb. Cortex , 127 (2000); O. Shefi et al. , Phys. Rev. E , 021905(2002); O. Sporns, Biosystems , 55 (2006); C. J. Honey et al. , Proc. Natl. Acad. Sci. U.S.A. , 10240 (2007); O.Sporns, C. J. Honey, and R. Kötter, PLoS ONE , e1049(2007).[3] V. Latora and M. Marchiori, Phys. Rev. Lett. , 198701(2001).[4] T. Nishikawa et al. , Phys. Rev. Lett. , 014101 (2003).[5] B. Haider et al. , J. Neurosci. , 4535 (2006).[6] R. Monasson, Eur. Phys. J. B , 555 (1999); M. E. J.Newman and D. J. Watts, Phys. Lett. A , 341 (1999).[7] R. Tönjes, N. Masuda, and H. Kori, Chaos , 033108(2010).IS03[8] W. Singer, Neuron , 49 (1999).[9] P. J. Uhlhaas and W. Singer, Neuron , 155 (2006).[10] C. Hauptmann and P. A. Tass, Biosystems , 173 (2007).[11] R. FitzHugh, Biophys. J. , 445 (1961); J. Nagumo, S.Arimoto, and S. Yoshizawa., Proc. IRE , 2061 (1962).[12] J. D. Murray, Mathematical Biology , Vol. 19 of
Biomathe-matics Texts , 2nd ed. (Springer, Berlin Heidelberg, 1993);A. S. Mikhailov,
Foundations of Synergetics Vol. I , 2 ed.(Springer, Berlin, 1994); J. P. Keener and J. Sneyd,
Math-ematical physiology (Springer, New York, Berlin, 1998);C. Koch,
Biophysics of Computation: Information Pro-cessing in Single Neurons (Oxford University Press, NewYork, 1999); H. J. Wünsche et al. , Phys. Rev. Lett. ,023901 (2001); A. Ganopolski and S. Rahmstorf, Phys.Rev. Lett. , 038501 (2002).[13] E. M. Izhikevich, Int. J. Bifurc. Chaos , 1171 (2000).[14] E. Schöll et al. , Phil. Trans. R. Soc. A , 1079 (2009);M. A. Dahlem et al. , Int. J. Bifur. Chaos , 745 (2009).[15] L. M. Pecora and T. L. Carroll, Phys. Rev. Lett. , 2109(1998).[16] V. Flunkert et al. , Phys. Rev. Lett. , 254101 (2010).[17] E.g., for C = 0 . and τ = 1 the radius is 0.303 in thedirection of the negative α -axis. p-5. Lehnert et al. [18] M. G. Earl and S. H. Strogatz, Phys. Rev. E , 036204(2003); C.-U. Choe et al. , Phys. Rev. E , 025205(R)(2010).[19] F. Sorrentino and E. Ott, Phys. Rev. E , 056114 (2007);I. Kanter et al. , Europhys. Lett. , 66001 (2011).[20] A. L. Hodgkin, J. Physiol. , 165 (1948).[21] G. Hu et al. , Phys. Rev. Lett. , 807 (1993).[22] W. Just et al. , Phys. Rev. Lett. , 203 (1997).[23] E. Schöll, Nonequilibrium Phase Transitions in Semicon-ductors (Springer, Berlin, 1987).(Springer, Berlin, 1987).