Connecting empirical phenomena and theoretical models of biological coordination across scales

Coordination in living systems—from cells to people—must be understood at multiple levels of description. Analyses and modelling of empirically observed patterns of biological coordination often focus either on ensemble-level statistics in large-scale systems with many components, or on detailed dynamics in small-scale systems with few components. The two approaches have proceeded largely independent of each other. To bridge this gap between levels and scales, we have recently conducted a human experiment of mid-scale social coordination specifically designed to reveal coordination at multiple levels (ensemble, subgroups and dyads) simultaneously. Based on this experiment, the present work shows that, surprisingly, a single system of equations captures key observations at all relevant levels. It also connects empirically validated models of large- and small-scale biological coordination—the Kuramoto and extended Haken–Kelso–Bunz (HKB) models—and the hallmark phenomena that each is known to capture. For example, it exhibits both multistability and metastability observed in small-scale empirical research (via the second-order coupling and symmetry breaking in extended HKB) and the growth of biological complexity as a function of scale (via the scalability of the Kuramoto model). Only by incorporating both of these features simultaneously can we reproduce the essential coordination behaviour observed in our experiment.


Introduction
Coordination is central to living systems and biological complexity at large, where the whole can be more than and different from the sum of its parts. Rhythmic coordination [1][2][3] is of particular interest for understanding the formation and change of spatio-temporal patterns in living systems, including, e.g. slime mould [4,5], fireflies [6,7], social groups [8,9] and the brain [10][11][12][13][14]. Theoretical descriptions of biological coordination are often in terms of coupled oscillators, whose behaviour is constrained by their phase relations with each other [2,[15][16][17][18]. Existing studies of phase coordination often focus on systems of either very few (small-scale, mostly N = 2) [13,19,20], or very many oscillators (large-scale, N → ∞) [21][22][23]. Here we inquire how the two might be connected and applied to midscale systems with neither too many nor too few components. The present work answers this question by modelling empirically observed coordinative behaviour in midscale systems (N = 8), based on data collected in a specially designed human experiment [24]. The resultant model that captures all key experimental observations happens to also connect previous theories of small-and large-scale biological coordination in a single mathematical formulation.
But first, how are small-and large-scale models different? Small-scale models were usually developed to capture empirically observed coordination patterns, as in animal gaits [25][26][27], bimanual movement coordination [28,29], neuronal coordination [30], interpersonal coordination [31,32], human-animal coordination [33] and human-machine coordination [34,35]. They describe multiple stable coordination patterns (multistability) and the transitions between them (order-to-order transitions), e.g. from a trot to a gallop for a horse [36]. In humans, dyadic coordination patterns like inphase and antiphase (synchronization, syncopation) were found across neural, sensorimotor and social levels (see [13,14] for reviews), well captured by the extended Haken-Kelso-Bunz (HKB) model [29,37,38]. However, the extended HKB was restricted to describing coordination phenomena at N = 2 (i.e. not directly applicable to higher-dimensional coordination phenomena). By contrast, large-scale models are concerned more about statistical features like the overall level of synchrony and disorder-to-order transitions, but not so much about patterns at finer levels. As a representative, the classical Kuramoto model [2] is applicable to describing a wide range of large-scale coordination between, e.g. people [23,39], fish [40] and neural processes [22], and is often studied analytically for its incoherence-to-coherence transition (at the statistical level, for N → ∞; see [41,42] for reviews).
Although the extended HKB and the classical Kuramoto model emerged separately, they connect to each other by an interesting difference: the Kuramoto model with N = 2 is almost the extended HKB model except that the former lacks the term responsible for antiphase coordination in the latter (more accurately, the bistability of inphase and antiphase). Bistability of inphase and antiphase coordination, with associated order-to-order transitions and hysteresis, happens to be a key observation in small-scale human experiments [28,43]. This begs the question of whether there is a fundamental difference between large-and small-scale coordination phenomena. Does the existence of antiphase, multistability, and order-to-order transitions depend on scale N? With these questions in mind, we recently conducted a human experiment [24] at an intermediate scale (N = 8), such that the system is large enough for studying its macro-level properties, yet small enough for examining patterns at finer levels, ideal for theories and empirical data to meet at multiple levels of description. In the following sections, we demonstrate how the marriage between the two models (not either one alone) is sufficient for capturing empirical observations at multiple levels of description and we discuss its empirical and theoretical implications for biological coordination.

Human coordination at intermediate scales
Before getting into the model, we briefly review the mid-scale experiment and key results [24]. In the experiment (dubbed the 'Human Firefly' experiment), ensembles of eight people (N = 8, total 120 subjects) spontaneously coordinated rhythmic movements in an all-to-all network (via eight touchpads and eight ring-shaped arrays of eight LEDs as in figure 1; see Material and methods for details), even though they were not explicitly instructed to coordinate with each other. To induce different grouping behaviour, subjects were paced with different metronomes prior to interaction such that each ensemble was split into two frequency groups of equal size with intergroup difference δf = 0, 0.3 or 0.6 Hz (referred to as levels of 'diversity'), and were asked to maintain that frequency during interaction after the metronome was turned off. Subjects' instantaneous tapping frequencies from three example trials (figure 2a-c) show intuitively the consequences of frequency manipulations: from (a) to (c) a supergroup of eight gradually split into two frequency groups of four as diversity increased from δf = 0 to 0.6 Hz.
Key results involve multiple levels of description, in terms of intergroup, intragroup and interpersonal relations. The level of intergroup integration is defined as the relationship between intragroup and intergroup coordination (β 1 , slope of regression lines in figure 3a. Here intragroup coordination is measured by the average pair-wise phase-locking value over all intragroup dyads, and likewise, intergroup coordination over all intergroup dyads. Phase-locking value per se is a measure of stability of a relative phase pattern within a period of time, which equals to one minus the circular variance. See section 'Phase-locking value and level of integration' for technical details). Intuitively, we say that two groups are integrated if intragroup and intergroup coordination facilitate each other ( positive relation between respective phase-locking values, β 1 > 0), and segregated if intragroup and intergroup coordination undermine each other (negative relation between respective phase-locking values, β 1 < 0. We will see later in figure 4 how this measure meaningfully captures coordination dynamics). In the experimental result, two frequency groups were integrated when diversity is low or moderate (δf = 0, 0.3 Hz, blue and red lines, slope β 1 > 0) and segregated when diversity is high (δf = 0.6 Hz, yellow line, slope β 1 < 0). A critical level of diversity demarcating the regime of intergroup integration and segregation was estimated to be δf* = 0.5 Hz. Within the frequency groups, coordination was also reduced by the presence of intergroup difference (figure 3b, left, red and yellow bars shorter than blue bar). At the interpersonal level, inphase and antiphase were preferred phase relations (inphase much stronger than antiphase; distributions in figure 2d-f ), especially when the diversity was very low (figure 2d, peaks around ϕ = 0, π, in radians throughout this paper), but both were weakened by increasing diversity (figure 2e,f; in episodes of strong coordination, antiphase is greatly amplified and much more susceptible to diversity than inphase, see [24]). Notice that subjects did not remain locked into these phase relations but rather engaged and disengaged intermittently (two persons dwell near and escape from preferred phase relations recurrently, a sign of metastability [13]; see figure 6a red trajectory for example), reflected also as 'kissing' and 'splitting' of frequency trajectories (e.g. in figure 2b). In the following sections, we present a model that captures these key experimental observations at their respective levels of description.

A minimal experiment-based model of multiagent coordination
Our model of coordination is based on a family of N oscillators, each represented by a single phase angle w i . We will show that a pair-wise phase coupling [2,25,29] of the form The coefficients a ij > 0 and b ij > 0 are parameters that govern the coupling. Equations (2.1) include a number of well-studied models as special cases. For instance, setting ϕ := w 1 − w 2 , δω := ω 1 − ω 2 ,ã :¼ a 12 þ a 21 and 2b :¼ b 12 þ b 21 for N = 2, the difference of the two resulting equations (2.1) yields the relative phase equation of the extended HKB model [37]. The HKB model [29] was originally designed to describe the dynamics of human bimanual coordination, corresponding to equation (2.2) with δω = 0 (i.e. describing the coordination between two identical components). The extended HKB introduces the symmetry breaking term δω to capture empirically observed coordinative behaviour between asymmetric as well as symmetric components (i.e. the HKB model is included in the extended HKB model, which is further included in equations (2.1)). It has since been shown to apply to a broad variety of dyadic coordination phenomena in living systems, e.g. [13,14,19,43,44]. Equations (2.1) can be considered a generalization of the extended HKB model from 2 to N oscillators.
It is remarkable that such a direct generalization can reproduce key features of the collective rhythmic coordination in ensembles of human subjects at multiple levels of description. Another well-studied special case of equations (2.1) is the Kuramoto model [2], which has b ij = 0 (and typically a ij = a, independent of i and j). We will see below, however, that the Kuramoto model cannot exhibit at least one feature of the experimental data. Namely, the data show a secondary peak in the pairwise relative phase of experimental subjects at antiphase (see figure 2d-f ). Simulations using the Kuramoto model do not reproduce this effect, while simulations of equations (2.1) model do (compare figure 5d-f and g-i). We give additional analytical support for this point by studying relevant fixed points of both models in the electronic supplementary materials (section 'Multistability of the present model').

Weak coupling captures human behaviour
Given the spatially symmetric set-up of the 'Human Firefly' experiment (all-to-all network, visual presentation at equal distance to fixation point), it is reasonable to further simplify equations (2.1) by letting a ij = a and b ij = b (a, b > 0), where ϕ ij = w i − w j is the relative phase between oscillators i and j (henceforth we use the notation ϕ ij instead of the subtraction, since relative phase is the crucial variable for coordination [28,29]). The behaviour of the model itself clearly depends on the coupling strength (a, b) and frequency diversity (distribution of ω i 's). While the latter was explicitly manipulated in the human experiment [24], the former was unknown. A qualitative look at simulated dynamics (see examples figure 4a-c for δf = 0.6 Hz) indicates that weak coupling better captures human behaviour (members of the same group do not collapse to a single trajectory in figure 4a as in figure 2). By contrast, stronger coupling (figure 4b,c) deprives the system of much of the metastability. Quantitatively, we fitted the coupling strength (assuming a = b) to the human data based on the level of intergroup integration (β 1 ) (see distribution of model β 1 in figure 4d) particularly for diversity condition δf = 0.3 Hz (i.e. using only one-third of the data to prevent overfitting). We show below how the model captures human behaviour across all diversity conditions and levels of description under the best-fit coupling strength (a = b = 0.105; see section 'Choosing the appropriate coupling strength' in the electronic supplementary materials for more details).
Here we did not estimate the critical diversity δf* the same way as for the human data (by linear interpolation), since we found theoretically that the level of integration depends touchpad LED array (top view) LED array (front view) group A group B Figure 1. Experimental set-up for multiagent coordination. In the Human Firefly experiment [24], eight subjects interacted simultaneously with each other via a set of touch pads and LED arrays. Each subject's movements were recorded with a dedicated touchpad. Taps of each subject were reflected as the flashes of a corresponding LED on the array presented in front of each subject. In each trial, each subject was paced with a metronome prior to interaction. The metronome assignment split the ensemble of eight into two frequency groups of four (group A and B, coloured red and blue, respectively, for illustrative purposes; the actual LEDs are all white). The frequency difference δf between group A and B was systematically manipulated to induce different grouping behaviour. See text and [24] for details. (Online version in colour.) royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 16: 20190360 nonlinearly on diversity δf, and as a result the theoretical δf* is 0.4 Hz (figure 4d). This prediction can be tested in future experiments by making finer divisions between δf = 0.3 and 0.6 Hz.
In the human experiment, not only did we uncover the effect of diversity on intergroup relations, but also, nontrivially, on intragroup coordination (outside affects within, a sign of complexity). Statistically, this is shown in figure 3b (three bars on the left): with the presence of intergroup difference (δf > 0), intragroup coordination was reduced (red, yellow bars significantly shorter than blue bar). This is well captured by the model as shown in figure 3d (two-way ANOVA interaction effect, F 2,19194 = 3416, p < 0.001; the simulated data also capture the rapid decline of intergroup coordination with increasing δf in human data, shown in figure 3b,d, right). In addition to capturing this statistical reduction of intragroup coordination due to intergroup difference, the model, more importantly, provides a window to the dynamical mechanism underlying such statistical phenomena. For example, comparing three simulated trials with identical intragroup properties but different levels of intergroup difference (figure 5a-c), we see that the presence of intergroup difference (figure 5b,c for δf = 0.3, 0.6 Hz) dramatically elevates metastability in the system (compare intermittently converging-diverging dynamics in figure 5b,c to the rather constant behaviour in figure 5a for δf = 0). This suggests that the decrease of intragroup coordination in a statistical sense reflects the increase of metastability in a dynamical sense (see section 'Examples of dynamics with intergroup coupling removed' in the electronic supplementary materials for baseline dynamics when intergroup coupling is removed). Indeed, if we remove intragroup metastability from all simulations (by reducing intragroup frequency variability), they no longer capture the empirically observed statistical result (see section 'Effect of reduced intragroup variability in natural frequency' in the electronic supplementary materials).
At the interpersonal level, human subjects tended to coordinate with each other around inphase and antiphase, especially when the diversity is low (δf = 0 Hz; figure 2d, peaks around ϕ = 0, π separated by a trough near ϕ = π/2); and the preference for inphase and antiphase both diminishes as diversity increases (δf = 0.3, 0.6, figure 2e,f ). Both aspects are well reproduced in simulations of the model (figure 5d-f ). Note that these model-based distributions are overall less dispersed than the more variable human-produced distributions (figure 2d-f ), likely due to the deterministic nature of the model (i.e. no stochastic terms). Yet as demonstrated above, a deterministic model is sufficient for capturing key empirical results at all three levels of description, i.e. coexistence of inphase and antiphase tendencies and their reduction with diversity, reduction of intragroup coordination with the   . A critical coupling ratio κ c = 1 demarcates the regimes of monostability (only all-inphase is stable for κ < 1) and multistability (any combination of inphase and antiphase is stable for κ > 1). This critical ratio (for equation (2.3)) is identical to the critical coupling of the HKB model [29], where the transition between monostability (inphase) and multistability (inphase and antiphase) occurs (equation (2.2); parameters in the two equations map to each other by a ¼ã=2 and b ¼b). This shows how equation (2.3) is a natural N-dimensional generalization of the extended HKB model, in terms of multistability and order-to-order transitions.

The effect of non-uniform coupling
So far, our model has captured very well experimental observations with the simple assumption of uniform coupling. However, loosening this assumption is necessary for understanding detailed dynamics. Here is an example from [24] (figure 6a), where coordination among three agents (1, 3 and 4, labels of locations on LED arrays) is visualized as the dynamics of two relative phases (ϕ 13 where each oscillator can have its own coupling style (oscillator specific coupling strength a i and b i ). In the present case, we are interested in what happens when a 3 ≠ a 4 for i ∈ {1, 3, 4}. Two simulated trials are shown in figure 6b and c with non-uniform versus uniform coupling (same initial conditions and natural frequencies across trials, estimated from the human data). The bumps in ϕ 34 , accompanying dwells in ϕ 13 , are reproduced when a 3 ≫ a 4 (figure 6b) but not when a 3 = a 4 (figure 6c; see section 'Additional triadic dynamics' in the electronic supplementary materials for more analyses). This example shows that to understand interesting dynamic patterns in specific trials, non-uniform coupling strength is important.  and all other parameters identical (members of the slower group, in warm colours, spread evenly within the interval 1.2 ± 0.08 Hz, similarly for members of the faster group, in cold colours, in the interval 1.8 ± 0.08 Hz; initial phases are random across oscillators but the same across trials). When the coupling is too strong (c), all oscillators lock to the same steady frequency. When the coupling is moderate (b), oscillators split into two frequency groups, phase-locked within themselves, interacting metastably with each other (dwell when trajectories are close, escape when trajectories are far apart). When the coupling is weak (a), intragroup coordination also becomes metastable seen as episodes of convergence (black triangles) and divergence. (d ) Level of intergroup integration quantitatively (β 1 , colour of each pixel) for each combination of frequency diversity δf and coupling strength a = b. White curve indicates the critical boundary between segregation (blue area on the left, β 1 < 0, minβ 1 = −0.2) and integration (red and yellow area on the right, β 1 > 0). Within the regime of integration, the yellow area indicates complete integration (β 1 ≈ 1) where there is a high level of phase locking, and the red area indicates partial integration (0 < β 1 ≪ 1) suggesting metastability. Dashed grey lines label δf 's that appeared in the human experiment. Solid grey line labels the empirically estimated critical diversity. (Online version in colour.) royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 16: 20190360

Discussion
The present model successfully captures key features of multiagent coordination in mid-scale ensembles at multiple levels of description [24]. Similar to the HKB model [29], secondorder coupling is demanded by the experimental observation of antiphase (and associated multistability) but now in eightperson coordination; and similar to the extended HKB [37], the model captures how increasing frequency difference δf weakens inphase and antiphase patterns, leading to segregation but now between two groups instead of two persons. This cross-scale consistency of experimental observations may be explained by the scale-invariant nature of the critical coupling ratio κ c = 1, the transition point between monostability (only an all-inphase state) and multistability (states containing any number of antiphase relations). The scale invariance suggests that experimental methods and conclusions for small-scale coordination dynamics have implications for multistability, phase transitions and metastability at larger scales, and enables a unified approach to biological coordination that meshes statistical mechanics and nonlinear dynamics.
Another generalization of the classical Kuramoto model by Hong & Strogatz [45] also allows for antiphase-containing patterns (π-state) by letting the sign of the first-order coupling (a) be positive for some oscillators (the conformists) and negative for others (the contrarians). However, in contrast to our model, antiphase induced this way does not come with multistability, nor the associated order-to-order transitions observed in human rhythmic coordination [13,46]. The second-order coupling in our model allows each individual to be both a conformist and a contrarian but possibly to different degrees [47]. The simple addition of a second stable state may not seem like a big plus at N = 2 (2 stable states), but it rapidly expands the system's behavioural repertoire as the system becomes larger (2 N−1 stable states for N oscillators; with only first-order coupling, the system always has 1 N−1 = 1 stable state, and therefore does not benefit from scaling up). This benefit of scale may be how micro-level multistability contributes to the functional complexity of biological systems [43,48].
Outside of the mathematical context of stability analysis, we have to recall that spontaneous social coordination is highly metastable (e.g. figure 2a) [24], captured by the model when frequency diversity is combined with weak coupling (e.g. figure 4a, in contrast to b,c under stronger coupling). Individuals did not become phase-locked in the long term, but coordinated temporarily when passing by a preferred state (inphase and antiphase) [14,43] (e.g. red trajectory in figure 6a). For N > 2, an ensemble can visit different spatial  (a-c) An example of how intergroup difference may affect intragroup coordination using frequency dynamics of three simulated trials (a = b = 0.105; note that frequency is the time derivative of phase divided by 2π, and consequently the distance between two frequency trajectories reflects the rate of change of the corresponding relative phase, which increases and decreases intermittently during metastable coordination). These three trials share the same initial phases and intragroup frequency dispersion but different intergroup difference i.e. δf = 0, 0.3, 0.6 Hz, respectively. When intergroup differences are introduced (b,c), not only is intergroup interaction altered but intragroup coordination also loses stability and becomes metastable (within-group trajectories converge at black triangles and diverge afterwards). The timescale of metastable coordination also changes with δf, i.e. the inter-convergence interval is shorter for (b) than (c).  [24]), forming patterns that extend in both space and time (electronic supplementary material, figure S4 for intragroup patterns), which further expands the repertoire of coordinative behaviour (see section 'A note on metastability' in the electronic supplementary materials). By allowing complex patterns to be elaborated over time, metastability makes a viable mechanism for encoding complex information as real-world complex living systems do (e.g. the brain) [13,14,22,[49][50][51][52]. By contrast, highly coherent patterns like collective synchronization can be less functional and even pathological [53,54]. Our results call for more attention to these not-quite coherent but empirically relevant patterns of coordination. Besides the multistability or multi-clustering in micro patterns (a general feature endowed by higher-order coupling, e.g. [55][56][57]), existing mathematical studies suggest that the presence of second-order coupling should also manifest at the macro level in large-scale coordination. Naturally, second-order coupling induces multistability of the order parameter in the thermodynamic limit [58][59][60][61]. It also alters the critical scaling of macroscopic order (see [41] for a summary), i.e. for coupling strength K > K c near K c , the order parameter ∥H∥ (norm of the order function [62]) is proportional to (K À K c ) b , with β = 1/2 for the classical Kuramoto model and β = 1 when second-order coupling is added [63,64]. For complex biological systems like the brain which appears to operate near criticality [65], these two types of scaling behaviour may have very different functional implications. When modelling empirical data of biological coordination, one may want to have a closer examination or re-examination of the data for multistability and critical scaling of the order parameter, especially if finer level details are not available.
Key experimental observations are captured by our model under the assumptions of uniform coupling (everyone coordinates with others in the same way) and constant natural frequency, but these assumptions may be loosened to reflect detailed dynamics. For example, introducing individual differences in coupling style (equation (2.5)) gives more room to explain how one metastable phase relation may exert strong influence on another (figure 6a). Long timescale dynamics observed in the experiment (see section 'Additional triadic dynamics' in the electronic supplementary materials) may also be explained by frequency adaptation, which has been observed in dyadic social coordination [66]. A systematic study of the consequences of asymmetric coupling and frequency adaptation on coordination among multiple agents seems worthy of further experimental and theoretical exploration.
To conclude, we proposed a model that captured key features of human social coordination in mid-sized ensembles royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 16: 20190360 [24], and at the same time connects empirically validated large-scale and small-scale models of biological coordination. The model provides mechanistic explanations of the statistics and dynamics already observed, as well as a road map for future empirical exploration. As an experimental-theoretical platform for understanding biological coordination, the value of the middle scale should not be underestimated, nor the importance of examining coordination phenomena at multiple levels of description.

Methods of the human experiment
A complete description of the methods of the 'Human Firefly' experiment can be found in [24]. Here we provide as many details as necessary for understanding the present paper. A total of 120 subjects participated in the experiment, making up 15 independent ensembles of eight people. The protocol was approved by Florida Atlantic University Institutional Review Board and is in agreement with the Declaration of Helsinki. Informed consent was obtained from all participants prior to the experiment.
For an ensemble of eight people, each subject was equipped with a touchpad that recorded his/her tapping behaviour as a series of zeros and ones at 250 Hz (1 = touch, 0 = detach; green rectangles in figure 1), and an array of eight LEDs arranged in a ring (yellow in figure 1), each of which flashed when a particular subject tapped. For each trial, subjects were first paced with metronomes for 10 s, later interacting with each other for 50 s (instructed to maintain metronome frequency while looking at others' taps as flashes of the LEDs). Between the pacing and interaction period, there was a 3 s transient, during which subjects tapped by themselves. Participants were instructed to match their own tapping frequency to the metronome frequency during the 10 s pacing period, and remain tapping at that frequency throughout the rest of the trial even after the metronome disappeared.
During pacing, four subjects received the same metronome (same frequency, random initial phase), and the other four another metronome. The metronome assignments created two frequency groups (say, group A and B) with intergroup difference δf = |f A − f B | = 0, 0.3 or 0.6 Hz (same average ( f A + f B )/2 = 1.5 Hz). This gives rise to three conditions: (1) 1.5 Hz versus 1.5 Hz, (2) 1.65 Hz versus 1.35 Hz, and (3) 1.8 Hz versus 1.2 Hz. Each ensemble completed six trials per condition (a total of 18 trials in random order). From a single subject's perspective, the LED array looks like the legend of figure 2a (all LEDs emit white light; colour-coding only for labelling locations): a subject always saw his/her own taps as the flashes of LED 1, members of his/her own frequency group LED 2-4, and members of the other group LED 5-8 (members from two groups were interleaved to preserve spatial symmetry).
From the tapping data (rectangular waves of zeros and ones), we obtained the onset of each tap, from which we calculated instantaneous frequency and phase. Instantaneous frequency is the reciprocal of the interval between two consecutive taps. Phase (w) is calculated by assigning the onset of the nth tap phase 2π(n − 1), then interpolating the phase between onsets with a cubic spline. The relative phase between the ith and jth subject at time t is ϕ ij (t) = w i (t) − w j (t).

Estimating the distribution of natural frequencies
Human subjects have variable capability to match the metronome frequency and maintain it, which in turn affects how they coordinate. To reflect this kind of variability in the simulations, the oscillators' natural frequencies were drawn from a probability distribution around the 'metronome frequency' (central frequencies f A and f B for groups A and B). To estimate this distribution from human data, we first approximated the 'natural frequency' of each subject in each trial with the average tapping frequency during the transient between pacing and interaction periods (see Methods of the human experiment), and subtracted from it the metronome frequency (see blue histogram in electronic supplementary material, figure S3 from the 'Human Firefly' experiment [24]). We then estimated the distribution non-parametrically, with a kernel density estimator in the form ofP e Ày 2 =2 . Here n = 2072 (259 trials × 8 subjects) from the experiment. We choose the bandwidth h = 0.0219, which is optimal for a normal density function according to [67], where σ is the measure of dispersion, estimated bỹ where y i 's are samples [68]. The result of the estimation is shown in electronic supplementary material, figure S3 (red curve).

Phase-locking value and level of integration
The (short-windowed) phase-locking value (PLV) between two oscillators (say x and y) during a trial is defined as where ϕ xy = w x − w y , W is the number of windows which each ϕ trajectory is split into, and M is number of samples in each window (in the present study, W = 16 and M = 750, same as [24]). Intragroup PLV (PLV intra ) is defined as where A and B are two frequency groups of four oscillators, corresponding to the design of the 'Human Firefly' experiment [24], A = {1, 2, 3, 4}, B = {5, 6, 7, 8} and |A| = |B| = 4. Intergroup PLV (PLV inter ) is defined as (4:6) In both the human and simulated data, comparisons of PLV intra and PLV inter for different levels of δf were done using two-way ANOVA with Type III sums of squares, and Tukey honest significant difference tests for post-hoc comparisons (shown in figure 3b,d ).
The level of integration between two frequency groups is defined based on the relationship between intragroup coordination (measured by PLV intra ) and intergroup coordination (measured by PLV inter ). The groups are said to be integrated if intragroup coordination is positively related to intergroup coordination, and segregated if negatively related. Quantitatively, for each combination of intergroup difference δf and coupling strength a (assuming a = b for our model, assuming b = 0 for the classical Kuramoto model), we use linear regression PLV (df,a) inter,k ¼ b where PLV (df ,a) Á,k is the inter/intra-group PLV for the kth trial simulated with the parameter pair (δf, a), and the slope of the regression line b (df,a) 1 is defined as the measure of the level of integration between two frequency groups. If β 1 > 0, the groups may be said to be integrated; if β 1 < 0, segregated. The set {(df, a)jb (df ,a) 1 ¼ 0} is the critical boundary between the domains of intergroup integration and segregation. Theoretical analyses (section 'Choosing the appropriate coupling strength' in the electronic supplementary materials) show that this measure is meaningful (i.e. reflecting qualitative differences between dynamics; figure 4a-c).

Method of simulation
All simulations were done using the Runge-Kutta 4th-order integration scheme, with a fixed time step Δt = 0.004 for duration T = 50 (matching the sampling interval and the duration of interaction period of the human experiment [24]; second may be used as unit), i.e. for system _ X ¼ f(X), with initial condition X(0) = X 0 , the (n + 1)th sample of the numeric solution can be solved recursively The solver was implemented in CUDA C++, ran on a NVIDIA graphics processing unit, solving every 200 trials in parallel for each parameter pair (δf, a). For each trial, initial phases (of eight oscillators) were drawn randomly from a uniform distribution between 0 and 2π, and natural frequencies from distributions defined by equation (4.2) (reflecting the design of and variability observed in the human experiment [24]). Here 200 trials are used per condition, greater than that of the human experiment (see [24] and section 'Design of the human experiment' in the electronic supplementary materials for details) to obtain a more accurate estimate of the mean.
Ethics. The protocol of the human experiment was approved by Funding. The authors of this work were supported by NIMH grant no. MH080838, NIBIB grant no. EB025819, the FAU Foundation, FAU I-SENSE and FAU Brain Institute.