Effects of time-delay in a model of intra- and inter-personal motor coordination

Motor coordination is an important feature of intra- and inter-personal interactions, and several scenarios — from finger tapping to human-computer interfaces — have been investigated experimentally. In the 1980s, Haken, Kelso and Bunz formulated a coupled nonlinear two-oscillator model, which has been shown to describe many observed aspects of coordination tasks. We present here a bifurcation study of this model, where we consider a delay in the coupling. The delay is shown to have a significant effect on the observed dynamics. In particular, we find a much larger degree of bistablility between in-phase and anti-phase oscillations in the presence of a frequency detuning.


Introduction
Many joint-action tasks demand some degree of movement coordination.Moreover, the degree of movement coordination plays an important role in inter-personal interactions; e.g., it can affect the level of affiliation between interacting people [21].In the case of (near) periodic movements the collective patterns of coordination are well captured by the properties of the relative phase, φ rel , between the individual coupled oscillating subsystems [10].In the case of two coupled oscillators the simplest coordination pattern is observed when the phase of the two oscillators coincide to give in-phase monostable coordination (where φ rel = 0).Monostable anti-phase coordination (where φ rel = π) can also occur, and an example of such monostable behaviour is observed in competitive games [3,6].In many real systems stable anti-phase coordination coexists with stable in-phase coordination [8,10,20,33].The development of the well-known Haken-Kelso-Bunz (HKB) [8] model (see (1) below) has been inspired by the in-phase and anti-phase coordination dynamics observed in bimanual coordination experiment [26].In the past 30 years the HBK model has been widely applied as a paradigm for studying dynamics of intra-and inter-personal motor coordination.It was found to be representative of a wide range of human movement experiments [4], suggesting that the dynamics observed in the HBK model are somehow fundamental [10].
The importance of perceptual-motor delays have been acknowledged in the original HKB model paper [8] and later discussed in [22].However, the possible role of time delays in the HKB model has not been analysed in a systematic way.More generally, there have been very few studies concerning models of motor coordination that incorporate time delays.A study of an excitator model with time delay has been presented in [2], and a model of relative phase dynamics with time delay can be found in [28,29].More recently, simulations of the HKB model with delays have been used to interpret interpersonal coordination patterns observed in the interaction differences between schizophrenic patients and healthy participants [31].
The full HKB model [8] with delays in the coupling function can be written as: In this system of four first-order differential equations with delays the variables x 1 (t) and x 2 (t) represent the positions and y 1 (t) and y 2 (t) the velocities of the two individual oscillators at time t, where time is measured in seconds.The parameters τ 1 and τ 2 denote time delays arising due to cognitive processes and/or physiological properties of the neuromuscular system.Hence, x 2 (t − τ 1 ), y 2 (t − τ 1 ) and x 1 (t − τ 2 ), y 1 (t − τ 2 ) represent the positions and velocities at times t − τ 1 and t − τ 2 , respectively.Further, ω ∈ R + is the frequency of the oscillations (either the natural or eigenfrequency or the external pacing) and ∆ is the detuning of the second oscillator with respect to ω.The parameters α, β, γ ∈ R in (1) govern the intrinsic dynamics of the single HKB oscillator; these include the Rayleigh oscillator term βy(t) 3 , as well as linear and Van der Pol-type nonlinear damping terms γy(t) and αx(t) 2 y(t), respectivley.Finally, the parameters a 1 , a 2 , b 1 and b 2 control the strength of coupling between the two HBK oscillators.Indeed, for τ 1 = τ 2 = 0 system (1) is exactly the four-dimensional HBK ODE.On the other hand, the HBK model with τ 1 , τ 2 = 0 is a delay differential equation (DDE); hence, it has as its phase space the infinite-dimensional space of continuous functions from the (maximal) delay interval into the four-dimensional (x 1 , x 2 , y 1 , y 2 )-space [5].
In this paper, we present a bifurcation study, by means of numerical continuation, of the HKB model with time delays that investigates further and explains simulation results in [31].Since delay differential equations can exhibit richer dynamics than ordinary differential equations, this type of study constitutes an important step towards a deeper understanding of more realistic models for motor coordination.More specifically, we find that time delay has a large effect with regard to regions of multistability between in-phase and anti-phase solutions.Our results provide a better understanding of the experiments and, importantly, allow us to make experimentally testable predictions.All computations were performed with the latest version of the continuation package DDE-Biftool v3.1 [27] under Matlab.
More specifically, we consider (1) with equal coupling strengths a = a 1 = a 2 and b = b 1 = b 2 .This is usually regarded as a model of intra-personal coordination, such as coordination of the left and right hand of the same person; therefore, we also assume an equal delay τ = τ 1 = τ 2 .Throughout, the intrinsic parameters of the two oscillators are fixed to values estimated directly from experimental data for wrist movements [9], namely to α = 12.457, β = 0.007905 and γ = 0.641; these values were also used in a relatively recent study of a virtual partner interaction [11].

Influence of the pacing frequency
We first investigate how the solutions of (1) change with increasing frequency ω of the oscillations, which is commonly interpreted as an increase in pacing frequency during an experiment.Here we assume that there is no detuning and, hence, set ∆ = 0. To evaluate the influence of the delay τ , we compare one-parameter bifurcation diagrams in ω of (1) for τ = 0 and for τ = 0.14.The latter value of the time delay was chosen because it is well within the range 70 − 150 ms that has been reported for responses to continuous stimuli [14]; it is also close to the value of 170 ms that was measured for refractory reactions [17].
The two bifurcation diagrams are shown in Fig. 1 in terms of the amplitude |x 1 | of the position of the first oscillator; note that |x 1 | = |x 2 | since the two coupled oscillators are identical.In-phase solutions are depicted as blue curves and anti-phase solutions as red curves, while their stability is indicated by the thickness of the respective curve.In Fig. 1 we set the coupling coefficients to a = −0.2 and b = 0.2 as in [8].
Figure 1(a) shows that for τ = 0 the in-phase solution is the only stable solution over the entire ω-range shown.Hence, without delay we do not find any bistability for the chosen experimentally validated values of α, β and γ. (We remark that these are very different from 0.3 0.9 Stability of the periodic solutions is gained or lost via branch points (where a single real Floquet multiplier crosses through 1; note that there is always an additional trivial Floquet multiplier 1).We remark that the stability properties of the periodic solutions of (1) depend also on the strength of coupling, but this is not explored further here.Panels (c)-(f) of Fig. 1 show phase portraits and time series of representative in-phase and anti-phase solutions for ω = 1.3 (indicated by the vertical line in the bifurcation diagrams in panels (a) and (b)).This value has been chosen because it lies within the range of frequencies at which bistability between in-phase and anti-phase solutions has been observed in experiments [9,19].Panel (c1) shows the stable in-phase periodic orbit for τ = 0 in the (x 1 , y 1 )-plane, and panel (c2) shows its positions and velocities over one period.The unstable solution that exists for the same value of ω = 1.3 is shown in the same way in panels (d1) and (d2), and note that now the two positions and velocities are indeed exactly in anti-phase.Panels (e) and (f) also show the in-phase and anti-phase solutions at ω = 1.3 but now for τ = 0.14.These solutions are quite similar to those for τ = 0, but they are both stable.

Regions of multistability in the (a, τ )plane
In order to gain insight into the effect of the time delay τ on the stability of the in-phase and anti-phase solutions we compute the twoparameter bifurcation diagram of (1) in the (a, τ )-plane of coupling coefficient and time delay.We investigate the dependence of the solutions on the coupling parameter a, where b = 1 is fixed.This choice is based on the observation that the bistability between in-phase and anti-phase solutions in the HKB model without delay depends on a for a range of values for b; see [1].In this bifurcation analysis we fix the frequency to ω = 1.3, as in Fig. 1(c)-(f), and explore a quite large range of τ .Namely, we not only consider values of τ that are within the physiologically relevant range of time delays observed in human physiology and cognition τ < 0.8 [7,13,18,30], but also larger time delays up to τ = 2.The latter range is relevant in coordination games with virtual partners [11] or coordination between two people who are communicating (with lags) over long distances.Similarly, we explore a relatively large range of the coupling coefficient a because there are no reliable estimates from experimental data [19,23,25].
Figure 2(a) shows the two-parameter bifurcation diagram of (1) in the (a, τ )-plane for ω = 1.3, with regions of stable in-phase and anti-phase solutions coloured accordingly; labels indicate the number   of stable in-phase (subscript i) and stable anti-phase (subscript a) solutions.In the white regions there are neither stable in-phase nor stable anti-phase solutions.The different regions are bounded by curves of Hopf bifurcations, branch points and torus bifurcations.It is important to keep in mind that we show only those parts of the respective bifurcation curves that bound regions of stable solutions.Moreover, we emphasise that we are only considering here solutions that appear via Hopf bifurcations of the trivial steady state (0,0,0,0).Additional periodic solutions (for example, those associated with symmetry-breaking or locking on tori) may exist in the the blue, red and white regions of the (a, τ )-plane; however, these are beyond the scope of this study.In Fig. 2(a) there are Hopf bifurcations only for a < 0, while we find that periodic solutions lose or gain stability at branch points and torus bifurcations.The curves of Hopf bifurcations associated with the in-phase and anti-phase solutions intersect at double Hopf points, which give rise to two branches of torus bifurcations.There are many parameter values where regions of stable in-phase and anti-phase solutions overlap, giving rise to bi-or multistability.Figure 2(a) clearly shows that, as τ increases, regions of stable in-phase and anti-phase solutions alternate.This can be understood intuitively by considering the relationship between the time delay τ and the frequency ω: stable in-phase or anti-phase can be found separated by a time shift equal to half of the period of the oscillations; for ω = 1.3 this corresponds to a distance in τ of 0.3846.
To illustrate the multistability further, Fig. 2(b) presents the oneparameter bifurcation diagram of (1) for increasing values of a and fixed τ = 0.954, which is a value (grey horizontal line in panel (a)) for which there is a complicated structure of different periodic solutions.There are four branches of periodic solutions that appear via Hopf bifurcations of the steady state: two of them extend to the left and two to the right.We further found four isolated branches: three of in-phase and one of anti-phase periodic solutions.The three branches of inphase solutions are actually connected by fold bifurcations at different values of τ with the branch of the in-phase solutions that emanates from a Hopf bifurcation and extends to the right; the same applies to the separate branch of the anti-phase solutions.Taken together, panels (a) and (b) of Fig. 2 demonstrate that there is a large degree of multistability, which increases with τ and is much more pronounced for a positive coupling strength a.
Figure 2(c)-(f) are examples of the different types of stable periodic solutions that can be found along the solution branches in Fig. 2(b); they are shown again as periodic orbits in the (x 1 , y 1 )-plane and time series of positions and velocities over one period.The periodic solutions in panels (c) and (d), along the solution branches emanating from Hopf bifurcations, are very similar to those from Fig. 1(e) and (f).The periodic solutions in Fig. 2(e) and (f), on the other hand, show an interesting new feature, namely an increasing number of smaller maxima, which correspond to little loops in projection onto the (x 1 , y 1 )-plane.We remark this phenomenon is not related to changes in the stability of the solutions; it is similar to spiking observed in model of coupled neural populations presented in [16].A detailed analysis of the emergence of maxima in ( 1) is beyond the scope of this paper.

Frequency detuning and relative phase
We now consider the case that there is a difference between the frequencies of the two oscillators, which is expressed in (1) as the detuning ∆ of the second oscillator with respect to the frequency ω of the first oscillator.A detuning between the oscillators is one of the main parameters that can be controlled in experiments.For example, it is a common experimental practice in studies of inter-personal coordination to detune the intrinsic frequencies between the pendula driven by the wrist movement by |∆| < 0.4 Hz [24,31,32].Furthermore, even for the case of intra-personal coordination of wrists movements (or if the pendula in the experimental set-up are identical) one should expect a small detuning due to small differences between the left and the right hand (or due to individual differences between participants) [9,24].As a result of detuning (∆ = 0) the relative phase φ rel between the two oscillators may take values other than 0 and π.In other words, When ∆ is increased from 0, the in-phase and anti-phase solutions generalise to solutions with a phase difference φ rel near 0 and near π, respectively; these intermediate-phase solutions are still referred to as (generalised) in-phase and anti-phase solutions for simplicity, provided the change in φ rel is reasonably small.To determine the phase of an oscillator we consider the Hilbert transform of a periodic orbit as a proto-phase, which we then transform into an observable-independent phase φ that grows linearly in time; see [12] for more details of this technique.
Figure 3 illustrates the two-parameter bifurcation diagram of (1) in the (a, ∆)-plane for τ = 0 and for τ = 0.14; also shown are examples of periodic solutions for τ = 0.14.In the bifurcation diagrams in panel (a) and (b) the central case ∆ = 0 signifies that the two oscillators have the same frequency, as was the case in previous sections.The blue regions indicate stable (generalised) in-phase periodic solutions that originate from the in-phase oscillation for ∆ = 0 with φ rel = 0, while the red regions indicate stable (generalised) anti-phase periodic solutions that originate from the anti-phase oscillation for ∆ = 0 with φ rel = π.The curves are contours of equal relative phase φ rel , which is given and labelled in term of the phase angle.As before, in the white regions of Fig. 3(a) and (b) there are no stable in-phase nor anti-phase solutions that appear via Hopf bifurcations of the trivial steady state (0,0,0,0).Figure 3(a) for τ = 0 shows a region of stable periodic solutions with phases near φ rel = 0 • , which exists for a < 0 and features φ rel up to about ±40 • , and a region of stable periodic solutions with phases near φ rel = 180 • , which exists for a > 0 and features φ rel in the range 180 • ± 50 • .Notice that the contours of φ rel are very regular and practically symmetric with respect to ∆ = 0.The blue region is bounded mostly by loci of Hopf bifurcations and the red region is bounded entirely by  Figure 3(b) for τ = 0.14 shows that the introduction of delay has a significant effect on the two typs of stability regions in the (a, ∆)-plane.First of all, the blue region of stable in-phase periodic solutions for a < 0 almost disappears.On the other hand, there is now a large region of stable in-phase periodic solutions for a > 0, which overlaps to a large extend with the red region of stable anti-phase periodic solutions.This results in a considerable and experimentally significant region of bistability between in-phase and anti-phase coordination regimes.Again, these regions are bounded by loci of branch points as well as Hopf and torus bifurcations.Notice further that the stability regions for a > 0 and the associated contours of relative phase are no longer symmetric with respect to ∆ = 0.In particular, the variation in φ rel of the stable anti-phase solutions near ∆ = 0 is larger than that of the in-phase solutions.These results are consistent with experimental findings demonstrating that intra-as well as inter-personal anti-phase coordination is less stable than the in-phase coordination [9,24,31,32].
Finally, panels (c)-(f) of Fig. 3 are examples of periodic solutions for τ = 0.14 that can be found in panel (b) for a = 9.1685 (indicated by the grey vertical line through the region of bistability).For a quite small frequency detuning of ∆ = −0.25 one finds the stable in-phase periodic solution with φ rel = 5.4 • and the stable anti-phase periodic solution with φ rel = 205.2• shown in Fig. 3(c) and (d), respectively.For a larger detuning of ∆ = −1.14one finds the quite similar stable periodic solutions with φ rel = 17.9 • in (e) and with φ rel = 240.3• in (f).This demonstrates that for sufficiently strong coupling it is possible to achieve coordination even in the face of a large frequency detuning between the two oscillators.Notice from the projections in (c1)-(f1) that the traces of the two oscillators in the (x 1 , y 1 )-plane and the (x 2 , y 2 )-plane are no longer identical, and compare with the time series in (c1)-(f1).Moreover, the stable periodic solutions originating from the anti-phase solutions for ∆ = 0 exhibit behaviour similar to the in-phase solutions for large a in Fig. 2(e)-(f).Interestingly, the additional maxima are much clearer in the time series of the velocities (black curves) than in the time series of positions.This means that they might be difficult to detect in experimental (and, hence, noisy) position data.

Conclusions
We presented a bifurcation study of the HKB model for physically relevant parameter settings that are consistent with experimental data, and with delay in the coupling.The focus was on the effect of the delay on stable in-phase and anti-phase periodic solutions, where we considered also the case when a frequency detuning is present.More specifically, the comparison of two-parameters bifurcation diagrams, without and with delay, clearly demonstrates that time delays may have a significant effect, especially on the size and features of regions of bistability between in-phase and anti-phase oscillations.Furthermore, our analysis of the HKB model with time delay and frequency detuning provides a possible dynamical explanation why the anti-phase solutions are less stable in practice.Namely, already small differences between the intrinsic frequencies of the two oscillators result in a much larger deviations from the desired coordination pattern for the anti-phase compared to the in-phase solutions.
Our study is merely a first attempt at a bifurcation analysis of the full four-dimensional HKB model with time delays.Clearly, there are many directions for future research.Indeed, a more comprehensive bifurcation anlysis with regard to the different parameters of the system is a next step and subject of our ongoing work.Given the experimental relevance of the time delays involved in human movement coordination dynamics, it would be of particular importance to study further the effects of heterogeneity in the intrinsic properties and coupling strengths of the two oscillators in the presence of delays.Another interesting future direction to explore would be the analysis of a two-tier model (different type of coupling) [2], which has some neurological support; see e.g.[15].

Figure 2 :
Figure 2: Two-parameter bifurcation diagram of (1) in the (a, τ )-plane for ω = 1.3 (a), associated one-parameter bifurcation diagram in a for τ = 0.95, and periodic solutions for τ = 0.95 and a = 5 (c)-(d), a = 8 (e) and a = 11.5 (f).Blue regions in panel (a) indicate stable in-phase solutions and red regions stable anti-phase solutions; multistability is indicated by different shades of red, blue and purple.Thick black curves are loci of Hopf bifurcations, thin black curves are loci of branch points and dashed black curves are loci of torus bifurcations; in panel (b) these bifurcations are indicated by dots (•), crosses (x) and asterisks (*), respectively.Here α = 12.457, β = 0.007095, γ = 0.641, b = 1 and ∆ = 0.