Desynchronization transitions in adaptive networks
DDesynchronization transitions in adaptive networks
Rico Berner , , ∗ Simon Vock , Eckehard Sch¨oll , , and Serhiy Yanchuk Institut f¨ur Theoretische Physik, Technische Universit¨at Berlin, Hardenbergstr. 36, 10623 Berlin, Germany Institut f¨ur Mathematik, Technische Universit¨at Berlin,Straße des 17. Juni 136, 10623 Berlin, Germany and Bernstein Center for Computational Neuroscience Berlin,Humboldt-Universit¨at, Philippstraße 13, 10115 Berlin, Germany (Dated: August 26, 2020)Adaptive networks change their connectivity with time, depending on their dynamical state.While synchronization in structurally static networks has been studied extensively, this problemis much more challenging for adaptive networks. In this Letter, we develop the master stabilityapproach for a large class of adaptive networks. This approach allows for reducing the synchro-nization problem for adaptive networks to a low-dimensional system, by decoupling topological anddynamical properties. We show how the interplay between adaptivity and network structure givesrise to the formation of stability islands. Moreover, we report a desynchronization transition andthe emergence of complex partial synchronization patterns induced by an increasing overall cou-pling strength. We illustrate our findings using adaptive networks of coupled phase oscillators andFitzHugh-Nagumo neurons with synaptic plasticity.
In nature and technology, complex networks serve asa ubiquitous paradigm with a broad range of appli-cations from physics, chemistry, biology, neuroscience,socio-economic and other systems [1]. Dynamical net-works are composed of interacting dynamical units, suchas, e.g., neurons or lasers. Collective behavior in dy-namical networks has attracted much attention over thelast decades. Depending on the network and the spe-cific dynamical system, various synchronization patternsof increasing complexity were explored [2–5]. Even insimple models of coupled oscillators, patterns such ascomplete synchronization [6], cluster synchronization [7–11], and various forms of partial synchronization havebeen found, such as frequency clusters [3], solitary [13]or chimera states [14–22]. In brain networks, particu-larly, synchronization is believed to play a crucial role:for instance, under normal conditions in the context ofcognition and learning [23, 24], and under pathologicalconditions, such as Parkinson’s disease [25], epilepsy [26–30], tinnitus [31, 32], schizophrenia, to name a few [33].Also in power grid networks, synchronization is essentialfor the stable operation [34–37].The powerful methodology of the master stability func-tion [38] has been a milestone for the analysis of syn-chronization phenomena. This method allows for sep-arating dynamical from structural features for a givendynamical network. It drastically simplifies the problemby reducing the dimension and unifying the synchroniza-tion study for different networks. Since its introduction,the master stability approach has been extended and re-fined for multilayer [39], multiplex [40, 41] and hyper-networks [42, 43]; to account for single and distributeddelays [44–49]; and to describe the stability of clusteredstates [50–53]. The master stability function has beenused to understand effects in temporal [54, 55] as well asadaptive networks [56] within a static formalism. Beyondthe local stability described by the master stability func- tion, Belykh et. al. have developed the connection graphstability method to provide analytic bounds for the globalasymptotic stability of synchronized states [57–60]. De-spite the apparent vivid interest in the stability featuresof synchronous states on complex networks, only little isknown about the effects induced by an adaptive networkstructure. This lack of knowledge is even more surprisingregarding how important adaptive networks are for themodeling of real-world systems.Adaptive networks are commonly used models forsynaptic plasticity [61–66] which determines learning,memory, and development in neural circuits. More-over, adaptive networks have been reported for chemi-cal [67, 68], epidemic [69], biological [70], transport [71],and social systems [72, 73]. A paradigmatic exampleof adaptively coupled phase oscillators has recently at-tracted much attention [1–5, 41, 74, 75, 77, 81], and it ap-pears to be useful for predicting and describing phenom-ena in more realistic and detailed models [7, 82, 83, 85].Systems of phase oscillators are important for under-standing synchronization phenomena in a wide range ofapplications [86–88].In this Letter, we report on a surprising desynchroniza-tion transition induced by an adaptive network structure.We find various parameter regimes of partial synchro-nization during the transition from the synchronized toan incoherent state. The partial synchronization phe-nomena include multi-frequency-cluster and chimera-likestates. By going beyond the static network paradigm,we develop a master stability approach for networks withadaptive coupling. We show how the adaptivity of thenetwork gives rise to the emergence of stability islandsin the master stability function that result in the desyn-chronization transition. With this, we establish a generalframework to study those transitions for a wide range ofdynamical systems. In order to provide analytic insights,we use the generalized Kuramoto-Sakaguchi system on a r X i v : . [ n li n . AO ] A ug an adaptive and complex network. Finally, we show thatour findings also hold for a more realistic neuronal set-up of coupled FitzHugh-Nagumo neurons with synapticplasticity.We consider the following general class of N adaptivelycoupled systems [1–5, 41, 74, 75, 77, 89]˙ x i = f ( x i ) − σ N (cid:88) j =1 a ij κ ij g ( x i , x j ) , (1)˙ κ ij = − (cid:15) ( κ ij + a ij h ( x i − x j )) , (2)where x i ∈ R d , i = 1 , . . . , N , is the d -dimensional dy-namical variable of the i th node, f ( x i ) describes the lo-cal dynamics of each node, and g ( x i , x j ) is the couplingfunction. The coupling is weighted by scalar variables κ ij which are adapted dynamically according to Eq. (2) withthe nonlinear adaptation function h ( x i − x j ). We assumethat the adaptation depends on the difference of the cor-responding dynamical variables, similar to the neuronalspike timing-dependent plasticity [62, 63, 90, 91]. Thebase connectivity structure is given by the matrix ele-ments a ij ∈ { , } of the N × N adjacency matrix A whichpossesses a constant row sum r , i.e., r = (cid:80) Nj =1 a ij for all i = 1 , . . . , N . The assumption of the constant row sumis necessary to allow for synchronization. The Laplacianmatrix is L = r I N − A where I N is the N -dimensionalidentity matrix. The eigenvalues of L are called Lapla-cian eigenvalues of the network. The parameter σ > (cid:15) > (cid:15) is small.Complete synchronization is defined by the N − x = x = · · · = x N . Denoting the synchroniza-tion state by x i ( t ) = s ( t ) and κ ij = κ sij , we obtain fromEqs. (1)–(2) the following equations for s ( t ) and κ sij ˙ s = f ( s ) + σrh (0) g ( s , s ) , (3) κ sij = − a ij h (0) . (4)In particular, we see that s ( t ) satisfies the dynamicalequation (3), and κ sij are either − h (0) or zero, if the cor-responding link in the base connectivity structure exists( a ij = 1) or not ( a ij = 0), respectively.In order to describe the local stability of the syn-chronous state, we introduce the variations ξ i = x i − s and χ ij = κ ij − κ sij . The linearized equations for thesevariations read˙ ξ i = D f ( s ) ξ i − σg ( s , s ) N (cid:88) j =1 a ij χ ij (5)+ σh (0) N (cid:88) j =1 a ij (D g ( s , s ) ξ i + D g ( s , s ) ξ j ) , ˙ χ ij = − (cid:15) ( χ ij + a ij D h (0)( ξ i − ξ j )) , (6) where D f and D h are the Jacobians ( d × d matrix and1 × d matrix, respectively), and D g and D g are the Ja-cobians with respect to the first and the second variable,respectively.The system (5)–(6) is used to calculate the Lyapunovexponents of the synchronous state; it possesses veryhigh dimension N + N d . However, one can introducea new coordinate frame which separates an N ( d + 1)-dimensional master from an N ( N − N ( d + 1)-dimensional mas-ter system into blocks of d + 1 dimensions. Hence, thedynamics in each block is described by the new coordi-nates ζ and κ which are d - and one-dimensional dynam-ical variables, respectively. Our analysis shows that thecoupling structure enters just as a complex parameter µ ,the network’s Laplacian eigenvalue. For all details andthe proof of the master stability function, we refer to theSupplemental Material [105].As a result, the stability problem is reduced to thelargest Lyapunov exponent Λ( µ ), depending on a com-plex parameter µ , for the following system˙ ζ = (cid:18) D f ( s ) + σrh (0) (cid:2) D g ( s , s )+ (1 − µr )D g ( s , s ) (cid:3)(cid:19) ζ − σg ( s , s ) κ, (7)˙ κ = − (cid:15) ( µ D h (0) ζ + κ ) , . (8)The function Λ( µ ) is called master stability function.Note that the first bracketed term in ζ of (7) resemblesthe master stability approach for static networks, which,in this case, is equipped by an additional interaction rep-resenting the adaptation.To obtain analytic insights into the stability featuresof synchronous states that are induced by an adaptivecoupling structure, we consider the following model of N adaptively coupled phase oscillators [1, 3]˙ φ i = ω + σ N (cid:88) j =1 a ij κ ij sin( φ i − φ j + α ) , (9)˙ κ ij = − (cid:15) ( κ ij + a ij sin( φ i − φ j + β )) , (10)where φ i represents the phase of the i th oscillator, ω isits natural frequency which we set to zero in a rotatingframe.The synchronous state of (9)–(10) is given by s ( t ) =( σr sin α sin β ) t and κ sij = − a ij sin β . Using (7)–(8), thestability of the synchronous state is described by the (d) (e)(a) (b)Λ I m ( σ µ ) I m ( σ µ ) Re( σµ ) Re( σµ ) (c)(f) Re( σµ )FIG. 1. Master stability function Λ( µ ) for the adaptive phaseoscillator network (9)–(10). Regions belonging to negativeLyapunov exponents Λ are colored blue. The curve whereΛ( µ ) = 0 is given as a black solid line. In panels (a) and(b) the case without adaptation ( (cid:15) = 0) is presented for β = − . π and β = 0 . π , respectively. Other panels: (cid:15) = 0 . β = − . π , (d) β = − . π , (e) β = 0 . π , and (f) β = 0 . π . In all panels α = 0 . π . quadratic characteristic polynomial λ + ( (cid:15) − σµ cos( α ) sin( β )) λ − (cid:15)σµ sin( α + β ) = 0 . (11)The master stability function for the synchronous solu-tion is given as the maximum real part Λ = max Re( λ , )of the solutions λ , of the polynomial (11). These solu-tions λ , should be considered as functions of the com-plex parameter µ determining the network structure. Itis convenient, however, to use the parameter σµ in ourcase.Figure 1 displays the master stability function deter-mined for different adaptation rules controlled by β . Theblue-colored areas correspond to regions that lead to sta-ble dynamics. By changing the control parameter β , var-ious shapes of the stable regions are visible. For someparameters, e.g., Fig. 1(c,d,e), almost a whole half-spaceeither left or right of the imaginary axis belongs to thestable regime. This resembles the case of no adaptationwhere the stability of the synchronous state is solely de-scribed by the sign of the real part of σµ sin β cos α , seeFig. 1(a,b). We also find parameters where most values σµ correspond to unstable dynamics, except for an is-land, i.e., a bounded region in σµ parameter space, seeFig. 1(f).To understand the emergence of the stability islands,we analyze the boundary that separates the stable (Λ <
0) from the unstable region (Λ > λ = 0, or, equivalently, λ = i γ . Substituting this into Eq. (11), we obtain a pa-rameterized expression for the boundary as a function of γ that has the form σµ = Z ( γ ) with Z ( γ ) given ex-plicitly in the Supplemental material [105]. The latterparametrization of the boundary is displayed in Fig. 1 asthe solid black line. It is straightforward to show that astability island exists if sin( α + β ) / (cos α sin β ) <
0. Thelatter condition indicates a certain balance between thecoupling and adaptation function. We emphasize thatthe emergence of stability islands is a direct consequenceof adaptation. Without adaptation, the boundary sim-plifies to the axis Re µ = 0, see Figs. 1(a,b).In the following, we analyze the behavior of the adap-tive network of phase oscillators (9)–(10) in the presenceof a stability island, and show how such an island in-troduces a desynchronization transition with increasingoverall coupling σ . To measure the coherence, we use thecluster parameter R C [1, 2], which is given by the numberof pairwise coherent oscillators normalized by the totalnumber of pairs N . In the case of complete synchro-nization, frequency clustering, or incoherence, the clusterparameter values are R C = 1, 1 < R C <
0, or R C = 0,respectively, see Supplemental Material for details [105].The top panel in Fig. 2 shows the cluster parameter R C for different values of the overall coupling constant σ . We observe that for small σ , the synchronous stateis stable, see Fig. 2(a,d,g). This stability follows directlyfrom the master stability function since all values σµ i forall Laplacian eigenvalues lie within the stability island,see Fig. 2(a).By increasing the coupling strength σ , the values σµ i move out of the stability island ( µ i remain the same), andthe synchronous state becomes unstable, see Fig. 2(b,c).For intermediate values of σ , multiclusters with hierar-chical structure in the cluster size emerge, see Fig. 2(e,h)for a three-cluster state. Increasing the coupling con-stant further leads to the emergence of incoherence. InFig. 2(f,i), the coexistence of a coherent and an inco-herent cluster is presented. Such chimera-like stateshave been numerically studied for adaptive networks in[1, 2, 5].In the following, we show how our findings are trans-ferred to a more realistic set-up of coupled neurons withsynaptic plasticity. For this, we consider a networkof FitzHugh-Nagumo neurons [92–95] coupled throughchemical excitatory synapses [96–98] equipped with plas-ticity: τ ˙ u i = u i − u i − v i − σu i N (cid:88) j =1 a ij κ ij I j , (12)˙ v i = u i + a − bv i , (13)˙ I i = α ( u i )(1 − I i ) − I i /τ syn , (14)˙ κ ij = − (cid:15) (cid:16) κ ij + a ij e − β ( u i − u j + β ) (cid:17) . (15)Here u i denotes the membrane potential and v i summa-rizes the recovery processes for each neuron; I i describes (g) (h) (i)Index i (d) (e) (f)Re( σµ )(b)(a) (c) σ (a,d,g) (b,e,h) (c,f,i) φ i h ˙ φ i i R C I m ( σ µ ) FIG. 2. Dynamics in the network of 200 oscillators (9)–(10) with random adjacency matrix A c [105], and differentvalues of overall coupling strength σ . Adiabatic continuationfor increasing σ with the stepsize of 0 . φ i = 0, κ ij = − a ij sin β . The top panelshows the cluster parameter R C vs σ . For the three valuesof σ : (a,d,g) σ = 0 . σ = 0 . σ =0 . σµ i , where µ i arethe N Laplacian eigenvalues of A c ; in (d,e,f) snapshots for φ i at t = 30000; and in (g,h,i) the temporal average of thephase velocities (cid:104) ˙ φ i (cid:105) over the last 5000 time units. Otherparameters: α = 0 . π , β = 0 . π , (cid:15) = 0 . the synaptic output for each neuron; the parameters a = 0 . b = 0 . τ = 0 .
08 and (cid:15) = 0 .
01 are fixed time scaleseparation parameters between the fast activation andslow inhibitory processes in each neuron, and betweenthe fast oscillatory dynamics and the slow adaptation ofthe coupling weights, respectively. The synaptic recoveryfunction is given by α ( u ) = 2 / (0 . − u/ . τ syn = 5 /
6. For more details onthe model, we refer to [98, 105]. The form of the synapticplasticity is similar to the rules used in [6, 7]. We con-sider β and β as control parameters of the adaptation function. Note that β and β are uniquely determinedby the values of h (0) and D h (0) of the plasticity rule, andthese are the only essential parameters of the plasticityfunction, regarding the stability of the synchronous state,see Eqs. (7)–(8).The synchronous state of the network of FitzHugh-Nagumo neurons (12)–(15) satisfies Eqs. (3)–(4), and itis periodic for the chosen parameter values. Using our ex-tended master stability approach, we determine numeri-cally the master stability function which is the maximumLyapunov exponent of Eqs. (7)–(8).In Fig. 3(a,b,c), we show the master stability functionin dependence on the parameter µ/r for different valuesof the overall coupling constant σ . We observe a sta-bility island for the chosen set of parameters, see theSupplemental material for other parameter values [105].In contrast to the phase oscillator network in Fig. 2,the master stability function does not scale linearly with σ . This is due to the non-diffusive coupling function inEq. (12). Moreover, with increasing σ , the size of thestability island shrinks. Since all Laplacian eigenvalues µ i are independent of σ , we observe that µ i /r move outof the stability island with increasing σ . For the globallycoupled network, in particular, we have either µ i /r = 0or µ i /r = 1. Therefore, with increasing σ , we find atransition from complete coherence, see Fig. 3(a,d,g) topartial synchronization and incoherence. We further ob-serve that closely after destabilization, a large frequencycluster remains visible, see Fig. 3(b,e,h). For higher over-all coupling, the cluster sizes shrink, and the number ofsmall clusters increases, see Fig. 3(c,f,i).In summary, we have developed a master stability ap-proach for a general class of adaptive networks. Thisapproach allows for studying the subtle interplay be-tween nodal dynamics, adaptivity, and a complex net-work structure. The master stability approach has beenfirst applied to a paradigmatic model of adaptively cou-pled phase oscillators. We have presented several typ-ical forms of the master stability function for differentadaptation rules, and observed adaptivity-induced sta-bility islands. Besides, we have shown that stability is-lands give rise to the emergence of multicluster states andchimera-like states in the desynchronization transitionfor an increasing overall coupling strength. Qualitativelythe same phenomena have been shown for a more realis-tic network of non-diffusively coupled FitzHugh-Nagumoneurons with synaptic plasticity. In this set-up, the emer-gence of a stability island and a desynchronization tran-sition have been found as well.The theoretical approach introduced in this Letter pro-vides a powerful tool to study collective effects in more re-alistic neuronal network models, including synaptic plas-ticity [32, 82]. While our approach is presented fordifferentiable models, it might be generalized to non-continuous models of spiking neurons equipped withspike timing-dependent plasticity [90, 91]. Our findings h f i i (f)(e)(d) I m ( µ / r ) (a) (c)(b) u i (g) (h) (i)Index i Re( µ/r )FIG. 3. Dynamics of globally coupled network of 200FitzHugh-Nagumo neurons with plasticity Eqs. (12)–(15).Adiabatic continuation for an increasing overall couplingstrength σ with the step size 0 . σ : (a,d,g) σ = 0 . σ = 0 . σ = 0 . µ i /r ,where µ i are the Laplacian eigenvalues (color code as inFig. 1), in (d,e,f) the average frequency (cid:104) f i (cid:105) , and in (g,h,i)snapshots for u i at t = 10000. Here (cid:104) f i (cid:105) = M i / M i is the number of rotations (spikes) of neuron i during thetime interval of length 1000. The control parameters for theadaptation rule β and β are chosen such that h (0) = 0 . h (0) = (80 , , on the transition from coherence to incoherence revealthe role adaptivity plays for the formation of partiallysynchronized patterns which are important for under-standing the functioning of neuronal systems [100]. Be-yond neuronal networks, adaptation is a well-known con-trol paradigm [101–104]. Our extended master stabilityapproach provides a generalized framework to study var-ious adaptive control schemes for a wide range of dynam-ical systems.This work was supported by the German Re-search Foundation DFG, Project Nos. 411803875 and440145547. ∗ [email protected][1] M. E. J. Newman, SIAM Review , 167 (2003).[2] A. Pikovsky, M. Rosenblum, and J. Kurths, Synchroniza-tion: a universal concept in nonlinear sciences , 1st ed.(Cambridge University Press, Cambridge, 2001).[3] S. H. Strogatz, Nature , 268 (2001).[4] A. Arenas, A. D´ıaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Phys. Rep. , 93 (2008).[5] S. Boccaletti, A. N. Pisarchik, C. I. del Genio, andA. Amann,
Synchronization: From Coupled Systems toComplex Networks (Cambridge University Press, Cam-bridge, 2018).[6] Y. Kuramoto,
Chemical Oscillations, Waves and Turbu-lence (Springer-Verlag, Berlin, 1984).[7] S. Yanchuk, Y. Maistrenko, and E. Mosekilde, Math.Comp. Simul. , 491 (2001).[8] F. Sorrentino and E. Ott, Phys. Rev. E , 056114 (2007).[9] I. Belykh and M. Hasler, Chaos , 016106 (2011).[10] M. Golubitsky and I. Stewart, Chaos , 094803 (2016).[11] Y. Zhang and A. E. Motter, arXiv:2003.05461 (2020).[3] R. Berner, E. Sch¨oll, and S. Yanchuk, SIAM J. Appl.Dyn. Syst. , 2227 (2019).[13] P. Jaros, S. Brezetsky, R. Levchenko, D. Dudkowski,T. Kapitaniak, and Y. Maistrenko, Chaos , 011103(2018).[14] Y. Kuramoto and D. Battogtokh, Nonlin. Phen. in Com-plex Sys. , 380 (2002).[15] D. M. Abrams and S. H. Strogatz, Phys. Rev. Lett. ,174102 (2004).[16] A. E. Motter, Nat. Phys. , 164 (2010).[17] M. J. Panaggio and D. M. Abrams, Nonlinearity , R67(2015).[18] E. Sch¨oll, Eur. Phys. J. Spec. Top. , 891 (2016).[19] S. Majhi, B. K. Bera, D. Ghosh, and M. Perc, Phys. LifeRev. , 100 (2019).[20] O. E. Omel’chenko and E. Knobloch, New J. Phys. ,093034 (2019).[21] E. Sch¨oll, A. Zakharova, and R. G. Andrze-jak, Front. Appl. Math. Stat. , 62 (2019), doi:10.3389/fams.2019.00062.[22] A. Zakharova, Chimera Patterns in Networks: Interplaybetween Dynamics, Structure, Noise, and Delay , Under-standing Complex Systems (Springer, 2020).[23] W. Singer, Neuron , 49 (1999).[24] J. Fell and N. Axmacher, Nat. Rev. Neurosci. , 105(2011).[25] C. Hammond, H. Bergman, and P. Brown, Trends Neu-rosci. , 357 (2007).[26] P. Jiruska, M. de Curtis, J. G. R. Jefferys, C. A. Schevon,S. J. Schiff, and K. Schindler, J. Physiol. , 787(2013).[27] V. K. Jirsa, W. C. Stacey, P. P. Quilichini, A. I. Ivanov,and C. Bernard, Brain , 2210 (2014).[28] A. Rothkegel and K. Lehnertz, New J. Phys. , 055006(2014).[29] R. G. Andrzejak, C. Rummel, F. Mormann, andK. Schindler, Sci. Rep. , 23000 (2016).[30] M. Gerster, R. Berner, J. Sawicki, A. Zakharova,A. Skoch, J. Hlinka, K. Lehnertz, and E. Sch¨oll,arXiv:2007.05497 (2020).[31] P. A. Tass, I. Adamchic, H. J. Freund, T. von Stackel-berg, and C. Hauptmann, Restor. Neurol. Neurosci. ,137 (2012).[32] P. A. Tass and O. V. Popovych, Biol. Cybern. , 27(2012).[33] P. Uhlhaas, G. Pipa, B. Lima, L. Melloni, S. Neuen-schwander, D. Nikolic, and W. Singer, Front. Integr. Neu-rosci. , 17 (2009).[34] M. Rohden, A. Sorge, M. Timme, and D. Witthaut,Phys. Rev. Lett. , 064101 (2012).[35] A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa, Nat. Phys. , 191 (2013).[36] P. J. Menck, J. Heitzig, J. Kurths, and H. J. Schellnhu-ber, Nat. Commun. , 3969 (2014).[37] H. Taher, S. Olmi, and E. Sch¨oll, Phys. Rev. E ,062306 (2019).[38] L. M. Pecora and T. L. Carroll, Phys. Rev. Lett. , 2109(1998).[39] A. Brechtel, P. Gramlich, D. Ritterskamp, B. Drossel,and T. Gross, Phys. Rev. E , 032307 (2018).[40] L. Tang, X. Wu, J. L¨u, J. Lu, and R. M. D’Souza, Phys.Rev. E (2019), 012304.[41] R. Berner, J. Sawicki, and E. Sch¨oll, Phys. Rev. Lett. , 088301 (2020).[42] F. Sorrentino, New J. Phys. , 033035 (2012).[43] R. Mulas, C. Kuehn, and J. Jost, Phys. Rev. E ,062313 (2020).[44] C. U. Choe, T. Dahms, P. H¨ovel, and E. Sch¨oll, Phys.Rev. E , 025205(R) (2010).[45] V. Flunkert, S. Yanchuk, T. Dahms, and E. Sch¨oll, Phys.Rev. Lett. , 254101 (2010).[46] A. Keane, T. Dahms, J. Lehnert, S. A. Suryanarayana,P. H¨ovel, and E. Sch¨oll, Eur. Phys. J. B , 407 (2012).[47] Y. N. Kyrychko, K. B. Blyuss, and E. Sch¨oll, Chaos ,043117 (2014).[48] J. Lehnert, Controlling synchronization patterns in com-plex networks , Springer Theses (Springer, Heidelberg,2016).[49] R. B¨orner, P. Schultz, B. ¨Unzelmann, D. Wang, F. Hell-mann, and J. Kurths, arXiv:1911.09730 (2020).[50] T. Dahms, J. Lehnert, and E. Sch¨oll, Phys. Rev. E ,016202 (2012).[51] L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E.Murphy, and R. Roy, Nat. Commun. , 4079 (2014).[52] C. I. del Genio, J. G´omez-Garde˜nes, I. Bonamassa, andS. Boccaletti, Sci. Adv. , e1601679 (2016).[53] K. A. Blaha, K. Huang, F. Della Rossa, L. M. Pecora,M. Hossein-Zadeh, and F. Sorrentino, Phys. Rev. Lett. , 014101 (2019).[54] D. J. Stilwell, E. M. Bollt, and D. G. Roberson, SIAMJ. Appl. Dyn. Syst. , 140 (2006).[55] V. Kohar, P. Ji, A. Choudhary, S. Sinha, and J. Kurths,Phys. Rev. E , 022812 (2014).[56] C. Zhou and J. Kurths, Phys. Rev. Lett. , 164102(2006).[57] V. N. Belykh, I. Belykh, and M. Hasler, Physica D ,159 (2004).[58] I. Belykh, E. de Lange, and M. Hasler, Phys. Rev. Lett. , 188101 (2005).[59] I. Belykh, V. N. Belykh, and M. Hasler, Physica D ,42 (2006).[60] I. Belykh, V. N. Belykh, and M. Hasler, Chaos ,015102 (2006).[61] H. Markram, J. L¨ubke, M. Frotscher, and B. Sakmann,Science , 213 (1997).[62] L. F. Abbott and S. Nelson, Nat. Neurosci. , 1178(2000).[63] N. Caporale and Y. Dan, Annu. Rev. Neurosci. , 25(2008).[64] C. Meisel and T. Gross, Phys. Rev. E , 061917 (2009).[65] K. Mikkelsen, A. Imparato, and A. Torcini, Phys. Rev.Lett. , 208101 (2013).[66] K. Mikkelsen, A. Imparato, and A. Torcini, Phys. Rev.E , 062701 (2014). [67] S. Jain and S. Krishna, Proc. Natl. Acad. Sci. , 543(2001).[68] C. Kuehn, Math. Model. Nat. Phenom. , 402 (2019).[69] T. Gross, C. J. D. D’Lima, and B. Blasius, Phys. Rev.Lett. , 208701 (2006).[70] S. R. Proulx, D. E. L. Promislow, and P. C. Phillips,Trends Ecol. Evol. , 345 (2005).[71] E. A. Martens and K. Klemm, Front. Phys. , 62 (2017).[72] T. Gross and B. Blasius, J. R. Soc. Interface , 259(2008).[73] L. Horstmeyer and C. Kuehn, Phys. Rev. E , 022305(2020).[74] R. Guti´errez, A. Amann, S. Assenza, J. G´omez-Garde˜nes, V. Latora, and S. Boccaletti, Phys. Rev. Lett. , 234103 (2011).[75] X. Zhang, S. Boccaletti, S. Guan, and Z. Liu, Phys. Rev.Lett. , 038701 (2015).[1] D. V. Kasatkin, S. Yanchuk, E. Sch¨oll, and V. I. Nekorkin,Phys. Rev. E , 062211 (2017).[77] M. M. Asl, A. Valizadeh, and P. A. Tass, Front. Physiol. , 1849 (2018).[5] D. V. Kasatkin and V. I. Nekorkin, Chaos , 093115(2018).[2] D. V. Kasatkin and V. I. Nekorkin, Eur. Phys. J. Spec.Top. , 1051 (2018).[4] R. Berner, J. Fialkowski, D. V. Kasatkin, V. I. Nekorkin,S. Yanchuk, and E. Sch¨oll, Chaos , 103134 (2019).[81] P. Feketa, A. Schaum, and T. Meurer, IEEE Trans. Au-tom. Control (2019), 10.1109/tac.2020.3012528.[82] O. V. Popovych, M. N. Xenakis, and P. A. Tass, PLoSONE , e0117205 (2015).[83] L. L¨ucken, O. V. Popovych, P. Tass, and S. Yanchuk,Phys. Rev. E , 032210 (2016).[7] S. Chakravartula, P. Indic, B. Sundaram, and T. Killing-back, PLoS ONE , e0178975 (2017).[85] V. R¨ohr, R. Berner, E. L. Lameu, O. V. Popovych, andS. Yanchuk, PLoS ONE , e0225094 (2019).[86] M. Breakspear, S. Heitmann, and A. Daffertshofer,Front. Hum. Neurosci. , 190 (2010).[87] A. Nabi and J. Moehlis, J. Neural Eng. , 065008 (2011).[88] C. Bick, M. Goodfellow, C. R. Laing, and E. A. Martens,J. Math. Neurosci. , 9 (2020).[89] Y. Maistrenko, B. Lysyansky, C. Hauptmann,O. Burylko, and P. A. Tass, Phys. Rev. E , 066207(2007).[90] J. Ladenbauer, J. Lehnert, H. Rankoohi, T. Dahms,E. Sch¨oll, and K. Obermayer, Phys. Rev. E , 042713(2013).[91] S. Coombes and R. Thul, Eur. J. Appl. Math. , 904(2016).[92] R. A. Stefanescu and V. K. Jirsa, PLoS Comput Biol ,e1000219 (2008).[93] I. Omelchenko, O. E. Omel’chenko, P. H¨ovel, andE. Sch¨oll, Phys. Rev. Lett. , 224101 (2013).[94] W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski, Neuronal Dynamics: From single neurons to networks andmodels of cognition (Cambridge University Press, 2014).[95] D. S. Bassett, P. Zurn, and J. I. Gold, Nat. Rev. Neu-rosci. , 566 (2018).[96] J. Drover, J. Rubin, J. Su, and G. B. Ermentrout, SIAMJ. Appl. Dyn. Syst. , 69 (2004).[97] M. Wechselberger, SIAM J. Appl. Dyn. Syst. , 101(2005).[98] X. Li, J. Wang, and W. Hu, Phys. Rev. E , 041902 (2007).[6] W. J. Yuan and C. Zhou, Phys. Rev. E , 016116 (2011).[100] E. Tang and D. S. Bassett, Rev. Mod. Phys. , 031003(2018).[101] P. De Lellis, M. di Bernardo, F. Garofalo, and M. Por-firi, IEEE Trans. Circuits Syst. I , 2132 (2010).[102] W. Yu, P. DeLellis, G. Chen, M. di Bernardo, andJ. Kurths, IEEE Trans. Autom. Control , 2153 (2012).[103] J. Lehnert, P. H¨ovel, A. A. Selivanov, A. L. Fradkov,and E. Sch¨oll, Phys. Rev. E , 042914 (2014).[104] E. Sch¨oll, S. H. L. Klapp, and P. H¨ovel (eds.), Control ofself-organizing nonlinear systems , (Springer, Berlin, 2016).[105] See Supplemental Material at http... for the proof anddetails on the master stability function for adaptive net-works of non-diffusively coupled oscillators, for detailson the master stability function for the phase oscillatormodel, for the definition of the cluster parameter, for de-tails on the desynchronization transition in the phase oscil-lator model, for details on the network model of coupledFitzHugh-Nagumo neurons with synaptic plasticity, andfor details on the master stability function and the desyn-chronization transition in the model of adaptively coupledFitzHugh-Nagumo neurons.
Supplemental Material on:Desynchronization transitions in adaptive networks
Rico Berner , , ∗ , Simon Vock , Eckehard Sch¨oll , , and Serhiy Yanchuk Institut f¨ur Theoretische Physik, Technische Universit¨at Berlin, Hardenbergstr. 36, 10623 Berlin, Germany Institut f¨ur Mathematik, Technische Universit¨at Berlin, Straße des 17. Juni 136, 10623 Berlin, Germany Bernstein Center for Computational Neuroscience Berlin, Humboldt-Universit¨at, Philippstraße 13, 10115 Berlin,Germany
I. DERIVATION OF THE MASTER STABILITYFUNCTION FOR ADAPTIVE COMPLEXNETWORKS
In this section, we derive the master stability functionfor system (1)–(2) from the main text. For convenience,we repeat these equations here:˙ x i = f ( x i ) − σ N (cid:88) j =1 a ij κ ij g ( x i , x j ) , (S1)˙ κ ij = − (cid:15) ( κ ij + a ij h ( x i − x j )) , (S2)where the adjacency matrix has constant row sum r = (cid:80) Nj =1 a ij .Let ( s ( t ) , κ sij ) be the synchronous state, i.e., x i = s ( t )and κ ij = κ sij for all i, j = 1 . . . , N . This state solves theset of differential Eqs. (3)–(4) of the main text.In order to describe the local stability of the syn-chronous state, we derive the variational equation forsmall perturbations close to this state. For this, we in-troduce the following vector variables denoting the devi-ations from the synchronized state: ξ = x − I N ⊗ s , and χ = κ − κ s with x = ( x , · · · , x N ) T , κ = ( κ , · · · , κ N , κ , · · · , κ NN ) T , where ⊗ denotes the Kronecker product. Using the fol-lowing notations a i = ( a i , . . . , a iN ) , diag( a i ) = a i . . . a iN , and the N × N , N × N , and N × N matrices B = a . . . a N ,C = B T − D,D = diag( a )...diag( a N ) , respectively, the variational equation reads (cid:18) ˙ ξ ˙ χ (cid:19) = (cid:18) S − σB ⊗ g ( s , s ) − (cid:15)C ⊗ D h (0) − (cid:15) I N (cid:19) (cid:18) ξχ (cid:19) , (S3)where S = I N ⊗ D f ( s )+ σh (0) ( r I N ⊗ D g ( s , s ) + A ⊗ D g ( s , s )) . We note that matrices
B, C , and D satisfy the relations B · B T = r I N , B · D = A , and B · C = L , which can beobtained by straightforward calculation.Due to the structure of the variational equation (S3),there exist N − N eigenvalues λ = − (cid:15) . The correspond-ing time-independent eigenspace can be found from (cid:18) S − (cid:15) I Nd − σB ⊗ g ( s , s ) − (cid:15)C ⊗ D h (0) 0 (cid:19) (cid:18) ξχ (cid:19) = 0 . One can see that ( ξ , χ ) such that ξ = 0 and B χ = 0are the time-independent eigenvectors. Moreover, therelation B χ = 0 defines N − N linearly independenteigenvectors spanning the eigenspace corresponding tothe eigenvalues λ = − (cid:15) . This follows from the fact that χ is N -dimensional and rank( B ) = N if the row sum r of A is non-zero.With these prerequisites we are now able to simplifythe local stability analysis on adaptive networks and finda master stability function. Let (S1) – (S2) possess a synchronous solution ( s , κ sij ) .Further, let (S3) be the variational equations around thissynchronous solution and assume that the Laplacian ma-trix L is diagonalizable. Then, the synchronous solutionis locally stable if and only if for all eigenvalues µ ∈ C ofthe Laplacian matrix, the largest Lyapunov exponent (ifit exists), i.e., the master stability function Λ( µ ) , of thefollowing system is negative d ζ d t = (cid:18) D f ( s ) + σrh (0) (cid:0) D g ( s , s )+ (1 − µr )D g ( s , s ) (cid:1)(cid:19) ζ − σg ( s , s ) κ, (S4)d κ d t = − (cid:15) ( µ D h (0) ζ + κ ) . (S5)Here, ζ ∈ C d and κ ∈ C .In the following we present the derivation of (S4)–(S5). As it is shown above, there are N − N in-dependent vectors w l ( l = 1 , . . . , N − N ) spanningthe kernel of B , i.e. B w l = 0. Using the Gram-Schmidt procedure we find an orthonormal basis forker( B ) = span { v , . . . , v N − N } . With this, we define the N × ( N − N ) matrix Q = ( v , . . . , v N − N ). Considernow the ( N + N d ) × ( N + N d ) matrix R = (cid:18) I Nd /r ) B T Q (cid:19) with left inverse R − = I Nd B Q T , i.e., R − R = I N + Nd . Introduce the new coordinatesgiven by R (cid:18) ˆ ξ ˆ χ (cid:19) = (cid:18) ξχ (cid:19) for which the variational equa-tion then readsdd t (cid:18) ˆ ξ ˆ χ (cid:19) = R − (cid:18) S − σB ⊗ g ( s , s ) − (cid:15)C ⊗ D h (0) − (cid:15) I N (cid:19) R (cid:18) ˆ ξ ˆ χ (cid:19) . We further obtain R − (cid:18) S − σB ⊗ g ( s , s ) − (cid:15)C ⊗ D h (0) − (cid:15) I N (cid:19) R = R − (cid:18) S − σ I N ⊗ g ( s , s ) 0 − (cid:15)C ⊗ D h (0) − (cid:15)/rB T − (cid:15)Q (cid:19) = S − σ I N ⊗ g ( s , s ) 0 − (cid:15)L ⊗ D h (0) − (cid:15) I N − (cid:15)Q T C ⊗ D h (0) 0 − (cid:15) I N − N . These equations yield that there are
N d + N coupleddifferential equations leftdd t (cid:18) ˆ ξ ˜ χ (cid:19) = (cid:18) S − σ I N ⊗ g ( s , s ) − (cid:15)L ⊗ D h (0) − (cid:15) I N (cid:19) (cid:18) ˆ ξ ˜ χ (cid:19) (S6)with ˜ χ = ˆ χ that determine the stability for the syn-chronous state, and N − N slave equationsdd t ¯ χ = (cid:0) − (cid:15)Q T C ⊗ D h (0) 0 − (cid:15) I N − N (cid:1) ˆ ξ ˜ χ ¯ χ with ¯ χ = ( ˆ χ T , . . . , ˆ χ TN ) T which are driven by the vari-ables ˆ ξ and ˜ χ and, hence, can be solved explicitly oncethe latter once are known. By assumption, there is aunitary matrix D L = U H LU where D L is the diago-nalization of the Laplacian matrix L . Transforming thedifferential equation (S6) by using the unitary transfor-mation U , we getdd t (cid:18) ζκ (cid:19) = (cid:18) I N ⊗ Df ( s ) + σh (0) ( r I N ⊗ D g ( s , s ) + ( r I N − D L ) ⊗ D g ( s , s )) − σ I N ⊗ g ( s , s ) − (cid:15)D L ⊗ D h (0) − (cid:15) I N (cid:19) (cid:18) ζκ (cid:19) where (cid:18) U ⊗ I d U (cid:19) (cid:18) ˆ ξ ˜ χ (cid:19) = (cid:18) ζκ (cid:19) .Remarkably, the master stability function Λ dependsexplicitly on the row sum r . Moreover, the master sta-bility function seems to depend on σ , r , and µ indepen-dently. The time scale separation parameter (cid:15) is alwayskept fixed. However, in any case, one parameter can bedisregarded. To see this, we note that the solution to theEq. (S5) is explicitly solvable and the solution reads κ = κ e − (cid:15) ( t − t ) − (cid:15)µ D h (0) (cid:90) tt e − (cid:15) ( t − t (cid:48) ) ζ ( t (cid:48) ) dt (cid:48) , where the first term vanishes for t → ∞ and hence canbe neglected (when studying asymptotic stability for t →∞ ). We use this and rewrite the asymptotic dynamics of (S4)–(S5) in its integro-differential formd ζ d t = (D f ( s ) + σrh (0) (D g ( s , s )+(1 − µr )D g ( s , s ) (cid:17)(cid:17) ζ + (cid:15)σr µr g ( s , s )D h (0) (cid:90) tt e − (cid:15) ( t − t (cid:48) ) ζ ( t (cid:48) ) dt (cid:48) . Hence, the master stability function can be regarded as afunction of two parameters, i.e., Λ( σ, µ, r ) = Λ( σr, µ/r ).Furthermore, in case of diffusive coupling, i.e., g ( x , y ) = g ( x − y ), the master stability function can be regardedas a function of only one parameter Λ( σ, µ, r ) = Λ( σµ ). II. MASTER STABILITY FUNCTION FORADAPTIVE PHASE OSCILLATOR NETWORKS
In this section, we provide a brief analysis of the masterstability function for the adaptive Kuramoto-Sakaguchinetwork (9)–(10) of the main text. Using the result ofSection I, the stability of the synchronous state of sys-tem (9)–(10) of the main text is governed by the twodifferential equationsdd t (cid:18) ζκ (cid:19) = (cid:18) µσ cos( α ) sin( β ) − σ sin( α ) − (cid:15)µ cos( β ) − (cid:15) (cid:19) (cid:18) ζκ (cid:19) , where µ ∈ C stands for all eigenvalues of the Laplacianmatrix L corresponding to the base network described bythe adjacency matrix A . The characteristic polynomialin λ of the latter system is of degree two and reads λ + ( (cid:15) − σµ cos( α ) sin( β )) λ − (cid:15)σµ sin( α + β ) = 0 . (S7)The master stability function is given as Λ( σµ ) =max(Re( λ ) , Re( λ )) where λ and λ are the two so-lutions of the quadratic polynomial (S7). Figure 1 ofthe main text displays the master stability function fordifferent parameters.The boundary of the region in σµ parameter space thatcorresponds to stable local dynamics, is given by λ = i γ with γ ∈ R . Plugging this into Eq. (11) of the main text,we obtain σµ = Z ( γ ) = a ( γ ) + i b ( γ )with a ( γ ) = (cid:15) γ (cos α sin β − sin( α + β )) γ cos α sin β + (cid:15) sin ( α + β ) ,b ( γ ) = γ cos α sin β + (cid:15) γ sin( α + β ) γ cos α sin β + (cid:15) sin ( α + β ) . Due to the symmetry of the master stability function,a necessary condition to observe a stability island isthat the curve σµ ( γ ) possesses two crossings with thereal axis, i.e., two real solutions for b ( γ ) = 0. Thethree crossings are given by γ = 0 and as real solu-tions γ and γ of γ cos α sin β = − (cid:15) sin( α + β ). Fromthis we deduce the existence condition for stability is-lands: sin( α + β ) / (cos α sin β ) < (cid:15) > a ( γ ) = a ( γ ). III. THE CLUSTER PARAMETER
In this section, we introduce the cluster parameter R C as a measure for coherence in a system of coupled phaseoscillators. A measure that can be used in order to detectfrequency synchronization between two oscillators relies on the mean phase velocity (average frequency) of eachphase oscillatorΩ i = lim T →∞ T ( φ i ( t + T ) − φ i ( t )) . (S8)The frequency synchronization measure between nodes isgiven by Ω ij = (cid:40) , if Ω i − Ω j = 0 , , otherwise . (S9)Numerically the limit is approximated by a very longaveraging window. In addition, we use a sufficiently smallthreshold (cid:36) in order to detect frequency synchronizationnumerically, i.e., Ω ij = 1 if Ω i − Ω j < (cid:36) . For the analysispresented here and in the main text, we use (cid:36) = 0 . ij , we define the cluster parameter R C = 1 N N (cid:88) i,j =1 Ω ij . (S10)The cluster parameter measures the following. First, foreach frequency cluster, the total number of pairwise syn-chronized nodes is computed. Second, all pairs of twonodes from the same cluster are summed up and normal-ized by the number of all possible pairs of nodes N . Incase of full synchronization, frequency clustering, or inco-herence the values of the cluster parameter are R C = 1,1 < R C <
0, or R C = 0, respectively. A similar measurecan be found in Refs. [1, 2]. IV. DESYNCHRONIZATION TRANSITIONAND THE FORMATION OF PARTIALSYNCHRONIZATION PATTERNS IN ADAPTIVEPHASE OSCILLATOR NETWORKS
In this section, we provide further details on the desyn-chronization transition in a network of adaptively cou-pled phase oscillators (9)–(10).Figure S1 shows the cluster parameter R C for differ-ent values of the coupling constant σ . In the adiabaticcontinuation, we increase σ step-wise after an integrationtime of t = 10000. For each simulation, the final stateof the previous simulations is taken as the initial condi-tion with an additional small perturbation. Note that R C = 1 refers to full in-phase synchrony of the oscilla-tors. We observe that, for small σ , the synchronous stateis stable, see Fig. S1(d,g,j). Here, the stability of the syn-chronous state is directly implied by the master stabilityfunction. We note that all Laplacian eigenvalues µ i of aglobally coupled network are given by either µ i = 0 or µ i = N . In Figure S1(a), all master function parameters σµ lie within the stability island.By increasing the coupling constant, the values σµ i move out of the stability regions and the synchronous (a) (b) (c)(f)(e)(d)(g) (h) (i)(l)(k)(j) I m ( σ µ ) Re( σµ ) φ j h ˙ φ j i I nd e x i R C (a,d,g,j)(b,e,h,k) (c,f,i,l) σ Index jκ ij −| sin β | | sin β | σ than in Fig.2 of the maintext. Adiabatic continuation for increasing σ with the step-size of 0 . φ i = 0, κ ij = − a ij sin β . The top panel shows the cluster parame-ter R C vs σ . For the three values of σ : (a,d,g,j) σ = 0 . σ = 0 . σ = 0 . σµ i , where µ i are the N Laplacian eigenvalues of A ; in (d,e,f) snapshots for φ i at t = 30000; in (g,h,i) the temporal average of the phase veloci-ties (cid:104) ˙ φ i (cid:105) over the last 5000 time units; and in (j,k,l) snapshotsfor the coupling matrix κ ij at t = 30000. Other parameters: α = 0 . π , β = 0 . π , (cid:15) = 0 . state becomes unstable. For intermediate values of σ the emergence of multiclusters with hierarchical struc-ture in the cluster size are observed. In Figure S1(e,h,k)a multicluster states is shown with three clusters. Note u i membrane potential/activator v i recovery/inhibitor variable I i synaptic output variable κ ij variable coupling weights N number of oscillators a ij entries of adjacency matrix, a ij ∈ { , } σ overall coupling strength r row sum, i.e., r = (cid:80) Nj =1 a ij a = 0 . , b = 0 . τ = 0 .
08 controls time separation between fastactivation and slow inhibition (cid:15) = 0 .
01 controls time separation between fastoscillation and slow adaptation τ syn = 5 / u shp = 0 .
05 coupling shape parameter β , β adaption control parametersTABLE S1. The table provides the meaning for each variableand parameter used in (S11)–(S14). that for the system (9)–(10) of the main text, in-phasesynchronous and antipodal clusters have the same prop-erties [3, 4]. In Refs. [3, 4] the role of the hierarchicalstructure of the cluster sizes have been discussed. In-creasing the coupling constant further shows the emer-gence of incoherence. In Figure S1(f,i,l), we show the co-existence of a coherent and an incoherent cluster. Thesestates, also called chimera-like states, have been numeri-cally analyzed in Refs. [1, 2, 5]. V. NETWORK OF COUPLEDFITZHUGH-NAGUMO NEURONS WITHSYNAPTIC PLASTICITY
In this section, we describe the model of coupledFitzHugh-Nagumo neurons with synaptic plasticity andpresent the synchronous state used in the main text. Themodel is given by τ ˙ u i = u i − u i − v i − σ N (cid:88) j =1 a ij κ ij u i I j , (S11)˙ v i = u i + a − bv i , (S12)˙ I i = α ( u i )(1 − I i ) − I i /τ syn , (S13)˙ κ ij = − (cid:15) (cid:16) κ ij + a ij e − β ( u i − u j + β ) (cid:17) , (S14)where α ( u i ) = 2 / (0 . − u i / . β and β as control parameters. In par-ticular, we have β = − h (0) / (2 Dh (0) β ) and β =(2D h (0) ln(D h (0) )) /h (0) where D h (0) denotes thefirst component of D h (0).The synchronous state of the equations (S11)–(S14) is u v I FIG. S2. Limit cycle in Eqs. (S15)–(S17) as solid line andthe projection onto the u - v -plane as dashed line. Parameters: σ = 0 . r = 200, h (0) = 0 . h (0) = (80 , , given by a solution of τ ˙ u s = u s − u i − v s + σru s I s e − β β , (S15)˙ v s = u s + a − bv s , (S16)˙ I s = α ( u s )(1 − I s ) − I s /τ syn , (S17) κ sij = − a ij e − β β , (S18)where ( u i , v i , I i ) = s = ( u s , v s , I s ) for all i = 1 , . . . , N .In Fig. S2, we display a limit cycle as a stable numericalsolution of (S15)–(S18) for the set of parameters used inthe main text. VI. THE MASTER STABILITY FUNCTIONAND DESYNCHRONIZATION TRANSITION INADAPTIVE NETWORKS OFFITZHUGH-NAGUMO NEURONS
In this section, we consider the model of adaptivelycoupled FitzHugh-Nagumo neurons (S11)–(S14). Wegive insights into the derivation of the system’s mas-ter stability function as well as on the desynchronizationtransition induced by the adaptivity.In order to investigate the local stability of the syn-chronous states that solves Eqs. (S15)–(S18), see Fig. S2,we linearize Eqs. (S11)–(S14) around these states. Usingthe results of Section I, the stability of the synchronoussolution is governed by the set of equationsd ζ d t = (cid:18) D f ( s ) + σrh (0) (cid:0) D g ( s , s )+ (1 − µr )D g ( s , s ) (cid:1)(cid:19) ζ − σg ( s , s ) κ, d κ d t = − (cid:15) ( µ D h (0) ζ + κ ) . I m ( µ / r ) Re( µ/r ) I m ( µ / r ) Re( µ/r )(c) (d)(a) (b)FIG. S3. The master stability functions for the synchronoussolution of (S11)–(S14) and different plasticity rules are dis-played (color code as in Fig. 1 of the main text). Regionsbelonging to negative Lyapunov exponents are colored blue.Parameters: the control parameters β and β are chosensuch that (a) h (0) = 0 .
8, D h (0) = (50 , ,
0) (b) h (0) = − . h (0) = (0 , , h (0) = 0 .
8, D h (0) = (10 , , h (0) = 0 .
4, D h (0) = (50 , , σ = 0 . Here, the derivatives of the functions f , g , and h areD f ( s ) = τ (cid:0) − u s (cid:1) − τ − b τ ( α ( u s )) (1 − I s ) α u shp exp( usu shp ) − α ( u s ) − τ syn , D g ( s , s ) = I s , D g ( s , s ) = u s , D h (0) = (cid:0) − β β exp( − β β ) 0 0 (cid:1) . Using this, we are able to determine numerically themaximum Lyapunov exponents and hence the stabilityof the periodic orbit displayed in Fig. S2. In Fig. S3, weshow different shapes of the master stability function de-pending on the form of the plasticity rule, i.e., dependingon h (0) and Dh (0). We observe that for certain parame-ters almost complete half spaces in the µ/r -plane refer tostable or unstable local dynamics, see Fig. S3(a,b). Thisis similar to Fig. 1(d,e) of the main text where we dis-play the master stability function of the phase oscillatormodel. Most remarkably, similar to the phase oscillatormodel (9)–(10) we find parameters for which stability is-lands exist, see Fig. S3(d).As we know from the example of phase oscillators, thepresence of a stability island may induce a desynchroniza- κ ij h f j i (f)(e)(d) I m ( µ / r ) u j (g) (h) (i)Index j (l)(k)(j) I nd e x i Re( µ/r )(a) (c)(b)FIG. S4. Dynamics of globally coupled network of 200FitzHugh-Nagumo neurons with plasticity Eqs. (12)–(15).Adiabatic continuation for an increasing overall couplingstrength σ with the step size 0 . σ : (a,d,g,j) σ = 0 . σ = 0 . σ = 0 . µ i /r , where µ i are the N Laplacianeigenvalues (color code as in Fig. 1 of the main text), in(d,e,f) the average frequency (cid:104) f i (cid:105) , in (g,h,i) snapshots for u i at t = 10000, and in (j,k,l) snapshots for the coupling ma-trices κ ij at t = 10000. Here (cid:104) f i (cid:105) = M i / M i isthe number of rotations (spikes) of neuron i during the timeinterval of length 1000. The control parameters for the adap-tation rule β and β are chosen such that h (0) = 0 . h (0) = (80 , , tion transition for an increasing overall coupling strength σ . In order to show this transition, we follow the sameapproach already presented in Fig. S1. The results ofthe adiabatic continuation on a globally coupled networkare shown in Fig. S4. We note that in contrast to thecase of phase oscillators, here, the shape of the master stability function depends explicitly on σ . The desyn-chronization is described in the main text. Additionallyto the figure given in the main text, we provide plotsfor the coupling matrices in Fig. S4(j,k,l). The couplingmatrices show very nicely the emergence of partial syn-chronization structures in the transition from coherenceto incoherence which is induced by the stability island. Index j I nd e x i FIG. S5. Adjacency matrix A c of a connected, directed ran-dom network of N = 200 nodes with constant row sum r = 50.The illustration shows the adjacency matrix where black andwhite refer to whether a link between two nodes exist or not,respectively. ∗ [email protected][S1] D. V. Kasatkin, S. Yanchuk, E. Sch¨oll, and V. I. Neko-rkin, Phys. Rev. E , 062211 (2017).[S2] D. V. Kasatkin and V. I. Nekorkin, Eur. Phys. J. Spec.Top. , 1051 (2018).[S3] R. Berner, E. Sch¨oll, and S. Yanchuk, SIAM J. Appl.Dyn. Syst. , 2227 (2019).[S4] R. Berner, J. Fialkowski, D. V. Kasatkin, V. I. Nekorkin,S. Yanchuk, and E. Sch¨oll, Chaos , 103134 (2019).[S5] D. V. Kasatkin and V. I. Nekorkin, Chaos , 093115(2018).[S6] W. J. Yuan and C. Zhou, Phys. Rev. E , 016116 (2011).[S7] S. Chakravartula, P. Indic, B. Sundaram, and T. Killing-back, PLoS ONE12