Complex Rotating Waves and Long Transients in a Ring-Network of Electrochemical Oscillators with Sparse Random Cross-Connections

4.1 MB·

by  Michael Sebek

Complex Rotating Waves and Long Transients in a Ring-Network of Electrochemical
Oscillators with Sparse Random Cross-Connections
1 2 1
Michael Sebek, Ralf Tönjes, and István Z. Kiss
1Department of Chemistry, Saint Louis University,
3501 Laclede Ave., St. Louis, Missouri 63103, USA
2Institute of Physics and Astronomy, Potsdam University, 14476 Potsdam-Golm, Germany
(Dated: January 13, 2016)
6 We perform experiments and phase model simulations with a ring network of oscillatory elec-
1 trochemical reactions to explore the effect of random connections and non-isochronocity of the
0 interactions on the pattern formation. A few additional links facilitate the emergence of the fully
2 synchronized state. With larger non-isochronicity, complex rotating waves or persistent irregular
phase dynamics can derail the convergence to global synchronization. The observed long tran-
n
sientsofirregular phasedynamicsexemplifythepossibility ofasuddenonset of hypersynchronous
a
behaviorwithout any external stimulus or network reorganization.
J

2
PACSnumbers: 89.75.Hc05.45.Xt82.47.-a
1

]
Wave propagation of activity of oscillatory units in by phase model calculations that predict the presence of
n
n rings or linear chains is a fundamental type of pattern complexrotatingwavesandlongtransientsinsmallran-
- formation, that occurs in many biological systems, e.g., dom networks with sufficiently non-isochronous oscilla-
s
i motion of leach [1], the segmentation clock [2], or brain tions. The experimental conditions allow the analysis of
d
waveactivitiesinthecortex[3]. Arotatingpinwheelwas the dependence of pattern formation on the randomness
.
t oneofthefirsttypeofchemicalpatternformationidenti- ofthenetworktopologyandthelevelofnon-isochronicity
a
m fied in the BZ reaction on a ring [4]. As a ring geometry of the interactions among the units.
is often used in chemistry, rotating phase wave patterns To study the properties of complex rotating waves on
-
d havebeenobservedina largenumberofsystems,e.g.,in networks we consider weakly coupled, identical limit cy-
n electrochemicalreactions[5], heterogeneouscatalysis [6], cleoscillatorswithaKuramototypephasemodel[17]for
o
coupled BZ reactors [7] and micro droplets [8]. Mathe- phase differences in a co-rotating frame of reference
c
[ matical analysis using phase models interpreted the ex-
istence and local stability of rotating waves in ring net- N
1 works [9–11]. It was found that the fully synchronized, ϑ˙n = Anmg(ϑm−ϑn). (1)
v
zero phase lag, non-rotating state is the most attracting mX=1
0
7 solution, locally and globally. However, with increasing
whereAnm representsacouplingmatrixandg(∆ϑ)isthe
8 system size, rotating waves with higher winding number
average effect of the coupling for oscillators with phase
2 become more probable in the aggregate [12].
0 difference ∆ϑ. A phase attractive coupling is assumed
Complex engineered and biological systems can often be
. with interaction function g(∆ϑ)=sin(∆ϑ−α)+sin(α).
1 described as networks of discrete, interacting units [13].
The phase shift parameter α is an important system
0
Considering the prevalence of phase waves on rings and
6 property determined by the average shear flow near the
chains,afundamentalquestionishowtherotatingwaves
1 limitcycleinthedirectionofperturbationcausedbycou-
v: manifestinnetworksthatarecomposedofaregularring pling[17],i.e. αquantifiesthenon-isochronicityoftheos-
backbonewithafewadditionalrandomconnections. Nu-
i cillationsinducedbyinteractions. Wenotethatsincethe
X mericalsimulationswithphasemodelsonsparsedirected
dynamics of the model equations (Eq.1) is invariant un-
r networkswithrandominitialphaseshaveshownthatfor
dera changeofα→−α, ϑ→−ϑandt→−t,the phase
a sufficiently non-isochronous oscillations, while the fully
differences and the frequency shift are inverted when α
synchronized state is locally stable, persistent irregular
changes sign, such that sources become sinks of rotating
phase dynamics is the typically observed behavior [14].
waves and vice versa. We assume non-normalized, bidi-
Such prolonged transient behavior can severely impact
rectionalcoupling Amn =Anm ∈{0,1}onringnetworks
systemresponsewhenrobustsynchronizationisrequired
ofN =500oscillatorswithNsc =σN additionalrandom
as it was demonstrated with power grid models [15] or
bidirectional links. The initial conditions for the simula-
when synchronization is undesirable, e.g., in hyper syn-
tionsandthe experimentsis a rotatingwaveonthe ring.
chronous neuronal discharges during seizures [16].
(Random initial conditions give comparable results.)
Inthispaper,weexplorethetypeofspatiotemporalpat- ′
Due to phase attractive coupling g (0) = cosα > 0 and
terns that can be obtained with oscillatory chemical re-
complete connectedness of the network, the fully syn-
actions on bidirectional ring networks with random long
chronized, one-cluster state is always a linearly stable
range connections. The experimental work is motivated
solution of (1) [14]. However, the typical behavior of
2
thenetworkstartingfromgloballydesynchronizedinitial a b
4
conditions, is far more complex than the intuitively ex- 1.5 irregular phase dynamics 0.9 1.5
pectedrelaxationto the one-clusterstate. We character- 0.8 3
ize the stateofthe systembydifferentorderparameters. α1 partial synch. R100..67 α1 ˙varϑ
The k-cluster order parameters Rk =hexp(ikϑn)i where 0.5 2
frozen complex waves
the average is taken instantaneously over all oscillators, 0.4
0.5 0.3 0.5
measurethecoherenceofthedistributionofphasesintok 1
0.2
evenlyspacedclusters. Thevariancevarϑ˙ ofthephaseve- s y n ccohmropnleiztaetion 0.1
locitiesisameasureforfrequencysynchronization. These 00 0.2 0.4σ 0.6 0.8 1 00 0.2 0.4 σ0.6 0.8 1 0
ensembleaveragedmeasuresareshownintheσ vs. αpa- c 1 d 1.5
rameter plane in Figs. 1a,b. 0.8 ααα===000...058 0.15 RR76
With smallσ, the originalring is divided into linear seg- varϑ˙ 1
0.6 0.1 ˙ϑ
ments between the end points of shortcuts, which can R10.4 R67, var
support traveling phase waves. At the interfaces where 0.05 0.5
two such traveling waves meet, the phase differences in 0.2
a state of stable synchronization are restricted. At low 00 0.2 0.4 σ0.6 0.8 1 000 00..55 α 11 11..550
non-isochronicity and low shortcut density these inter-
faces can be frozen when all topological boundary con-
ditions can be met simultaneously. We refer to such a FIG.1: (Coloronline)Orderparametersinmodel(1)attime
pattern, which does not change in time, as frozen com- t=500asfunctionsofshortcutdensityσandnonisochronic-
plexrotatingwavepattern. Whenthetopologicalbound- ity parameter α averaged over ten random network realiza-
tions with N =500 oscillators and rotating wave initial con-
ary conditions are not met (this is likely to occur with
ditions. (a) Color coded Kuramoto order parameter R1 and
a largenumber of oscillators),slowly changinginterfaces ˙
(b)variancevarϑofphasevelocities. Theblackandthewhite
are obtained, reminiscent of vortex glasses in 2d oscilla- linein(a)and(b),respectively,markthecontourofR1 =0.9.
torymedia[18]. BoththeKuramotoorderparameterR1 The light (green) line marks the contour line of varϑ˙ = 0.2.
and the variance of the phase velocities are small in this (c) Kuramoto order parameter as a function of σ for three
regime. When the shortcut density is increased, there different values of α (cf.3c). (d) Mean cluster order parame-
existsatopologicalcross-overtoarandomnetworkwith- ters R6,R7 and variance of phase velocities as functions of α
outlinearchainsegments. Therefore,whenσisincreased at σ=0.15.
the system cannot maintain rotating wavesand the one-
cluster state becomes globally attractive with R1 ≈ 1. plex rotating wave pattern becomes frozen, very long
When α is increasedfrom zero,higher shortcutdensities transient dynamics can be observed. Figure 2 shows
are required for complete synchronization(Figs. 1a-c). an example of such transient dynamics, from random
At large values of non-isochronicity (α & 1) a qualita- initial conditions and with negative non- isochronicity.
tively different type of behavior exists. The topological The time evolution of next-neighbor phase differences
boundary conditions for rotating phase waves along the demonstrates the competition between different phase
ring segments with stationary phase differences are very patterns,withstationarysinksofthephasewaveslocated
difficult tosatisfy simultaneously. Instead,the dominant atthenetworkheterogeneitiesanddynamicallyrearrang-
behavior is persistent irregular dynamics with nonzero ing sources which may form or annihilate upon collision
varϑ˙. The distribution of phase differences during the with a sink or at phase slip events. Figures 2d,e show a
transient and in frozen complex rotating patterns be- transientandastationaryphaseprofile,respectivelyand
comes bimodal, suggesting a preferred phase difference Fig.2b illustrates the time evolution of the phase differ-
that depends on α (Fig. 2b). The emergence of persis- encedistribution. Thephasedifferencesinthestationary
tent irregular dynamics is demonstrated in Fig. 1d by phase pattern are peaked sharply around 2π/7 result-
fixing the shortcutdensity σ, andincreasingthe value of ing in a very precise wavelength and the formation of 7
α. At the transition between frozen and unfrozen com- globalphaseclusters. As showninFig.1a,therealsoex-
plex rotating patterns global clustering can arise result- ists a narrowregime of partialsynchronizationfor larger
inginasharpincreaseintheorderparametersR6 orR7. shortcutdensities whichis replacedin asharpdiscontin-
This global order is mediated by the end points of the uous transition by persistent incoherent phase dynamics
shortcuts in the network (Fig. 2f) : Due to the narrow at values of α approaching π/2 and which may be ana-
distribution of phase differences in a phase locked state, lyzed in a mean field approach [19].
oscillatorsatthesamedistancetoacrosslinkednodehave To confirm the modeling results, experiments were per-
the samephase. Globalclusteringisonlyobservableina formedwithanarrayofN =20,1.00mmdiameternickel
narrow parameter region at criticality and after a long, wires on which an oscillatory metal dissolution reaction
system size dependent transient. takesplacemeasuredbycurrents. Numericalsimulations
In addition, around the transition point, before a com- indicate [23] that regions of frozen rotating patterns,
3
complete synchronization, and irregular phase dynamics a b
500 6
can be clearly distinguished even in such a small setup. 5
The electrodes are coupled into a ring topology with ad- t340000 ∆ϑ)4
ditional random cross-connectionvia resistances and ca- 200 p(23
pacitances. We report the conductance accross the cou- 100 1
pling resistance as coupling strength K. Capacitance is − 3 −2 −1 0∆ϑ1 2 3 0
usedtointroducenon-isochronicitythroughaphaseshift c
π
ionftthheeecxopueprliimngenctuirsreanrto[t2a0t]in[g23w]a.veTahseshinoiwtinalincoFnidgi.ti3oan. 450000 -ϕ)/+1n0.5
Firstwedescribetheresultswithα=0,i.e.,resistivecou- t300 (ϕn0
pling. When a random cross connection was added, one =
200 ϕ
of two scenarios occurred. If the random shortcut con- ∆−0.5
100
nectedtwoelementsatadistancelargerthan4units,the
system quickly converged to a fully synchronized state 0 50 100 150 200 250 n 300 350 400 450 500
3.14
similar to that shown in Fig. 3b. When the distance
between the shortcut elements was smaller, the rotating d 0 ϑn
150 160 170 180 190 200 n 210 220 230 240 −3.14
waves jumped across the connection. Figure 3d shows 3.14
e
such pattern with two cross connections. 0 ϑn
Wehaveperformed16independenttrialsaddingcross- 150 160 170 180 190 200 n 210 220 230 240 −3.14
500
connections successively to a ring configuration. The f 300n
mean order parameter as a function of the added num- 100
−3.1415 −2 −1 0 ϑn 1 2 3.1415
berofshortcutsisshowninFig.3c. Onlythreeshortcuts
were required for the average order parameter to exceed FIG.2: (Coloronline)Longtransienttoacomplexfrozenro-
0.90. Therefore,we canconclude that with α=0 a rela- tating wave pattern from random initial phases in a network
tively small number of shortcuts (σ ≈ 0.15) induces full ofN= 500 phaseoscillators with α=-1.15 andshortcut den-
synchrony. When α was changed to -0.97, with a paral- sity σ= 0.05. (a) Ring network with 26 additional random
shortcuts (b) Time evolution of the density of next neighbor
lel RC coupling, we still observed rotating waves jump-
phase differences. The solid (red) line marks the average of
ing over the connection when the first cross-connection
|∆ϑ|. (c) Time evolution of next neighbor phase differences
was placed up to a distance of 5 units between the ele-
modulus π (color coded). Dark (blue) colors indicate phase
ments. However, when the distance was larger, instead wavestotherightandlight(red)colorsindicatephasewaves
of full synchronization, we observed frozen complex ro- totheleft. Thewhitelinesindicatephaseprofilesϑnbetween
tatingpatternviatheformationofasourceandsinkpair 150≤n≤250showninsubfigures(d)att=300and(e)att=
(Fig. 3e). The order parameter vs. number of shortcuts 500. Large (red) circles and solid (red) lines in panels (a,d-
f) indicate nodes with shortcut connections. Black squares
graphinFig.3cshowsthatwithα=−0.97themeanor-
in panels (d,e) indicate dynamically realized centers of phase
derparameterstartstoincreaseforn>3,anditrequires
waves. Panel(f)showsclusteringofthephasesatt=500with
relatively large number (n > 6) of random shortcuts to R7 ≈0.9 and cross-links (solid lines) connecting neighboring
achieve Kuramoto order larger than 0.90. The presence clusters.
of jumping waves and complex rotating wave patterns
thus contributes to resisting complete synchronization. for over 700 oscillations before the system settled into
When α was further changed to -1.3, the trend of resist- theone-clusterstate. Atα=0thesamenetworkrelaxes
ing the fully synchronized state continued. Figures 3e,f exponentially to synchronization in 230s [23] . Numeri-
demonstrate the inversion of the stationary phase pro- cal simulations with the phase model confirmed that by
file, andthus the directionof the rotating wavesand the changing α from 0 to -1.3 the lifetime of the transient
reversal of source and sink, upon switching the sign of increases about 10 times and diverges as α approaches
the non-isochronicity by adding a capacitance to the in- |π/2| [23]. A density plotofnext neighborphasediffer-
dividual current instead of the coupling current [20] [23] encesduringandafterthetransientisdepictedinFig.4b
. All the patterns in Fig. 3 were reproduced by phase and shows an asymmetric distribution with a preferred
model simulations in [23] . Long transients to both wavelength,skewedtowardsthe initialleft-handedrota-
frozencomplex wavepatterns and identical synchroniza- tional state. In addition, a wide distribution of peak to
tion were observed in experiment. A long transient over peakperiodscanbeobservedduringthetransientasseen
1700s in a network with five random shortcuts near α=- in Fig. 4c, which marks the presence of irregular phase
1.3 is shownin Fig. 4. The time evolutionof next neigh- dynamics. A typical snapshot of a transient complex ro-
bor phase differences during the transient is shown in tating waveis shownin Fig.4d. The arrowsindicate the
Fig. 4a. At least two competing wave patterns, which direction of the rotations as well as the sources (oscil-
aremetastableoverthecourseoftensofoscillationsand lators 9 and 18) and sinks (oscillators 2 and 11) of the
transformvia intermittent phase slips could be observed unstable rotational wave pattern. The order parameter
4
a π b
/ transient
1400 ϑ)n0.8 final
1200 -+1n0.6
(ϑ0.4 1
(s)1000 ∆ϑ=0.2 ∆ϑ)
t800 0 (
p
600 −0.2 0.5
−0.4
400
−0.6
200
−0.8
0
0 10n 20 −3.14−2−1∆0ϑ1 23.14
e1
0.8
R
0.6
0.4
0.2
0
0 500 1000 t(s) 1500 2000 2500
FIG. 4: Long transient to synchronization at α = −1.3
(V=1110 mV, Cc=82 µF, K=0.033 mS). (a) Time evolution
of next neighbor phase differences. (b) Histogram of next
neighbor phase differences during the transient and the final
synchronized state. (c) Histogram of the period per cycle
of all electrodes (white bars are the individual periods and
FIG. 3: The impact of shortcuts and α on the formation of
dark bars the transient). (d) Network topology and a typi-
complexrotatingwavesandonsynchronization. (a)Initialro-
cal complex rotating wave pattern during the transient. Ar-
tatingwave(V=1105mV,K=0.20mS,α=0). (b)Synchrony
rows indicate the wave direction from sources to sinks. (e)
induced by one long distance shortcut. (c) The mean or-
Time evolution of the Kuramoto order parameter where the
der parameter with increasing number of shortcuts at α= 0
arrows indicate transitions between wave patterns through
(circles), -0.97 (squares) and -1.3 (triangles). V=1110 mV,
phase slips.
Cc=82 µF,K=0.10 mSat α=-0.97 andK=0.025 mSat α=-
1.3. (d)Shortdistanceshortcutsatα=0yieldjumpingwaves
(V=1110 mV). (e) Long distance shortcut at α=-0.97 yields occur without apparent external perturbation or change
a frozen complex rotating pattern with a source (triangle) of network topology. The experimentally recorded, ir-
awayfromtheheterogeneity(V=1110mV).(f)Longdistance regulartransientdynamicscontributestothefewexperi-
shortcut at α=0.76 yields a frozen complex rotating pattern
mentalexamplesofhighdimensionaltransientchaos[21],
withasource(fullsquare)ontheheterogeneity(V=1245mV,
where system size effects on the lifetime of the transient
Cind=1 mF, K=0.40 mS).
irregular state could be studied.
This material is based upon work supported by the Na-
tional Science Foundation under Grant Number CHE-
(Fig. 4e) clearly exhibits the transient behavior as its
1465013. Themanuscripthassupplementalmaterial [23]
value changes between approximately 0.2 and 0.6 irreg-
.
ularly throughout the transient until synchronization is
achieved.
In conclusion, we have observed frozen complex rotat-
ing patterns and long transients to synchronization in
electrochemical oscillations and numerical simulations References
of phase oscillators on a ring network topology with
sparse random shortcuts. Depending on the sign of non-
isochronicityα,eitherthesinksorthesourcesarepinned
to endpoints of cross-connections in the network. At
increased non-isochronicity the variability in the phase
differences in a phase locked state decreases until syn- [1] T.Iwasaki, J.Chen,andW.O.Friesen,PNAS111,978
(2014).
chronization is no longer possible and persistent or very
[2] V. M. Lauschke, C. D. Tsiairis, P. Francois, and
long transient phase dynamics occurs. The presence of
A. Aulehla, Nature493, 101 (2012).
long transients of irregular phase dynamics could have [3] G. Ermentrout and D. Kleinfeld, Neuron 29, 33 (2001).
relevanceinthefunctioningofbiologicalsystems,e.g.,in [4] Z.Noszticzius,W.Horsthemke,W.D.McCormick,H.L.
neurondynamicswherepathologicalsynchronizationcan Swinney,and W.Y. Tam, Nature 329, 619 (1987).
1
[5] H.Varela,C.Beta,A.Bonnefont,andK.Krischer,Phys.
Chem. Chem. Phys.7, 2429 (2005).
[6] D. Luss and M. Sheintuch, Catalysis Today 105, 254
(2005).
[7] J. P. Laplante and T. Erneux, J. Phys. Chem. 96, 4931
(1992).
[8] N.Tompkins et al., PNAS 111, 4397 (2014).
[9] G. Ermentrout,J. Math. Biol. 23, 55 (1985).
[10] G.B.Ermentrout,SIAMJ.Appl.Math.52,1665(1992).
[11] N.KopellandG.Ermentrout,Comm.PureAppl.Math.
39, 623 (1986).
[12] D. A. Wiley, S. H. Strogatz, and M. Girvan, Chaos 16,
015103 (2006).
[13] R. Albert and A. L. Barabasi, Rev. Mod. Phys. 74, 47
(2002).
[14] R. Toenjes, N. Masuda, and H. Kori, Chaos 20, 033108
(2010).
[15] P.J. Menck et al., Nat.Commun. 5, 3969 (2014).
[16] P.J. Uhlhaas and W. Singer, Neuron 52, 155 (2006).
[17] Y. Kuramoto, Chemical oscillations, waves and turbu-
lence (Springer, Berlin, 1984).
[18] C. Brito, I. S. Aranson, and H. Chaté, Phys. Rev. Lett.
90, 068301 (2003).
[19] T.-W. Ko and G. B. Ermentrout, Phys. Rev. E 78,
016203 (2008).
[20] M. Wickramasinghe and I. Z. Kiss, Phys. Rev. E. 88,
062911 (2013).
[21] T. Tamás and Y.-C. Lai, Phys. Rep.460, 245 (2008).
[22] M. Sebek, R.Tönjes, and I. Z. Kiss, submitted to Phys- FIG. S1: Schematic of the experiments where the electrodes
ical Review (2015). (blackcircles) areconnectedtotheworkingpoint (W)ofthe
[23] SupplementalMaterial potentiostat. RcareillustratedasemptysquaresandRindare
illustrated by the filled squares. (a) Diagram of the twenty
Supplemental Material electrode ring network with the location of Rc and Rind dis-
played. (b)Simplifieddiagram showingthreeelectrodes with
theadditionofCc. (c)Simplifieddiagramshowingthreeelec-
trodes with theaddition of Cind.
EXPERIMENTAL SETUP
studied in great detail in a previous publication [20].
Experiments areperformedby a standardelectrochemi-
cal cell with 3 M H2SO4 as the electrolyte. The counter
electrodeisaplatinumcoatedtitaniumrod,thereference
METHOD
electrode is Hg/Hg2SO4/sat. K2SO4, and the working
electrodeisanarrayof20nickel,1.00mmdiameter,wire
cross sections embedded in epoxy. The temperature is Once coupled into the ring topology, global negative
o
held at 10 C. Each nickel electrode in the array is con- feedback is applied to linearly increase the phase differ-
nected to the potentiostat (ACM Instruments GillAC) ence between neighboring oscillators according to ϑn =
through an individual resistance Rind. The frequencies 2πn/N. Thephasepatternisstablewhenfeedbackisre-
oftheelectrodesarecontrolledbyanincreaseordecrease moved (see Fig. S2a). Additional connections are added
of Rind and regulated to 0.4Hz ±2.5mHz at Rind =1kΩ randomlyintothesystemwiththesameRc andCc value
in the beginning of the experiment. The system is cou- as the existing connections in the network. If the sys-
pled locally into a ring topology by coupling resistance tem did not synchronize in-phase during the addition
(Rc) as shown in Fig. S1a. A constant potential (V) is of the connection, a second connection is chosen by the
applied and the oscillatory current of each electrode in samemethodas the firstandaddedtothe network. The
the array is measured at 200Hz by a potential drop over process is repeated until in-phase synchronization is ob-
Rind. ThecouplingstrengthK isreportedastheinverse served. Further additions of shortcuts are assumed to
1/Rc of the coupling resistance. Non-isochronicity is in- have no desynchronizing effect. The non-rotating, syn-
troducedintothenetworkbytheadditionofcapacitance chronizedstateisshownasaspace-timeplotinFig.S2b.
C.Thenon-isochronticityisnegativeforacapacitanceCc Transient dynamics to a synchronized state may be ex-
paralleltoRc andpositiveforacapacitanceCind parallel ponentialrelaxation(Fig.S3)orirregularphasedynamics
toRind asseeninFig.S1b,c. The levelofnonisochronic- withperiodsandamplitudesfluctatingintime(Fig.S2c).
ity induced by the capacitancein a pair of ocillatorswas Insynchrony,theperiodisidenticalforalloscillatorsand
2
b
FIG. S2: Current oscillation of electrodes. (a) Space-time
plotoftheelectrodecurrentsontheringnetworkinaphase-
locked rotational wavestate with no shortcus (V= 1105 mV,
K= 0.20 mS, α= 0). (b) Space-time plot of the electrode
currents on the ring network with in-phase synchronization
with a1-11 cross connection. (c) Measured currentsof cross-
FIG. S3: Short transient to synchronization at α = 0 (V=
cut oscillators (1,2,3,6,10,11,13,17,19) showing typical irreg-
1105mV,K=0.033mS).(a)TimeevolutionoftheKuramoto
ular transient behavior (V= 1110 mV, Cc= 82µF K= 0.033
order parameter. (b) Space-time plot of the current during
mS, α= -1.3). (d) Measured currents of cross-cut oscillators
thetransientfromtheinitialrotationalstatetoin-phasesyn-
(1,2,3,6,10,11,13,17,19) exhibiting in-phasesynchronization.
chronization.
late the non-isochronicity (or phase shift parameter):
the amplitudes aremuchmorehomogeneousasshownin
Fig. S2d.
A1
α=arctan (S1)
(cid:18)B1(cid:19)
where A1 is the sine coefficient and B1 is the cosine co-
efficient of the first harmonic of the interaction function
DETERMINATION OF ALPHA
obtained by fast fourier transform. The same procedure
is used for positive α with the product of Rind and Cind
The experimental value of non-isochronicity is deter- held constant.
mined from the interaction function of two coupled os-
cillators with a frequency difference of approximately 20
20 PHASE OSCILLATORS
mHz,achievedbyanRind differenceof100Ω. Inorderto
examine the interaction function, the coupling strength
is lowered from the value of the experiments such that
phaseslipsarepresent. TheproductofRc andCc iskept While the main text [22] provides simulations for sys-
constant to maintain the same non-isochronicity as in tems with N= 500 phase oscillators, simulations with
theexperiments. Theinteractionfunctionisobtainedby N= 20 phase oscillators were performed to demonstrate
plottingthederivativeofthephaseofanoscillator(minus the qualitative agreement between experiments and the
its natural frequency) as a function of the phase differ- phase model. Figure 3 of the main text [22] is repro-
encebetweenthecoupledoscillators. Moredetailsonthe duced with phase oscillators using the same method as
procedureisgiveninapreviouspublication[20]. Theex- intheexperiments. Astablerotationalstate(Fig.S5a)is
perimentallyobtainedinteractionfunctions areshownin found when there are no shortcuts and α=0. The stable
Fig. S4a,b,c. The value of the sine and cosine harmonic state of complete synchronization exists with a connec-
coefficients of the interaction function are used to calcu- tion between oscillator 1 and 11 (Fig. S5b). However,
3
FIG. S4: Experimentally measured interaction functions
(V=1105 mV, Rind1=1 kΩ, Rind2=900 Ω). (a) α=0 (Rc=
20kΩ). (b)α=-0.97(Rc=23kΩ,Cc=35.65µF). (c)α=-1.3
(Rc= 69 kΩ, Cc= 47.54 µF).
in contrast to the experiments this state is not reached
FIG. S5: Dynamicsof N= 20 phaseoscillators. (a) Rotating
from an initial rotating wave. Instead a slight break of
wavephasesimulation,(α=0). (b)Synchronizedstateinthe
the symmetry in the initial rotatingwavewill leadto in-
presenceofa1-11crossconnection. (c)Meanorderparameter
phase synchronization on one side of the ring (jumping
after transient as a function of sequentially added shortcuts.
wave) and a rotating wave on the other part of the ring (d)Jumpingwavephasesimulation, (α=0). (e)Phasesimu-
(Fig. S6). If α is increased, the number of shortcuts re- lation of frozen complex rotating wave pattern with a source
quired to achieve synchrony increases as seen in the plot (triangle) away from the shortcut heterogeneity, (α= -0.97).
ofthemeanorderparameterasafunctionofthenumber (f)Phasesimulation offrozencomplexrotatingwavepattern
with a sink (triangle) away from the shortcut heterogeneity,
of shortcuts shown in Fig. S5c and Fig. S7a. A larger
(α= 0.76).
number of short cuts than in the experiments is needed
to achieve complete synchronizationin the phase model.
We attribute this to the lack of the amplitude degree of transient time as α is increased. At low α in the regime
freedominthe phasemodelwhichcanincreasethe basin of frozen complex wave patterns the transient time me-
of attraction of frozen complex rotating wave patterns. dian is approximately constant but increases by a fac-
Jumpingwaves(Fig S5d)formwiththeadditionofcross- tor of 10 at α = 1.3. This behavior is consistent with
connections over small distances at α= 0. If α is suffi- the experimentally observedlong and irregular transient
ciently negative a wave source and a sink form with the at α = −1.3 which with Ttrans ≈ 1700s ([22] Fig.4e)
source position away from the cross-connection network is about 7 times longer than the exponential relaxation
heterogeneity whereas for sufficiently positive α a wave with Ttrans ≈ 230s for the same network and the same
source and sink form in opposing positions of negative α coupling strength K=0.033 mS at α=0 (Fig. S3a).
(Fig.S5e,f). Thetypicalsystemstatesoffrozencomplex
wave patterns, complete synchronization and irregular
CLUSTERING
phase dynamics shown in Figure 1 in the main text [22]
are also observedwith as few as N=20 phase oscillators.
Fig. S7a,b show the ensemble averaged Kuramoto order Clusteringisobservedfornetworksofidenticalphaseos-
parameterandthe varianceof phasevelocities asa func- cillators with shortcut densities between σ = 0.05 and
tion of the number of short cuts and non-isochronicity 0.2 in a narrow parameter region along the critical line
α for 200 random network realizations. The transient betweentheregionsoffrozencomplexwavepatternsand
timemediantoasynchronizedstatewherethedifference irregularphase dynamics. We havenotobservedcluster-
between the maximal and the minimal phase velocity is ing in the experiments. The reasons for that might be
less than 10−3 in ring networks of N=20 phase oscilla- the small number of oscillators, small heterogeneities in
tors with Nsc=5 short cuts as a function of α is shown the frequencies, amplitude effects or simply not enough
in Fig. S8. Our simulations show a divergence of this observations. ClusteringinthesimulationswithN =500
4
α=0.0 phase oscillators is evident in the cluster order parame-
3
ters shown in in Fig. S9. Figure 1d in the main text [22]
2 showsR6 and R7 ofthe same data with restrictedshort-
cutdensityσ =0.05. InoursimulationsoftheKuramoto
ϑ 1 phase equations on our random network model we have
n observed clustering in networks as large as N = 105 os-
0
cillators and Nsc = 5000 shortcuts. In the clustering
−1 region the cluster number and strength of the clustering
a system may display is distributed with a dependence
−2
on the position on the critical line. Shown in Fig. S9
arethe averageclusterorderparameters. Individualsys-
−3
1 5 n 10 15 20 tems may demonstrate stronger or weaker clustering, a
different cluster number or no clustering at all.
FIG. S6: At α = 0 with a single connection between oscil-
lator 1 and 11 in a ring of N = 20 phase oscillators, and a
slightly perturbed rotating wave initial condition a jumping
waveaccrossonepartoftheringandarotatingwaveaccross
the other part of the ring forms. This configuration is not
observed in theexperimentswhich synchronize completely.
a N=20 b N=20
1.5 1 1.5 4.5 0.3
00..89 4 1.5 R6 1.5 R7 1.5 R8
0.7 ˙ϑ3.5
α1 R 00..56 α1 var23.5 1 1 1 0.25
0.4 2
0.5 0.3 0.5 1.5 α
0.2 1
0.1 0.5
0 0 5 10 15 20 0 0 5 10 15 20 0 0.5 0.5 0.5 0.2
shortcuts shortcuts
Rk
FIG.S7: (a)Orderparameterand(b)varianceofphaseveloc- 00 0.1σ0.2 0.300 0.1σ0.2 0.300 0.1σ0.2 0.3 0.15
iNti=es2a0vperhaagseedosocviellrat2o0r0sraesalaizfautniocntisonoforfaαndaonmdntheetwnourmksbewritohf 1.5 R1 1.5 R5 1.5 R9
shortcuts. Initial conditions are a rotating wave of winding
number one. The data was evaluated after reaching a phase 0.1
locked state or at most t=500 time units. 1 1 1
α
N =20,N =5
sc
1000
0.5 0.5 0.5 0.05
800
nc 600 00 0.1σ0.2 0.300 0.1σ0.2 0.300 0.1σ0.2 0.3
y
s
T 400 Tsync≈300
200 1.3 FIG. S9: k-cluster order parameters Rk - R6, R7, R8, R1,
Tsync≈30 α= R5 and R9 as functions of shortcut density, σ, and non-
0 isochronicity α,forN=500oscillators, averagedover10net-
0 0.5 1 1.5
α work realizations at t= 500 time units.
FIG. S8: Median time to phase locking as a function of α
averagedover200realizations,foranN=20systemwithfive
random shortcuts.