Two dimensionless parameters and a mechanical analogue for the HKB model of motor coordination.

The widely cited Haken-Kelso-Bunz (HKB) model of motor coordination is used in an enormous range of applications. In this paper, we show analytically that the weakly damped, weakly coupled HKB model of two oscillators depends on only two dimensionless parameters; the ratio of the linear damping coefficient and the linear coupling coefficient and the ratio of the combined nonlinear damping coefficients and the combined nonlinear coupling coefficients. We illustrate our results with a mechanical analogue. We use our analytic results to predict behaviours in arbitrary parameter regimes and show how this led us to explain and extend recent numerical continuation results of the full HKB model. The key finding is that the HKB model contains a significant amount of behaviour in biologically relevant parameter regimes not yet observed in experiments or numerical simulations. This observation has implications for the development of virtual partner interaction and the human dynamic clamp, and potentially for the HKB model itself.

Introduction

Coordinated human and animal movement is achieved through complex vascular, skeletal, muscular and neural interactions. The degrees of freedom problem (first recognised by Bernstein (1967)) ask how the high dimensionality of such systems is reduced, enabling an organism to perform macroscopic functional tasks. Coordination dynamics (Fuchs and Jirsa 2007; Kelso 2009) has approached this question through the study of the simplest observable quantity that characterizes coordination: the relative phase of two oscillators in periodic motion (Kelso 1995). These oscillators may, for example, represent the fingers or limbs of one or more experimental subjects, teams in group sports or neural oscillators. Stable in-phase () synchronisation is often found to be the simplest to maintain (Buchanan and Ryu 2006; Kelso 1981), while stable anti-phase () motions are also commonly observed (Bourbousson et al. 2010). Stable phase-lagged (intermediate values of ) states are also found in a variety of situations (Collins and Stewart 1993; Duarte et al. 2012).

In 1985, Haken et al. (1985) proposed a model, which we shall call the basic HKB model, for the potential function of the relative phase of two oscillators, which has equilibria at and . Variation of the model’s single parameter induces changes in the one-dimensional dynamics derived from the potential function, in particular a loss of stability of the anti-phase mode causing an abrupt shift to the (always stable) in-phase mode. These dynamics have been used to model phase transitions, analogous to the switching of an animal’s gait, originally seen in bimanual finger experiments (Kelso 1981), but also later in interpersonal (Schmidt et al. 1990) and sensorimotor (Kelso et al. 1990) contexts. The suitability of the model for these applications (and its stochastic extension Schöner et al. 1986) has been well established through measurements of critical fluctuations and critical slowing down around the transition (Kelso et al. 1986; Scholz et al. 1987).

The natural next step was to consider the dynamics of the component oscillators. A self-sustaining so-called hybrid oscillator (consisting of Rayleigh and Van-der-Pol style damping terms) was found to fit an observed (inverse) monotonic relationship between amplitude and frequency (Kay et al. 1987). This oscillator contains four intrinsic parameters—one linear and two nonlinear damping coefficients and a natural frequency (considered to be a control parameter set by the experimenter, for example a pacing metronome). Many authors have commented on the difficulty in ascribing physical significance to these parameters (Peper et al. 2004). When coupled nonlinearly, two hybrid oscillators were shown (Haken et al. 1985) to produce the relative phase dynamics of the potential in the basic HKB model (under slowly varying amplitude and rotating wave approximations, and considering equal and constant amplitudes). The coupling contains additional parameters multiplying its linear and nonlinear terms, whose physical significance has again been difficult to determine. We call this model the full HKB model.

To make analytical progress, an approximate HKB (aHKB) model was obtained by averaging the system of coupled oscillators (Leise and Cohen 2007), under the assumption of weakly nonlinear damping and weakly nonlinear coupling.1 Such analysis can then be used to verify any numerical work and to shed light on the processes that give rise to the observed behaviour.

A bifurcation analysis of the full HKB model has been carried out (Avitabile et al. 2016). This approach considers arbitrary coupling strengths, together with no constraints on any other parameter. In this way, the full range of qualitative behaviour in the model can be explored, as well as uncovering the transitions (bifurcations) between such behaviours as parameters are varied. This is part of the general approach to studying such systems, known as nonlinear dynamics (Strogatz 2018, which has had an extensive and prolonged impact in many fields. Such an approach has already been used to create a dynamical framework for motor behaviour, using functional architectures and structured flows on manifolds (Huys et al. 2014).

As well as describing and predicting the dynamics of social coordination between two people (Kay et al. 1987), the HKB model has also been used to study real-time interaction between one human and a machine, via virtual partner interaction (VPI) (Kelso et al. 2009), where novel behaviours were found, and its extension to the human dynamic clamp (HDC) (Dumas et al. 2014). Multiple human players modelled as networks of HKB oscillators have also been studied (Alderisio et al. 2016).

Hence, for VPI and HDC to be effective, or to aid in understanding the more complex situation of many coupled oscillators, it is important that all the different types of behaviour contained within the two-oscillator HKB model are explored and catalogued and the overall dynamics understood. In particular, it is important to understand the parameter regimes where the approximate HKB model is valid and how results differ when the full HKB model is considered.

In this paper, we extend work on the approximate HKB model (Leise and Cohen 2007) and compare it to our own bifurcation analysis of the full HKB system. We find that the dynamics of the approximate HKB model are governed by just two dimensionless parameters: (17) the ratio of the linear damping to linear coupling and (18) the ratio of nonlinear damping to nonlinear coupling. This nondimensionalisation clarifies the weakly nonlinear regimes where we see the emergence of two new synchronisation behaviours; bistability of in-phase, anti-phase and phase-lagged solutions, without the need for extensive computations. We find excellent agreement between these analytical solutions and numerical computations (for moderate amplitudes), and through this process uncover extra solution branches not present in the analysis of Avitabile et al. (2016), leading us to provide, for the first time, the complete picture of the dynamics in the parameter range chosen by these authors.

The new results in this paper are:Our paper is organised as follows. In Sect. 2, we present the full HKB model (Haken et al. 1985). Then, we introduce the linear HKB model and a corresponding mechanical analogue, which does not appear to have been considered before. This section also contains a nondimensionalisation that reduces the aHKB model to a system with just two dimensionless parameters.

the discovery of normal modes (5) in the linear HKB model (4) and the mechanical analogy to which it corresponds (Fig. 2),

Fig. 2

Mechanical analogue of the linear HKB model (4). Both masses are connected to rigid surfaces by a linear spring, with spring constant , and a linear dashpot, with damping coefficient . The masses themselves are connected by another linear dashpot, with linear damping coefficient . The governing equations of this system are precisely the same as those of the linear HKB model (4), when

the nondimensionalisation (15) that reduces the aHKB model (13) to an Eq. (16) with just two dimensionless parameters,

determination of the existence and stability of steady states of this nondimensional aHKB model (16) in terms of arbitrary parameter values,

the use of these steady states (Fig. 4) to predict the dynamics in the full HKB model in arbitrary parameter ranges,

Fig. 4

Stability regions of steady states of (50). The vertical dashed lines correspond to parameter sweeps in the corresponding quadrant of (Avitabile et al. 2016, Figure 4). All curves are labelled according to (37)–(41). This Figure shows how much varied behaviour there is in the approximate aHKB model (16), and how to find it in -parameter space, for any value of a, d. We expect the full HKB model (2) to contain even more solutions as the nonlinear terms induce bifurcations. Bear in mind that each point in corresponds to many different parameter choices, thus allowing efficient design and development of VPI and HDC

an alternative nondimensionalisation more suited to comparison with numerical computations (Sect. 4.1),

our own numerical computations that provide the full picture (Sect. 4) of the dynamics in the parameter ranges used in Avitabile et al. (2016),

a mechanical analogue (Fig. 5) of the nondimensional aHKB model (16),

Fig. 5

A weakly nonlinear mechanical analogue of the HKB system. This is the same as Fig. 2, except that we have added three weakly nonlinear dampers (shown in red). The extra damping is not directly proportional to the velocity. Instead the coefficients represent an averaged effect. But since the averaging is the same for these new elements, it is only the ratio that matters. This observation is independent of the precise form of HKB model, which is why the weakly nonlinear general HKB model (88) can also be written in terms of two dimensionless parameters. Note that d and have opposite signs, since the nonlinear damping term is softening for positive parameter values, whereas the nonlinear coupling term is hardening

the reduction in a very general form of aHKB model (88) to one with just two dimensionless parameters.

In Sect. 3, we determine general existence and stability criteria of the steady states of the nondimensionalised aHKB model in terms of and that correspond to synchronisation. In Sect. 4, we use these results to explain and extend recent numerical continuation results (Avitabile et al. 2016) of the full HKB model. This section contains Fig. 4, possibly the most important figure in the paper. It can be used to predict the behaviour of the weakly coupled, weakly damping aHKB system in an arbitrary range of parameters. Much of this behaviour has not yet been seen in experiments or numerical simulations to date. In Sect. 5, we discuss our results and present our conclusions. The paper ends with a number of appendices, including Appendix E where we show that a large family of weakly damped, weakly coupled HKB oscillators can be reduced to our nondimensional aHKB system, suggesting that its dynamics are more universal than the particular form of HKB model that is the current paradigm.

The Haken–Kelso–Bunz model

In its most general form, the Haken–Kelso–Bunz (HKB) model (Haken et al. 1985) of two coupled nonlinear oscillators is given by the ordinary differential equationswhere differentiation with respect to time is denoted by a dot, are the oscillator amplitudes, frequency (the control or pacing parameter), is a nonlinear damping term and is a nonlinear coupling term.

Much work has been done in establishing the correct form of and to use when seeking to explain human-human or human-virtual player interactions. The following form covers all the different types of full HKB model (Haken et al. 1985) studied in the subsequent literature:where is a linear damping coefficient, is the Van der Pol damping coefficient, is the Rayleigh damping coefficient, and a is a linear coupling coefficient, b and c are nonlinear coupling coefficients. The dimensions of these coefficients are given byThe rationale for the form of (1) and (2) is fully explained in Haken et al. (1985), to which the interested reader is referred. Terms in these equations are needed to describe oscillations; these are given by the left hand side of both (1) and (2). Then, since “movement has a more or less stable amplitude, the equations must be nonlinear” (Haken et al. 1985, p. 351) and that experimentally this amplitude “drops with increasing ” (Haken et al. 1985, p. 352). One form of that meets these criteria is given by . The coupling term J was a subject of great discussion in Haken et al. (1985), with the form chosen “in the sense of a minimal model” (Haken et al. 1985, p. 353) that produced “the correct phase relationship between the two oscillators” (Haken et al. 1985, p. 352).

In (2), the nonlinear damping term is softening for positive parameter values , whereas the nonlinear coupling term is hardening for positive parameter values . The values of are either chosen by fitting with experimental data or selected in parameter sweeps. In Table 1, we give representative values of the other parameters selected from the literature. All authors except (Leise and Cohen 2007) set . In most papers, the pacing frequency , although values as high as have been used. Theoretical analysis is simplified under the assumption (Haken et al. 1985).

Table 1

Representative values of parameters used in the full HKB model (2) selected from the literature

References	γ	α	β	a	b
Avitabile et al. (2016)	∈[-2,8]	1	1	±0.5	±0.5
Fink et al. (2000)	1	1	1	−0.2	0.5
Haken et al. (1985)	1	0	1	−0.2	0.2
Jirsa et al. (2000)	1	1	1	−0.2	0.2
Leise and Cohen (2007)	0.5	0.38	0.001	−0.05	0.036
Słowiński et al. (2016)	0.641	12.457	0.008	∈[-15,15]	1
Varlet et al. (2012)	0.7	1	1	−3.2	3.2

All authors except (Leise and Cohen 2007) set

The linear HKB model

We begin by studying the linear HKB model, obtained by neglecting the nonlinear terms in (2) to getThis model does not seem to have been considered before, within the context of the HKB system. But it sheds important light on both the origin of in-phase and anti-phase synchronisation and the nature of the coupling.

The linear HKB model (4) consists of two normal modes, given by and , which satisfy the equationscorresponding to in-phase and anti-phase motion, respectively. The resulting instability chart in parameter space is shown in Fig. 1 (left pane), with the analytic details given in Appendix A. Our analytic results should be compared with (Avitabile et al. 2016, Fig. 5a), reproduced here in Fig. 1 (right pane) which was obtained by fixing the value of the nonlinear coupling at and performing a numerical continuation in .

Fig. 1

(Left pane) Regions of parameter space in which in-phase and anti-phase normal modes are unstable, where we find in-phase and anti-phase synchronisation. In region S, we find stable steady states only and hence no limit cycles. We find that and , both independent of the nonlinear coefficients. (Right pane) Numerical continuation results from (Avitabile et al. 2016, Figure 5a), with nonlinear coupling term . The lines , , , and in this pane are nonlinear effects. Most studies of the full HKB model (2) are carried out in the fourth quadrant of these figures (see Table 1)

These two panes show agreement in the lines and (where unstable modes become saturated limit cycles). We have shown that these lines, which are a fundamental part of the model originating in the linear HKB equations (4), do not change as the nonlinear coefficients , , b, c vary.

The lines , , , and in the right pane are all nonlinear effects, which do vary as the nonlinear coefficients vary and so are absent from the left pane. In Sect. 3 below, we will obtain analytic expressions for the weakly nonlinear versions of these lines (except which occurs at very large amplitudes only, outside the range of validity of that analysis).

Mechanical analogue of the linear HKB model

The form of (4) leads us to propose a mechanical analogue of the linear HKB model, shown in Fig. 2.

Mass is connected to a rigid surface by a linear spring, with spring constant , and a linear dashpot, with damping coefficient2. Another mass is connected in an identical way to a separate rigid surface. The two masses are themselves connected by a linear dashpot whose damping coefficient is . When , the governing equations of the system in Fig. 2 are precisely those of the linear HKB model (4).

In this analogy, the coupling term can be seen to be a form of damping. Clearly, when both , there is negative damping in the whole system and the equilibrium is impossible. Similarly, if both , we have positive damping, and all oscillations die out as energy is taken out of the system. We see both these cases in Fig. 1 (in the first and third quadrants, respectively).

It is clear why the in-phase motion is independent of the coupling coefficient a and dependent on the sign of . As move in phase, the coupling has no influence and the only damping that can affect the motion is .

On the other hand, anti-phase motion is clearly dependent on the relative size of the coefficients . As we shall see in Sect. 2.4, the ratio of these two coefficients plays a crucial role in the dynamics of the full HKB model (2).

Hopf bifurcations

We now consider how in-phase and anti-phase motions manifest themselves in the full (nonlinear) HKB model (2). For in-phase motion, we set in (2) to findUnder the assumption of weak nonlinearity, it can be shown (Haken et al. 1985; Leise and Cohen 2007) that a Hopf bifurcation occurs in (6) at , to produce equal amplitude in-phase synchronisation with amplitude given byprovided . is supercritical, degenerate, subcritical (Avitabile et al. 2016) for , respectively.

For anti-phase motion, we set in (2) to findSo, under the assumption of weak nonlinearity, by relabelling terms in (6), it can be seen that a Hopf bifurcation occurs in (8) at , to produce equal amplitude anti-phase synchronisation with amplitude given byprovided . Hence, must be supercritical, degenerate, subcritical for , respectively. Both and are shown in Fig. 1.

The nondimensional approximate HKB (aHKB) model

In this section, we present the approximate HKB (aHKB) model, as derived in Haken et al. (1985) and show that it can be represented by a set of equations involving only two dimensionless parameters.

Let us rewrite (2) in terms of a small parameter and suitably redefined damping and coupling coefficients:We look for solutions of the formwhere and represent two time scales, and . Then, we takecorresponding to a limit cycle of slowly varying amplitude r(T) and phase .

Using averaging (Leise and Cohen 2007, eq. (8)–(10)) or two-timing (Cass 2019), we obtain the approximate HKB (aHKB) model:whereis the relative phase between the two oscillators.

The aHKB model (13) depends on six parameters: . We have found a scaling that greatly simplifies (13), leading to a model with just two dimensionless parameters.

Assume, in line with experimental data, that each of is positive (see Table 1). One solution of (13) is given (Leise and Cohen 2007) by , where is given by (7), corresponding to equal amplitude in-phase synchronisation.

Let us introduce nondimensionalised amplitudes and time s as follows:Then, equations (13) becomewhere differentiation with respect to s is (still) denoted by a dot and the dimensionless parameters are given byWe see that is the ratio of the linear damping coefficient to the linear coupling coefficient a. We have already seen the importance of in the mechanical analogue (Fig. 2) of the linear HKB model (4). We can think of as being the ratio of the combined nonlinear coupling coefficient d, whereto the combined nonlinear damping coefficient , where3To the best of our knowledge, the derivation of (15)–(18) has not been reported before in the literature.4

In this section, we have considered the linear HKB model (4). We have shown the presence of normal modes (5), how their loss of stability corresponds to the generation of stable limit cycles and the central role they play in understanding the stability structure of the full HKB model (Fig. 1). We have also shown that the linear HKB model has a mechanical analogue, Fig. 2, which sheds light on the fundamental mechanisms behind the HKB model. We conclude this section with a note of the importance of dimensionless parameters in the HKB model. This significant development means (for example) that we will get the same results for as we would for for any value of . Hence, we can search parameter space far more efficiently to find dynamics relevant to the application of the HKB model under consideration.

Steady states of the nondimensional aHKB model

To understand the dynamics of (16) as dimensionless parameters vary, we consider the existence and stability of steady states of these equations. Stable steady states of the dimensionless aHKB model correspond to stable limit cycles of the full HKB model. We revisit and adapt the approach in Leise and Cohen (2007) to obtain results in an arbitrary range of parameters. We will see in the following analysis the emergence of two new synchronisation behaviours not possible in the linear HKB model; bistability of in-phase and anti-phase, and phase-lagged synchronisation.

Existence of steady states

Let us denote steady states of (16) bycorresponding to limit cycles of constant, possibly different, amplitudes separated by a constant phase difference . Since therefore , we must have from (16) eitherorIf we assume (21) holds, then . Let , corresponding to in-phase motion. From (16), we must haveWe solve these equations to find three possibilities for given by:Steady state I has in-phase equal amplitudes, exists and corresponds to in (9). Steady states have in-phase unequal amplitudes and steady state is degenerate.

Still assuming (21), but now with , this case corresponds to anti-phase motion. Again there are three possibilities for given by :Steady state A has anti-phase equal amplitudes. Steady states have anti-phase unequal amplitudes, and steady state is degenerate.

Finally, if (22) holds, there are three possibilities for corresponding to phase-lagged steady states given by:Both of have equal oscillation amplitudes, whereas the have unequal oscillation amplitudes. Steady state is degenerate.

Table 2 summarises the steady states of (16), together with those regions of parameter space in which they exist.

Table 2

Steady states of (16), corresponding to limit cycles of (2), grouped according to the relative size of

	Name	R1∗	R2∗	ϕ^∗	Existence regions
R1∗=R2∗≠0
1	I	1	1	0	Arbitrary κ,μ
2	A	1+2μ1-8κ	1+2μ1-8κ	π	κ<18,    μ>-12 or κ>18,    μ<-12
3	L^±	1	1	±cos^-11+μ2κ	0<μ<-4κ or -4κ<μ<0
R1∗≠R2∗(≠0)
4	N0±	121+3μ+4κ1+κ±121-μ1+κ	121+3μ+4κ1+κ∓121-μ1+κ	0	κ<-1,μ>1,μ+κ<0 or κ>-1,μ<1,μ+κ>0
5	Nπ±	121-μ1+κ±121+3μ+4κ1+κ	121-μ1+κ∓121+3μ+4κ1+κ	π	κ<-1,3μ+4κ+1<0,μ+κ>0 or κ>-1,3μ+4κ+1>0,μ+κ<0
6	N±π21,2	1/0	0/1	±π2	μ+κ=0
R1∗=R2∗=0
7	Z_0	0	0	0	Arbitrary κ,μ
8	Z_π	0	0	π	Arbitrary κ,μ
9	Z_ϕ	0	0	Arbitrary	μ=0

Rows 1–5 correspond to (Leise and Cohen 2007, Table 1, rows 2–6); rows 7–9 cover (Leise and Cohen 2007, Table 1, row 1). But (Leise and Cohen 2007, Table 1) has no equivalent of row 6 and (Leise and Cohen 2007, Table 1, row 7) is not a solution of (16)

A similar table was given in (Leise and Cohen (2007), Table 1). However, there are some differences5 and our nondimensionalisation was not carried out.

Stability of steady states

In this section, we consider the stability properties in the plane of the steady states of (16) found in the previous section. The regions of stability/instability we identify in this section are given in Fig. 3. As before, we adapt the approach of Leise and Cohen (2007).

Fig. 3

Stability regions (on the left) and instability regions (on the right) of steady states of (16), where we assume that . The labelled lines correspond to , , , , . Stable I and A co-exist and are stable in the fourth quadrant on the left in the purple region between the lines

First, we consider the three cases of equal amplitude synchronisation (the first three rows of Table 2) given by (25), (28) and (31), respectively.

For I, the case of equal amplitude in-phase synchronisation, the eigenvalues are given by Cass (2019)Hence, I is stable for and unstable for .

For A, the case of equal amplitude anti-phase synchronisation, the eigenvalues are given by Cass (2019)Hence, A is stable when between the lines and and unstable either when , or when between the lines and . From (17), we see that corresponds to the line in Fig. 1.

For , the case of equal amplitude phase-lagged synchronisation, the eigenvalues are given by Cass (2019)where . It is straightforward to show that is stable between the lines and when and unstable between the same lines when .

The cases , with unequal amplitudes are all unstable where they exist (see Cass 2019 for details).

We plot existence and stability regions for solutions in Fig. 3. Regions of stable solutions are shown on the left and regions of unstable solutions on the right. Several lines have been labelled, as follows:A full explanation of the importance of these lines is given in Appendix B.

We see that bistability and phase-lagging are possible only when and take different signs. Note the co-existence of stable equal amplitude in-phase I and anti-phase A synchronisation between the lines in the plane. Since we assume each of is positive, this region in the plane corresponds to the following region in the (a, b) coupling plane: Furthermore, the lines and mark the boundary of stable phase-lagged solutions. The above simple characterisation of these two important behaviours shows the power of our dimensionless analytical approach.

Steady states for arbitrary parameter values

If we drop the assumption that each of ?, a, is positive, for example, when designing an HDC (Dumas et al. 2014), then it can be shown that (16) becomeswhere . Equations (43) have been nondimensionalised usingThere are four cases to consider, according to the sign of , the linear damping coefficient, and the sign of , the combined nonlinear damping coefficient (20). Each case leads to a separate stability diagram in the plane. Results are given in Appendix C.

In this section, we have extended older work (Leise and Cohen 2007) on the steady states of the aHKB model (16), corresponding to constant amplitude limit cycles, by presenting results in terms of our dimensionless parameters (see Table 2). We have also given expressions for bifurcation curves lines , , and found numerically in recent work (Avitabile et al. 2016).

Comparison with numerical continuation methods

A detailed bifurcation analysis of (2) using numerical continuation methods. Doedel et al. (1997) was carried out recently (Avitabile et al. 2016). In this section, we compare our analytical results for the aHKB model from Sect. 3 with our own numerical continuation results using AUTO,6 for the full HKB model, in order to validate that analysis and determine its region of validity. We will then show how to use our results to predict behaviours in arbitrary parameter regimes. We extend the results in Avitabile et al. (2016) to provide the full picture of the dynamics in the parameter regime considered by these authors, by showing extra solution branches not seen in Avitabile et al. (2016).

Alternative nondimensionalisation

In Avitabile et al. (2016), the bifurcation analysis was performed with (see Table 1). But with our current nondimensionalisation (15)–(18), corresponds to infinite. So we propose an alternative nondimensionalisation that is better suited for comparison with Avitabile et al. (2016). We define new dimensionless parameters given byFrom (7), we have thatNow, we define and as follows:and redefine , we find that (43) becomeswhere we have dropped the tildes over the , and a dot now denotes differentiation with respect to .

Existence and stability of solutions

We present results for the existence and stability of solutions to (50), corresponding to limit cycles of the aHKB model, as dimensionless parameters vary. The methodology is identical to that in Sect. 3. Our aim is to use these results to predict behaviour in an arbitrary range of parameters for the full HKB model (2).

Our results are shown in Fig. 4, arranged in the four quadrants of (a, d) space.

In Avitabile et al. (2016), results are presented of numerical continuation of the full HKB model (2) for damping parameters , coupling parameters and pacing frequency7. From (45), (46), this corresponds to . Hence, in the first and second quadrants of Fig. 4, we have included a dashed vertical line at , and in the third and fourth quadrants, a dashed vertical line at . These lines correspond to parameter values used in (Avitabile et al. 2016, Fig. 4) when . We emphasise here that numerical studies of the parameter space are restricted to studying one such line at a time. Along these lines, we use the nondimensional aHKB model (50) to predict what we might find for weakly nonlinear amplitudes in numerical computations of the full HKB model (2).

Figure 4 is possibly the most important figure in this paper. It can be used to predict dynamics along any line in space for arbitrary values of the coupling parameters a and d.

Comparison between analysis and numerics

Our analytical results from Sect. 4.2 are valid for weakly nonlinear amplitudes. We should not expect these results to be valid for large amplitudes. To test this, we carried out our own numerical continuations of the full HKB model (2) and plotted these numerical results against our theoretical results from Sect. 3.1, in dimensional form, given in (93)–(97) of Appendix F. Detailed comparisons are given in Appendix D. In summary, we find that

for small amplitudes, our theoretical results are in complete agreement with our numerical results;

for intermediate amplitudes, our theoretical results are qualitatively similar to our numerical results, which in turn reveal additional solutions not seen in Avitabile et al. (2016), Fig. 4;

Discussion and conclusions

The original HKB paper Haken et al. (1985) has attracted considerable scientific attention, owing to its wide applicability in the area of human coordination (see reference in Avitabile et al. 2016).

In this paper, we have provided a fundamental understanding of the HKB model that has significant potential for the development of VPI and HDC.

Our first contribution was to drop all the nonlinear terms from the most widely used form of the HKB model (2) of two coupled oscillators to give the linear HKB model (4). A straightforward analysis yields the presence of two system normal modes; in-phase and anti-phase motion with differing stability properties, see Fig. 1.

In turn, the linear HKB model led us to a mechanical analogue (Fig. 2) where the coupling term is seen to be a form of damping between the two subsystems. This analogue leads to the notion that the ratio of linear damping to linear coupling a must be important for the dynamics of the full HKB model and that the in-phase normal mode must be independent of the coupling, whatever form takes.

To solve the full HKB model (2), we make use of approximation techniques, part of the method of multiple scales Nayfeh (2008), in which a) any limit cycle has a slowly varying amplitude and b) the linear coefficients (the rotating wave approximation). The resulting equations of the approximate HKB (aHKB) model (16) were shown to depend on just two dimensionless parameters: , the ratio of linear coupling to linear damping (17) and , the ratio of combined nonlinear coupling to combined nonlinear damping (18).

The discovery that aHKB dynamics is governed by just two parameters has far-reaching consequences. As mentioned above, we know for example that dynamics at will be the same as those at for any value of . Hence, we can model a variety of different experimental observations using the same model. We could even use results in one experiment to predict behaviour in another.

This applies in particular to phase-lagged results. It is well-known Avitabile et al. (2016) that the relative phase in many real-world applications can take values different from both (in-phase) and (anti-phase). For example, a stable relative phase of is seen in both the amble-to-walk gait of quadrupeds Collins and Stewart (1993) and unsuccessful defences in football Duarte et al. (2012). From Table 2, we see that the general value of the phase-lag is given by . So when , we know that such solutions lie on the line . We can often estimate , so we have a value of where we can find these solutions. That value then corresponds (18) to different nonlinear coefficients that we can select, based on the experiment under consideration.

The existence and stability of solutions to the aHKB model (16) were considered in Sect. 3. Some, but not all, of the existence results in that section have been derived earlier Leise and Cohen (2007). But these authors made the restrictive assumption that , and , and only worked with the dimensional form of the governing equations.

In Sect. 4, we used the dimensionless approximate HKB model (16) to make predictions about the dynamics of the full HKB model (2). The agreement was excellent for small to moderate amplitudes, as expected, even when . In addition, we discovered additional dynamic behaviour missing from Avitabile et al. (2016). This section contains Fig. 4, possibly the most important figure in the paper.

Considering dynamics in space can be thought of as weakly nonlinear vs. linear space. Hence, we propose a weakly nonlinear mechanical analogue, shown in Fig. 5.

This is the same as Fig. 2, except that we have added three weakly nonlinear dampers (shown in red). Just as the ratio is important in the linear mechanical analogue, so the ratio in the weakly nonlinear analogue is also important. Thus, the relative strength of the linear dampers is complemented by the relative strength of the weakly nonlinear dampers.

One feature that can explained using Fig. 5 is the point , , visible in Fig. 3 (right hand side) and Fig. 4. This occurs when , , and all intersect. In dimensional terms, this point corresponds to (the linear damping and coupling coefficients being equal) and (this happens when the combined nonlinear corrections to damping and coupling are equal). At this point, the weakly nonlinear damping and coupling are identical.

As mentioned in Sect. 2, much work has been done in establishing the correct form of and in (1). In their original paper Haken et al. (1985), the authors described establishing the form of the coupling term as “...the central problem, namely to derive a suitable coupling between ... and .”

We can now consider their two other choices within the framework of our mechanical analogue. The first choice (Haken et al. 1985, equation (3.11a)) was . When , this corresponds to replacing the central dashpot in Fig. 2 with a spring of stiffness a. When , this corresponds to replacing both the central dashpot with a spring of stiffness a and the central red weakly nonlinear dashpot with a weakly nonlinear spring of stiffness . This choice of was rejected because both linear and weakly nonlinear forms failed to produce the correct type of equation for . The second choice of (Haken et al. 1985, equation (3.11b)) added a time delay to the first choice, in an integral sense that became under approximation. So we can modify our analogue accordingly to include mechanical elements that integrate resistance over time. This second form of did produce the correct relationship between and but was not considered further by the authors.

In Appendix E, we show that a very general form of the aHKB equation (88) can be represented by one set of nondimensionalised equations (16) involving only dimensionless parameters , already given in (17), and a redefined , in (90). This can be expected from Fig. 5, where the precise form of the HKB equation is not important. These more general models may be of use to the VPI/HDC community in choosing models of coordination for experiments.

We end with some general observations about the form of the HKB model (1). The methods of nonlinear dynamics (including numerical continuation) are far more powerful than the weakly nonlinear techniques used to study the HKB equations to date. But the new solutions found by these methods are not just of passing interest. They are present in the model and so should be seen in experiments. We hope that experiments will soon be performed to test out these predictions.

But, absent any experimental evidence, another possibility is that the HKB model itself needs revision. The original authors themselves already suggested another acceptable form of , as we have just mentioned. Does that contain any of the solutions we see here? Then, there are the extensions to the HKB model that explicitly include time delays Banerjee and Jirsa (2006), Słowiński et al. (2016), Słowiński et al. (2020). But these seem to include even more behaviour than seen in experiments to date. Perhaps one place to start is if we ask the question posed by the linear HKB model (4) and its mechanical analogue, Fig. 2. What is the physiological meaning of the damping that connects the two oscillators?

Table 3

Analytic expressions for features seen in Figs. 3–10 and (Avitabile et al. 2016, Fig. 4), in terms of both pairs of dimensionless parameters and , and in terms of the system parameters . The expressions involving system parameters are obtained by undoing the scaling and setting

Label	(κ,μ)	(σ,ν)	System parameters of (2) with c=0
HBI	μ infinite	ν=0	γ=0
HBA	μ=-12	ν=-2	2a+γ=0
BPII	μ=1	ν=1	a-γ=0
BPAA	3μ+4κ+1=0	ν(σ+4)+3σ=0	3a(α+3βω^2)+4γb+γ(α+3βω^2)=0
BPAL	μ+4κ=0	4ν+σ=0	a(α+3βω^2)+4γb=0
BPN	μ+κ=0	ν+σ=0	a(α+3βω^2)+bγ=0
BPIL	μ=0	ν infinite	a=0
FH	(κ,μ)=(18,-12)	(σ,ν)=(8,-2)	(a,b)=(-12γ,18(α+3βω^2))