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.

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 φ = 0 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 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) ), 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: -the discovery of normal modes (5) in the linear HKB model (4) and the mechanical analogy to which it corresponds ( Fig. 2), -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, -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), -the reduction in a very general form of aHKB model (88) to one with just two dimensionless parameters.
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.
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 equations (1) where differentiation with respect to time is denoted by a dot, x 1 , x 2 are the oscillator amplitudes, frequency ω > 0 (the control or pacing parameter), h(x,ẋ) is a nonlinear damping term and J (x 1 − x 2 ,ẋ 1 −ẋ 2 ) is a nonlinear coupling term. Much work has been done in establishing the correct form of h(x,ẋ) and J (x 1 − x 2 ,ẋ 1 −ẋ 2 ) 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 by The 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 h(x,ẋ) that meets these criteria is given by h(x,ẋ) = (γ − αx 2 − βẋ 2 )ẋ. 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 h(x,ẋ) is softening for positive parameter values α, β, whereas the nonlinear coupling term J (x 1 − x 2 ,ẋ 1 −ẋ 2 ) is hardening for positive parameter values b, c. The values of γ, α, β, a, b, c 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 c = 0. In most papers, the pacing frequency ω ∈ [1, 10], although values as high as ω = 12π have been used. Theoretical analysis is simplified under the assumption γ ω (Haken et al. 1985).

The linear HKB model
We begin by studying the linear HKB model, obtained by neglecting the nonlinear terms in (2) to geẗ x 1 + ω 2 x 1 = γẋ 1 + a(ẋ 1 −ẋ 2 ), This 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 η I ≡ x 1 + x 2 and η A ≡ x 1 − x 2 , which satisfy the equationsη corresponding to in-phase and anti-phase motion, respectively. The resulting instability chart in (γ , a) 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   Table 1) of the nonlinear coupling at b = 0.5 and performing a numerical continuation in (γ , a). These two panes show agreement in the lines H B I and H B A (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 B P I I , B P I L , S N A , B P AL and B P A A 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 S N A 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 m 1 is connected to a rigid surface by a linear spring, with spring constant ω 2 , and a linear dashpot, with damp-ing coefficient 2 −γ . Another mass m 2 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 −a. When m 1 = m 2 = 1, the governing equations of the system in Fig. 2 are precisely those of the linear HKB model (4).
In this analogy, the coupling term a(ẋ 1 −ẋ 2 ) can be seen to be a form of damping. Clearly, when both a, γ > 0, there is negative damping in the whole system and the equilibrium x 1 = x 2 = 0 is impossible. Similarly, if both a, γ < 0, 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 x 1,2 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 a, γ . As we shall see Fig. 2 Mechanical analogue of the linear HKB model (4). Both masses m 1,2 are connected to rigid surfaces by a linear spring, with spring constant ω 2 , and a linear dashpot, with damping coefficient −γ . The masses themselves are connected by another linear dashpot, with linear damping coefficient −a. The governing equations of this system are precisely the same as those of the linear HKB model (4), when m 1 = m 2 = 1 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 x 1 = x 2 = x I in (2) to find Under the assumption of weak nonlinearity, it can be shown (Haken et al. 1985;Leise and Cohen 2007) that a Hopf bifurcation H B I occurs in (6) at γ = 0, to produce equal amplitude in-phase synchronisation with amplitude x I = r I given by provided γ α+3βω 2 > 0. H B I is supercritical, degenerate, subcritical (Avitabile et al. 2016) for α + 3βω 2 > 0, = 0, < 0, respectively.
For anti-phase motion, we set (2) to find So, under the assumption of weak nonlinearity, by relabelling terms in (6), it can be seen that a Hopf bifurcation H B A occurs in (8) at γ + 2a = 0, to produce equal amplitude anti-phase synchronisation with amplitude x A = r A given by provided γ +2a

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 1 and suitably redefined damping and coupling coefficients: We look for solutions of the form where τ = t and T = t represent two time scales, and x i 1,2 = O(1), (i = 0, 1). Then, we take corresponding to a limit cycle of slowly varying amplitude r (T ) and phase φ(T ).
Using averaging (Leise and Cohen 2007, eq. (8)-(10)) or two-timing (Cass 2019), we obtain the approximate HKB (aHKB) model: where is the relative phase between the two oscillators. The aHKB model (13) depends on six parameters: a, b, c, α, β, γ . 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 r 1 = r 2 = r I , φ = 0, where r I is given by (7), corresponding to equal amplitude in-phase synchronisation.
Let us introduce nondimensionalised amplitudes R 1 , R 2 and time s as follows: Then, equations (13) becomė where differentiation with respect to s is (still) denoted by a dot and the dimensionless parameters μ, κ are given by We 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, where to the combined nonlinear damping coefficient δ, where 3 To 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 a = a 0 , γ = γ 0 as we would for a = ka 0 , γ = kγ 0 for any value of k = 0. 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) by corresponding to limit cycles of constant, possibly different, amplitudes R * 1 , R * 2 separated by a constant phase difference φ * . Since thereforeφ = 0, we must have from (16) either If we assume (21) holds, then φ * = 0, π. Let φ * = 0, corresponding to in-phase motion. From (16), we must have We solve these equations to find three possibilities for (R * 1 , R * 2 , φ * ) given by: Steady state I has in-phase equal amplitudes, exists ∀μ, κ and corresponds to r I in (9). Steady states N ± 0 have in-phase unequal amplitudes and steady state Z 0 is degenerate.
Still assuming (21), but now with φ * = π , this case corresponds to anti-phase motion. Again there are three possibilities for (R * 1 , R * 2 , φ * ) given by : Steady state A has anti-phase equal amplitudes. Steady states N ± π have anti-phase unequal amplitudes, and steady state Z π is degenerate.
Finally, if (22) holds, there are three possibilities for (R * 1 , R * 2 , φ * ) corresponding to phase-lagged steady states given by: Both of L ± have equal oscillation amplitudes, whereas the N (1,2) ± π 2 have unequal oscillation amplitudes. Steady state Z φ is degenerate. Table 2 summarises the steady states of (16), together with those regions of (κ, μ) parameter space in which they exist.
A similar table was given in (Leise and Cohen (2007), Table 1). However, there are some differences 5 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).
First, we consider the three cases of equal amplitude synchronisation I , A, L ± (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 λ i , (i = 1, 2, 3) are given by Cass (2019) Hence, I is stable for μ < 0 and unstable for μ > 0. For A, the case of equal amplitude anti-phase synchronisation, the eigenvalues are given by Cass (2019)
For L ± , the case of equal amplitude phase-lagged synchronisation, the eigenvalues are given by Cass (2019) where . It is straightforward to show that L ± is stable between the lines μ = 0 and μ+4κ = 0 when κ < 0 and unstable between the same lines when κ > 0. The , with unequal amplitudes are all unstable where they exist (see Cass 2019 for details).
We plot existence and stability regions for solutions 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 μ = 0, κ = 1 8 , μ + 4κ = 0 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 μ = 0 and B P AL 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) becomeṡ where s γ = sgn(γ ), s δ = sgn(δ). Equations (43) have been nondimensionalised using There 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 B P AL , B P I I , B P A A and B P N 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). 6 AUTO is a widely used computer package that detects bifurcations and other nonlinear phenomena, by treating a set of parameterized ordinary differential equations and their initial conditions as a boundary value problem. There are currently more than a dozen similar packages available.

Alternative nondimensionalisation
In Avitabile et al. (2016), the bifurcation analysis was performed with γ ∈ [−2, 8] (see Table 1). But with our current nondimensionalisation (15)-(18), γ = 0 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 by From (7), we have that Now, we definer I ands as follows: and redefineR i = r ĩ r I , (i = 1, 2), we find that (43) becomeṡ where we have dropped the tildes over the R i , s a = sgn(a), s d = sgn(d) and a dot now denotes differentiation with respect tos.

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 γ ∈ [−2, 8], α = β = 1, coupling parameters a = ±0.5, b = ±0.5, c = 0 and pacing frequency 7 ω = 2. From (45), (46), this corresponds to ν = ±2γ, σ = ±26. Hence, in the first and second quadrants of Fig. 4, we have included a dashed vertical line at σ = 26, and in the third and fourth quadrants, a dashed vertical line at σ = −26. These lines correspond to parameter values used in (Avitabile et al. 2016, Fig. 4) when γ ∈ [−2, 8]. 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)

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.  (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 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 J (x 1 − x 2 ,ẋ 1 −ẋ 2 ) 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 vary-ing amplitude and b) the linear coefficients γ, a ω (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 a = a 0 , γ = γ 0 will be the same as those at a = ka 0 , γ = kγ 0 for any value of k = 0. 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 wellknown Avitabile et al. (2016) that the relative phase in many real-world applications can take values different from both 0 • (in-phase) and 180 • (anti-phase). For example, a stable relative phase of 90 • 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 φ * = ± cos −1 1 + μ 2κ . So when φ * = 90 • , we know that such solutions lie on the line μ + 2κ = 0. We can often estimate μ, so we have a value of κ where we can find these 90 • 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 γ + 2a > 0, β > 8c and α > 8b, 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 μ = a γ is important in the linear mechanical analogue, so the ratio κ = d δ 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 (μ, κ) = (1, −1), [(ν, σ ) = (1, −1)], visible in Fig. 3 (right hand side) and Fig. 4. This occurs when B P N , B P I I , and B P A A all intersect. In dimensional terms, this point corresponds to γ = a (the linear damping and coupling coefficients being equal) and d = −δ (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 h(x,ẋ) and J (x 1 − x 2 ,ẋ 1 − x 2 ) in (1). In their original paper Haken et al. (1985), the authors described establishing the form of the coupling term J (x 1 − x 2 ,ẋ 1 −ẋ 2 ) as "…the central problem, namely to derive a suitable coupling between …x 1 and x 2 ." 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) . When b = 0, this corresponds to replacing the central dashpot in Fig. 2 with a spring of stiffness a. When b = 0, 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 a + 3bω 2 . This choice of J (x 1 − x 2 ,ẋ 1 −ẋ 2 ) 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 becamė under approximation. So we can modify our analogue accordingly to include mechanical elements that integrate resistance over time. This second form of J (x 1 − x 2 ,ẋ 1 −ẋ 2 ) did produce the correct relationship between x 1 and x 2 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

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 d, δ represent an averaged effect. But since the averaging is the same for these new elements, it is only the ratio κ = d δ = 1 σ 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 authors themselves already suggested another acceptable form of J (x 1 − x 2 ,ẋ 1 −ẋ 2 ), 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?
Funding (information that explains whether and by whom the research was supported): Not applicable

Conflicts of interest/Competing interests (include appropriate disclosures): Not applicable
Ethics approval (include appropriate approvals or waivers): Not applicable

Consent for publication (include appropriate statements): Not applicable
Code availability (software application or custom code): Not applicable Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Normal modes
For in-phase motion x 1 = x 2 , and we have η A = 0; for anti-phase motion x 1 = −x 2 , we have η I = 0. Note the asymmetry between the dynamics of η I ,A in (5). The normal modes η I ,A are identical to the symmetric and anti-symmetric coordinates ψ ± of Fuchs and Jirsa (2000). But the connection with normal modes was not made in that paper. From (5), if γ > 0, the in-phase normal mode η I is always unstable. Then, (i) if γ +2a > 0, both normal modes η I ,A are unstable; (ii) if 2a + γ < 0, the anti-phase normal mode η A is stable. If γ < 0, the in-phase normal mode η I in (5) is always stable. Then, (iii) if γ + 2a > 0, η A is unstable-this happens because the coupling coefficient a is strong enough to overcome the damping, leading to instability; (iii) if 2a + γ < 0, then η A is stable and we would not expect to see any form of limit cycle in the full HKB model (2) in this region of parameter space. These four cases are shown in Fig. 1. From Table 1, we observe that nearly all studies of the full HKB model (2) are carried out in the fourth quadrant, where γ > 0, a < 0. The dynamics within the regions of Fig. 1 are as follows.
From (5), we have where C ± I ,A are constants and where we have assumed that γ, a ω. Since x 1 = 1 2 (η I + η A ) and x 2 = 1 2 (η I − η A ), we see that x 1,2 are both stable in S, but at least one of them is unstable elsewhere.

B Bifurcations
In Figs. 3-10, we labelled certain curves, following the convention in Avitabile et al. (2016). We can now explain their role in our approximate analysis and the relationship with the numerical solutions. Expressions for them in terms of both pairs of dimensionless parameters (κ, μ) and (σ, ν) and in terms of the system parameters of (2) with c = 0 are given in Table 3.
Curves HB I ,A , BP I I ,A A appear in all four quadrants in Fig. 4  Curve BP AL : 4ν + σ = 0 appears in Fig. 4, where it has two sets of four different roles in the second and fourth quadrants. We have seen it, labelled B P (a) AL , in one role in Fig. 8 and in another role in Fig. 10. The numerical equivalent of these roles is labelled B P ( f ) AL in the same figures. Curve BP N : ν + σ = 0 separates N ± 0 solutions from N ± π solutions. It features in Fig. 4, where it has a number of roles in the second and fourth quadrants. We do not see it in Figs. 7-10 since it occurs at very large values of γ .
Curve BP I L appears in (Avitabile et al. 2016, Figure 5(a)), as the line a = 0. It corresponds to a phase bifurcation between equal amplitude in-phase synchronisation I and equal amplitude phase-lagged synchronisation L ± .
Finally, we mention the point FH (a fold-Hopf bifurcation) that appears in (Avitabile et al. 2016, Figure 5(b)).It occurs around γ = −1, b ≈ 1.625 when a = 0.5 (Avitabile et al. 2016, p. 209). In our analysis, this corresponds to the point (κ, μ) When a = 0.5, we have γ = −1, and substituting in the relevant parameter values, we predict that FH occurs at b = 13 8 = 1.625, very close to the stated numerical value Avitabile et al. (2016). In Fig. 3, FH can be seen as the intersection of the three lines H B A , B P A A and B P AL . FH also appears in the second and fourth quadrants of Fig. 4. In our numerical computations, FH would correspond to varying parameters a, b in Fig. 8

so that the bifurcation points H B A , B P A A and B P (a)
AL all coincide.

C Different cases
If we drop the assumption that each of γ, α, β is positive, then it can be shown that the aHKB model (13) takes on a new form given by (43). There are four cases to consider, according to the signs of γ and δ (20). The results for stable solutions are shown in Fig. 6, where only the three cases of equal amplitude synchronisation I , A, L ± appear.
The plot in the first quadrant of Fig. 6 is identical to the left hand figure in Fig. 3, since we assumed γ > 0 and taking α, β > 0 ensures that δ = 1 4 (α + 3βω 2 ) > 0. In the second and third quadrants of Fig. 6, we have γ < 0. We find regions of stable equal amplitude anti-phase synchronisation A for all nonzero values of δ, as can be expected from the analysis in Sect. 2.1.
Note that there is no stable equal amplitude synchronisation I , A, L ± in the fourth quadrant of Fig. 6, where γ > 0, δ < 0. So all stable synchronisation for γ > 0 are captured in the first quadrant, corresponding to parameter values used in most experimental studies of the HKB model (see Table 1).

D Comparison
In what follows, we organise our results according to the quadrants of (Avitabile et al. 2016, Fig. 4), which correspond directly to those in our own Fig. 4. Throughout this comparison, we set α = 1, β = 1, ω = 2.  Figure 7 is a bifurcation diagram of max(x 1 ), r 1 vs. γ , computed for coupling strengths a = 0.5, b = 0.5, using AUTO Doedel et al. (1997), where r 1 is defined in (12).  Table 3.  (Avitabile et al. 2016, Fig. 4), in terms of both pairs of dimensionless parameters (κ, μ) and (σ, ν), and in terms of the system parameters γ, α, β, a, b, ω. The expressions involving system parameters are obtained by undoing the scaling and setting c = 0 Label (κ, μ) (σ, ν) System parameters of (2) with c = 0  Table 3. Although details are slightly different, both our theory and our numerics predict that stable A loses stability and gives rise to stable phase-lagged motion. AL at γ = 13 4 = 3.25 for these parameter values (see Table 3). At B P (a) AL , stable I and unstable A collide to give rise to stable I , A and unstable equal amplitude phaselagged solutions L ± (see inset for more detail) Note that I , L ± have the same amplitude at this level of approximation and so the solution branches are indistinguishable

E Generalised HKB model
We show that the nondimensional version of the aHKB equation (16) holds for quite general forms of damping and coupling, subject only to redefinition of r I (89) and κ (90). The rescaled time s (15) and parameter μ (17) remain the same.

E.1 Generalised nonlinear damping
To begin with, let us consider a general nonlinear damping term of the form for p ∈ [0, 1]. The full HKB model (2) contains two such terms: the Rayleigh and Van der Pol terms corresponding to p = 1 and p = 0, respectively. Before considering two coupled oscillators, let us demonstrate our approach by considering a single uncoupled oscillator, with one general damping term, of the form We examine the dynamics of (54), using two-timing Cass (2019). We give some details of our working, since the elim-ination of secular terms is more involved than outlined in Sect. 2.4. Taking x(t, ) = x 0 (τ, T ) + x 1 (τ, T ) + O( 2 ) as in (11), we have Hence from (55), we have corresponding to a limit cycle of slowly varying amplitude r (T ) and phase φ(T ), where is the complex amplitude. The next step is to substitute (57) into the right hand side of (56). Let f (τ ) ≡ −β|x 0 | 2−2 p |∂ τ x 0 | 2 p . Then using (57), we find f (τ ) = βω 2 p+1 r 3 | cos(ωτ + φ)| 2−2 p To find the secular terms in f (τ ) at O( ), we calculate the complex Fourier series expansion of f (τ ). The coefficient c 1 of e iωτ in this expansion is given by Setting θ = ωτ + φ and substituting (60) into (61), we find The expansion cos(θ − φ) = cos θ cos φ + sin θ sin φ then leads to two integrals on the right hand side of (62). The first integral vanishes, since the integrand is an odd function. The second integral can be simplified to become a 1 = 4βω 2 p+1 r 3 sin φ π π 2 0 cos 2−2 p (θ ) sin 2+2 p (θ ) dθ. (63) Then, using the identity π 2 0 cos 2s−1 θ sin 2z−1 θ dθ = 1 2 we have where Similarly, from (61) we find and hence where the complex amplitude A is defined in (59).
We then obtain an evolution equation for A(T ) by setting to zero the secular terms in (56) to find and hencė The amplitude of the resulting limit cycle is therefore We can compare this result with (7), since the equation for a single oscillator is identical to the equation for equal AL at γ = 13 4 = 3.25 for these parameter values (see Table 3). At B P (a) AL , unstable I and stable A collide to give rise to unstable I , A and stable equal amplitude phase-lagged synchronisation L ± (see inset for more detail). Note that I , L ± have the same amplitude at this level of approximation and so the solution branches are indistinguishable amplitude in-phase oscillations (see Sect. 2.3). For the Van der Pol oscillator, we have p = 0 and β = α. Hence, from (66), g(α, 0) = 1 π Γ 2 ( 3 2 )α = 1 4 α since Γ ( 3 2 ) = 1 2 √ π . So from (72), we have r I = γ g(α,0) = 2 γ α in agreement with (7) when β = 0.

E.2 Generalised nonlinear coupling
Now, we consider generalised nonlinear coupling of the form where q ∈ [0, 1]. The full HKB model (2) contains two such terms, one with q = 0 and the other with q = 1.
As with the general damping terms, any additional general coupling terms (76) are additive.

E.3 The general class of HKB models
We now present the most general class of HKB model, for coupled oscillators with n general nonlinear damping terms and m general nonlinear coupling terms, of the form: