Table Of ContentSynchronization as Aggregation: Cluster Kinetics of Pulse-Coupled Oscillators
Kevin P. O’Keeffe,1 Pavel L. Krapivsky,2 and Steven H. Strogatz3
1Department of Physics, Cornell University, Ithaca, NY 14853, USA
2Center for BioDynamics, Center for Polymer Studies,
and Department of Physics, Boston University, Boston, MA 02215
3Department of Mathematics, Cornell University, Ithaca, NY 14853, USA
(ΩDated: January 20, 2015)
We consider models of identical pulse-coupled oscillators with global interactions. Previous work
showedthatundercertainconditionssuchsystemsalwaysendupinsync,butdidnotquantifyhow
small clusters of synchronized oscillators progressively coalesce into larger ones. Using tools from
thestudyofaggregationphenomena,weobtainexactresultsforthetime-dependentdistributionof
cluster sizes as the system evolves from disorder to synchrony.
5
PACSnumbers: 05.45.Xt,05.70.Ln
1
0
2
In one of the first experiments on firefly synchroniza- transient dynamics leading up to synchrony—the open-
n tion, the biologists John and Elizabeth Buck captured ing and middle game, as opposed to the end game. Ag-
a
hundredsofmalefirefliesalongatidalrivernearBangkok gregation theory offers a new set of tools to explore this
J
and then released them at night, fifty at a time, in their prelude to synchrony.
6
1 darkened hotel room [1]. As they looked on in wonder, We introduce a toy model, which we call scrambler
they observed that “centers of synchrony began to build oscillators. It consists of N (cid:29)1 identical integrate-and-
] up slowly among the fireflies on the wall. In one area fire oscillators coupled all to all. Each oscillator has a
O
we would notice that a pair had begun to pulse in uni- voltage-like state variable x that increases linearly ac-
A
son; in another part of the room a group of three would cording to x˙ = 1, rising from a baseline value of 0 to
n. be flashing together, and so on.” Synchronized groups a threshold value of 1. Whenever any oscillator reaches
i continued to emerge and grow, until as many as a dozen threshold, it fires and does three things. (i) It kicks ev-
l
n fireflies were blinking on and off in concert. The Bucks ery other oscillator (and every cluster of oscillators) to
[ realized that the fireflies were phase shifting each other a new random voltage, independently and uniformly—
1 with their flashes, driving themselves into sync. in this sense, it scrambles the other oscillators. (ii) It
v Here we study stylized models of oscillators akin to then “absorbs” any scrambled oscillators that lie within
5 the fireflies, in which synchrony builds up stepwise, in adistance1/N ofthreshold, bybringingthemtothresh-
1
expandingclusters. Byborrowingtechniquesusedtoan- old and thereby synchronizing with them. To avoid the
1
alyze aggregation phenomena in polymer physics, mate- complicationsthatwouldbecausedbychainreactionsof
4
0 rials science, and related subjects [2, 3], we give the first firings, we assume that the oscillators being brought to
. analytical description of how these synchronized clusters thresholddonotgettofireuntilthenext timetheyreach
1
0 emerge, coalesce, and grow. We hasten to add, however, threshold. (iii) The oscillator that fired resets to x = 0
5 that the models we discuss are not even remotely real- along with the oscillators it absorbed.
1 istic descriptions of fireflies; they are merely intended as If a cluster of j oscillators does the firing, the same
:
v tractable first steps toward understanding how clusters rules apply, except that now any oscillators within a
i evolve en route to synchrony. distance j/N of threshold get absorbed. The assumed
X
Our work is part of a broader interdisciplinary ef- proportionality to j is natural, if each member of the
r
a fort [4, 5]. Oscillators coupled by sudden pulses have cluster contributes to the pulse strength. We study
been used to model sensor networks [6–10], earthquakes otherplausiblecouplingrulesintheSupplementalMate-
[11, 12], economic booms and busts [13], firing neurons rial [26].)
[14, 15], and cardiac pacemaker cells [16]. Diverse forms The motivation for this scrambler model is that it
of collective behavior can occur in these pulse-coupled leads to the simplest possible mean-field approximation.
systems, depending on how the oscillators are connected Intheinfinite-N limit,wewouldlikeclustersofeverysize
inspace. Systemswithlocalcouplingoftendisplaywaves to be uniformly distributed in voltage at all times. This
[17,18]orself-organizedcriticality[11,19,20],withpos- convenient property would greatly ease the derivation of
sible relevance to neural computation [15] and epilepsy therate equationsfor thecluster kinetics. Aswe willsee
[21]. Incontrast,systemswithglobalcoupling,whereev- below, the predictions that follow from this approxima-
ery oscillator interacts equally with every other, tend to tion agree reasonably well with simulations.
fall into perfect synchrony. Rigorous convergence results Assume the initial voltages x , for i = 1,...,N, are
i
have been proven for this case [20, 22–25]. But the tech- independentanduniformlydistributed. Atfirst, nothing
niquesusedpreviouslyhavenotrevealedmuchaboutthe interesting happens. The oscillators increase their volt-
2
ageswithoutinteracting. Butthenoneoscillatorreaches and since the cluster density is c , the corresponding
i
thresholdandfires. Theremainingoscillatorsgetscram- background rate of firing is R0 =c v =c . Next, to find
i i i i
bled, and perhaps some get absorbed. Then another os- theactualR ,wemustsubtractfromR0therateatwhich
i i
cillator fires, and so on. After this process has gone on clusters of size i are being absorbed and hence deprived
for a while, the system has formed synchronized clusters of their chance at firing. Call this absorption rate Ra.
i
of various sizes. Clustersofsizeiareabsorbedwhenclustersofsizej fire,
Let N (t) denote the number of clusters of size j at for j =1,...,N, taking a fraction j/N of the uniformly
j
time t. Thus there are N (t) singleton oscillators, N (t) distributed i-clusters along with them. Since there are
1 2
pairs of synchronized oscillators, N (t) triplets, and so N clustersofsizeiandthej-clustersfireatrateR ,the
3 i j
on. The N are correlated random quantities. They are total rate at which i-clusters are being absorbed is given
j
correlatedbecauseoscillatorsbelongingtoclustersofone by Ra =(cid:80) (j/N)N R =(cid:80) jc R =c (cid:80)jR .
i j i j j i j i j
size are unavailable to clusters of another size, and they Putting all this together gives R = R0 − Ra =
i i i
arerandombecauseoftherandomnessintheinitialcon- c −c (cid:80) jR = c (1−(cid:80) jR ). Let β = 1−(cid:80) jR .
i i j j i j j j j
ditions and the scrambling procedure. It does not seem Note that β is the same for all i, which enables it to be
feasibletounderstandthetime-evolutionoftheNj unless determined self-consistently, as follows. From Ri = βci
theyaresolargethattheirfluctuationsfromonerandom we obtain β = 1−(cid:80) jR = 1−(cid:80)j(βc ). Now invoke
j j j
realization to another are negligible. theidentity(cid:80) jc =j(N /N)=1,whichexpressescon-
j j j
So assume from now on that Nj (cid:29) 1 for all j and servation of oscillators. Solving for β then gives β =1/2
replace these random quantities by their ensemble av- and therefore R =c /2.
i i
erages. Let c = N−1(cid:104)N (cid:105) denote the average clus-
j j Next, plug the expressions derived for R and L into
i i
ter densities. One hopes that relative fluctuations are (cid:80)
the rate equation c˙ = − R L . The result is c˙ =
small; more precisely, N−1N = c + O(N−1/2). An (cid:80) (cid:80) i i i
j j − (c /2)(ic) = −(c/2) ic = −c/2. Recalling that
even stronger assumption is that the densities of dif- i i i i
c(0)=1, we conclude that
ferent sub-populations are asymptotically uncorrelated:
N−2N N =c c +O(N−1/2).
i j i j c(t)=exp(−t/2). (1)
These c allow us to define a natural disorder param-
j
(cid:80)
eter, given by the total density c(t) = jcj(t). It mea- Figure 1 shows this result matches simulations.
sures the extent of the system’s fragmentation. To see
this,notethatatt=0eachoscillatorisalone; onlyclus-
ters of size 1 exist. Accordingly c (0) = 1 and all other
1
c (0) = 0 for j > 1. Hence c(0) = 1, correctly indicat-
j
ingthatthesystemstartsoutmaximallyfragmented. At
the opposite extreme, as t→∞ only one giant cluster of
synchronized oscillators exists. The system is then mini-
mally fragmented: c(t)=1/N →0 as N →∞.
Toderivearateequationforthedeclineofc(t),letR
i
betherateatwhichclustersofsizeifire,fori=1,...,N,
and let L be the number of clusters lost to absorption
i
(cid:80)
in each such firing. Then c˙=− R L .
i i i
TofindL ,recallthatwhenaclusterofsizeifires,all
i
theotherclustersgetassignedanewvoltageuniformlyat
random. Moreover, any clusters assigned to the interval FIG.1: Theoreticalandsimulatedc(t)andc1(t). Solidlines
[1−i/N,1)getbroughttothresholdandabsorbed. Since showtheoreticalcurvesobtainedanalytically(seetext). Data
points show simulation results for N =104 oscillators.
the voltages of these other clusters are uniformly dis-
tributedon[0,1],afractioni/N ofthemwillbeabsorbed.
Thereare(cid:80)jNj clustersintotal. Hencethenumberab- Howdotheindividualclusterdensitiesci behave? To
sorbed is L =(i/N)(cid:80) N =i(cid:80)c =ic. derive their rate equations, note that since the voltage
i j j j
TherateR takesmoreworktocalculate. Sincesome space is the interval [0,1], a segment of length N−1 con-
i
clustersgetabsorbed,noteveryclustergetsthechanceto tainsonaverageNc×N−1 =cclusters. Infact,theprob-
fire. Wemustaccountforthisdepletionwhencalculating ability that it contains n clusters (of any sizes) is given
Ri. Firstconsiderthebackgroundrateoffiringofclusters by the Poisson distribution: Πn =cne−c/n!. This is the
of size i in the absence of absorptions. In other words, mathematical expression of the assumption that clusters
pretendforamomentthatwhenani-clusterfires,itsim- are distributed randomly without correlations.
plyscrambleseveryotherclusterandrestartsitsowncy- With this in mind, let us solve for c (t), the density
1
cle without absorbing anyone. Call this background rate of singletons. It is the easiest c (t) to analyze, since it
j
R0. Since all oscillators move with velocity v = x˙ = 1, can only decrease. Two mechanisms cause the loss of
i i i
3
singletons: (i)asingletonfiresandabsorbsacluster,and Weseefromthesumthattheequations(6)arerecursive.
(ii) a cluster fires and absorbs p≥1 singletons. They can be solved one by one, though not analytically,
Consider mechanism (i). From the expression for R so we resort to numerical integration. Figure 2 shows
i
obtained above, singletons fire at a rate R = c /2. that the theoretical and simulated c agree.
1 1 i
Whentheyfire,theyabsorbanyclusterlyinginthevolt-
age segment [1−1/N,1). According to the Poisson dis-
tribution, the probability that this segment contains one
ormoreclustersisgivenby1−e−c. Hencesingletonsare
lostbymechanism(i)atarate(c /2)(cid:2)1−e−c(t)(cid:3).
1
The loss due to mechanism (ii) is handled similarly.
Suppose p singletons lie the interval [1−j/N,1) at the
moment when a cluster of size j fires, for j = 1,...,N.
This event happens with probability e−jc1(jc1)p/p!, and
when it does, it consumes p singletons. The j-clusters
fire at a rate R = c /2. Hence singletons are lost by
j j
mechanism (ii) at a rate
(cid:88)cj ×(cid:88)p(jc1)pe−jc1 =c /2. (2)
2 p! 1 FIG. 2: Theoretical and simulated cluster densities c2(t)
j≥1 p≥1 through c (t). Solid lines show theoretical predictions com-
5
Summing the loss rates from (i) and (ii) gives putedfromnumericalintegrationofEq.(6). Datapointsshow
simulation results for N =5×104 oscillators.
dc c
1 =− 1(2−e−c(t)). (3)
dt 2 Although we cannot find all the c (t) explicitly, we
i
This equation has a closed-form solution in terms of ex- can get their moments M (t) = (cid:80) jnc (t). The first
n j j
(cid:80)
ponential integrals: two moments are trivial: M = c (t) = c(t) and
0 j j
(cid:80)
M (t) = jc (t) = 1, from conservation of oscillators.
c (t)=exp(−t+Ei(−1)−Ei(−e−t/2)), (4) 1 j j
1 The higher moments can be obtained from a generating
where we have used the initial condition c (0)=1. Fig- function. Let G(z,t)=(cid:80) c (t)ekz. Then the infinite
1 k≥1 k
ure 1 shows good agreement between the theoretical and set of differential equations (6) transforms into
numerical c (t).
1
∂G(z,t) 1
Fori>1,therateequationforc includesgainterms
i +G(z,t)= G[z−c(t)+G(z,t),t]. (7)
aswellaslossterms. Clustersofsizei>1canbecreated ∂t 2
when two or more smaller clusters coalesce, or destroyed
Thisequationlooksneat,butitisfarfromtrivial,asthe
when they themselves coalesce with at least one other
right-hand side involves G in a very nonlinear manner.
cluster. Thelosstermisastraightforwardgeneralization Using the identity M (t) = ∂nG| , we can however
of that for c , and is given by (c /2)(cid:2)2−e−ic(t)(cid:3). n ∂zn z=0
1 i derive the following equations for the moments:
To find the gain term, imagine that a cluster of size
k fires. The segment [1−k/N,1) may contain a1 clus- M˙ = 3M
ters of size 1, a clusters of size 2, etc. This event hap- 2 2 2
2
pareensuswinitghtphreoabsasbuimlitpyti(oknca11t)!ah1at×cl(ukcsa2t2)e!ar2s×of.d.i.ffe(wrehnetresizwees M˙3 = 27M3+3(M2)2 (8)
15 3
are independent as well as Poisson distributed). If the M˙ = M +16M M + (M )3.
segment contains a combination of clusters such that 4 4 4 2 3 2 2
k + a + 2a + 3a + ··· = i, then a cluster of size i
1 2 3 Like the c equations (6), these moment equations are
i
will form. We sum over all such combinations for a clus-
recursiveandcanbesolvedinsuccession. Wefind
ter of size k firing, and then sum over all k, to get the
rate at which clusters of size i are created: M (t)=e3t/2
2
(cid:88)i−1 cke−kc (cid:88) (cid:89) (kcp)ap. (5) M3(t)=74e478t/2−6e3t 217 (9)
k=1 2 a1+2a2+···=i−k p≥1 ap! M4(t)= 5 e5t−128e9t/2+ 5 e15t/4−4e27t/8.
Combining the loss and gain terms, and transferring
Figure 3 plots theoretical and simulated values of the
c e−ic into the gain term, we finally obtain
i M . TheagreementisgoodforM butworseforM and
i 2 3
c˙i =−ci+k(cid:88)=i1c2ke−kc(cid:80)p(cid:88)ap=i−kp(cid:89)≥1(kacpp)!ap. (6) Mpforro4m.ceasTtsih,oinssuobisfjebtcoitgbtcoeluflesutxecprtesuc.attSeidion.ncseEdMaocmnh(itnc)ia(=tte)d(cid:80)isbjyajtnshctejo(ccthh)a,asntthcieec
4
systems with local coupling, network structure, hetero-
geneity, delays, and other realistic features.
[1] J. Buck and E. Buck, Science 159, 1319 (1968).
[2] F. Leyvraz, Physics Reports 383, 95 (2003).
[3] P.L.Krapivsky,S.Redner,andE.Ben-Naim,AKinetic
ViewofStatisticalPhysics (CambridgeUniversityPress,
Cambridge, UK, 2010).
[4] A. Pikovsky, M. Rosenblum, and J. Kurths, Synchro-
nization (Cambridge University Press, Cambridge, UK,
2003).
FIG. 3: (Color online) Log plot of the first three nontrivial [5] S. Strogatz, Sync (Hyperion, New York, 2003).
moments M ,M ,M . Black curves, theoretical predictions [6] Y.-W.HongandA.Scaglione,in2003 IEEE Conference
2 3 4
obtained from (8); red dots, average simulation results for on Ultra Wideband Systems and Technologies (IEEE,
100 realizations of N =104 oscillators. 2003), pp. 190–194.
[7] X.Y.WangandA.B.Apsel,in50thMidwestSymposium
on Circuits and Systems, 2007 (IEEE, 2007), pp. 1524–
higher moments amplify these fluctuations more and are 1527.
therefore noisier themselves. [8] G. Werner-Allen, G. Tewari, A. Patel, M. Welsh, and
Variants of the scrambler model can also be solved. R. Nagpal, in Proceedings of the 3rd International Con-
ference on Embedded Networked Sensor Systems (ACM,
For example, suppose that when a cluster of size j fires,
2005), pp. 142–153.
itabsorbsalloscillatorswithinadistancekj/N ofthresh-
[9] J. Nishimura and E. J. Friedman, Phys. Rev. Lett. 106,
old, where k >0 is a tunable coupling strength. Or sup-
194101 (2011).
pose that the pulse strength is independent of the size [10] J. Nishimura and E. J. Friedman, Phys. Rev. E 86,
j of the firing cluster. Both these cases are discussed in 025201 (2012).
[26]. Nothing qualitatively new happens, but there are [11] A. V. M. Herz and J. J. Hopfield, Phys. Rev. Lett. 75,
some interesting quantitative differences, e.g., in the lat- 1222 (1995).
[12] S.BottaniandB.Delamotte,PhysicaD103,430(1997).
ter case the disorder parameter c(t) decays algebraically
[13] S. Gualdi, J.-P. Bouchaud, G. Cencetti, M. Tarzia, and
rather than exponentially.
F. Zamponi, ArXiv e-prints (2014), 1409.3296.
Amoresubstantialmodificationwouldbetomakethe
[14] C.Kirst,T.Geisel,andM.Timme,Phys.Rev.Lett.102,
model deterministic. Scrambler oscillators respond ran- 068101 (2009).
domly, but real oscillators respond deterministically, ac- [15] J. J. Hopfield and A. V. M. Herz, Proc. Natl. Acad. Sci.
cording to a phase-response curve [27]. This determinis- 92, 6655 (1995).
ticcharacterisincorporatedinRefs.[16,22]byassuming [16] C. S. Peskin, Mathematical Aspects of Heart Physiology
(CourantInstituteofMathematicalSciences,NewYork,
that when a cluster of size j fires, it adds a voltage pulse
1975), pp. 268–278.
j(cid:15) to every other oscillator, or pulls it up to threshold,
[17] P.C.Bressloff,S.Coombes,andB.deSouza,Phys.Rev.
whichever is less. For the case where (cid:15) = 1/N we show
Lett. 79, 2791 (1997).
in[26]thatthisdeterministicsystemcanalsobewellap- [18] P.C.BressloffandS.Coombes,Phys.Rev.Lett.80,4815
proximatedbytheaggregationmethodsusedabove. The (1998).
main new feature is that c(t) and the cluster densities [19] A. Corral, C. J. P´erez, A. Diaz-Guilera, and A. Arenas,
c (t) become piecewise linear. But their overall shapes Phys. Rev. Lett. 74, 118 (1995).
i
[20] S. Bottani, Phys. Rev. Lett. 74, 4189 (1995).
still resemble those seen in the scrambler model.
[21] G. A. Worrell, S. D. Cranstoun, J. Echauz, and B. Litt,
Therearemanyavenuestoexploreinfuturework. For
NeuroReport 13, 2017 (2002).
example, throughout our analysis each oscillator obeyed
[22] R. E. Mirollo and S. H. Strogatz, SIAM J. Appl. Math.
x˙i = 1. Such linear sawtooth waveforms are reasonable 50, 1645 (1990).
fortheoscillatorsusedinsensornetworks[8],butnotfor [23] W. Gerstner and J. L. van Hemmen, Phys. Rev. Lett.
neurons or cardiac pacemaker cells. Changing the con- 71, 312 (1993).
cavityorshapeofthewaveformwouldmaketheanalysis [24] U. Ernst, K. Pawelzik, and T. Geisel, Phys. Rev. Lett.
74, 1570 (1995).
more difficult, and might even affect the cluster kinetics,
[25] S. Bottani, Phys. Rev. E 54, 2334 (1996).
becauseclusterswouldnolongerbeuniformlydistributed
[26] See Supplemental Material at
aswehaveassumedthroughout. Coulddifferenttypesof
http://link.aps.org/supplemental/ for results for other
pulses be engineered to preserve a uniform distribution? coupling rules.
Or perhaps different pulses would distribute the clusters [27] A. T. Winfree, The Geometry of Biological Time
in a nonuniform but still tractable manner. It would (Springer, New York, 2001).
also be interesting to study cluster kinetics in oscillator
