Entrainment and synchronization in networks of Rayleigh–van der Pol oscillators with diffusive and Haken–Kelso–Bunz couplings

We analyze a network of non-identical Rayleigh–van der Pol (RvdP) oscillators interconnected through either diffusive or nonlinear coupling functions. The work presented here extends existing results on the case of two nonlinearly coupled RvdP oscillators to the problem of considering a network of three or more of them. Specifically, we study synchronization and entrainment in networks of heterogeneous RvdP oscillators and contrast the effects of diffusive linear coupling strategies with the nonlinear Haken–Kelso–Bunz coupling, originally introduced to study human bimanual experiments. We show how convergence of the error among the nodes’ trajectories toward a bounded region is possible with both linear and nonlinear coupling functions. Under the assumption that the network is connected, simple, and undirected, analytical results are obtained to prove boundedness of the error when the oscillators are coupled diffusively. All results are illustrated by way of numerical examples and compared with the experimental findings available in the literature on synchronization of people rocking chairs, confirming the effectiveness of the model we propose to capture some of the features of human group synchronization observed experimentally in the previous literature.


Introduction
Interpersonal coordination and synchronization between the motion of two individuals have been extensively studied over the past few decades (Schmidt and Turvey 1994;Varlet et al. 2011). Synergetic movements of two or more people mirroring each other frequently occur in many activities such as handling objects, manipulating a common workpiece, dancing, choir singing, and movement therapy (Himberg and Thompson 2009;Valdesolo et al. 2010;Mörtl et al. 2012;Lorenz et al. 2013;Repp and Su 2013). It is of great importance to reveal not only the effects of mirroring movements among people on human physiological and mental functions, but also to understand the link between intrapersonal and interpersonal coordination. In social psychology, it has been shown that people prefer to team up with others possessing similar morphological and behavioral features and that they tend to coordinate their movement unconsciously (Folkes 1982;Lakens and Stel 2011). Moreover, much evidence suggests that interpersonal motor coordination is strictly related to social attachment, meaning that synchronous activities between individuals, that occur even more when the kinematic features of their movements share similar patterns (Słowiński et al. 2016), may produce positive emotions (Wiltermuth and Heath 2009).
In order to explain the experimental observations of human interpersonal coordination, mathematical models are usually derived to capture the key features of the observed behavior. A classical example is the nonlinear RvdP oscillator with HKB coupling, which was introduced in Haken et al. (1985) to explain the transition from phase to antiphase synchronization in bimanual coordination experiments (for more details, see Kelso et al. 1987;Jirsa et al. 1998). Such model was shown to be able to capture many features of human coordination even beyond the bimanual synchronization experiments it was derived to explain. For example, RvdP oscillators were used in Zhai et al. (2014aZhai et al. ( ,b, 2015a; Alderisio et al. (2016) to design virtual players in the context of the mirror game (Noy et al. 2011), presented as an important paradigmatic case of study to investigate the onset of social motor coordination between two players imitating each other's hand movements. Furthermore, in Kelso et al. (2009) the authors take inspiration from the dynamic clamp of cellular and computational neuroscience in order to probe essential properties of human social coordination by reciprocally coupling human subjects to a computationally implemented model of themselves (RvdP oscillator with HKB coupling), referred to as virtual player. Such concept, namely the human dynamic clamp, was further investigated and developed in Dumas et al. (2014) in order to cover a broader repertoire of human behavior, including rhythmic and discrete movements, adaptation to changes of pacing, and behavioral skill learning as specified by the virtual player. Besides, RvdP oscillators were also used in Schmidt and Turvey (1994) in order to capture the rhythmic coordination between two individuals swinging handheld pendulums, in Varlet et al. (2011) in order to model spontaneous interpersonal postural coordination between two human people and account for the competition between the coupling to a visual target to track and the coupling to the partner, in Richardson et al. (2007) in order to qualitatively explain interpersonal movement synchronization between two human beings involved in rhythmic paradigms, and in Amazeen et al. (1995) in order to account for the frequency detuning of the phase entrainment dynamics of two people involved in interlimb coordination.
While coordination of two human players has been studied in numerous previous investigations, the case of multiple human players has been seldom studied in the existing literature, due to a combination of practical problems in running the experiments and lack of a formal method able not only to model the considered scenario but also to quantify and characterize the synchronization level of the ensemble. Multiplayer games involve a group of three or more people engaged in a communal coordination task. The variety of scenarios that can be considered is vast due to the countless activities the players might be involved in (limb movements, finger movements, head movements, walking in a crowd, or more in general music and sport activities), the many ways in which the participants can interact and communicate with each other and the different ways all the players can be physically located with respect to each other while performing the specified task.
Some of the existing works on coordination of multiple human players include studies on choir singers during a concert (Himberg and Thompson 2009), rhythmic activities as, for example, "the cup game" and marching tasks (Iqbal and Riek 2015), rocking chairs (Frank and Richardson 2010;Richardson et al. 2012), and coordination of rowers' movements during a race (Wing and Woodburn 1995). In these papers, the authors provide several experimental results in order to analyze the behavior of a group of people performing some coordinated activities, but a rigorous mathematical model capable of capturing the observed results and explaining the features of the movement coordination among them is still missing. In particular, in Frank and Richardson (2010) the authors study both unintentional and intentional coordination by asking the players to try and synchronize the oscillations of the rocking chairs they are sitting on with their eyes shut or open. Synchronization is observed to spontaneously emerge when players observe each other's movements. Another study in which multiplayer activities are analyzed but a mathematical model is missing is carried out in Yokoyama and Yamamoto (2011) where the theory of symmetric Hopf bifurcations in coupled oscillators is used to investigate the synchronized patterns of three people during sport activities.
Further results about multiplayer activities deal with spontaneous group synchronization of arm movements and respiratory rhythms. For example, in Codrons et al. (2014) the authors test whether pre-assigned arm movements performed in a group setting spontaneously synchronize and whether synchronization extends to heart and respiratory rhythms. In their study, no explicit directions are given on whether or how the arm swingings are to be synchronized among participants, and experiments are repeated with and without external cues. Interestingly, when an external auditory rhythm is present, both motor and respiratory synchronization is found to be enhanced among the group. Also, the overall coordination level is observed to increase when compared to that detected when the same experiments are again carried out in the absence of the external cue.
The main objective of this paper is to propose and analyze a model able to account for movement synchronization in multiplayer scenarios and explain some of the features observed experimentally in the existing literature. Specifically, we consider networks of heterogeneous nonlinear RvdP oscillators as a good model of multiplayer coordination and, as already done in Mörtl et al. (2012) for the case of two agents only, we regard it as a synchronization problem. Each equation is used to model the movement of a different player and is therefore characterized by a different set of parameters to account for human-to-human variability. The effects of different interaction models, linear and nonlinear, are investigated to understand under what conditions synchronization is observed to emerge. Our analysis suggests that bounded synchronization is indeed a common emergent property in these networks whose occurrence can also be accounted for analytically in a number of different scenarios. Also, as expected from existing theoretical results, we find that the structure of the interactions among players has an effect on the coordination level detected in the network.
Furthermore, the effects of adding an external sinusoidal signal are studied in order to understand whether synchronization can be improved by means of an entrainment signal (Russo et al. 2010). Our analysis suggests that the coordination level of the ensemble can indeed increase when the oscillation frequency of the external signal is similar to the natural angular velocity of the agents in the network. However, in all the other cases, the external signal acts as a disturbance and leads to a decrease in the coordination among the agents.
We wish to emphasize that the study reported in this paper will form the basis of future experimental investigations which are currently being planned.
The rest of the paper is organized as follows. In Sect. 2, some notation that will be used in later sections is introduced. In Sect. 3, the equation that describes the network is presented, in terms of both internal dynamics of each agent and coupling function thanks to which they can interact with each other. In Sect. 4, some metrics are introduced to characterize the quality and the level of coordination in human groups. In Sect. 5, a testbed scenario of multiplayer coordination in networks of human people is presented, while in Sect. 5.1, the key synchronization features experimentally observed are reproduced by considering a network of heterogeneous RvdP oscillators, and the effects of three different coupling strategies thanks to which they are interconnected are explored. In Sect. 6, the effects of adding an external entrainment signal are analyzed with respect to the overall coordination level of the network. In Sect. 7, bounded synchronization of the network when its nodes are connected through a linear diffusive coupling function is analytically proven to be achieved, and some numerical examples are provided in order to both illustrate the effectiveness of our analysis and to show that bounded synchronization can be achieved also when considering different couplings. Finally, in Sect. 8, a summary of our results and some possible future developments are presented.

Preliminaries and background
We denote with ⊗ the Kronecker product between two matrices. The operator λ k (·) defined over a matrix indicates the kth eigenvalue of the matrix itself, and λ M (·) indicates its maximum eigenvalue when the matrix is real and symmetric and as a consequence all the eigenvalue are real as well.
A graph is a tuple G = {V, E} defined by a set of nodes V = {1, . . . , N } and a set of edges E ⊆ V × V. A graph is said to be undirected if (i, j) ∈ E ⇐⇒ ( j, i) ∈ E. In an undirected graph, two nodes i and j are said to be neighbors is called adjacency matrix, and a i j ≥ 0 is called strength of the interaction between the pair (i, j). In particular, a graph is said to be unweighted if the interaction between two neighbors is equal to 1. A path between nodes h and k is a sequence of nodes, with h and k as endpoints, such that every two consecutive nodes are neighbors. A graph is said to be simple if a ii = 0 ∀i ∈ V, and it is said to be connected if there exists a path between any two of its nodes. The matrix L = {l i j } ∈ R N ×N defined as is called Laplacian matrix of the graph (or simply Laplacian).
The Laplacian of any simple undirected graph is symmetric with zero row sum and is a positive semidefinite matrix with as many null eigenvalues as there are components in the graph. In particular, a connected graph has only one null eigenvalue.
Throughout the paper, we will consider a connected simple undirected network of N agents assuming that any two players interact symmetrically with one another.
Before analyzing a multiplayer scenario, it is worth considering the simpler case of only two human players interacting with each other. The system that can be used to model the interaction between them is described in terms of two coupled RvdP oscillators and can be given as follows Fuchs and Jirsa 2008): where x i ∈ R denotes the position of the ith player, with i = 1, 2. The right-hand side of both equations represents the coupling term between the two players: in particular, using the model proposed by Haken, Kelso and Bunz in Haken et al. (1985), it is given by The term αx 2 i + βẋ 2 i − γ ẋ i represents the nonlinear damping of the oscillatory movement of player i. Specifically, the sign of γ determines whether, in the absence of coupling, the oscillation is persistent (γ > 0) or vanishes (vice versa) as time goes by: it is trivial to verify this by studying the stability of the origin and checking the sign of the eigenvalues of the Jacobian of the system (Avitabile et al. 2015). Moreover, α and β determine the amplitude of such oscillation, while ω i is related to its frequency. It has been proven that this model of two nonlinearly coupled oscillators accounts for the features observed during experimental data in bimanual experiments (see Haken et al. (1985) for further details).

Human-to-human coordination as a synchronization problem
We have pointed out that the dynamics of two coupled RvdP oscillators has been used to describe different kinds of interpersonal coordination tasks between two people, including bimanual coordination experiments, mirror game, social postural coordination, and rocking chairs. According to the particular scenario considered, the state vector of each oscillator is used to represent position and velocity of the particular body part of interest of either of the players (finger, hand, head, and so forth). Following the same approach, we can consider a scenario in which more than two human beings are performing a multiplayer coordination task involving some oscillatory motion, as, for example, arm or hand rhythmic movements, rocking chairs, and head tracking of a visual target. In these cases, the state vector of each node represents position and velocity of the particular body part of interest of each player. Therefore, the dynamics of each player when moving in isolation will be described by the following nonlinear system: where x i = [x i 1 x i 2 ] T ∈ R 2 is the state vector, with x i 1 , x i 2 representing position and velocity of the ith human player, respectively. To model the interaction between different players, we assume that the dynamics of each of them is affected by some coupling function u i which depends on the difference between the state of the ith player and that of his/her neighbors. In what follows, we will explore the effects of three possible selections for such a function. We are interested in analyzing which one leads to synchronization features which are the closest to those observed in previous experimental work about human ensembles involved in a joint coordination task, e.g. (Richardson et al. 2012).
1. Full state coupling. With this kind of coupling, we assume that players adjust both their velocities and accelerations proportionally to the average mismatch between their own position and velocity and those of their neighbors. Mathematically, we have: In particular, N i > 0 is the number of neighbors of node i, while c > 0 is the coupling strength among the agents. 2. Partial state coupling. Next, we explore the case where players only adjust their accelerations according to the position and velocity mismatches from their neighbors: In particular, N i > 0 is the number of neighbors of node i, while c 1 , c 2 > 0 represent the position and the velocity coupling strengths, respectively. 3. HKB coupling. Finally, we consider an interaction model which is the direct extension to multiplayer coordination problems of the interaction function used to describe the bimanual experiments (Haken et al. 1985;Fuchs et al. 1996). Specifically, we choose the following nonlinear function: Once again, N i > 0 is the number of neighbors of node i, while c > 0 represents the coupling strength among the agents.
The resulting network model describing the interaction of a group of N players can then be written aṡ where the coupling function u i can be chosen as one of those listed above. We now explore under what conditions coordination, and hence synchronization, emerges for each of the three scenarios of interest. We wish to emphasize that, since the node parameters are heterogeneous, complete synchronization as defined in Li et al. (2010) cannot be achieved. We will consider instead the case where bounded synchronization, as defined in Hill and Zhao (2008) and below, emerges. Namely, we define the average trajectory as and its distance from the state of each node i as We also define the parameters vector for each node i as ϑ i := [α i β i γ i ω i ] T ∈ R 4 , and we introduce the stack vectors Definition 1 We say that a network of non-identical RvdP oscillators achieves bounded synchronization if and only if there exists some time instantt such that for any initial condition x i,0 and parameter vector ϑ i of the nodes in the network.
Definition 2 If a network of non-identical RvdP oscillators achieves bounded synchronization, we define the relative synchronization error bound χ as the upper bound of the ratio between the error norm η(t) when the oscillators are coupled, and its maximum valueη when they are uncoupled, that is According to Definition 1, we will say that the network has achieved bounded synchronization if the error norm reaches and remains into a bounded region for all t >t. We will then use the relative synchronization error bound χ , as described in Definition 2, in order to evaluate the improvement in synchronization among the oscillators when coupling through the network is present. Assuming that when the nodes are coupled η(t) ≤η, we have that χ ∈ [0, 1]. In particular, χ = 0 represents the ideal case of complete synchronization, while χ = 1 represents the worst-case scenario. The relative synchronization error bound can be lowered, in general, by increasing the strength of the coupling function (e.g., see Theorem 2 and Sect. 7.3).

Coordination metrics
In order to quantify and analyze the coordination level in a network of more than two agents, we use the metrics introduced in Richardson et al. (2012) to characterize the quality of synchronization in human groups.
Let x k (t) ∈ R ∀t ∈ [0, T ] be the continuous time series representing the motion of each agent, with k ∈ [1, N ], where N is the number of individuals and T is the duration of the experiment. Let x k (t i ) ∈ R, with k ∈ [1, N ] and i ∈ [1, N T ], be the respective discrete time series of the kth agent, obtained after sampling x k (t), where N T is the number of time steps and T := T N T is the sampling period. Let θ k (t) ∈ [−π, π] be the phase of the kth agent, which can be estimated by making use of the Hilbert transform of the signal x k (t) (Kralemann et al. 2008). We define the cluster phase or Kuramoto order parameter, both in its complex form q (t) ∈ C and in its real form q(t) ∈ [−π, π] as which can be regarded as the average phase of the group at time t.
be the relative phase between the kth participant and the group phase at time t. We can define the relative phase between the kth participant and the group averaged over the time interval [t 1 , t N T ], both in its complex formφ k ∈ C and in its real formφ k ∈ [−π, π] as In order to quantify the degree of synchronization for the kth agent within the group, we define the following order parameter which simply gives information on how much the kth agent is synchronized with the average trend of the group. The closer ρ k is to 1, the better the synchronization of the kth agent itself.
In order to quantify the coordination level of the entire group at time t, following Richardson et al. (2012), we define the following parameter which simply represents the group synchronization, with is to 1, the better the coordination level of the group at time t. Its value can be averaged over the whole time interval [0, T ] in order to have an estimate of the mean coordination level of the group during the total duration of the performance: Besides, if we denote with φ d k,h (t) := θ k (t) − θ h (t) the relative phase between two participants in the group at time t, it is possible to estimate their dyadic synchronization, that is the coordination level between participants k and h over the total duration T of the trial (Richardson et al. 2012): Note that high dyadic coordination levels can coexist with low group synchronization values. Also note that, in general, the tighter bounded synchronization is (i.e., the smaller in Eq. 11 or equivalently χ in Eq. 12), the higher is the group synchronization index defined in Eq. 18. Intuitively, the closer the trajectories x k (t) of all the nodes are, the closer are their phases θ k (t), and so are their differences with the phase of the group φ k (t) and their respective values averaged over timeφ k . A numerical example of such relationship is reported in Sect. 7.3.
Remark 1 For strong coupling, the HKB model entails bistable coordination as an important feature; hence, it is possible that the network converges toward regimes where some of the oscillators evolve in anti-phase. According to our definitions, these scenarios correspond to unwanted situations. Indeed, we focus on the case where all the oscillators exhibit relatively small phase differences corresponding to large values of the group synchronization index (typically over 90 %).

Testbed example
As a testbed scenario, we consider the synchronization of rocking chairs motion studied in Richardson et al. (2012). In particular, participants sit on six identical wooden rocking chairs disposed as a circle and are supposed to rock them in two different conditions: 1. Eyes closed: participants are required to rock at their own preferred frequency while keeping their eyes closed; 2. Eyes open: participants are required to rock at their own preferred frequency while trying to synchronize their rocking chair movements as a group.
In the eyes closed condition, the participants are not visually coupled, meaning that the oscillation frequency of each of them is not influenced by the motion of the others, while in the eyes open condition each player is asked to look at the middle of the circle in order to try and synchronize their In Fig. 1, the typical trend of the group synchronization index ρ g (t) and its mean value and standard deviation are depicted for each of the three aforementioned trials. In particular, in Fig. 1a it is possible to appreciate that the mean value ρ g of the group synchronization index is around 0.4 in the eyes closed condition, while it is around 0.85 in the eyes open condition. This means that, when the participants are not visually coupled, as expected synchronization does not emerge, while when the participants are visually coupled and explicitly told to rock their own chair movements as a group, the coordination level significantly increases. In Fig. 1b, it is possible to appreciate that in the eyes closed condition the amplitude of the oscillations of the group synchronization is higher than that observed during the eyes open trials.
In Table 1, we show typical values of the degree of synchronization ρ k of the participants involved in the rocking chairs experiments, for both the eyes closed and the eyes open condition. It is possible to appreciate how, as expected, the value of ρ k is much higher for almost all the participants when they are visually coupled. Interestingly enough, agent 6 does not undergo an improvement of ρ 6 with respect to the eyes closed condition, meaning that such participant has In Table 2, we show typical values of the dyadic synchronization ρ d k,h between participants involved in the rocking chairs experiments for the eyes open condition. As expected, the lowest values are those corresponding to the participant that struggled the most to synchronize with the rest of the group, that is participant 6.

Modeling results
In this section, we uncover the synchronization features that the three different coupling functions introduced earlier lead to, with respect to the rocking chairs experiments introduced earlier as a testbed scenario (Richardson et al. 2012). We will explore whether and how the model of coupled RvdP oscillators we propose in this paper can reproduce the key features of the observed experimental results. In so doing, we will explore: • the effects of choosing different coupling functions; • how varying the coupling strength affects the coordination level of the agents.
We simulate a network of N = 6 heterogeneous RvdP oscillators whose parameters and initial values are heuristically set as described in Table 3 and we set T = 200 s. We suppose that the network is simple, connected, unweighted, and undirected and we assume that each node is connected to all the others (complete graph), which we believe well represents the topology implemented in the rocking chairs experiments of Richardson et al. (2012) for the eyes open condition.
In particular, since we are interested in replicating the key features of the rocking chairs experiments for both conditions (eyes open and eyes closed), in Fig. 2 we show the group synchronization obtained with and without coupling between the agents. In particular, in Fig. 2a we show the mean value and standard deviation of ρ g (t) for the case in which the nodes are uncoupled, and for the cases in which they are connected through each of the three coupling functions presented earlier. For each of the cases being considered, after a simulation of duration T , the standard deviation σ ρ g is computed as where ρ g is the mean value of the group synchronization index defined in Eq. 19. From Fig. 2a, it is possible to appreciate that in the absence of connections among the nodes, which corresponds to u i = 0 ∀i ∈ [1, N ], the group synchronization has a mean value approximately equal to 0.4, while it significantly increases (approximately 0.9) when connecting the nodes with any of the three coupling functions introduced above. These results confirm the observations previously made for a network of six human people involved in rocking chairs experiments (Fig. 1a). In particular, we heuristically chose c = 0.15 for the full state coupling, c 1 = c 2 = 0.15 for the partial state coupling, and a = b = −1, c = 0.15 for the HKB coupling in order to match the experimental observations. In Fig. 2b, we show the time evolution of the group syn- of brevity, we do not show the evolutions of ρ g (t) obtained with partial state and HKB coupling as they are similar to that obtained with full state coupling. Our model is able to reproduce another key feature observed in Richardson et al. (2012): when the nodes are uncoupled, which corresponds to the eyes closed condition, the amplitude of the oscillations of ρ g (t) is higher than that obtained when the nodes are coupled, which corresponds to the eyes open condition (Fig. 1b). However, the way a network of oscillators is coupled might in general affect the overall coordination level of the agents for certain values of the coupling strength, especially because of the different nature (linear or nonlinear) of the coupling. Indeed, when increasing the value of the coupling strength, the group synchronization ρ g (t) turns out to be qualitatively different in each of the three different cases (Fig. 2c). The oscillations observed in the rocking chairs experiments at steady state are preserved only when using partial state coupling or HKB coupling, with the latter yielding results which are more similar to those observed experimentally in Richardson et al. (2012) in terms of the oscillations amplitude.
In Table 4, we show the degree of synchronization ρ k obtained for each node of the network, both in the absence of coupling among the agents and in its presence. It is possible to appreciate how, for each node k in the network, ρ k has much higher values when some coupling is present, confirming what observed in Richardson et al. (2012) when contrasting synchronization levels when the participants have their eyes open or closed. Moreover, we observe that, as in the experiments, despite the group synchronization index taking higher values when the nodes are coupled, there are some agents that struggle to keep up with the general trend of the group, therefore showing lower values in terms of ρ k (node 5 in our simulations).
The fact that the degree of synchronization of node 5 is lower than that of the others can be explained by noting that its parameters are the furthest from the average values, with respect to all the other nodes in the network. In particular, if we define ω * i = ω 2 i and the average values of all the parameters as follows we can evaluate their percentile error from the mean values through the following quantities: As a measure of the overall distance with respect to the average values of the parameters, we use the norm of the vector Table 5, it is possible to appreciate that ||θ 5 || = max i ∈[1,N ]θi and therefore player 5 is the furthest away from the rest of the group in parameter space.
In Table 6, we show the dyadic synchronization ρ d k,h for all the possible couples of nodes in the network: again, our simulation results confirm what observed for the rocking chairs experiments. Indeed, the pairs of nodes with lower dyadic synchronization correspond to pairs in which at least one of the two nodes had trouble synchronizing with the general trend of the group (node 5 in our simulations). For the sake of clarity, we show ρ d k,h only when connecting the nodes through full state coupling. Analogous results can be obtained also for the other two strategies introduced earlier in this paper. It is easy to foresee that, regardless of the coupling function the nodes are connected through, the value of the coupling strength has a direct impact on the group synchronization in terms of its mean value and its standard deviation. We now show how ρ g varies quantitatively as the coupling strength varies for all the three coupling functions introduced earlier in this paper, when considering an unweighted complete graph of N = 6 heterogeneous RvdP oscillators whose parameters and initial values are those given in Table 3. Moreover, we once again set T = 200 s. In Figs. 3, 4 and 5, we show the mean value and standard deviation of the group synchronization index ρ g (t) obtained for different values of the coupling strength when considering full state coupling, partial state coupling and HKB coupling.
From Fig. 3a, it is clear that, when considering full state coupling, the group synchronization index increases as the coupling strength c increases: in particular, a relatively small value of the coupling strength is sufficient for the network to synchronize within an acceptable bound (c 0.15 leads to ρ g 0.9, see Fig. 3b). This means that the stronger the influence of each player on the others is, the better the overall synchronization of the human participants.
In Fig. 4a, we show how, when considering a partial state coupling, the group synchronization index varies for increasing values of c 1 while keeping c 2 constant and equal to 0, and vice versa in Fig. 4b. It is possible to appreciate how the influence that c 2 has on the group synchronization is stronger than that provided by c 1 . A possible interpretation of this effect is that human players react better to changes in the velocity of their neighbors rather than to changes in their position. This result is also confirmed in terms of the mean value of the group synchronization index as c 1 and c 2 are simultaneously varied (Fig. 6a).
Finally from Fig. 5a it is clear that, when considering HKB coupling while keeping a and b constant and equal to −1, the group synchronization increases as the coupling strength c increases. In particular, like in the case of full state coupling, in order for the network to well synchronize it is sufficient to choose a relatively small value for the coupling strength (c 0.15 leads to ρ g 0.9, see Fig. 5b). A possible interpretation of this effect is that the stronger the influence that each player has on the others, the better the overall synchronization of the human participants. This result is confirmed also in Fig. 6b in which we show how the mean value of the group synchronization changes as a and b are simultaneously varied while keeping c constant and equal to 1. It is possible to appreciate how, as the values of |a| and |b| increase, then so does the average of the group synchronization index.

Entrainment of the network
In this section, we analyze the effects on the group synchronization of an entrainment signal affecting all nodes, modeled Our main objective is understanding whether, and possibly under what conditions, such entrainment signal leads to better coordination level in a network of heterogeneous RvdP oscillators when compared to the case in which the signal is absent. This will help us uncover whether the presence of an external auditory or visual stimulus can improve synchronization in a group of individuals performing some joint coordination task.
Following the approach of Russo et al. (2010), we model such a scenario in the following way: where ζ(t) = A ζ sin ω ζ t represents the entrainment signal and u i (t) one of the coupling functions introduced earlier in this paper. We introduce the entrainment index ρ E ∈ [0, 1] in order to quantify the overall coordination level between the network and the external signal ζ(t): is the phase of the kth node, θ ζ (t) is the phase of ζ(t), T is the duration of the experiment, and N is the number of nodes in the network. The closer ρ E is to 1, the stronger the entertainment of the group with the external signal.
The simulation scenario is the same as that described in Sect. 5.1, with the addition of the entrainment signal.
In Fig. 7, we show the values of the entrainment index for different values of the frequency ω ζ and the amplitude A ζ of the entrainment signal ζ(t) when considering full state coupling with c = 0.15. It is possible to appreciate how, for each value of ω ζ , the entrainment index increases as A ζ increases, meaning that the network better synchronizes with ζ(t) for increasing values of its amplitude. Moreover, for a given value of A ζ , the highest values of ρ E are achieved when the frequency of the entrainment signal is close to the average value Ω of the natural frequencies ω i of the nodes (in this case Ω 0.5 ). These results confirm the findings of Schmidt et al. (2007); Varlet et al. (2015), in which it is shown that spontaneous unintentional synchronization between the oscillation of a handheld pendulum swung by an individual and an external sinusoidal stimulus (which corresponds to our external entrainment signal) emerges only when the frequency of the signal itself is similar to the preferred frequency of the player. For the sake of brevity, we do not show how ρ E varies as ω ζ and A ζ vary when considering partial state coupling and HKB coupling, since we obtain results which are similar to those shown in Fig. 7 for full state coupling.
In Fig. 8, we show the mean value and standard deviation of the group synchronization index ρ g (t) when considering different parameters of the entrainment signal for all the three coupling functions we have presented. As we are interested in understanding whether an additive external sinusoidal signal can improve the coordination level of the network, the values of the coupling strengths chosen in these simulations for all the three coupling functions are the same as those previously used in the absence of entrainment signal (Fig. 2a). From Fig. 8, it is easy to observe that, for all coupling functions, the group synchronization index of the network improves only when the entrainment index ρ E has high values (see black line compared to the green one). In the other two cases (blue line and red line), the entrainment signal acts as a disturbance for the dynamics of the nodes and the group synchronization index decreases. This confirms that it is possible to further enhance the coordination level of participants only when the entrainment signal has an oscillation frequency which is close to the average of the natural oscillation frequencies of the individuals involved and its amplitude is sufficiently high.

Convergence analysis
As anticipated earlier, we are also interested in understanding under what conditions synchronization is theoretically observed to emerge. In particular, in this section we are going to show that bounded synchronization can be analytically guaranteed for a network of N diffusively coupled non-identical RvdP oscillators by making use of two dif-ferent approaches, namely contraction theory and Lyapunov theory.

Contraction theory
Let | · | be a norm defined on a vector w ∈ R n with induced matrix norm || · ||. As stated in Russo et al. (2013), given a matrix P ∈ R n×n , the induced matrix measure is defined as μ(P) := lim h→0 + (||I +h P||−1) h . Definition 3 Let us consider the systemẇ = F(t, w) defined ∀t ≥ 0, ∀w ∈ C ⊂ R n . We say that such system is contracting with respect to a norm | · | with associated matrix measure μ(·) iff where J is the Jacobian of the system.
The key stage in the application of contraction theory to synchronization of networks of oscillators is the construction of the so-called virtual system (Jouffroy and Slotine 2004).
Definition 4 Let us consider a network of identical systems. The virtual system is defined as a dynamical system expressed in some auxiliary variable set whose particular solutions are the trajectories of each of the nodes in the network.
Formally, the virtual system depends on the state variables of the agents in the network and on some virtual state variables. Substitution of the state variables of a certain node i with the virtual ones returns the dynamics of the ith node of the network itself. Then, if the virtual system is contracting, complete synchronization is achieved (Russo and di Bernardo 2009;Wang and Slotine 2005 ;di Bernardo et al. 2016).
As noted in Russo and di Bernardo (2009), virtual systems can also be used to study networks of oscillators with parameter mismatches. Indeed, it is shown that if a network of identical systems achieves complete synchronization, then it achieves bounded synchronization when parameter mismatches among the nodes are present (see Russo and di Bernardo 2009 for further details and example of application to networks of heterogeneous biological oscillators).
In Russo and di Bernardo (2009), a simple algorithm is provided that allows to check whether the virtual system of a certain network of N non-identical agents is contracting, which leads to bounded synchronization of the network itself. In particular, rather than verifying Eq. 24 in order to check whether the virtual system is contracting, the algorithm consists in checking the truth of some statements regarding the single elements of the Jacobian of the virtual system and imposing some conditions: 1. build the Jacobian J of the virtual system; 2. check whether the following statements are true or false 3. generate a set of conditions for synchronization (CFS) according to the truth or the falsity of the previous statements.
In particular, denoting with n 0 i the number of 0 elements in the ith row of the Jacobian of the virtual system, the CFS are generated in the following way: Note that if statement S1 is not true, it is not possible for the virtual system to be contracting.

Theorem 1 Suppose to have a network of N heterogeneous
RvdP oscillators interconnected via full state coupling as described in Eq. 5. Let us also assume that the network topology is a connected, simple, undirected, and unweighted complete graph. If the following hypothesis of inequalities is satisfied whereα,ω,γ are the average values of parameters α i , ω i , γ i , respectively, and z 1 max , z 2 max are the bounds of the two virtual state variables introduced below in Eq. 27, then bounded synchronization is achieved by the network.
Proof Let us consider an unweighted complete graph of N heterogeneous RvdP oscillators interconnected via full state coupling, that isẋ where x i ∈ R 2 is the state variable of node i andĉ := c N −1 since each node in a connected complete graph has N − 1 neighbors. The virtual system corresponding to this network can be written aṡ where z ∈ R 2 is the vector of auxiliary variables of the virtual system andα : values of all the different parameters of the oscillators in the heterogeneous network. The Jacobian of the virtual system is: In order to prove bounded synchronization of the network, we need the virtual system to be contracting. In order to do so, we apply the algorithm presented in Russo and di Bernardo (2009) to Eq. 28. When i = 1, j = 2, it is immediate to see that statement S1 is true, while S2 and S3 are false (c might be in general time varying), leading to |J 12 | > |J 11 | and |J 21 | < |J 22 |. When i = 2, j = 1 instead, inequalities to satisfy and the associated CFS depend on the sign ofα andβ. Supposing without loss of generality thatα,β > 0 as usually done in the literature (Fuchs and Jirsa 2008;Kelso et al. 2009), it is immediate to see that an inequality corresponding to the fulfillment of S1 needs to be added to the list of CFS generated by the algorithm (a worst-case scenario is −ĉN +γ < 0) and that both S2 and S3 are again false, leading to the two same conditions. This means that the network achieves bounded synchronization when the following system of inequalities is satisfied: Supposing that the dynamics of the virtual system is bounded, meaning that |z 1 (t)| ≤ z 1 max , |z 2 (t)| ≤ z 2 max ∀t ≥ 0, we can consider the following worst-case scenario which, sinceĉ = c N −1 , yields Eq. 25. So we can conclude that if the coupling strength c fulfills Eq. 25, the network of heterogeneous RvdP oscillators overlying a complete graph achieves bounded synchronization.
Remark 2 Note that when the number of nodes in the network N tends to infinity, then N −1 N → 1. This means that bounded synchronization is achieved when: 2αz 1 max z 2 max +ω 2 +γ < c < 1. (31)

Lyapunov theory
Let D be the set of diagonal matrices, D + the set of positive definite diagonal matrices, and D − the set of negative definite diagonal matrices. Let us now define QUAD and QUAD Affine vector fields (DeLellis et al. 2015).
Definition 5 Given n × n matrices P ∈ D + , W i ∈ D, the vector field f i is said to be QUAD(P, W i ) iff for any z, w ∈ R n and for any t ≥ 0.
Definition 6 Given n × n matrices P ∈ D + , W i ∈ D, the vector field f i is said to be QUAD (P, W i Let us consider a network of N non-identical agents interconnected via linear coupling: where Γ ∈ R n×n . Note that this is a generalization of the full state coupling previously introduced, which can be obtained by setting Γ = I n . As reported in DeLellis et al. (2015) in details, in order to prove bounded synchronization for a network of N non-identical QUAD Affine systems interconnected via a linear coupling function, we need h i (t, x i ) to be QUAD(P, W i ) with W i ∈ D − for all the nodes in the network at any time instant. However, in a network of N heterogeneous RvdP oscillators with vector fields described by Eq. 4, regardless of the way we define h i and g i it is never possible to satisfy the following condition with definite negative matrices W i . Indeed, the right-hand side is always negative, while the left-hand side can be positive for any value of P > 0. On the other hand, in order to avoid conditions on the sign of the matrices W i , it is necessary to write the dynamics of the nodes in the following way , ∀z ∈ R n , and with all the terms g i being bounded, at any time instant. In particular, in DeLellis et al. (2015) the authors prove the following theorem.

Theorem 2 Let us consider a network of N non-identical agents interconnected via linear coupling as described in
Eq. 33. Let us suppose that f i (t, and that: 1. the network is made up of N QUAD(P, W ) Affine systems, with P ∈ D + and W ∈ D; 2. Γ is a positive semidefinite diagonal matrix; 3. if W is made up of l ∈ [0, n] nonnegative elements, which without loss of generality can be collected in its l × l upper-left block, then Γ is made up ofl positive elements, where l ≤l ≤ n, which again without loss of generality can be collected in itsl ×l upper-left block; 4. ∃ 0 <M < ∞ such that ||g i (t, x i )|| 2 <M ∀i = 1, 2, . . . , N , ∀t ≥ 0.

Then, bounded synchronization is achieved by the network.
In particular, if we define the matrix L N = {l N i j } as we can state that ∃ 0 <c < ∞, > 0 such that lim t→∞ η(t) ≤ ∀c >c, wherē with W l , P l , Γ l representing the l×l upper-left block of matrices W, P, Γ , respectively, and where for a given value of c >c = min with W n−l representing the (n − l) × (n − l) lowerright block of matrix W and with the assumption that Proof See DeLellis et al. (2015).
We can thus derive the following corollary.
Proof First of all we need to write the dynamics of each node in the network as f i (t, . This is possible if we suppose that γ i =γ ∀i ∈ [1, N ] and define: Then we need to verify whether the nodes in the network are QUAD(P, W ) Affine systems. In particular, this means that we need h to be QUAD(P, W ), with P ∈ D + and W ∈ D. Therefore, if we define P = diag{P 11 , P 22 } with P 11 , P 22 > 0, W = diag{W 11 , W 22 } and h(t, z) = [0γ z 2 ] T ∀z ∈ R 2 , we have to satisfy: Choosing W 22 =γ P 22 , it possible to reduce Eq. 42 to which is true for any W 11 > 0. This means that the first hypothesis of Theorem 2 simply reduces to choosing any P ∈ D + and W = diag{W 11 ,γ P 22 } for any W 11 > 0. Since W ∈ R 2×2 is made up of 2 nonnegative elements, we have that l =l = 2. Therefore, in order to satisfy the second and the third hypotheses of Theorem 2, we need Γ to be a diagonal positive definite matrix, that is Γ ∈ D + (note that this is true since the nodes are connected through full state coupling, and hence Γ = I 2 ).
Finally, the last hypothesis left to satisfy is related to the boundedness of the terms g i at any time instant. As already shown before, we have chosen: Since the dynamics of each RvdP oscillator is bounded (Zhai et al. 2015c), we can define Therefore, from Eq. 44 we get: Besides, we have that This means that the fourth hypothesis of Theorem 2 is always satisfied in the case of RvdP oscillators, and the bound M is defined in Eq. 46.
In order to find an easier expression of the minimum value required for the coupling strength and of the upper bound for the error norm, we can take advantage of the particular form of matrices P and W : Therefore, from Eq. 37 we have that the minimum valuē c for the coupling strength that guarantees bounded synchronization of the network is given by Eq. 39, while from Eq. 38 we have that the upper bound of the error norm is given by Eq. 40 for a given c >c. So we can conclude that if c >c > 0, wherec is defined in Eq. 39, bounded synchronization is achieved.
Remark 3 Note that for increasing values of the coupling strength c, the estimated error bound decreases (see Eq. 40 and Eq. 41), hence so does the relative synchronization error bound χ from Definition 2. As shown qualitatively earlier in Sect. 4 and numerically later in Sect. 7.3, this also implies that the coordination level of the network, captured by the group synchronization index ρ g , increases.

Numerical validation
As previously shown for a connected simple undirected network of N heterogeneous RvdP oscillators, by making use of contraction theory it is possible to guarantee bounded synchronization if we suppose that the underlying topology is an unweighted complete graph (all-to-all network). On the other hand, by making use of Lyapunov theory, bounded synchronization can be achieved regardless of the topology and the weights of the interconnections, although an assumption has to be made on one of the parameters of the nodes (γ i =γ ∀i ∈ [1, N ]).
In order to be able to study the most general case, we consider a weighted random graph of N = 5 heterogeneous RvdP oscillators characterized by γ i =γ ∀i ∈ [1, N ]. The graph is generated by assuming an edge formation probability of 60%, and edge weights are randomly picked between 0 and 2 (Fig. 9). The parameters and the initial conditions of the nodes are selected as in Table 7. With this parameter choice, the variables needed to prove convergence according to Corollary 1 can be computed as p M = 2.6, v M = 0.96, Fig. 9 Underlying topology-simple connected weighted graph  Fig. 10, we show x 1 (t) for all the nodes in the network when they are connected through full state coupling with c = 1.45. We also show the corresponding error norm η(t). We observe that the nodes quickly converge toward each other with a residual error norm which is asymptotically bounded by 0.25, as opposed to a maximum upper bound when the oscillators are uncoupledη 3.89 (corresponding to a relative synchronization error bound χ 0.065).
In Fig. 11a, we show that, when considering full state coupling, bounded synchronization can actually be achieved for smaller values of the coupling strength (c = 0.085) and that it can also be achieved with the two other coupling functions presented earlier in this paper, yielding 1.7 and χ 0.44. Furthermore, in Fig. 11b we show how the upper bound , and as a consequence, the relative synchronization error bound χ can be lowered at will by increasing the coupling strength. For instance, by setting c = 0.15 in the full state coupling, the error norm bound reduces to 1.16, corresponding to χ 0.30.
For the sake of completeness, in Fig. 12a we show the trend of the group synchronization obtained in each of three cases of coupling functions considered here: it is possible to appreciate how, after an initial transient, ρ g (t) reaches a much higher value, confirming what observed in Richardson et al. (2012). In particular, the trend obtained when considering an HKB coupling resembles the most that obtained in the real experiments involving human people. Indeed, from Fig. 12b it is possible to appreciate how, during the transient before synchronization reaches a high value, ρ g (t) exhibits oscillations as observed in the rocking chairs experiments only when using an HKB coupling, while its trend is mostly exponential with the other two couplings.
Finally, in Table 8 we show the link between the coupling strength c, the relative synchronization error bound χ as defined in Eq. 12, and the mean value of the group synchronization ρ g as defined in Eq. 19, when considering full state coupling (analogous results can be obtained when considering partial state coupling or HKB coupling). As expected the interaction among all the nodes in the network, the smaller the distance among their trajectories, and the higher their coordination level.

Conclusion
We proposed a mathematical model based on a network of coupled heterogeneous RvdP oscillators to explain movement synchronization in a group of three or more people performing a joint oscillatory task. Each oscillator in the model is characterized by a different set of parameters to account for human-to-human variability. Three different coupling models, both linear and nonlinear, were investigated, and the effects of adding an external entrainment signal were analyzed. Also, we found analytical conditions for a  Note that, for all the different values of c, the group synchronization index has been averaged over only the first 150 s of the simulation, despite its total duration being T = 500 s. This has been necessary in order to appreciate the quantitative difference of the coordination level during the transient before ρ g (t) reaches a value almost equal to unity at steady-state connected simple undirected network to achieve bounded synchronization when considering full state coupling among the nodes, while we numerically verified that bounded synchronization can be achieved also when considering partial state coupling or HKB coupling. In particular, we observed that it is possible to replicate some of the synchronization features observed experimentally in the rocking chairs set-up described in Richardson et al. (2012) with all the three coupling functions proposed in this paper; the closest matching with the experimental data being achieved when connecting the nodes through a nonlinear HKB coupling function. The link between coupling strength, bounded synchronization, and coordination level of the network was also investigated, showing that, as expected from the theory of synchronization in complex networks, increasing the coupling reduces the synchronization error bound and improves group coordination.
Ongoing work is focused on exploiting the model analyzed in this paper to capture further experimental observations related to the coordination of hand movements in groups of more than two players. Also, research effort is being spent to obtain analytical conditions for synchronization of networks of RvdP oscillators with more general coupling functions.