A Coupled Oscillator Model for the Origin of Bimodality and Multimodality
AA coupled oscillator model for the origin of bimodality and multimodality
J.D. Johnson and D.M. Abrams Department of Engineering Sciences and Applied Mathematics, Northwestern UniversityMcCormick School of Engineering and Applied Science2145 Sheridan Road Evanston, IL 60208 (Dated: 25 February 2020)
Perhaps because of the elegance of the central limit theorem, it is often assumed that distributions in nature will approachsingly-peaked, unimodal shapes reminiscent of the Gaussian normal distribution. However, many systems behavedifferently, with variables following apparently bimodal or multimodal distributions. Here we argue that multimodalitymay emerge naturally as a result of repulsive or inhibitory coupling dynamics, and we show rigorously how it emergesfor a broad class of coupling functions in variants of the paradigmatic Kuramoto model.
In this paper we employ oscillators as a test system forunderstanding how bimodality—the splitting of oscillatorsinto two rather than one cluster—may emerge as a resultof coupling between interacting units. We present numer-ical and analytical results showing that repulsive couplingcan lead to bimodality (or multimodality) for a wide rangeof detailed interaction dynamics.
I. INTRODUCTION
Synchronization is a widespread phenomenon observedin biological, chemical, physical, and socialsettings.
A paradigmatic mathematical model that canexplain synchronization in many contexts is the Kuramotomodel.
Much work has been done on understanding thecomplex and surprising dynamics of the Kuramoto model andits variants, but the vast majority of that research focuses onthe case of attractive coupling; here we are interested in thecase where the coupling is repulsive.Repulsive (or inhibitory) coupling is of physical interest asit arises frequently in the context of neuronal networks (e.g.,see refs. 2 and 20), chemical interactions (e.g., refs. 4, 21, and22), and many other systems (see refs. 23–28). Some coupledoscillator models have examined repulsive coupling: Giveret al. developed a local variant of the Kuramoto model withrepulsive coupling based on the interaction between watermicro-droplets with reactants of the Belousov-Zhabotinskyreaction. Hong and Strogatz developed two variants of theKuramoto model that involved mixes of positive and negativecoupling.
The relationship between network structure and repulsivecoupling has also been analyzed, with Levnaji´c show-ing that, given the network coupling structure, many differentphase configurations can arise. Recently, it has been shownthat synchronization can arise in both repulsive and attrac-tive coupling scenarios subject to common noise.
Gonget al., inspired by the work of Gil et al., studied instanceswhere common noise can lead to clustering in the phase dis-tribution of oscillators for repulsive coupling.Nakamura et al. investigated the effect of time-delayednearest-neighbor coupling in the Kuramoto model and foundthat it could lead to the development of clustered states forboth attractive and repulsive coupling. Mishra et al. demon- Beetle Horn Size P D F Galaxy Color P D F
600 700 800 900 1000 1100
Appearance Time of Colonies
Atlantic Salmon Weight (a)(c) (d)(b)
FIG. 1.
Selected examples of bimodality.
Histograms (normalized)for (a) size of beetle horns [mm], (b) Atlantic salmon body mass[g] (c) color of galaxies at redshift 0.1 (d) inverse growthrates of bacteria [min − ]. strated that “chimeralike” states could arise with globally cou-pled Liénard systems incorporating both attractive and re-pulsive mean-field feedback. Yeldesbay et al. establishedthat chimeralike states can arise in the Kuramoto-Sakaguchimodel. They also considered a model with oscillators thatcould be synchronous (attractive coupling) or asynchronous(repulsive coupling) depending on their natural frequencies.They found that in this case a chimera state arises.Golomb et al. showed that clustering is possible in a cou-pled oscillator model with repulsive coupling that is suited forstrong interactions between the limit-cycle oscillators. Theyfurther provided theory for when a frequency locked station-ary phase distribution and when a nonperiodic attractor canarise.Tsimring et al. showed that heterogeneous globally cou-pled oscillators obeying the standard Kuramoto model cancluster with all configurations having a zero order parameter,but this clustering breaks down as the number of oscillatorsincreases. They also showed that, with local coupling, clus-tering can occur for nonidentical oscillators given sufficientlylarge coupling strength.Closest to the work we present here, Okuda looked at theeffect that an arbitrary coupling function may have on oscilla-tors and developed theory as to when an n -cluster state, withall clusters being the same size, can arise. He found that har-monics in the coupling function are necessary for clusters toarise. a r X i v : . [ n li n . AO ] F e b The central limit theorem may influence us to ex-pect that distributions in nature should tend to a singly-peaked, unimodal shape akin to the Gaussian normal distri-bution. Yet, bimodality and multimodality can be observedin biological, social, and chemical contexts andbeyond (see Fig. 1 for selected examples). In this paperwe demonstrate that multimodality may arise as a result ofrepulsive or inhibitory coupling dynamics and we give an in-depth explanation of how it can arise for a range of couplingfunctions. II. MODEL WITH ANTISYMMETRIC REPULSIVECOUPLING
We begin by considering a system of N phase oscillatorscharacterized by natural frequencies ω i , i = . . . N . The oscil-lators are globally coupled with coupling strength K throughan interaction function f that depends only on the phase dif-ference between each pair of oscillators:˙ θ i = ω i + KN N ∑ i = f ( θ j − θ i ) , i = , . . . , N . (1)Here K > K < f ( u ) , u ∈ ( − π , π ] , thatsatisfy the following conditions: f ( ) = f (cid:48) ( ) > f ( u ) = − f ( − u ) (2c) f (cid:48) ( u ) continuous (2d) f ( π ) = lim u →− π + f ( u ) . (2e)These conditions impose: (2a) no coupling effects betweenoscillators in sync; (2b) locally attractive (repulsive) couplingnear sync state for K > K < f (cid:48) ( u ) ; (2e) 2 π -periodic interac-tion function on ( − π , π ] domain. We point out that conditions ( c ) and ( e ) lead to f ( π ) = lim u →− π + f ( u ) = A. Identical Oscillators
We assume that oscillators frequencies are drawn from aknown frequency distribution g ( ω ) . For simplicity we firstconsider the case of identical oscillators, i.e., we set the distri-bution to be g ( ω ) = δ ( ω − ω ) , so the system becomes˙ θ i = ω + KN N ∑ j = f ( θ j − θ i ) , i = , . . . , N . (3) B. Bimodal equilibria
We assume that the number of oscillators is large, N (cid:29) − π − π π π − π − π π π f ( ψ ) ψ FIG. 2.
Sample interaction functions.
Two cases of coupling func-tions that we consider. Case 1 (red, dashed) is an odd, 2 π -periodicfunction with a continuous derivative, no zeros in between 0 and π ,and has a positive slope at 0. Case 2 (blue, solid) is similar to case 1but has a zero of order 1 in between 0 and π . of an oscillator phase distribution h ( θ ) = x δ ( θ − θ ) + ( − x ) δ ( θ − θ ) , where 0 < x < θ = ω + KN (cid:32) xN ∑ i = f ( θ − θ ) + N ∑ i = xN + f ( θ − θ ) (cid:33) = ω + K ( − x ) f ( θ − θ ) (4)˙ θ = ω + KN (cid:32) xN ∑ i = f ( θ − θ ) + N ∑ i = xN + f ( θ − θ ) (cid:33) = ω − Kx f ( θ − θ ) . (5)We define a new phase-difference variable ψ = θ − θ andwrite its dynamical system by subtracting Eq. (4) from Eq. (5):˙ ψ = − K f ( ψ ) . (6)We observe that the fixed points of the system for ψ are fullydetermined by the zeros of f ( ψ ) . From the assumptions above f ( ψ ) must have zeros at ψ = ψ = π . Furthermore, ifconditions (2a–2e) hold and f ( ψ ) has no other zeros (as in thecase of the red dashed curve from Fig. 2), then it is impliedthat f (cid:48) ( π ) ≤
0. Hence, within the bimodal manifold, the fixedpoint at ψ = π should be stable with ψ = ψ = π corresponds to a bimodal equilibrium with two clustersof oscillators separated by 180 ◦ of phase.If additional roots of f ( ψ ) exist between 0 and π , thesewill also correspond to bimodal fixed points with alternatingstability (again restricted to the bimodal manifold). We focuson the cases where there are no other fixed points or there isexactly one other fixed point ψ in ( , π ) ; other cases are sim-ilarly tractable. Figure 2 illustrates the typical general shapesof the interaction functions that we consider.
1. Stability of bimodal equilibrium
To investigate the broader stability of solutions to pertur-bations outside the bimodal manifold, we consider the pertur-bation of a single oscillator by a small amount ε . Because N (cid:29)
1, we approximate the dynamics of the two clusters asunaffected by this perturbation. We examine the evolution ofdistance between the perturbed oscillator and the group fromwhich it was perturbed, ε ( t ) , to evaluate whether the systemreturns to its initial state.For convenience, we move into a rotating frame by re-defining θ i → θ i + ω t , which is equivalent to setting ω = N fromthe θ cluster for the perturbation and assume θ =
0, andthus θ = ψ ≤ π (assuming for now that our interactionfunction has only one or zero fixed points in ( , π ) ). Then θ N = θ − ε = ψ − ε , and˙ ε = − KN (cid:34) xN ∑ i = f ( θ − ψ + ε ) + N − ∑ i = xN + f ( θ − ψ + ε ) (cid:35) = − Kx f ( − ψ + ε ) − K ( − x ) f ( ε ) . We expand the functions f in a Taylor series to linear order:˙ ε ≈ − ε K (cid:2) x f (cid:48) ( ψ ) + ( − x ) f (cid:48) ( ) (cid:3) . Assuming that K < x f (cid:48) ( ψ ) + ( − x ) f (cid:48) ( ) < . (7)A nearly identical calculation starting with the perturbationof a single oscillator from the θ (zero phase) cluster leads toa similar equation, ( − x ) f (cid:48) ( ψ ) + x f (cid:48) ( ) < . (8)Since Eqs. (7) and (8) must be simultaneously satisfied for sta-bility of the full bimodal distribution, the following inequalitymust hold: f (cid:48) ( ) < ( − x )[ f (cid:48) ( ) − f (cid:48) ( ψ )] < − f (cid:48) ( ψ ) . (9)Interestingly, this implies that the slope of the interactionfunction f ( ψ ) must be steeper at ψ = ψ compared to ψ = x in inequality ( ) : f (cid:48) ( ) f (cid:48) ( ) − f (cid:48) ( ψ ) < x < − f (cid:48) ( ψ ) f (cid:48) ( ) − f (cid:48) ( ψ ) . (10) III. CONCRETE EXAMPLE
As a concrete example, we consider a simple class of inter-action functions f ( u ; a ) = π a u (cid:0) π − u (cid:1) (cid:0) a − u (cid:1) . (11) − π − π π π − − f ( ψ ) ψ FIG. 3.
Concrete interaction function.
The interaction functiondefined in Eq. (11) plotted for several different values of a : √ π / √ π / π / | a | approaches π the slope at zero stays fixed with slope 1and the slope at ± a decreases in magnitude. This relation between a and the slope values at ± a , combined with Eq. (10) leads to thethreshold for bimodality given by Eq. (12). These functions have roots on ( − π , π ] at 0, π , and ± a , andsatisfy all the conditions set forth earlier in section II. As longas 0 < | a | < π there are three roots in 0 ≤ u ≤ π , and one cancheck that f (cid:48) ( ) = a (see Fig. 3 for exampleplots). For inequality (10) to be satisfiable, we require π π − a < π − a π − a , which reduces to | a | < π / √ ≡ a crit . (12)We note that symmetry of the roots allows us to considerpositive a without loss of generality. Figure 4 shows the re-sults of numerical experiments where we test this predictedstability threshold. In each panel, Eq. (3) is implementedwith the interaction function from Eq. (11). We initialize xN oscillators at θ = ( − x ) N at θ = a , then add asmall random perturbation ξ i to each oscillator’s initial phase,where ξ i is drawn from the normal distribution N ( , δ ) ,with δ = . ψ = θ − θ → a with x final = x initial . We note that in theseexperiments we set coupling strength K = − N =
100 oscillators, ω = a = π /
2, consistent with the stability threshold fromEq. (12), a < a crit = π / √
2. The stable band of fractionationaccording to inequality (10) is then 2 / < x < /
5. In panel(a), we set x initial = .
55, below the band’s upper bound; inpanel (b), we set x initial = .
65, above the band’s upper bound.As expected, the bimodal equilibrium appears stable in panel(a), but unstable in panel (b), where eleven oscillators movebetween clusters to establish a different equilibrium within thestable fractionation band (2 / < x final = . < / -0 P ha s e Time -0 P ha s e Time(c) (d)(b)(a)
FIG. 4.
Numerical experiments with identical oscillators.
Usingexample from Eq. (11), top two panels show test for stability range offractionation x from Eq. (10); bottom two panels show test for criticalparameter a crit from Eq. (12). (a) When initial fractionation is in thestable range (here 0 . < x initial = . < .
6) perturbations shrink andthe solution returns to its initial state. (b) When initial fractionationis outside stable band (here x initial = . > .
6) perturbations growfor some oscillators until system evolves to a different fractionationstate. (c) When x initial = / a < a crit , perturbations shrink andthe solution returns to its initial state. (d) When x initial = / a > a crit , perturbations grow and the system moves away from theunstable bimodal state until it reaches a new trimodal equilibrium. In panels (c) and (d), we again use N =
100 oscillators and ω =
0, but here we examine the predicted stability thresh-old a crit = π / √ ψ ∗ = a to be unstable for all positive a > a crit (but note that this state ceases to exist when a > π ). Weset x initial = / a < a crit . In panel (c), we set a = a crit − .
1, just below the threshold for stability; in panel(d), we set a = a crit + .
1, just barely in the unstable domain.As expected, the bimodal equilibrium again appears stable inpanel (c), but it appears unstable in panel (d). Since no frac-tionation x will lead to a stable bimodal equilibrium, the sys-tem must move to an entirely different state, and it appears toconverge to a trimodal distribution of oscillator phases.We are able to understand why the system converges to atrimodal state by performing a similar analysis for the stabilityof three-cluster, or trimodal, oscillator distributions. One canshow that a necessary condition for stability is: f (cid:48) ( ) < − (cid:2) ( x + y ) f (cid:48) ( ψ ) + ( y + z ) f (cid:48) ( ψ )+( x + z ) f (cid:48) ( ψ + ψ ) (cid:3) (13)where ψ i is the angle separating clusters at θ i and θ i + ( θ identified with θ ), and x , y , and z are the fractionations ofthe three clusters at θ , θ , and θ respectively. With equalspacing between the clusters ψ = ψ = π − ψ − ψ = π / f (cid:48) ( ) < − f (cid:48) ( π ) . − π − π π π − π − π π π f ( ψ ) ψ FIG. 5.
Sample asymmetric interaction function.
This function(solid blue curve) does not satisfy f ( ψ ) = − f ( − ψ ) . Existence ofbimodal equilibria requires that it intersect its mirror reflection (dot-ted blue curve) or a scaled version of it (see Eq. (15)). The fixedpoints of the system for x = / For the example function shown in Eq. (11) this is a < (cid:114) π ≡ a tricrit ≈ . π . This implies that a trimodal state remains stable for all a < π .It stably coexists with the bimodal state for a < π / √
2, andmay coexist with other multimodal states for π / √ < a < π . In general different multimodal states may stably coexistover various parameter ranges. More details of the analysisfor trimodality can be found in the appendix. IV. GENERALIZATION TO ASYMMETRICINTERACTION FUNCTIONS
We can relax assumption (2c) of an antisymmetric cou-pling function and still find stability boundaries for multi-modal states. In place of Eq. (6) (which used oddness of thecoupling function), we find instead˙ ψ = Kx f ( − ψ ) − K ( − x ) f ( ψ ) . (14)Clearly ψ ∗ = ψ ∗ = π both remain fixed points. Otherfixed points exist if x f ( − ψ ∗ ) = ( − x ) f ( ψ ∗ ) (15)has a solution on − π < ψ ∗ ≤ π . Figure 5 shows an example ofan asymmetric interaction function. Geometrically this condi-tion can be understood as identifying intersections of f ( ψ ) and its reflection f ( − ψ ) when x = / x (cid:54) = / (a)(c) (d)(b) FIG. 6.
Numerical experiments with heterogeneous oscillators.
Here, N = N ( , ) and the perturbations, ξ i , i = ... N , are drawnfrom N ( , . ) . Using example from Eq. (11), panels (a) and (b)show the results for x initial = / a = π / < a crit (compare toFig. 4(a)). Panels (c) and (d) show the results for x initial = / a = π / √ + . > a crit (compare to Fig. 4(d)). V. GENERALIZATION TO NON-IDENTICALOSCILLATORS
We argue that real-world bimodal or multimodal distribu-tions may result from similar dynamics to those presented inthis paper. Of course, heterogeneity is inevitable in most real-world systems, yet we have focused thus far on the case ofidentical oscillators. While we leave the more general analy-sis for future work, we have conducted numerical experimentsthat appear to show that the predicted behavior occurs even inthe presence of oscillator heterogeneity.Again using the same example interaction function fromEq. (11), we now draw frequencies, ω i , from a normal distri-bution N ( , σ ) and set the initial phases of the oscillators to θ i = ξ i (fraction x ) or θ i = a + ξ i (fraction 1 − x ), where ξ i is asmall perturbation draw from the distribution N ( , δ ) . Fig-ure 6 shows the results of perturbation experiments analogousto those presented in Fig. 4, with analogous results except thatthe final phase distributions have phases that cluster about themodes rather than all converging to them precisely (right pan-els show histograms of final states).In Fig. 6 panels (a) and (b), we use N = a = π / < a crit and x initial = /
2. Even with perturbed ini-tial phases and heterogeneous natural frequencies, the oscilla-tors still remain in the bimodal state as predicted for a < a crit .Specifically, panel (b) shows that the steady state distributionof oscillators has finite-width clustering about the fixed pointpositions predicted from the identical-oscillator case. In pan-els (c) and (d), since a = π / √ + . > a crit , the bimodal statebreaks down (consistent with the prediction of the identical-oscillator theory) and the system appears to converge to a tri-modal equilibrium with three finite-width clusters. VI. DISCUSSION
Coupled oscillators are an excellent testbed for mod-els of synchronization or clustering. Even though real-world variables (e.g., sediment grain size, salmon bodysize, human communication frequency, dopamine re-ceptor density, neutron star mass, galaxy color, gammaray burst duration, tree height, animal ornament size )may not be oscillatory or confined to a periodic domain, bi-modality may emerge for qualitatively similar reasons. In ourmodel, the coupling of one unit’s dynamical behavior to thatof others is key to the phenomenon.For clarity of presentation we have focused on a single ex-ample of interaction function (Eq. (11)), but evaluation of twoother classes of interaction functions (triangle waves and an-tisymmetrized von-Mises kernels) also supports our analyti-cal results—see supplementary material for details. In sup-plementary material, we also present further results regardingdependence of bimodal equilibria on coupling strength K , aswell as some numerical evidence regarding sizes of basins ofattraction; each of these topics merits further in-depth study.The analysis we present here focuses exclusively on the caseof all-to-all coupling; we leave further investigation of the im-pact of network structure for future work.For real-world scenarios where bimodality or multimodal-ity is of interest, the interaction function may not be knownexactly. Nevertheless, we expect that it will often be possibleto assess whether the conditions expressed in Eqns. (2) holdin a particular case. It also seems plausible that functions de-scribing real-world interactions between coupled systems willhave no more than a handful of roots, making bimodality andtrimodality likely outcomes when repulsive or inhibitory cou-pling is imposed.One particularly important case occurs when the interactionfunction has only roots at zero and π , with the root at zerohaving larger or equal magnitude slope. That is the case in thestandard Kuramoto model with sinusoidal coupling. In such acase we expect that the incoherent splay state will be stable.In general, the splay state should be stable when the tendencyto cluster (due to long-distance interactions) cannot overcomethe oscillators’ locally repulsive interactions. VII. CONCLUSIONS
We have shown that, when coupling is repulsive, multi-modality of the oscillator distribution can be a stable con-figuration for a wide range of interaction functions. Weshowed that bimodality can be expected under repulsive cou-pling when the slope of the interaction function at the originis shallower than at the other root(s). We performed numeri-cal experiments for both identical and nonidentical oscillatorsand observed results consistent with theory.This demonstration that repulsive coupling can produceclustering under reasonable assumptions about the interac-tion dynamics is important as repulsive coupling is presentin many natural systems. Hence, the theory we present inthis paper provides an argument as to why one might expectmulti-modality instead of unimodality or incoherence in sys-tems known to have repulsive coupling.
SUPPLEMENTARY MATERIAL
See accompanying supplementary material for numericalexperiments using a selection of additional interaction func-tions, for discussion of basins of attraction for different states,and for discussion of coupling strength dependency.
Appendix: Trimodal equilibria
We again consider a function f ( u ) that satisfies condi-tions (2a)–(2e). We look for solutions with oscillators dis-tributed according to h ( θ ) = x δ ( θ − θ ) + y δ ( θ − θ ) + ( − x − y ) δ ( θ − θ ) , where x , y > x + y < θ , θ , and θ (we again assumethat the natural frequencies are identical):˙ θ = ω + K ( y f ( θ − θ ) + z f ( θ − θ )) (A.1)˙ θ = ω + K ( − x f ( θ − θ ) + z f ( θ − θ )) (A.2)˙ θ = ω − K ( x f ( θ − θ ) + y f ( θ − θ )) . (A.3)Here, z = − x − y . We define two variables ψ = θ − θ and ψ = θ − θ , so that the system reduces to˙ ψ = − K ( z [ f ( ψ + ψ ) − f ( ψ )] + ( x + y ) f ( ψ )) (A.4)˙ ψ = − K ( x [ f ( ψ + ψ ) − f ( ψ )] + ( y + z ) f ( ψ )) . (A.5)We set ˙ ψ i = i = , f ( ψ ) = x f ( ψ ) z (A.6) f ( ψ + ψ ) = − y f ( ψ ) z . (A.7)To set bounds on the fractionation of the clusters, we assumethat there exists points ψ , ψ ∈ ( − π , π ) such that Eqns. (A.6)and (A.7) are satisfied. Additionally, we put our system ofcoupled oscillators into a rotating frame so that θ i → θ i + ω t .In the rotating frame, we set θ = θ = ψ , and θ = ψ + ψ − π . As before we perturb an oscillator from one of thethree groups. We do this for all three groups and get a systemof inequalities x f (cid:48) ( ) + y f (cid:48) ( ψ ) + z f (cid:48) ( ψ + ψ ) < y f (cid:48) ( ) + x f (cid:48) ( ψ ) + z f (cid:48) ( ψ ) < z f (cid:48) ( ) + y f (cid:48) ( ψ ) + x f (cid:48) ( ψ + ψ ) < . (A.8c)All these must be simultaneously satisfied for stability of atrimodal state. Adding, we find f (cid:48) ( ) < − (cid:2) ( x + y ) f (cid:48) ( ψ ) + ( y + z ) f (cid:48) ( ψ )+( x + z ) f (cid:48) ( ψ + ψ ) (cid:3) . (A.9) This states that the weighted sum of the slopes of the couplingfunction at ψ = ψ i , where the weights are the proportions forthe groups separated by ψ i , is greater in magnitude than theslope at the origin. This condition reduces to f (cid:48) ( ) < − f (cid:48) (cid:18) π (cid:19) (A.10)if ψ = ψ = π − ψ − ψ = π /
3. As an example, we re-turn to the class of interaction functions that we introduced inSection II. We relax the assumption that | a | < π and considerthe case when ψ = ψ = π − ψ − ψ . To satisfy inequality ( A . ) , this means that1 < π + a a , (A.11)which reduces to | a | < (cid:114) π ≡ a tricrit ≈ . π . (A.12)Figure 7 shows the results of a numerical experiment wherewe test this threshold. In both panels we use N = x = y = z = / ψ = ψ = π − ψ − ψ = π /
3, and set ω = a > a tricrit .In panel (a) we set a = a tricrit + . ξ i , with values drawn from the distribution N ( , . ) . We can see that this perturbation leads to thesystem leaving the trimodal state and going to a bimodal statewith 180 ◦ phase difference.One might be interested in why the bimodal state is stable inpanel (a). Since there are only zeros at ψ = ψ = π , onemay check the stability by evaluating the derivative of f ( ψ ) at these points. One can show that if | a | > √ π (A.13)the 180 ◦ antiphase state is stable. Thus, when a = a tricrit + . > √ π the trimodal state becomes unstable and perturba-tions lead to the stable bimodal state.For the case, when √ π < a < a tricrit , both the trimodalstate and the bimodal state are stable configurations. Figure 8shows the result of the numerical experiment where we placethe parameter inside the previously stated interval and outsideof the interval. In all panels we use N = ω = ξ i from the predicted fixed points, whose values are drawnfrom the distribution N ( , . ) . In panels (a) and (b) weset a = . π ∈ ( √ π , a tricrit ) . In these cases we expect boththe bimodal state and the trimodal state to be stable for thisvalue of a . In panels (a) and (b), we set the fractionation to beequal in all groups, and we set the spacing between groups tobe equal. As expected, we see that the trimodal state and thebimodal state are stable under perturbation.In panel (c) we set a = . π − . < √ π < a crit . As ex-pected, we see that the bimodal state is unstable and the sys-tem goes in to trimodal state. Given the proximity of the clus-ters to ± π , we have added black dashed lines that at ± π , sothat one can see that the difference between the final state and -0 P ha s e Time -0 P ha s e (b)(a) FIG. 7.
Numerical experiments testing the threshold for tri-modality.
Panel (a): parameter value is a = a tricrit + .
1, and thetrimodal state appears to be unstable (as expected). Panel (b): pa-rameter va;ie is a = a tricrit − .
1, and the trimodal state appears to bestable (as expected). Both panels use the example interaction func-tion from Eq. (11), and both use equal fractionation ( x = y = z = / ψ = π /
3) in initial conditions. (a)(c) (d)(b)
FIG. 8.
Numerical experiments testing bistability.
Panel (a) and(b): we set a = . π ∈ ( √ π , a tricrit ) and both the bimodal stateand the trimodal state are stable (as predicted). Panel (c): we set a = . π − . < √ π < a crit and we see that the bimodal state isunstable (we have added black dashed lines so that one can see thatthe clusters away from the origin are not at ± π ). Panel (d): we set a = . π + . > a tricrit > √ π and the trimodal state is unstable(as predicted). In all panels N =
300 and the initial conditions areequally spaced and have equal fractionation with a random perturba-tion to all the phases of the oscillators. ± π . In panel (d), we set a = . π + . > a tricrit > √ π . Wealso observe an expected result, as trimodality appears to beunstable and the system converges to a bimodal equilibrium,which is stable given that a > √ π .In summary, we have a necessary condition for the stabilityof the trimodal equilibrium. Although, this condition is onlynecessary for stability, not sufficient, numerical experimentsseems to point to it being an accurate threshold in exampleswe have considered. Also, theory and numerical experimentsdemonstrate that multistability of different multimodal equi-libria is possible over parameter space. The theory for thestability of higher modes we leave for future work. REFERENCES T. Saigusa, A. Tero, T. Nakagaki, and Y. Kuramoto, Physical Review Let-ters , 018101 (2008). J. Myung, S. Hong, D. DeWoskin, E. De Schutter, D. B. Forger, andT. Takumi, Proceedings of the National Academy of Sciences , E3920(2015). A. F. Taylor, M. R. Tinsley, F. Wang, Z. Huang, and K. Showalter, Science , 614 (2009). M. Toiya, H. O. González-Ochoa, V. K. Vanag, S. Fraden, and I. R. Epstein,The Journal of Physical Chemistry Letters , 1241 (2010). S. Vaidyanathan, International Journal of ChemTech Research , 759(2015). Y.-N. Li, L. Chen, Z.-S. Cai, and X.-Z. Zhao, Chaos, Solitons & Fractals , 699 (2003). J. Pantaleone, American Journal of Physics , 992 (2002). R. Yoshida, M. Tanaka, S. Onodera, T. Yamaguchi, and E. Kokufuta, TheJournal of Physical Chemistry A , 7549 (2000). N. Henk and R.-a. Alejandro,
Synchronization of mechanical systems ,Vol. 46 (World Scientific, 2003). H. Ulrichs, A. Mann, and U. Parlitz, Chaos: An Interdisciplinary Journalof Nonlinear Science , 043120 (2009). B. H. Repp, Psychonomic bulletin & review , 969 (2005). S. Kirschner and M. Tomasello, Journal of Experimental Child Psychology , 299 (2009). A. De Paula, Journal of econometrics , 56 (2009). A. Pluchino, V. Latora, and A. Rapisarda, The European Physical JournalB-Condensed Matter and Complex Systems , 169 (2006). Y. Kuramoto,
Chemical oscillations, waves, and turbulence (Courier Cor-poration, 2003). Y. Kuramoto, in
International symposium on mathematical problems in the-oretical physics (Springer, 1975) pp. 420–422. S. H. Strogatz, Physica D: Nonlinear Phenomena , 1 (2000). J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler,Reviews of modern physics , 137 (2005). Y. Moreno and A. F. Pacheco, EPL (Europhysics Letters) , 603 (2004). Q. Wang, G. Chen, and M. Perc, PLoS one , e15851 (2011). I. R. Epstein, Physica D: Nonlinear Phenomena , 152 (1991). W. Hohmann, M. Kraus, and F. Schneider, The Journal of Physical Chem-istry A , 3103 (1998). A. Koseska, E. Volkov, A. Zaikin, and J. Kurths, Physical Review E ,031916 (2007). E. Ullner, A. Zaikin, E. I. Volkov, and J. García-Ojalvo, Physical reviewletters , 148103 (2007). E. Ullner, A. Koseska, J. Kurths, E. Volkov, H. Kantz, and J. García-Ojalvo,Physical Review E , 031904 (2008). S. A. Marvel, S. H. Strogatz, and J. M. Kleinberg, Physical review letters , 198701 (2009). R. Laje and G. B. Mindlin, Physical review letters , 288102 (2002). G. Balázsi, A. Cornell-Bell, A. B. Neiman, and F. Moss, Physical ReviewE , 041912 (2001). M. Giver, Z. Jabeen, and B. Chakraborty, Physical Review E , 046206(2011). H. Hong and S. H. Strogatz, Physical Review Letters , 054102 (2011). H. Hong and S. H. Strogatz, Physical Review E , 056210 (2012). D. J. Emlen, Evolution , 1219 (1996). S. M. Clifton, R. I. Braun, and D. M. Abrams, Proceedings of the RoyalSociety B: Biological Sciences , 20161970 (2016). S. Clifton, R. Braun, and D. Abrams, Dryad Digital Repository.(doi:10.5061/dryad. vb1pp) (2016). B. Damsgård, T. H. Evensen, Ø. Øverli, M. Gorissen, L. O. Ebbesson,S. Rey, and E. Höglund, Royal Society Open Science , 181859 (2019). B. Damsgård, T. Evensen, Ø. Øverli, M. Gorissen, L. Ebbesson, S. Ray,and E. H´’aglund, “Data from: Proactive avoidance behaviour and pace-of-life syndrome in atlantic salmon,” (2019). M. R. Blanton, D. W. Hogg, N. A. Bahcall, I. K. Baldry, J. Brinkmann,I. Csabai, D. Eisenstein, M. Fukugita, J. E. Gunn, Ž. Ivezi´c, et al. , TheAstrophysical Journal , 186 (2003). K. Abazajian, J. K. Adelman-McCarthy, M. A. Agüeros, S. S. Allam, S. F.Anderson, J. Annis, N. A. Bahcall, I. K. Baldry, S. Bastian, A. Berlind, et al. , The Astronomical Journal , 2081 (2003). I. K. Baldry, K. Glazebrook, J. Brinkmann, Ž. Ivezi´c, R. H. Lupton, R. C.Nichol, and A. S. Szalay, The Astrophysical Journal , 681 (2004). I. Ronin, N. Katsowich, I. Rosenshine, and N. Q. Balaban, Elife , e19599(2017). A. Goldberg, O. Fridman, I. Ronin, and N. Q. Balaban, Genome medicine , 112 (2014). Z. Levnaji´c, Physical Review E , 016231 (2011). Z. Levnaji´c, Scientific Reports , 967 (2012). A. V. Pimenova, D. S. Goldobin, M. Rosenblum, and A. Pikovsky, Scien-tific Reports , 38518 (2016). K. H. Nagai and H. Kori, Physical Review E , 065202 (2010). S. Gil, Y. Kuramoto, and A. S. Mikhailov, EPL (Europhysics Letters) ,60005 (2010). C. C. Gong, C. Zheng, R. Toenjes, and A. Pikovsky, Chaos: An Interdisci-plinary Journal of Nonlinear Science , 033127 (2019). Y. Nakamura, F. Tominaga, and T. Munakata, Physical Review E , 4849(1994). A. Mishra, C. Hens, M. Bose, P. K. Roy, and S. K. Dana, Physical ReviewE , 062920 (2015). A. Yeldesbay, A. Pikovsky, and M. Rosenblum, Physical review letters ,144103 (2014). D. Golomb, D. Hansel, B. Shraiman, and H. Sompolinsky, Physical ReviewA , 3516 (1992). L. Tsimring, N. Rulkov, M. Larsen, and M. Gabbay, Physical Review Let-ters , 014101 (2005). K. Okuda, Physica D: Nonlinear Phenomena , 424 (1993). P. Billingsley,
Probability and Measure. 1995 (John Wiley & Sons, NewYork, 1995) p. 357. C. Chaudhary, H. Saeedi, and M. J. Costello, Trends in Ecology & Evolu-tion , 670 (2016). P. Seeman, C. Ulpian, C. Bergeron, P. Riederer, K. Jellinger, E. Gabriel,G. Reynolds, and W. Tourtellotte, Science , 728 (1984). T. J. Lewis and J. Rinzel, Journal of Computational Neuroscience , 283(2003). Y. Wu, C. Zhou, J. Xiao, J. Kurths, and H. J. Schellnhuber, Proceedings ofthe National Academy of Sciences , 18803 (2010). P. DiMaggio, J. Evans, and B. Bryson, American Journal of Sociology ,690 (1996). B. R. Flay, Behavioral Science , 335 (1978). P. Bonifacio, E. Caffau, M. Spite, M. Limongi, A. Chieffi, R. Klessen,P. François, P. Molaro, H.-G. Ludwig, S. Zaggia, et al. , Astronomy & As-trophysics , A28 (2015). J. Reiter and I. R. Epstein, Journal of Physical Chemistry , 4813 (1987). Y. Ren-Chang, B. Jian-Chun, and F. Qiu-Jun, Atmospheric and OceanicScience Letters , 229 (2011). S.-S. Sun and R. W. Nesbitt, Earth and Planetary Science Letters , 429(1977). A. Dekel and Y. Birnboim, Monthly Notices of the Royal AstronomicalSociety , 2 (2006). C. Zhang, B. E. Mapes, and B. J. Soden, Quarterly Journal of the RoyalMeteorological Society: A Journal of the Atmospheric Sciences, AppliedMeteorology and Physical Oceanography , 2847 (2003). P. K. Kuhl and A. N. Meltzoff, Science , 1138 (1982). D. Sun, J. Bloemendal, D. K. Rea, Z. An, J. Vandenberghe, H. Lu, R. Su,and T. Liu, Catena , 325 (2004). J. Thorpe, Journal of Fish Biology , 175 (1977). Y. Wu, C. Zhou, J. Xiao, J. Kurths, and H. J. Schellnhuber, Proceedings ofthe National Academy of Sciences , 18803 (2010). P. Seeman, C. Ulpian, C. Bergeron, P. Riederer, K. Jellinger, E. Gabriel,G. Reynolds, and W. Tourtellotte, Science , 728 (1984). J. Schwab, P. Podsiadlowski, and S. Rappaport, The Astrophysical Journal , 722 (2010). S. Mao, R. Narayan, and T. Piran, The Astrophysical Journal , 171(1994). M. Eichhorn, Journal of forest research , 391 (2010). − π − π π π -2-1012 f ( ψ ) ψ FIG. S1.
Additional interaction functions.
Solid blue curve: tri-angle wave from Eq. (S1); solid red curve: antisymmetrized variantof the von Mises distribution from Eq. (S2) with κ <
0; dashed redcurve: antisymmetrized variant of the von Mises distribution fromEq. (S2) with κ >
0. Panels (a) and (b) of Fig. S2 use the trianglewave. Panels (c) and (d) use the antisymmetrized von Mises func-tion, with positive κ (dashed red) in panel (c) and negative κ (solidred) in panel (d). We note that for κ > ± π is neversteeper when compared to the origin and for κ < ± π . Supplementary Material: A cou-pled oscillator model for the originof bimodality and multimodality
S1. ADDITIONAL COUPLING FUNCTIONS
Figure S1 illustrates two additional coupling functions thatwe examined. We used a variant of the triangle wave (blue,solid) given by the equation f tri ( u ; c ) = (cid:40) uc | u | < c uc − π − sign ( u )( π c − π ) c ≤ | u | ≤ π , (S1)assuming that 0 < c < π , and an antisymmetrized variant ofthe von Mises distribution (red curves) given by f vM ( u ; µ , κ ) = sin ( u − µ ) e κ cos ( u − µ ) π I ( κ ) . (S2)We numerically probe the stability of the bimodal equi-librium using these interaction functions in Fig. S2. Here N = N ( , ) , the phase perturbation, ξ i , is drawn fromthe distribution N ( , . ) and we set K = − c = π /
4; this gives a stable fractionation threshold1 / < x < /
4. We test that threshold numerically by setting x initial = / < / x initial = / > / -0 P ha s e Time -0 P ha s e Time(a)(c) (d)(b)
FIG. S2.
Numerical experiments using additional interactionfunctions.
We test the stability of the bimodal equilibria for alterna-tive coupling functions shown in Fig. S1. (a) Triangle wave couplingwith initial fractionation in predicted stable range. (b) Triangle wavecoupling with initial fractionation outside predicted stable range. (c)Von Mises coupling with κ > κ > N =
100 and oscillators’ natural frequencies are drawn from the dis-tribution N ( , ) . Initial phases are bimodally distributed withmodes at 0 and π , with perturbations ξ i , i = ,..., N , are drawn from N ( , . ) . In panels (c) and (d) we use the antisymmetrized von-Misesfunction from Eq. (S2) with µ = x initial = /
2. In panel(c) we set κ =
10, and, as expected, we see that the bimodalequilibrium appears unstable; this is because there does notexist a range of x such that Eq. (10) can be satisfied giventhat the slope at the origin is far steeper than the slope at the ± π . We note that in (c) the system appears to tend to theincoherent state. In panel (d) we set κ = −
10 and observethat the bimodal state appears to be stable under perturbation,which is expected given that the slope at the ± π is steeperwhen compared to the origin. S2. BASINS OF ATTRACTION FOR MULTIMODALSTATES
We have conducted some preliminary numerical explo-ration of the sizes of basins of attraction for various equilib-ria for the example interaction function given in Eq. (11) ofthe main text. We simulated the system one hundred timeswith initial phases chosen independently at random from theuniform distribution over the circle, i.e. U ( − π , π ] , and eval-uated the fraction of the time that the system converged toeach distinct equilibrium state. Results are shown in Fig. S3,with N = K = − N ( , ) .Fig. S3 also shows the stability thresholds described inEqns. (12) (bimodal state), (A.13) (antiphase state), and tri-modal state (A.12) of the main text, visualized by the solidblack, and dot-dashed green, and magenta vertical lines re-spectively. In order to classify the observed equilibria, we usea k -means algorithm on the unit circle, with the number ofclusters, k , being decided by the gap statistic. We say that aequilibrium state is bimodal if k =
2, trimodal if k =
3, and so0 P r obab ili t y o f C on f i gu r a t i on FIG. S3.
Basins of attraction . We plot the fraction of uniform ran-dom initial conditions that end up in bimodal (blue circles), trimodal(orange asterisks), or higher order multimodal (purple xs) states forthe concrete system examined in the main text. Here N = K = − N ( , ) . We performed 100 unique simulationsfor each value of a . Final states (presumed equilibria) were identi-fied automatically via k -means clustering. Thresholds given in themain text for stability of bimodality and the antiphase state are givenby the solid black line and the dot-dashed green line, respectively.The threshold for the necessary condition for stability of the trimodalstate is given by the vertical dashed magenta line. on.We note that the results are consistent with our analysis inthat the probability of a configuration is always zero in rangesof a where it is excluded. Although, we have not analyzedequilibria with more than three modes, we observe that suchmodes are unlikely to be observed for most values of a , andthus have apparently small basins of attraction.Given that this experiment was conducted with heteroge-neous oscillators, this lends plausibility to the idea that thesystem will end up in a multimodal state for sufficiently largecoupling. More formal analysis of the basin size of the bi-modal and trimodal state will be left for future work. S3. CRITICAL COUPLING STRENGTH
In the standard Kuramoto model with attractive coupling,there exists a critical coupling strength K c at which the systembifurcates from an incoherent state to the ordered state. Tolook for K dependence in the system detailed in main text,we examine the simplest cases of N = N =
3, and alsoconduct several numerical experiments with results shown inFig. S4, though we leave more thorough exploration for futurework.Figure S4 shows how order varies as we increase couplingstrength among nonidentical oscillators with the concrete in-teraction function used in the main text. Here, we set N = N ( , σ ) .From here, we vary the quantity K / σ so that log ( K / σ ) runsfrom -2 to 4. Each curves shown above represents the result -3 -2 -1 0 1 2 3 400.51 O r de r P a r a m e t e r (a)(b) FIG. S4.
Critical coupling strength . We perform numerical experi-ments to demonstrate the existence of a critical coupling strength forour system and evaluate its dependence on parameter a using the in-teraction function defined in Eq. (11) of the main text. Here N = N ( , σ ) , and the ini-tial phase distribution is ρ ( θ ) = . δ ( θ ) + . δ ( θ − ψ ) , where ψ is the predicted phase separation given by the stable fixed points ofEq. (11). Here, each curve represents a different value of a (valuesindicated in legend). As in the standard Kuramoto model, the criticalcoupling strength is dependant on the size of the standard deviationof the distribution, but unlike the standard Kuramoto model, it ap-pears to also depend on a , which sets the shape of the interactionfunction. of an experiment for a given value of a . Here, the order pa-rameter is defined as follows: R = max (cid:40)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∑ j e i θ j N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∑ j e i θ j N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∑ j e π a i θ j N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:41) . (S1)Defining the order parameter in this fashion sets the value ofthe order parameter to be 1 whenever the final configurationis bimodal or an equally spaced trimodal solution. Just as inthe standard Kuramoto model, if the coupling strength K is notsufficiently large in magnitude, the system goes to the incoher-ent state due to intrinsic oscillator heterogeneity. We observethat the critical coupling strength appears to be proportionalto the standard deviation of the frequency distribution, similarto the result in the standard Kuramoto analysis, but we pointout that the critical coupling strength K c also appears to havedependence on the value of a . We believe that some insightinto this dependence can be gained from examining the simple N = N = N =
2, the system reduces to˙ ψ = ∆ ω − K f ( ψ ) (S2)where ∆ ω = ω − ω . Setting ˙ ψ =
0, we find that a fixed point ψ must satisfy the equation: ∆ ω K = f ( ψ ) . (S3)Note, this fixed point does not always exist, but if the couplingfunction f has zeros, a fixed point must arise as | K | → ∞ .1Even without explicitly defining ψ , we can observe scalingdependencies for the critical coupling strength K c , which isdefined such that f ( ψ max ) = ∆ ω K c (S4)where ψ max ∈ ( − π , π ] is the value such that f ( ψ max ) = max f ( ψ ) (the arg max). We observe that K c ∝ ∆ ω , which isexpected if K c ∝ σ as in the standard Kuramoto model (sincefor two oscillators σ ∝ ∆ ω ) and is observed in our numericalexperiments even for N (cid:29) K c scales with the maximum valueof the interaction function f , which in our numerical experi-ments depends on the parameter a . Similar dependence is alsoevident if we consider the N = N =
3, we take the natural frequencies (without lossof generality) to be 0 , − σ / , σ / ψ = θ − θ and ψ = θ − θ , and arrive at two conditions for existence of equilibria: σ K = f ( ψ − ψ ) − f ( ψ ) − f ( ψ ) (S5) σ K = f ( ψ − ψ ) + f ( ψ ) + f ( ψ ) , (S6) which simplify to σ K = f ( ψ − ψ ) + f ( ψ ) (S7) f ( ψ ) = − f ( ψ ) . (S8)Hence, a necessary condition K must satisfy for the existenceof equilibria is σ K ≤ f ( ψ max ) . (S9)So, just as in the N = K c is proportional to the oscillator heterogeneity σ and inversely proportional to the maximum of the interactionfunction f .We hypothesize that similar scaling laws hold for N (cid:29)(cid:29)