Asymmetric dynamic interaction shifts synchronized frequency of coupled oscillators

Interacting dynamic agents can often exhibit synchronization. It has been reported that the rhythm all agents agree on in the synchronized state could be different from the average of intrinsic rhythms of individual agents. Hinted by such a real-world behavior of the interaction-driven change of the average phase velocity, we propose a modified version of the Kuramoto model, in which the ith oscillator of the phase ϕi interacts with other oscillator j only when the phase difference \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\phi }}_{{j}}$$\end{document}ϕj − \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\phi }}_{i}$$\end{document}ϕi is in a limited range [−βπ, απ]. From extensive numerical investigations, we conclude that the asymmetric dynamic interaction characterized by β ≠ α leads to the shift of the synchronized frequency with respect to the original distribution of the intrinsic frequency. We also perform and report our computer-based synchronization experiment, which exhibits the expected shift of the synchronized frequency of human participants. In analogy to interacting runners, our result roughly suggests that agents tend to run faster if they are more influenced by runners ahead than behind. We verify the observation by using a simple model of interacting runners.

Synchronization is a widespread phenomenon and has been observed in various physical and biological systems like lasers, sperms, fireflies, and so on [1][2][3][4][5][6] . Winfree proposed a model of coupled oscillator system 7,8 and Kuramoto simplified the model later to make it analytically tractable 9,10 . Kuramoto model has also been studied in a variety of different interaction structures, including d-dimensional regular lattices and complex network structures 11-13 . There have been various extensions since the original Kuramoto model has been proposed. For example, the effect of the time delay in interaction between oscillators has been studied 14,15 , and the correlation between the characteristic frequency and the coupling strength to describe the neural activity in brain has been investigated [16][17][18] . However, most existing studies have assumed that the interaction structure is static and thus does not change in time. Furthermore, extension of the Kuramoto model in such a way that the dynamic states of the oscillators are closely coupled with interaction structure has not been tried. In reality, it is often observed that the interaction structure can be dynamically coupled with the internal state of individual agent, and thus the oscillator system can hardly be an exception. In our approach in this paper, the asymmetric dynamic interaction due to the difference of phases of oscillators is proposed, which has not been considered in existing studies. We believe that our proposed model is also realistic and perform computer-based experiments of human participants.
The synchronized rhythm of oscillators can be faster or slower than average frequency of isolated individual oscillators. For example, music play of finger tapping has been reported to show a faster paired synchronized rhythm when there is no timing cue by the director 19 . The hand clapping of students has also been studied and faster synchronized rhythm like for the finger tapping 19 has been reported 20 . For applause of audiences, however, slower synchronized rhythm has been observed 21 . Although a broad range of phenomena can be described through the use of the Kuramoto model, above mentioned behavior of the off-the-mean synchronized rhythm cannot be explained by the conventional model. One can easily understand this by a simple symmetry argument: The equation of motion for the conventional Kuramoto model is invariant under the phase reversal transformation φ φ → − i i for ∀i, and thus the shift of the average phase velocity φ 〈 〉 i  contradicts the symmetry. From this reason, we suggest that we need to break the phase reversal symmetry in the conventional Kuramoto model in order to explain the observed shift of the synchronized rhythm. In our model, the phase reversal symmetry is broken not explicitly, but dynamically, as will be explained below. Furthermore, we perform smartphone and

Model
We suggest a simple modified Kuramoto model of N coupled oscillators with dynamic interaction to mimic the shift of the average frequency: where > K( 0) is the coupling strength, and φ π π ∈ − [ , ) i and ω i are the phase variable and the time-fixed (quenched) intrinsic frequency of the ith oscillator, respectively. The intrinsic frequency ω i is chosen randomly from the normal distribution g( ) (1/ 2 )exp( /2) 2 ω π ω = − with zero mean and unit variance. In detail, we first generate N random numbers ω′ i from the normal distribution and use the shifted frequency i i N i 1 ω ω ω = ′ − ∑ ′ as the intrinsic frequency in Eq. (1) to make it sure that the average intrinsic frequency ω ∑ N i 1 is null. The null value of the mean frequency suggests that all phase variables and phase velocities are to be measured in the center of mass frame in which 0 does not mean that the oscillator is not oscillating but it only means that it is oscillating with the average intrinsic frequency. Our model differs from the conventional Kuramoto model in that the interaction structure dynamically depends on phase variables: t ( ) i Λ in Eq. (1) is the set of oscillators which the ith oscillator interact with and is written as Our key parameters are α and β in t ( ) i Λ , which determine the condition for interaction. It is to be noted that if α = β = 1, all oscillators interact with all other oscillators and thus our model becomes identical to the globally-coupled Kuramoto model. If α = β < 1, the interaction becomes limited but is still symmetric. If α > β (α < β), on the other hand, the interaction becomes asymmetric toward oscillators with advanced (lagged) phases. The two asymmetric cases are equivalent to each other under the phase reversal transformation, i i φ φ → − for ∀i and α β ↔ . Accordingly, we only focus on two cases: The symmetric interaction with α = β and forward-biased asymmetric interaction with α > β. In Fig. 1, we display a schematic diagram for the meaning of α and β in Eq. (2). numerical results. We numerically integrate Eq. (1) using the Heun's method, which is the second-order algorithm, with the discrete time step t 0 01 ∆ = . for 2 10 5 × time steps, corresponding to ∈ t [0, 2000]. After a sufficiently long time, the system arrives at the steady state, and the first half 10 5 time steps are discarded before we start to measure observables for the later half time steps. All the results in the present paper are obtained from the averages over 100 (200 for = N 25) independent and different random initial phase variables and intrinsic frequency realizations.
We measure the standard order parameter R for the synchronization transition defined as www.nature.com/scientificreports www.nature.com/scientificreports/ where  〈 〉 and  denote the ensemble and the temporal averages, respectively. We also measure the average phase velocity Ω in the steady state defined as We first present our result for the symmetric case that 1/2 α β = = in Fig In contrast, at any value of K, the average phase velocity Ω in Eq. (4) remains close to 0. A small bump in Fig. 2b around ≈ K K c does not depend much on the system size, but it becomes shallower as the number of samples used in the ensemble average is increased. We believe that the bump structure is a reminiscent of the large fluctuation of the phase velocity around the critically. We thus conclude that the symmetric interaction characterized by 0 α β = > does not change the average phase velocity from the average intrinsic frequency, which is set to null, irrespective of whether the system is in the asynchronous or in the synchronous phase.
Our observation of null average phase velocity for α β = can be understood through a simple symmetry argument: Due to the phase reversal symmetry under the transformation i i φ φ → − for ∀ i for the symmetric case of α β = , there is no reason that the phase velocity is biased toward a positive or a negative value. Consequently, the average phase velocity Ω must be identical to the average intrinsic frequency ω ∑ N i i 1 , which is set to null as described above. In contrast, the off-the-mean synchronized rhythm has been often reported in music play and hand clapping [19][20][21] . Accordingly, we suggest that such observation of the shift can be explained by the biased (or asymmetric) interaction.
For the asymmetric cases α β > ( ) , we first display our results for 1 α = and 0 1 β < < in Fig. 3 for = N 200. As shown in Fig. 3(a), the system again exhibits a well-defined synchronization transition and the order parameter R changes from as the coupling strength K is increased. The average phase velocity Ω in Fig. 3b, in contrast, shows a very interesting behavior, which sharply differs from what has been observed for the symmetric interaction in Fig. 2. As K is increased from null value, Ω first increases with K, and displays a peak, beyond which Ω decreases with K, as shown in Fig. 3b. The observed nonmonotonous behavior of Ω in Fig. 3b can be understood as follows: In the limit of the vanishing coupling strength → K 0, all oscillators are completely decoupled and the phase reversal symmetry is preserved to yield Ω = 0. As K is increased, the asymmetric interaction breaks the symmetry and Ω increases. In the opposite limit of the strong coupling → ∞ K ( ) , all oscillators are fully synchronized and interact with all other oscillators as in the globally-coupled Kuramoto model. Consequently, we expect 0 Ω → in both limits of K 0 → and K → ∞, which suggests that the existence of a peak of Ω in the middle is inevitable. It appears that the peak position of Ω is close to the synchronization transition point in Fig. 3a. Although not shown here, we also check the finite-size effect on the shift of Ω, only to find that the peak structure becomes more enhanced in a larger system. We thus conclude that our modified Kuramoto model with asymmetric interaction reveals the observed off-the-mean synchronized rhythm in previous studies 19,20 .
We next investigate the limiting case of the asymmetric interaction characterized by α > 0 and 0 β = . In words, this type of interaction means that each oscillator interacts with other oscillators only when their phases are ahead. Figure 4a displays the average phase velocity Ω for β = 0 and 1, 3/4, 1/2 α = , and 1/4 for = N 200. We www.nature.com/scientificreports www.nature.com/scientificreports/ find that the decrease of Ω in the strong-coupling limit observed in Fig. 3b disappears for β = 0, and Ω stays at a positive value even when K is increased to a very large value. We emphasize that the strong-coupling limit for β = 0 is different from the same limit for 0 β ≠ : The interaction cannot be of all-to-all type for the former case due to the strict bound 0 β = . We thus expect that the finite positive value of Ω should persist even when the coupling strength → ∞ K . Similarly to the peak structure in Fig. 3b for 0 β > , the peak of Ω is clearly seen in Fig. 4a for β = 0, near the synchronization transition point. We believe that the height of the peak of Ω can be related with maximum intrinsic frequency. The maximum value M N of N random numbers generated from the normal distribution is given by − and the Euler's constant γ 22,23 . In Fig. 4b, we display the peak height Ω max [see Fig. 4a], the measured value of max i i max ω ω = , and M N , versus N for ( , ) (1, 0) α β = . The latter two values, max ω and M N , almost coincide as expected, and max Ω follows the similar form of the logarithmic increase. However, the peak height Ω max is significantly larger than M N (e.g, Ω ≈ . M / 136 N max for = N 200), which means that the average phase velocity of oscillators is larger than the maximum intrinsic frequency. We believe that this result is particularly interesting, since the forward-only asymmetric interaction drives all oscillators to have angular velocity beyond the maximum intrinsic frequency. In contrast, the peak height in Fig. 3b for 1/10 β = is found to be less than the maximum intrinsic frequency (we find M / 068) N max Ω ≈ . . In order to examine details of the dynamics of the system, we also integrate Eq. (1) starting from a single random initial condition and a single configuration of ω { } i for α β = ( , ) (1,0) at coupling strengths = K 1 and 5. In Fig. 5a,b we display individual phase trajectories for N 25 = oscillators in the moving reference frame of the   Fig. 5b. We also display the individual phase velocity t ( ) Fig. 5c,d for K 1 = and 5, respectively. The behavior reported in Fig. 5d  experiments. We next perform smartphone and computer-based experiments with 31 individual human participants. The server computer plays the periodic sound of a metronome for 15 seconds and each participant is asked to press the touchpad(display) button of her/his notebook computer(android phone) in accord with the metronome sound. After 15 seconds, the metronome sound of the server is stopped, but each participant is asked to keep pressing the button for about 1 minute. We use two different programs in the client side to examine the effect of interaction among participants: (i) Without Interaction: The beeping sound of the client device is turned off. Accordingly, each player tries only to continue the rhythm of metronome (s)he listened for the first 15 seconds, without any interaction or interference with other participants. In this case, the interaction between players is unlikely since there is no sound cue in the later stage after 15 seconds and we ask players only to look at their computer or smartphone screens. (ii) With Interaction: The beeping sound is turned on, and all participants are asked to synchronize the beeping sound (s)he generates to the sound generated by other players. The advantage of using our platform is that the time instants of button pressing are automatically recorded in the server, which makes further analysis straightforward.
We use 64 and 120 BPM (beats per minute) for the metronome rhythm, which correspond to the metronome frequency 1.07 and 2.0 Hz, respectively. In Fig. 6a, we display partial results of our experiment for the 64 BPM metronome for the case of (ii) With Interaction, as an example. The time instants of button pressing are displayed in the form of the short vertical bars, as often used for the visualization of the neuron firing pattern in neuroscience 24 . In Fig. 6b, we show how the average frequency changes in time for the case of (i) Without and (ii) With Interaction. We compute the frequency of each individual by using the time window of 5 seconds, and results from four independent experiments are averaged. Note that the two curves almost coincide to the metronome www.nature.com/scientificreports www.nature.com/scientificreports/ frequency 1.07 Hz at around 10 second, before the metronome sound is stopped. As time goes on, it should be recognized that the average frequency for (ii) With Interaction case increases. Differently from previous study 19 , we have not observed downward shift of the synchronized frequency for (i) Without Interaction case. We believe that our experimental result in Fig. 6b supports our previous finding in Figs. 3 and 4, since the average frequency for (ii) With Interaction lies significantly higher than for (i) Without Interaction. However, for 120 BPM the difference between (i) Without and (ii) With Interaction case is found insignificant. Probably, 120 BPM could be too fast for some participants to comfortably follow.
interacting runners. We suggest that our main finding of the shift of the average phase velocity can also be tested by other more realistic phenomenon, like interacting runners. In this regard, it is interesting to note that researchers recorded two international marathon races and found the grouping behavior of the runners 25 . Furthermore, it has been found that the denser the group is, the faster the runners in the group are 25 , which indicates that the speed of a runner is not only an intrinsic property of the runner, but also can be an outcome from interaction with other runners.
Hinted by both the previous observation 25 and our results for the shift of the average phase velocity of oscillators, we propose a simple model to mimic interacting runners in a one-dimensional linear track and write the equation of motion for the position x i of the ith runner as where the characteristic speed v i is a Gaussian quenched random variable with the mean v 5 m/s 〈 〉 = and the standard deviation 1 m/s v σ = , = .
C( 0 1/s) is the competitive coupling strength, and Λ t ( ) i is the set of runners within the ith runner's sight limit ( 20 m) = : In Eq. (5), we try to mimic the situation that if one sees other runners in front in close distance, the runner tries to speed up within his limitation. However, if the distance to the front runner is too big, the runner may not try hard to catch up. Only for simplicity, we assume the runner's speed is limited by v i v σ + from that every human individual runner must have some physical and biological limitation. Note that if there is no front runner in the distance , the runner runs at his intrinsic speed of v i . Figure 7 shows our result for the forward-interacting runners for = N 50. We first observe that runners form groups of different speeds, and the runner's speed within each group is identical. It is to be noted that although our model of runners with forward-biased interaction is very simple, it produces the grouping behavior reported in previous study 25 . The average speed of all runners is measured 5.5 m/s, which is larger than the average characteristic speed v 5 m/s i 〈 〉= . The qualitative feature of the grouping behavior is observed robust if the value of C is not too small, nor too large. Although our interacting runner model does not include physiological factors which affect the performance of runners 26,27 , the shift of the average speed is observed as in the marathon race 28 . We emphasize that this finding of the shift of the average speed is in parallel to our observation for oscillators with asymmetric interaction in Figs. 3 and 4 as well as our experimental observation in Fig. 6. . When the players are interacting with each other, the synchronized frequency in later time is found to be significantly larger than the metronome frequency 64 BPM(=1.07 Hz).

Discussion
In this paper, we have suggested models of interacting oscillators and runners, in which the key ingredient is the forward-biased interaction. Through numerical investigations, we have observed that such asymmetric interaction shifts the average velocity in the two model systems, in accord with the findings in previous studies 19,20 . Our simple model of forward-interacting runners has also been shown to produce the grouping behavior reported in previous study 25 . We have also performed a computer-based synchronization experiment of interacting human individuals and observed again the shift of the average frequency. We believe that when we synchronize sound rhythm with others we tend to catch up the sound of others slightly ahead in time. Accordingly, we suggest that our modified Kuramoto model can have some relevance in explaining our experimental result of sound synchronization of human individuals. The asymmetric dynamic interaction has not been studied before in the research field of synchronization and our new computer-based experiments also show that our model study is realistic.  Runners start running at their own characteristic speed {v i } but as time goes on they form about 10 groups, in accord with the previous finding 25 . The size of each group is 4,8,14,8,9,5,1 and 1 in the order of group speed for t ≳ 2000. Runners within each group run at the same speed. (b) The average speed of all runners at later times is 5.5 m/s, which is larger than the average characteristic speed 5 m/s.