The uncoupling limit of identical Hopf bifurcations with an application to perceptual bistability

Abstract

We study the dynamics arising when two identical oscillators are coupled near a Hopf bifurcation where we assume a parameter ϵ uncouples the system at $\epsilon =0$. Using a normal form for $N=2$ identical systems undergoing Hopf bifurcation, we explore the dynamical properties. Matching the normal form coefficients to a coupled Wilson–Cowan oscillator network gives an understanding of different types of behaviour that arise in a model of perceptual bistability. Notably, we find bistability between in-phase and anti-phase solutions that demonstrates the feasibility for synchronisation to act as the mechanism by which periodic inputs can be segregated (rather than via strong inhibitory coupling, as in the existing models). Using numerical continuation we confirm our theoretical analysis for small coupling strength and explore the bifurcation diagrams for large coupling strength, where the normal form approximation breaks down.

1 Introduction

The Hopf bifurcation is a generic and well-characterised transition that a nonlinear system can undergo to create temporal patterns of behaviour on changing a parameter. At such a bifurcation, an equilibrium of an autonomous smooth dynamical system develops oscillatory instability and emits a small amplitude periodic orbit that, when followed, may be used to understand a wide variety of oscillatory phenomena. This includes many problems that appear in neuroscience applications [7].

For larger network systems composed of similar subsystems that undergo oscillatory instability, when coupled together, this can lead to the formation of non-trivial spatio-temporal patterns. Notably there is a large body of literature on coupled oscillators, viewed from a wide variety of theoretical view points and from the point of view of applications, e.g. [34]. Much of this theory either considers very specific models or makes an assumption of weak coupling which allows a reduction to a phase oscillator description such as that of Kuramoto [1], suitable for answering a lot of questions about synchronisation of system oscillations.

In this paper we consider identical subsystems undergoing a Hopf bifurcation that have an uncoupling limit. This approach gives a natural setting of two parameters that allows a thorough and generic analysis of the low-dimensional dynamics of coupled oscillator systems, by means of normal form theory. We use this analysis to understand the behaviour of a pair of Wilson–Cowan oscillators that arise in a model of perceptual bistability, which complements the results in [11].

The phenomenon of perceptual bistability motivates this study of oscillatory dynamics in a coupled dynamical system. For certain static but ambiguous sensory stimuli, two distinct perceptual interpretations (percepts) are possible, but only one can be held at a time. Not only can the initial percept be different from one short presentation of the stimulus to the next, but for extended presentations, the percept can switch dynamically. Perceptual bistability has been investigated in a number of different visual paradigms, e.g. ambiguous figures [33, 41], binocular rivalry [9, 10, 27], random-dot rotating spheres [47], motion plaids [19] and multistable barber-pole illusion [29]. Such ambiguous stimuli provide an opportunity to gain insights about the computations underlying perceptual competition in the brain. Whilst synchrony of oscillatory activity is known to play a role in the encoding of perceptually ambiguous stimuli [13], this mechanism has been widely overlooked.

Further background and motivation for the study of coupled oscillatory instabilities close to the uncoupling limit are given in Sect. 1.1, whilst further background and motivation for the study of oscillatory dynamics in the context of perceptual competition are given in Sect. 1.2 (not required reading if primarily interested in this paper’s mathematical results).

1.1 Coupled oscillatory instabilities

As noted by several authors, networks of oscillators near Hopf bifurcation allow one to explore not just the collective phase dynamics but also amplitude behaviour [16], and this allows one to use many of the tools of generic bifurcation theory with symmetry (in particular, the consequence of group actions on normal forms and the phase space) to understand the creating and properties of many oscillator patterns that may arise, biological applications including, for example, animal gaits and visual hallucination patterns [16].

A recent paper [8] explored coupled Hopf bifurcations in a two-parameter setting where one of the parameters results in uncoupling of the systems. In that setting, they found that it is possible to find not only a reduction to Kuramoto-like oscillators in a weak coupling close to threshold limit, but also to find the next order corrections that include multiple oscillator interactions. The setting also allows study of patterns where only part of the system is oscillating. More precisely, [8] considers N identical and identically interacting smooth ($C^{\infty }$) vector fields on $x_{i} \in \mathbb{R}^{d}\ (d \geq 2)$ and presents the normal form near a Hopf bifurcation.

In this paper we explore the dynamical properties of the special case $N=2$ with $d=2$. We give a dimension reduction via group-invariant coordinates in order to simplify dynamics. In the 2D normal form we look at the effects of coupling beyond the weak limit. A similar analysis was performed in [6] for the case of a linear coupling term, thus considering a particular sub-case of the normal form studied here. We then apply this theory to understand the appearance of a variety of oscillatory patterns in a model of perceptual bistability.

We emphasise that we explore a special case of two identical Hopf bifurcations that has symmetries and is close to $1:1$ resonance. The case of a double Hopf bifurcation without symmetries has been studied in [14] (see also [17]); by assuming non-resonant conditions on the Hopf bifurcation frequencies, the author provides a normal form for the bifurcation and performs a detailed study of the dynamics. Depending on the value of the coefficients, very rich dynamics can be found. However, our study examines different behaviours that are generic for systems with symmetries close to identical Hopf bifurcations, but not in the more general case.

1.2 Oscillatory models of perceptual bistability

Perceptual bistability can also arise with stimuli that change periodically. Apparent motion can be observed when a dot on a screen present at one location disappears and spontaneously reappears at a nearby location, as if travelling smoothly across the screen [3, 23]. Figure 1(A) shows two frames of such an apparent motion display,¹ where a black square to the left of a fixation point might reappear on the right of the fixation point. If two such frames alternate every, say, 200 ms as in the schematic Fig. 1(B), this can be perceived as a single square moving from side to side (“percept 1” in Fig. 1(C)). However, another interpretation is possible, of distinct squares blinking on and off either side of the fixation point (“percept 2” in Fig. 1(C)). Watching such a display, perception switches between percept 1 and percept 2 every few seconds; see [4, 36], references within and more recently [15, 32]. Perceptual bistability also occurs for the so-called auditory streaming paradigm [5, 35, 45].² The stimulus consists of interleaved sequences of tones A and B, separated by a difference in tone frequency Δf and repeating in an “ABABAB…” pattern (Fig. 1(D)). This can be perceived as one stream integrated into an alternating rhythm (“percept 1” in Fig. 1(E)) or as two segregated streams (“percept 2” in Fig. 1(E)); see recent reviews [30, 44]. There are commonalities between these visual and auditory paradigms: in percept 1 (Fig. 1(C) and (E)) the stimulus elements are linked into a single; in percept 2, the stimulus elements are separated into their distinct parts in space or in frequency. In both cases the stimulus alternates rapidly (in the range at 2–5 Hz for the visual stimulus [4]; in the range 5–10 Hz for the auditory stimulus [45]), whilst the perceptual interpretations are stable on the order of several seconds (over many cycles of the rapidly alternating stimuli).

Figure 1

Perceptually bistable stimuli that repeat periodically. (A): The visual stimulus alternates between two frames f1 and f2 with a black square flipping between a position to the left (L) or right (R) of a fixed fixation point (+). (B): Schematic of repeating stimulus, illustrating that the black square alternates between the L (f1) and R (f2) locations over time. (C): Two perceptual interpretations are possible for the same repeating stimulus: (1) apparent motion where a single square appears to travel between the L and R locations, (2) blinking squares, where the squares appear to separately blink on and off at two distinct locations. (D): The auditory stimulus features alternating pure tones at frequencies A and B separated by a difference in frequency Δf. (E): Two perceptual interpretations are possible for the same repeating stimulus: (1) an integrated percept where both the A and B tone sequences are heard in a single alternating stream “ABABAB…”, (2) a segregated percept where the A tone sequence “A_A_A_…” is heard separately from the B tone sequence “_B_B_B…” (“_” is a silent gap)

Models of perceptual bistability have successfully captured the dynamics of perceptual switching [24, 49, 50], the dependence of these dynamics on stimulus parameters [24, 31, 39, 42], mechanisms for attention [28], entrainment to slowly varying stimuli [22] and the effects of stimulus perturbations [38]. Generally models are based on competition between abstract, percept-based units [18, 28, 43, 49], but more recently models with a feature-based representation of competition have been developed [20, 24, 37, 39]. Some percept-based models have explored how rapidly alternating inputs (>2 Hz) can still give rise to stable perception over several seconds [28, 46, 50]. The models described above have considered competition directly between populations encoding different percepts, or between populations separated on a feature space. In general model studies of perceptual bistability have not explored how synchrony properties of oscillations entrained at the rate of a rapidly alternating stimulus could be the mechanism by which different perceptual interpretations emerge and coexist as bistable states (although see [48] for a large network approach to this problem). We hypothesise that oscillations play a key role in perceptual integration (such as “percept 1”) and perceptual segregation (such as “percept 2”). Towards exploring this hypothesis in future modelling studies of perceptual bistability, this paper lays the mathematical groundwork for studying the encoding perceptual states similar to those described above. An aim of the study is to identify regions of parameter space where such states coexist for a suitable neural oscillator model (but not transitions between these states).

Matching the normal form coefficients to a coupled Wilson–Cowan oscillator network allows for an understanding of the parameters in the model that govern different types of behaviour. Numerical continuation is used to confirm our theoretical analysis and to complete bifurcation diagrams for large coupling strength demonstrating where the normal form approximation breaks down. Finally, our analysis is extended with numerics to demonstrate that coexisting states akin to “percept 1” and “percept 2” persist in the presence of symmetrical periodic inputs. These coexisting states persist with low coupling strengths (down to the uncoupling limit), thus removing the need for the assumption of strong mutual inhibition between neural populations encoding different perceptual interpretations.

1.3 Outline

The structure of the paper is as follows: in Sect. 2 we use recent theoretical results in [8] to write the normal form of a system of two weakly coupled identical oscillators near a Hopf bifurcation. In Sect. 3 we perform a dynamical analysis of the system given by the dominant terms of the normal form. In particular, we study how the solutions for the uncoupled system persist for weak coupling. In Sect. 4 we identify different dynamical regimes depending on specific coefficients of the normal form and study the bifurcation diagrams. In Sect. 5 we write the equations for two mutually inhibiting Wilson–Cowan oscillators near a Hopf bifurcation, and we perform a change of coordinates to put the system in the normal form discussed in Sect. 2. For this example, we compare the theoretical predictions given by the normal form analysis with a bifurcation diagram computed numerically. Finally, we note that the results are of broad interest, extending beyond the study of neural oscillators and perceptual bistability to the study of any system involving two coupled oscillators.

2 Two identical Hopf bifurcations with an uncoupling limit

We will study systems consisting of two identically coupled oscillators of the form:

$$ \begin{gathered} \frac{d x_{1}}{dt} = H_{\lambda }(x_{1}) + \epsilon h_{\lambda , \epsilon }(x_{1}; x_{2}), \\ \frac{d x_{2}}{dt} = H_{\lambda }(x_{2}) + \epsilon h_{\lambda , \epsilon }(x_{2}; x_{1}), \end{gathered}\quad x_{1}, x_{2} \in \mathbb{R}^{2} \quad \epsilon , \lambda \in \mathbb{R} $$

(1)

having $S_{2}$ permutation symmetry. We assume that when system (1) is uncoupled ($\epsilon =0$), each system undergoes a Hopf bifurcation at the origin when the parameter λ crosses zero.

More concretely, we assume that the uncoupled system for $x \in \mathbb{R}^{2}$ given by

$$ \frac{dx}{dt} = H_{\lambda }(x) $$

has a stable focus at $x=0$ for $\lambda < 0$ that undergoes a supercritical Hopf bifurcation for $\lambda = 0$, which gives rise to a small amplitude stable limit cycle for $\lambda > 0$. For simplicity we assume that the eigenvalues of $DH_{\lambda }(0)$ are $\lambda \pm i\omega$ with $\omega \neq 0$. Moreover, without loss of generality, we assume that $(x_{1}, x_{2}) = (0,0)$ is an equilibrium point for $(\lambda , \epsilon )$ in some neighbourhood of $(0,0)$ for system (1).

2.1 Truncated normal form in complex coordinates

In [8], it is shown that systems as in (1), having $S_{2}$ symmetry and undergoing a supercritical Hopf bifurcation for $\lambda = 0$, can be written in the following normal form:

$$ \begin{gathered} \frac{dz_{1}}{dt} = U_{\lambda }(z_{1}) + \epsilon F_{N}(z_{1}, z_{2}, \epsilon ) + \mathcal{O}_{N+1}(z_{1}, z_{2}), \\ \frac{dz_{2}}{dt} = U_{\lambda }(z_{2}) + \epsilon F_{N}(z_{2}, z_{1}, \epsilon ) + \mathcal{O}_{N+1}(z_{2}, z_{1}), \end{gathered}\quad z_{1}, z_{2} \in \mathbb{C}, $$

(2)

where $F_{N}$ is an N-degree polynomial function that is equivariant under the rotational symmetries

$$ F_{N} \bigl(z_{1} e^{i\phi }, z_{2} e^{i\phi }, \epsilon \bigr) = e^{i \phi } F_{N}(z_{1}, z_{2}, \epsilon ). $$

If we consider the normal form up to order three and ignore the $\mathcal{O}_{4}(z)$ terms, we obtain the truncated normal form

$$\begin{aligned} \begin{aligned} \frac{dz_{1}}{dt}= {}& z_{1} \bigl(\lambda + i\omega + \alpha _{01} \vert z_{1} \vert ^{2} \bigr) + \epsilon \bigl[ z_{1} \bigl( \alpha _{\epsilon 0} + \alpha _{\epsilon 1} \vert z_{1} \vert ^{2} + \alpha _{\epsilon 2} \vert z_{2} \vert ^{2} + \alpha _{\epsilon 3} \bar{z}_{2}z_{1} \bigr) \\ &{}+ z_{2} \bigl(\beta _{\epsilon 0} + \beta _{\epsilon 1} \vert z_{1} \vert ^{2} + \beta _{\epsilon 2} \vert z_{2} \vert ^{2} + \beta _{\epsilon 3} \bar{z}_{1}z_{2} \bigr) \bigr] , \\ \frac{dz_{2}}{dt}={}& z_{2} \bigl(\lambda + i\omega + \alpha _{01} \vert z _{2} \vert ^{2} \bigr) + \epsilon \bigl[ z_{2} \bigl(\alpha _{\epsilon 0} + \alpha _{\epsilon 1} \vert z_{2} \vert ^{2} + \alpha _{\epsilon 2} \vert z_{1} \vert ^{2} + \alpha _{\epsilon 3}\bar{z}_{1}z_{2} \bigr) \\ &{}+ z_{1} \bigl(\beta _{\epsilon 0} + \beta _{\epsilon 1} \vert z_{2} \vert ^{2} + \beta _{\epsilon 2} \vert z_{1} \vert ^{2} + \beta _{\epsilon 3} \bar{z}_{2}z_{1} \bigr) \bigr] , \end{aligned} \end{aligned}$$

(3)

where the constants $\alpha _{01}, \alpha _{\epsilon i}, \beta _{\epsilon i} \in \mathbb{C}$ with the restriction $\operatorname{Re}(\alpha _{01}) < 0$ because the Hopf bifurcation is supercritical.

3 Dynamical analysis of the truncated normal form

3.1 Hopf bifurcations of the origin

It is straightforward to check that the origin

$$ \mathcal{S}_{0} = \{ z_{1} = z_{2} = 0 \} $$

(4)

is a fixed point of the normal form (3). Let us start by analysing its stability. The Jacobian matrix of system (3) evaluated at the origin is

$$ \begin{pmatrix} \lambda + i\omega + \epsilon \alpha _{\epsilon 0} & 0 & \epsilon \beta _{\epsilon 0} & 0 \\ 0 & \lambda - i\omega + \epsilon \bar{\alpha }_{\epsilon 0} & 0 & \epsilon \bar{\beta }_{\epsilon 0} \\ \epsilon \beta _{\epsilon 0} & 0 & \lambda + i\omega + \epsilon \alpha _{\epsilon 0} & 0 \\ 0 & \epsilon \bar{\beta }_{\epsilon 0} & 0 & \lambda - i\omega + \epsilon \bar{\alpha }_{\epsilon 0} \end{pmatrix}, $$

(5)

and their eigenvalues are given by

$$ \begin{aligned} \mu _{+} = \lambda + i \omega + \epsilon (\alpha _{\epsilon 0} + \beta _{\epsilon 0}),\qquad \mu _{-} = \lambda + i\omega + \epsilon (\alpha _{\epsilon 0} - \beta _{\epsilon 0}), \end{aligned} $$

(6)

and its complex conjugate pairs ($\bar{\mu }_{+}, \bar{\mu }_{-}$).

Clearly, when $\epsilon = 0$, the origin undergoes a double Hopf bifurcation at $\lambda = 0$. More interestingly, for $\epsilon \neq 0$, the origin undergoes two independent Hopf bifurcations, given by $\operatorname{Re} (\mu _{+})=0$ and $\operatorname{Re} (\mu _{-})=0$. These conditions define the following Hopf bifurcation curves $\mathcal{C}^{\pm }_{\mathrm{HB}}$ in the (λ, ϵ)-parameter space:

$$\begin{aligned} \begin{aligned} &C^{+}_{\mathrm{HB}}= \bigl\{ \operatorname{Re}(\mu _{+}) = 0 \quad \text{or equivalently} \quad \bar{ \alpha }^{+}:= \lambda + \epsilon (\alpha _{\epsilon 0R} + \beta _{\epsilon 0R}) = 0 \bigr\} , \\ &C^{-}_{\mathrm{HB}}= \bigl\{ \operatorname{Re}(\mu _{-}) = 0 \quad \text{or equivalently} \quad \bar{\alpha }^{-}:= \lambda + \epsilon (\alpha _{\epsilon 0R} - \beta _{\epsilon 0R}) = 0 \bigr\} . \end{aligned} \end{aligned}$$

(7)

At each curve $C^{\pm }_{\mathrm{HB}}$, a limit cycle, that will be denoted by $\mathcal{S}^{\pm }_{\mathrm{osc}}$, is born.

To study the stability of the origin of system (3), we analyse the sign of the real part of its eigenvalues $\mu ^{+}$ and $\mu ^{-}$ given in (6) at the Hopf bifurcation curves $C^{\pm }_{\mathrm{HB}}$ defined in (7). Thus,

$$ \begin{aligned} &\text{if}\quad (\lambda , \epsilon ) \in C^{+}_{\mathrm{HB}}\quad \rightarrow \quad\operatorname{Re}(\mu _{+}) = 0,\qquad \operatorname{Re}(\mu _{-}) = -2 \epsilon \beta _{\epsilon 0R}, \\ &\text{if}\quad (\lambda , \epsilon ) \in C^{-}_{\mathrm{HB}}\quad \rightarrow\quad \operatorname{Re}(\mu _{+}) = 2 \epsilon \beta _{\epsilon 0R},\qquad \operatorname{Re}(\mu _{-}) = 0. \end{aligned} $$

(8)

Therefore, we conclude that (see Fig. 2):

• If $\beta _{\epsilon 0R} > 0$, for $(\lambda , \epsilon ) \in C^{+} _{\mathrm{HB}}$, the solution $\mathcal{S}_{0}$ changes from a stable focus to a saddle-focus and a stable limit cycle $\mathcal{S}_{ \mathrm{osc}}^{+}$ emerges from $C^{+}_{\mathrm{HB}}$. Moreover, when $(\lambda , \epsilon ) \in C^{-}_{\mathrm{HB}}$, the solution $\mathcal{S}_{0}$ changes from a saddle-focus to an unstable focus and a saddle limit cycle $\mathcal{S}_{\mathrm{osc}}^{-}$ appears.

• If $\beta _{\epsilon 0R} < 0$, for $(\lambda , \epsilon ) \in C^{-} _{\mathrm{HB}}$, the solution $\mathcal{S}_{0}$ changes from a stable focus to a saddle-focus and a stable limit cycle $\mathcal{S}_{ \mathrm{osc}}^{-}$ emerges from $C^{-}_{\mathrm{HB}}$. Moreover, when $(\lambda , \epsilon ) \in C^{+}_{\mathrm{HB}}$, the solution $\mathcal{S}_{0}$ changes from a saddle-focus to an unstable focus and a saddle limit cycle $\mathcal{S}_{\mathrm{osc}}^{+}$ appears.

• If $\beta _{\epsilon 0R} = 0$, for $(\lambda , \epsilon ) \in C^{-} _{\mathrm{HB}} =C^{+}_{\mathrm{HB}}$, the solution $\mathcal{S}_{0}$ changes from a stable focus to an unstable focus and two stable limit cycles $\mathcal{S}_{\mathrm{osc}}^{+}$ and $\mathcal{S}_{ \mathrm{osc}}^{-}$ appear.

Figure 2

Sketch for the curves $C^{\pm }_{\mathrm{HB}}$ in (2). If $\beta _{\epsilon 0R} > 0$, a stable limit cycle emerges from $C^{+}_{\mathrm{HB}}$, whereas a saddle limit cycle emerges from $C^{-}_{\mathrm{HB}}$. The case $\beta _{\epsilon 0R} < 0$ is analogous just reversing ± by ∓. For the special case $\beta _{\epsilon 0R} = 0$, two stable limit cycles emerge at the coincident curves $C^{+}_{\mathrm{HB}}$ and $C^{-}_{\mathrm{HB}}$. For these plots, we assume $\beta _{\epsilon 0R} > \alpha _{\epsilon 0R} > 0$

In the next section we analyse the oscillatory solutions $\mathcal{S} _{\mathrm{osc}}^{\pm }$ that arise from the bifurcation curves $C^{\pm }_{\mathrm{HB}}$ of system (3).

3.2 Truncated normal form in polar coordinates

To perform the analysis of the oscillatory solutions $\mathcal{S}_{ \mathrm{osc}}^{\pm }$, we express the normal form in (3) in polar coordinates, that is, we write $z_{n} = r_{n} e^{i\varphi _{n}}$ with $r_{n} > 0$ and $\varphi _{n} \in \mathbb{T}$:

$$ \begin{aligned}& \dot{r}_{1}= r_{1} \bigl( \lambda + \alpha _{01R} r^{2}_{1} \bigr) + \epsilon f_{r}(r_{1}, r_{2}, \varDelta \varphi ), \\ &\dot{r}_{2}= r_{2} \bigl( \lambda + \alpha _{01R} r^{2}_{2} \bigr) + \epsilon f_{r}(r_{2}, r_{1}, -\varDelta \varphi ), \\ &r_{1} \dot{\varphi }_{1}= r_{1} \bigl( \omega + \alpha _{01I} r^{2} _{1} \bigr) + \epsilon f_{\varphi }(r_{1}, r_{2}, \varDelta \varphi ), \\ &r_{2} \dot{\varphi }_{2}= r_{2} \bigl( \omega + \alpha _{01I} r^{2} _{2} \bigr) + \epsilon f_{\varphi }(r_{2}, r_{1}, -\varDelta \varphi ), \end{aligned} $$

(9)

where $\varDelta \varphi = \varphi _{2} - \varphi _{1}$ and the subscript $X = {R, I}$ in $\alpha _{01}$ refers to its real and imaginary parts, respectively. The expression for the functions $f_{r}$ and $f_{\varphi }$ can be found in Eq. (53) in the Appendix. System (9) can be also written using the variable Δφ:

$$ \begin{aligned} &\dot{r}_{1}= r_{1} \bigl( \lambda + \alpha _{01R} r^{2}_{1} \bigr) + \epsilon f_{r}(r_{1}, r_{2}, \varDelta \varphi ), \\ &\dot{r}_{2}= r_{2} \bigl( \lambda + \alpha _{01R} r^{2}_{2} \bigr) + \epsilon f_{r}(r_{2}, r _{1}, -\varDelta \varphi ), \\ &\dot{\varDelta \varphi }= \alpha _{01I} \bigl(r ^{2}_{2} - r^{2}_{1} \bigr) + \epsilon f_{\varDelta \varphi }(r_{1}, r_{2}, \varDelta \varphi ), \\ &\dot{\varphi }_{1}= \omega + \alpha _{01I} r^{2} _{1} + \frac{\epsilon }{r_{1} } f_{\varphi }(r_{1}, r_{2}, \varDelta \varphi ), \end{aligned} $$

(10)

where the expression for the function $f_{\varDelta \varphi }$ can be found in Eq. (54) in the Appendix.

Remark 1

The general non-resonant case of the double Hopf bifurcation is discussed in [14] (see also [17]). The equations for the normal form in polar coordinates satisfy that the amplitudes $r_{1}, r_{2}$ decouple from the angles $\varphi _{1}, \varphi _{2}$. However, in our case (see system (10)) the equations for the amplitudes $r_{1}, r_{2}$ depend on $\varDelta \varphi = \varphi _{2} - \varphi _{1}$ leading to different generic dynamics than the one in [14], which we study in this paper.

Notice that the analysis of system (10) can be simplified by studying the system consisting of the first three equations, since they can be decoupled from the last one. Furthermore, we can further simplify the analysis by exploiting the $S_{2}$ permutation symmetry of the system. This symmetry acts on the phase space as

$$ K: (r_{1}, r_{2}, \varDelta \varphi ) \rightarrow (r_{2}, r_{1}, - \varDelta \varphi ) \quad \text{and}\quad K^{2} = \mathrm{Id}. $$

(11)

This action can be diagonalised using sum and difference variables $s = r_{1} + r_{2}$, $d = r_{1} - r_{2}$, with $s, d \in \mathbb{R} ^{+}\times \mathbb{R}$: in this case

$$ \tilde{K}(s, d, \varDelta \varphi ) \rightarrow (s, -d, -\varDelta \varphi ). $$

(12)

Thus, expressing the first three equations of system (10) in the variables $(s, d, \varDelta \varphi )$, we have

$$ \begin{aligned} &\dot{s}= s \biggl(\lambda + \frac{\alpha _{01R}}{4} \bigl(s^{2} + 3d^{2} \bigr) \biggr) + \epsilon g_{s}(s, d, \varDelta \varphi ), \\ &\dot{d}= d \biggl(\lambda + \frac{\alpha _{01R}}{4} \bigl(d^{2} + 3s^{2} \bigr) \biggr) + \epsilon g_{d}(s, d, \varDelta \varphi ), \\ &\dot{\varDelta \varphi }= -\alpha _{01I}sd + \epsilon g_{\varDelta \varphi }(s, d, \varDelta \varphi ), \end{aligned} $$

(13)

where the expressions for functions $g_{s}$, $g_{d}$ and $g_{\varDelta \varphi }$ are given in Eq. (55) in the Appendix.

System (13) will be the object of study for the rest of the section and will be referred to as the reduced system. As we will see in Sect. 3.2.2, working in the variables $s, d, \varDelta \varphi$ has the advantage that the linearised system about the solutions of interest becomes block diagonal.

3.2.1 Dynamical analysis of the reduced system in the uncoupled case ($\epsilon = 0$)

The general picture of the uncoupled case can be obtained straightforwardly from the original system (3) for $\epsilon =0$. Indeed, as we consider two identical systems having a supercritical Hopf bifurcation at $\lambda = 0$, the solutions of system (3) for $\lambda >0$ will correspond to the Cartesian product of solutions of each 2-dimensional system. In this section we show how the solutions for $\epsilon = 0$ are seen in the reduced system (13) so that we can explore how they evolve for $\epsilon \neq 0$. System (13) for $\epsilon = 0$ writes

$$ \begin{aligned} &\dot{s}= s \biggl(\lambda + \frac{\alpha _{01R}}{4} \bigl(s^{2} + 3d^{2} \bigr) \biggr), \\ &\dot{d}= d \biggl(\lambda + \frac{\alpha _{01R}}{4} \bigl(d^{2} + 3s^{2} \bigr) \biggr), \\ &\dot{\varDelta \varphi }= - \alpha _{01I}sd. \end{aligned} $$

(14)

Notice that in this case, the first two equations decouple from the third one and can be studied independently. As the variables $(s, d)$ are defined in $\mathbb{R}^{+}\times \mathbb{R}$, the fixed points of the first two equations of system (14) are given by

$$ (0, 0 ), \qquad \biggl(\sqrt{\frac{-4\lambda }{\alpha _{01R}}}, 0 \biggr),\qquad \biggl(+\sqrt{\frac{-\lambda }{\alpha _{01R}}},-\sqrt{ \frac{- \lambda }{\alpha _{01R}}} \biggr), \qquad \biggl(+\sqrt{\frac{-\lambda }{ \alpha _{01R}}},+ \sqrt{\frac{-\lambda }{\alpha _{01R}}} \biggr). $$

(15)

Then, as the Jacobian matrix for the two first equations of system (14) is given by

$$ \begin{pmatrix} \lambda + \frac{3\alpha _{01R}}{4}(s^{2} + d^{2}) & \frac{\alpha _{01R}}{4}6\,ds \\ \frac{\alpha _{01R}}{4}6\,ds & \lambda + \frac{3\alpha _{01R}}{4}(s^{2} + d^{2}) \end{pmatrix}, $$

(16)

it is straightforward to see that the eigenvalues of (16) for $(s, d) = (0, 0)$ are λ (double), for $(s, d) = (\sqrt{\frac{-4\lambda }{\alpha _{01R}}}, 0 )$ are $-2\lambda$ (double) and for $(s, d) = (\sqrt{\frac{- \lambda }{\alpha _{01R}}}, \pm \sqrt{\frac{-\lambda }{\alpha _{01R}}} )$ are λ and $-2 \lambda$.

Thus, as Fig. 3 shows, when $\lambda = 0$ the origin undergoes a bifurcation and changes from stable to unstable, while three new fixed points appear: one stable corresponding to $(s, d) = ( \sqrt{ \frac{-\lambda }{\alpha _{01R}}}, 0 )$ plus two unstable corresponding to $(s, d) = ( \sqrt{\frac{-\lambda }{\alpha _{01R}}}, \pm \sqrt{\frac{- \lambda }{\alpha _{01R}}} )$.

Figure 3

Bifurcation diagram of system (14) for $\epsilon =0$ as a function of λ. For the critical value $\lambda = 0$, the system undergoes a bifurcation

Now let us study the solutions of system (14) obtained from the fixed points (15) when considering the variable Δφ. The (singular) solution

$$ \bar{\mathcal{S}}_{0} = \{ s = d = 0, \quad \varDelta \varphi \in \mathbb{T} \} $$

(17)

corresponds to the origin $\mathcal{S}_{0}$ in (4) of system (3), which is a focus with eigenvalues $\lambda \pm i\omega$ (double) (see Sect. 3.1).

For any value $\varDelta \varphi _{0}$, the solution

$$ \bar{\mathcal{S}}_{1} (\varDelta \varphi _{0}) = \biggl\{ s = \sqrt{\frac{-4 \lambda }{\alpha _{01R}}},\quad d = 0,\quad \varDelta \varphi =\varDelta \varphi _{0} \biggr\} $$

(18)

is a fixed point of system (14) with eigenvalues $-2\lambda$ (double) and 0. These fixed points fill up the invariant curve

$$ \bar{\mathcal{T}}_{0} = \biggl\{ s = \sqrt{ \frac{-4\lambda }{\alpha _{01R}}},\quad d = 0,\quad \varDelta \varphi \in \mathbb{T} \biggr\} , $$

(19)

whose characteristic exponents are $-2\lambda$ (double). The fixed points $\bar{\mathcal{S}}_{1} (\varDelta \varphi _{0})$ and the invariant curve $\bar{\mathcal{T}}_{0}$ correspond in the original system (3) for $\epsilon = 0$ to the periodic orbits

$$ \begin{aligned}\mathcal{S}_{1} \bigl(\varphi _{2}^{0} \bigr) ={} &\biggl\{ z_{1}=\sqrt{ \frac{-\lambda }{\alpha _{01R}}}e^{i\varphi _{1}(t)}, \quad z_{2}=\sqrt{ \frac{- \lambda }{\alpha _{01R}}}e^{i(\varphi _{1}(t) + \varDelta \varphi _{0})}, \\ &\varDelta \varphi _{0} = \varphi ^{0}_{2} - \varphi ^{0}_{1}, \quad \varphi _{1}(t) = \varphi ^{0}_{1} + \biggl(\omega - \lambda \frac{ \alpha _{01I}}{\alpha _{01R}} \biggr)t ,\quad t \in \mathbb{R} \biggr\} \end{aligned} $$

(20)

and the 2-dimensional invariant torus $\mathcal{T}_{0}$

$$ \mathcal{T}_{0} = \bigcup _{\varphi _{2}^{0} \in \mathbb{T}} \mathcal{S}_{1} \bigl(\varphi _{2}^{0} \bigr)= \biggl\{ \vert z_{1} \vert = \vert z_{2} \vert = \sqrt{\frac{- \lambda }{\alpha _{01R}}}, \quad \varphi _{1}, \varphi _{2} \in \mathbb{T} \biggr\} , $$

(21)

respectively. Notice that the periodic orbits $\mathcal{S}_{1}(\varphi _{2}^{0})$ fill the torus $\mathcal{T}_{0}$. The characteristic exponents of $\mathcal{T}_{0}$ are the eigenvalues of the fixed point $(s, d) = (\sqrt{\frac{-4\lambda }{\alpha _{01R}}}, 0)$ of the first two equations of system (14) which are $-2\lambda$ (double).

The invariant 2-torus $\mathcal{T}_{0}$ is the product of two periodic orbits with the same period in the uncoupled case $\epsilon = 0$. Note that $\mathcal{T}_{0}$ is normally hyperbolic as each periodic orbit is linearly stable and the torus is foliated with periodic orbits; see for example [8]. We recall that roughly speaking an invariant manifold is normally hyperbolic if the dynamics in the normal directions expands or contracts at a stronger rate than the internal dynamics. In our case the normal dynamics near the torus behaves as $e^{-2\lambda t}$, whereas the internal dynamics is just a rotation. Therefore the torus $\mathcal{T}_{0}$ is normally hyperbolic.

The last two fixed points in (15) give rise to the following periodic orbits of system (14):

$$ \begin{aligned} &\bar{\mathcal{S}}^{2} = \biggl\{ s = d= \sqrt{\frac{-\lambda }{ \alpha _{01R}}}, \quad \varDelta \varphi = \varDelta \varphi _{0}- \frac{\alpha _{01I}}{ \alpha _{01R}}\lambda t, \quad t \in \mathbb{R} \biggr\} , \\ &\bar{ \mathcal{S}}^{3} = \biggl\{ s =- d= \sqrt{ \frac{-\lambda }{\alpha _{01R}}}, \quad \varDelta \varphi = \varDelta \varphi _{0}+ \frac{ \alpha _{01I}}{\alpha _{01R}}\lambda t,\quad t \in \mathbb{R} \biggr\} , \end{aligned} $$

(22)

whose characteristic exponents are λ and $-2\lambda$, so they are of saddle type. These solutions correspond to the periodic solutions

$$ \begin{aligned} &\mathcal{S}^{2} = \biggl\{ z_{1}= \sqrt{\frac{-\lambda }{\alpha _{01R}}}e^{i\varphi _{1}(t)},\quad z_{2} = 0,\quad \varphi _{1}(t) = \varphi ^{0}_{1} + \biggl(\omega - \lambda \frac{\alpha _{01I}}{\alpha _{01R}} \biggr)t, \quad t \in \mathbb{R} \biggr\} , \\ &\mathcal{S}^{3} = \biggl\{ z_{1}= 0,\quad z_{2} = \sqrt{\frac{- \lambda }{\alpha _{01R}}}e^{i\varphi _{2}(t)}, \quad \varphi_{2}(t) = \varphi ^{0}_{2} + \biggl(\omega - \lambda \frac{\alpha _{01I}}{\alpha _{01R}} \biggr)t, \quad t \in \mathbb{R} \biggr\} , \end{aligned} $$

(23)

of the original system (3) for $\epsilon = 0$ which have characteristic exponents $-2\lambda , \lambda \pm i\omega$. Therefore, they are hyperbolic periodic orbits of saddle type for $\lambda > 0$.

In conclusion, for $\epsilon = 0$, the 4D solutions $\mathcal{S}_{0}$, $\mathcal{T}_{0}$ and $\mathcal{S}^{2}$ and $\mathcal{S}^{3}$ arising from the union of solutions of each independent subsystem in (3) can be seen in the uncoupled reduced system (14) as two invariant curves filled with fixed points, $\bar{\mathcal{S}}_{0}$ and $\bar{\mathcal{T}}_{0}$, and two saddle periodic orbits, $\bar{\mathcal{S}}^{2}$ and $\bar{\mathcal{S}}^{3}$, respectively (see Fig. 4).

Figure 4

Phase space for the unperturbed system (14) for $\lambda > 0$. There are two invariant curves, $\bar{\mathcal{S}}_{0}$ (which is unstable) and $\bar{\mathcal{T}}_{0}$ (which is stable), filled with fixed points. Moreover, there exist two saddle periodic orbits $\bar{\mathcal{S}} _{2}$ and $\bar{\mathcal{S}}_{3}$

Solutions $\mathcal{S}^{2}$ and $\mathcal{S}^{3}$ are hyperbolic periodic orbits for $\lambda > 0$ and $\epsilon =0$. Therefore, for $\lambda > 0$ fixed and ϵ small enough, there exist periodic orbits $\mathcal{S}^{2}_{\epsilon }$ and $\mathcal{S}^{3}_{\epsilon }$ that are $C^{1}$-close to the unperturbed ones.

To ensure the persistence of the torus $\mathcal{T}_{0}$, we use the Fenichel theorem [12] which guarantees the persistence of normally hyperbolic invariant manifolds (with a certain degree of smoothness) for small enough perturbations.

Lemma 1

For a fixed value of $\lambda >0$, there exists $\epsilon _{0} = \epsilon _{0}(\lambda )$ such that, for any $0\le \epsilon \le \epsilon _{0}$, system (3) has a stable 2-dimensional torus $\mathcal{T}_{\epsilon }$ that is $\mathcal{C}^{1}$-close to $\mathcal{T}_{0}$.

The analytic continuation when ϵ increases of the periodic orbits $\mathcal{S}_{\epsilon }^{2}$, $\mathcal{S}_{\epsilon }^{3}$ and the invariant torus $\mathcal{T}_{\epsilon }$ provided by Lemma 1 is beyond the scope of this paper. We note that the periodic orbits $\mathcal{S}_{\epsilon }^{2}$ and $\mathcal{S}_{\epsilon }^{3}$ are limited only by hyperbolicity. Moreover, previous work [8] highlighted that continuation of the torus $\mathcal{T}_{\epsilon }$ with ϵ in Lemma 1 is only possible for $\epsilon =o(\lambda )$. Beyond this regime there will typically be loss of smoothness and breakup of the torus [2].

In Sect. 3.2.2 we are able to study the persistence, for $(\lambda ,\epsilon )$ small, of the periodic solutions $\mathcal{S} _{\mathrm{osc}}^{\pm }$ that are born at the bifurcation curves $C_{\mathrm{HB}}^{\pm }$ (see (7)). In Remark 2 we relate these periodic orbits $\mathcal{S} _{\mathrm{osc}}^{\pm }$ with the invariant torus $\mathcal{T}_{\epsilon }$ for λ fixed and ϵ small enough, i.e. where the existence of the invariant torus is guaranteed. Later, in Sect. 4 we give a detailed study of all the possible bifurcations of the solutions $\mathcal{S}_{\mathrm{osc}}^{\pm }$.

3.2.2 The oscillating solutions $\mathcal{S}_{ \mathrm{osc}}^{\pm }$ in the coupled case ($\epsilon > 0$)

We can take advantage of the $S_{2}$ symmetry of system (13) to look for solutions which remain invariant under the application of the permutation map K̃ in (12). Notice that by denoting $r_{1} = r_{2} = r ^{*}$, the curves $(r^{*}, r^{*}, 0)$ and $(r^{*}, r^{*}, \pi )$ are invariant for system (13). Then, if we write these curves in the ($s,d$) coordinates

$$ \varXi ^{+} = \bigl\{ (s, d, \varDelta \varphi ) = (s, 0, 0) \bigr\} ,\qquad \varXi ^{-} = \bigl\{ (s, d, \varDelta \varphi ) = (s, 0, \pi ) \bigr\} , $$

(24)

the dynamics for system (13) when restricted to $\varXi ^{\pm }$ reduces to

$$ \begin{aligned} &\dot{s} = \lambda s + \frac{s^{3}\alpha _{01R}}{4} + \epsilon \biggl[( \alpha _{\epsilon 0R} \pm \beta _{\epsilon 0R})s \\ &\phantom{\dot{s} =}{} +\frac{s^{3}}{4} \bigl(\overbrace{\alpha _{\epsilon 2R} + \alpha _{\epsilon 1R} + \beta _{\epsilon 3R} \pm (\beta _{\epsilon 2R} + \beta _{\epsilon 1R} + \alpha _{\epsilon 3R})}^{K_{\mathrm{stb}}^{\pm }} \bigr) \biggr], \\ &\dot{d} = 0, \\ &\dot{\varDelta \varphi } = 0, \end{aligned} $$

(25)

where the ± sign corresponds to $\varDelta \varphi = 0, \pi$, respectively.

It is straightforward to check that the equation for s in (25) has three steady solutions, namely $s=0$ (which corresponds to the solution $\bar{\mathcal{S}}_{0}$ studied before) and $s^{\pm }_{\mathrm{osc}}$ given by

$$ s^{\pm }_{\mathrm{osc}} = \sqrt{ \frac{-4 (\lambda + \epsilon (\alpha _{\epsilon 0R} \pm \beta _{\epsilon 0R}) )}{\alpha _{01R} + \epsilon K_{\mathrm{stb}}^{\pm }}}. $$

(26)

Notice that since $s \in \mathbb{R}^{+}$ we have discarded the negative solutions for the square root.

Taking into account that $\alpha _{01R}<0$, solutions $s^{\pm }_{ \mathrm{osc}}$ in (26) are only admissible when $\bar{\alpha }^{\pm }= \lambda + \epsilon (\alpha _{\epsilon 0R} \pm \beta _{\epsilon 0R}) > 0$. This restriction defines the following conditions for the bifurcation:

$$ \begin{aligned}\bar{\alpha }^{+} = \epsilon (\alpha _{\epsilon 0R} + \beta _{\epsilon 0R}) + \lambda = 0 \quad \text{for} \quad \varDelta \varphi = 0, \\ \bar{\alpha }^{-} = \epsilon (\alpha _{\epsilon 0R} - \beta _{\epsilon 0R}) + \lambda = 0 \quad \text{for}\quad \varDelta \varphi = \pi , \end{aligned} $$

(27)

which are exactly the conditions defining the curves $C^{\pm }_{ \mathrm{HB}}$ in (7) corresponding to the Hopf bifurcations of the origin.

Therefore, for ($\lambda , \epsilon$)-values on the right-hand side of curves $C^{\pm }_{\mathrm{HB}}$, we can define, respectively, the following fixed points of system (13):

$$ \begin{aligned} &\bar{\mathcal{S}}^{+}_{\mathrm{osc}} = (s, d, \varDelta \varphi ) = \bigl(s^{+}_{\mathrm{osc}}, 0, 0 \bigr), \\ &\bar{\mathcal{S}}^{-}_{\mathrm{osc}} = (s, d, \varDelta \varphi ) = \bigl(s^{-}_{\mathrm{osc}}, 0, \pi \bigr), \end{aligned} $$

(28)

which appear across a pitchfork bifurcation (whose character will be discussed below) of the origin in the s direction. Fixed points in (28) correspond to the periodic orbits $\mathcal{S} ^{\pm }_{\mathrm{osc}}$ of system (3) that appear at the Hopf bifurcation curves. Next, we will study its stability and possible bifurcations by using the reduced system (13).

The Jacobian matrix evaluated at the fixed points $\bar{\mathcal{S}} ^{\pm }_{\mathrm{osc}}$ is block diagonal

$$ \begin{pmatrix} c^{s}_{s} & 0 & 0 \\ 0 & c^{d}_{d} & c^{d}_{\varDelta \varphi } \\ 0 & c^{\varDelta \varphi }_{d} & c^{\varDelta \varphi }_{\varDelta \varphi } \end{pmatrix}, $$

(29)

where the terms $c^{s}_{s}$, $c^{d}_{d}$, $c^{d}_{\varDelta \varphi }$, $c^{\varDelta \varphi }_{d}$ and $c^{\varDelta \varphi }_{\varDelta \varphi }$ are different from zero, and their precise expressions are given in Eq. (56) in the Appendix.

Because of the block diagonal form of the Jacobian matrix, it is straightforward to check the stability in the s direction as it corresponds to the $1\times 1$ block. Thus, the eigenvalue $\bar{\mu }^{ \pm }_{1}$ takes the form

$$ \bar{\mu }^{\pm }_{1} = c_{s}^{s}=-2 \bigl(\epsilon (\alpha _{\epsilon 0R} \pm \beta _{\epsilon 0R}) + \lambda \bigr), $$

(30)

and therefore, the solutions $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$ are always stable in the s direction as they appear for $\bar{ \alpha }^{\pm }= \epsilon (\alpha _{\epsilon 0R} \pm \beta _{\epsilon 0R}) + \lambda > 0$. Therefore, the pitchfork bifurcations of the origin are supercritical (see Fig. 5).

Figure 5

Solutions $s^{\pm }_{\mathrm{osc}}$ appear through a supercritical pitchfork bifurcation of the origin in the s direction which takes place at the critical value $\bar{\alpha }^{\pm }= 0$ of the bifurcation parameter $\bar{\alpha }^{\pm }= \lambda + \epsilon ( \alpha _{\epsilon 0R} \pm \beta _{\epsilon 0R})$

As the solutions $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$ are always stable in the s direction, one has to consider the eigenvalues of the $2\times 2$ block, corresponding to the transverse directions in order to study possible bifurcations of the symmetric solutions $\bar{\mathcal{S}} ^{\pm }_{\mathrm{osc}}$. The trace (Tr^±) and the determinant (Det^±) of the $2\times 2$ block of (29) at $\bar{\mathcal{S}}^{\pm }_{ \mathrm{osc}}$ are given up to order two in $\lambda , \epsilon$ by

$$\begin{aligned} &\operatorname{Tr}^{\pm }(\lambda , \epsilon ) = c^{d}_{d} + c^{ \varDelta \varphi }_{\varDelta \varphi } = -2 \bigl(\lambda + \epsilon ( \alpha _{\epsilon 0R} \pm 3\beta _{\epsilon 0R}) \bigr), \end{aligned}$$

(31)

$$\begin{aligned} &\operatorname{Det}^{\pm }(\lambda , \epsilon ) = \pm 4\epsilon \bigl(\lambda + \epsilon (\alpha _{\epsilon 0R} \pm \beta _{\epsilon 0R}) \bigr) (C_{\mathrm{det}} + \beta _{\epsilon 0R}) + 4 \epsilon ^{2} \bigl(\beta ^{2}_{\epsilon 0I} + \beta ^{2}_{\epsilon 0R} \bigr), \end{aligned}$$

(32)

where

$$ C_{\mathrm{det}}:= \frac{\beta _{\epsilon 0I}\alpha _{01I}}{\alpha _{01R}}. $$

(33)

So, computing the discriminant

$$ \varDelta ^{\pm }= \bigl(\operatorname{Tr}^{\pm } \bigr) ^{2} - 4\operatorname{Det}^{\pm }= \bigl(\lambda + \epsilon (\alpha _{\epsilon 0R} \pm \beta _{\epsilon 0R}) \bigr) \bigl( \lambda + \epsilon (\alpha _{\epsilon 0R} \pm \beta _{\epsilon 0R}) \mp 4 \epsilon C_{\mathrm{det}} \bigr) - 4\epsilon ^{2}\beta ^{2}_{\epsilon 0I}, $$

(34)

we find that the eigenvalues of the $2\times 2$ block of the Jacobian matrix (29) write as

$$ \begin{aligned} &\bar{\mu }^{\pm }_{2} = - \bigl(\lambda + \epsilon (\alpha _{\epsilon 0R} \pm 3\beta _{\epsilon 0R}) \bigr) - \sqrt{\xi }, \\ &\bar{\mu }^{\pm }_{3}= - \bigl(\lambda + \epsilon ( \alpha _{\epsilon 0R} \pm 3\beta _{\epsilon 0R}) \bigr) + \sqrt{\xi }, \end{aligned} $$

(35)

where

$$ \xi = \bigl(\lambda + \epsilon (\alpha _{\epsilon 0R} \pm \beta _{\epsilon 0R}) \bigr) \bigl(\lambda + \epsilon (\alpha _{\epsilon 0R} \pm \beta _{\epsilon 0R}) \mp 4\epsilon C_{\mathrm{det}} \bigr) - 4\epsilon ^{2}\beta ^{2}_{\epsilon 0I}. $$

Next, we study the stability of the solutions $\bar{\mathcal{S}}^{ \pm }_{\mathrm{osc}}$ given in (28) when the parameters $\lambda ,\epsilon$ lie in the domain

$$ \mathcal{A}^{\pm }: = \bigl\{ (\lambda , \epsilon ) \in \mathbb{R}^{2}\quad | \quad \bar{\alpha }^{\pm }\geq 0,\quad \epsilon > 0 \bigr\} , $$

(36)

where $\bar{\alpha }^{\pm }$ are defined in (27). Notice that the domain $\mathcal{A}^{\pm }$ corresponds to the region on the right-hand side of curves $C^{\pm }_{\mathrm{HB}}$ and above the horizontal axis (see Fig. 2 left). Furthermore, as for the uncoupled case, we link the solutions for the reduced system (13) with the original system (3).

For $\bar{\alpha }^{\pm }= 0$, that is $(\lambda ,\epsilon ) \in C ^{\pm }_{\mathrm{HB}}$, the eigenvalues of the Jacobian matrix (29) at the fixed points $\bar{\mathcal{S}} ^{\pm }_{\mathrm{osc}}$ are given by

$$ \begin{aligned} &\bar{\mu }^{\pm }_{1} = 0, \\ &\bar{\mu }^{\pm }_{2}= \mp 2 \beta _{\epsilon 0R} - i 2\epsilon \beta _{\epsilon 0I}, \\ &\bar{\mu } ^{\pm }_{3}= \mp 2\beta _{\epsilon 0R} + i 2\epsilon \beta _{\epsilon 0I}. \end{aligned} $$

(37)

Therefore, when the parameters $(\lambda , \epsilon)$ cross the curves $C_{\mathrm{HB}}^{\pm }$ from left to right, if $\beta _{\epsilon 0R} > 0$, $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$ is a stable focus-node whereas $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ is a saddle-focus with a 1-dimensional stable manifold (corresponding to the s direction which is always stable) and vice versa if $\beta _{\epsilon 0R} < 0$. These results match exactly the results in Sect. 3.1: the 4D system has two periodic orbits that are born at different Hopf bifurcation curves $C^{+}_{\mathrm{HB}}$ and $C^{-}_{\mathrm{HB}}$ given in (7), and the stability of these periodic orbits depends on the sign of $\beta _{\epsilon 0R}$.

For ϵ small and $\bar{\alpha }^{\pm }\geq 0$, the eigenvalues of the Jacobian matrix (29) at the fixed points $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$ are given by

$$ \begin{aligned} &\bar{\mu }^{\pm }_{1} = -2\lambda + \mathcal{O}(\epsilon ) , \\ &\bar{ \mu }^{\pm }_{2} = -2\lambda + \mathcal{O}(\epsilon ) , \\ &\bar{ \mu }^{\pm }_{3} = \mp 2\epsilon (\beta _{\epsilon 0R} + C_{ \mathrm{det}}) + \mathcal{O} \bigl(\epsilon ^{2} \bigr), \end{aligned} $$

(38)

which are $\mathcal{O}(\epsilon )$-close to the ones of the uncoupled case, $-2\lambda$ (double) and 0. In particular, depending on the sign of $(\beta _{\epsilon 0R} + C_{\mathrm{det}})$, one fixed point is a stable node whereas the other is a saddle with a 1-dimensional unstable manifold.

We remark that, for $\lambda >0$ fixed and ϵ small enough, we know that there exists an invariant curve $\bar{\mathcal{T}}_{\epsilon }$ corresponding to the invariant torus $\mathcal{T}_{\epsilon }$ obtained in Lemma 1. Since this invariant curve is provided by Fenichel theory, it will contain the invariant points $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$. Consequently, if $\beta _{\epsilon 0R}+C_{\mathrm{det}}>0$, $\bar{\mathcal{T}}_{\epsilon }$ consists of the union of the saddle point $\bar{\mathcal{S}}^{-} _{\mathrm{osc}}$, its unstable 1-dimensional manifold and the stable node $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$ (and vice versa if $\beta _{\epsilon 0R}+C_{\mathrm{det}} < 0$) (see Fig. 6). In conclusion, for $\lambda >0$ fixed and ϵ small enough, the invariant torus $\mathcal{T}_{\epsilon }$ of system (3) contains the periodic orbits $\mathcal{S}^{+}_{\mathrm{osc}}$ and $\mathcal{S}^{-}_{\mathrm{osc}}$ with $\varDelta \varphi = 0$ and $\varDelta \varphi = \pi$, respectively, whose stability depends on the sign of $\beta _{\epsilon 0R} + C_{ \mathrm{det}}$.

Figure 6

Phase space of system (13) for $\beta _{\epsilon 0R} + C_{\mathrm{det}} > 0$, $\lambda > 0$, and $0 < \epsilon < \epsilon _{0}(\lambda )$ (in particular $\lambda + \epsilon (\alpha _{\epsilon 0R} \pm \beta _{\epsilon 0R}) > 0$). There exist two fixed points $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$, a stable node and a saddle point, respectively, which together with the unstable invariant manifold of the saddle point form the invariant curve $\bar{\mathcal{T}}_{\epsilon }$. Due to the coupling term, there are only two fixed points on $\bar{\mathcal{T}}_{\epsilon }$, whereas we had an infinite number in the unperturbed case. Notice that the dynamics on the s direction is always attracting

Remark 2

The existence of the invariant torus $\mathcal{T}_{\epsilon }$ is only guaranteed for $\lambda >0$ fixed and ϵ small enough by Lemma 1. The evolution and eventual breakdown of this torus $\mathcal{T}_{\epsilon }$ (or, equivalently, the invariant curve $\bar{\mathcal{T}}_{\epsilon }$) when ϵ increases is beyond the scope of this paper.

However, in Sect. 4, using system (13), we study the evolution and bifurcations of the periodic orbits $\mathcal{S}^{\pm }_{\mathrm{osc}}$ (corresponding to fixed points $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$) for $(\lambda ,\epsilon )$ small and no assumption on $\epsilon =o(\lambda )$. We leave as future work the exploration of the relationship between these bifurcations and the different mechanisms of destruction of the torus discussed in [2].

4 Bifurcation diagrams of the oscillating $S_{\mathrm{osc}}^{\pm }$ solutions

In the previous sections we have shown that when ϵ is small and $\bar{\alpha }^{\pm }\geq 0$ there exist two critical points $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$ of system (13) belonging to the curve $\bar{\mathcal{T}}_{\epsilon }$ which disappear at two independent curves $C^{\pm }_{\mathrm{HB}}$. Therefore, the points $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$ undergo several bifurcations in the domain $\mathcal{A}^{\pm }$ defined in (36). Table 1 shows the values of the trace Tr^± in (31), the determinant Det^± in (32) and the discriminant $\varDelta ^{\pm }$ in (34) of the Jacobian matrix of system (13) at $\bar{\mathcal{S}}^{ \pm }_{\mathrm{osc}}$ near the curves $C^{\pm }_{\mathrm{HB}}$ (given by the condition $\bar{\alpha }^{\pm }= 0$) and for $\bar{\alpha } ^{\pm }\geq 0$ and ϵ small. Notice that the sign of the constants $\beta _{\epsilon 0R}$ and $C_{\mathrm{det}} + \beta _{\epsilon 0R}$ is relevant to determine the local dynamics around the fixed points. In particular,

• $\beta _{\epsilon 0R}$ determines which of the two solutions $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$ can have a null trace. For $\beta _{\epsilon 0R} > 0$, it is $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$, whereas for $\beta _{\epsilon 0R} < 0$, it is $\bar{\mathcal{S}}^{-} _{\mathrm{osc}}$.

• The sign of $C_{\mathrm{det}} + \beta _{\epsilon 0R}$ determines which of the two solutions $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$ can have a null determinant. For $C_{\mathrm{det}} + \beta _{\epsilon 0R} > 0$, it is $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$, whereas for $C_{ \mathrm{det}} + \beta _{\epsilon 0R} < 0$, it is $\bar{\mathcal{S}} ^{+}_{\mathrm{osc}}$.

• Moreover, as we increase ϵ, the discriminant always changes from negative to positive. That is, consistently with the eigenvalues obtained in (37) and (38), the fixed points $\bar{\mathcal{S}} ^{\pm }_{\mathrm{osc}}$ change from a stable node and a saddle point to a stable focus and a saddle-focus.

Table 1 Values for the trace (Tr), the determinant (Det) and the discriminant (Δ) of the linearisation of system (13) at the fixed points $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$ near the curves $C^{\pm } _{\mathrm{HB}}$ ($\bar{\alpha }^{\pm }= 0$) and near to the uncoupled case ($\bar{\alpha }^{\pm }\geq 0$ and ϵ small)

Depending on the sign of $\beta _{\epsilon 0R}$ and $C_{\mathrm{det}} + \beta _{\epsilon 0R}$, we consider three different cases: (1) $\beta _{\epsilon 0R}>0$, $C_{\mathrm{det}} + \beta _{\epsilon 0R}>0$, (2) $\beta _{\epsilon 0R}<0$, $C_{\mathrm{det}} + \beta _{\epsilon 0R}>0$, and (3) $\beta _{\epsilon 0R}=0$, $C_{\mathrm{det}}>0$. The cases (i) $\beta _{\epsilon 0R}<0$, $C_{\mathrm{det}} + \beta _{\epsilon 0R}<0$, (ii) $\beta _{\epsilon 0R}>0$, $C_{\mathrm{det}} + \beta _{\epsilon 0R}<0$, and (iii) $\beta _{\epsilon 0R}=0$, $C_{ \mathrm{det}}<0$ are analogous to (1), (2) and (3), respectively, just replacing $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$ by $\bar{ \mathcal{S}}^{\mp }_{\mathrm{osc}}$. For each case, we study in detail the different bifurcations of the solutions $\bar{\mathcal{S}}^{ \pm }_{\mathrm{osc}}$ in the $(\lambda , \epsilon )$ parameter space, we link the results obtained for the 3D system (13) with the complete 4D system (3), and we discuss the regions of bistability.

4.1 Case $\beta _{\epsilon 0R} > 0$ and $C_{\mathrm{det}} + \beta _{\epsilon 0R} > 0$ (or $\beta _{\epsilon 0R} < 0$ and $C_{\mathrm{det}} + \beta _{\epsilon 0R} < 0$ )

4.1.1 Dynamics of $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$

For $\bar{\alpha }^{+} \geq 0$, λ fixed and ϵ small, the fixed point $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$ for system (13) is a stable node contained in the invariant curve $\bar{\mathcal{T}}_{\epsilon }$ (region B in Fig. 7), and as ϵ increases it becomes a stable focus at the curve $\varDelta ^{+} = 0$ (region A in Fig. 7). It disappears at a pitchfork bifurcation of the origin in the s-direction at $C^{+}_{\mathrm{HB}}$.

Figure 7

Bifurcation diagram for $\bar{\mathcal{S}}^{+}_{ \mathrm{osc}}$ in the case $\beta _{\epsilon 0R} > 0$ and $C_{ \mathrm{det}} + \beta _{\epsilon 0R} > 0$. The fixed point $\bar{ \mathcal{S}}^{+}_{\mathrm{osc}}$ appears at a supercritical pitchfork bifurcation of the origin occurring at the curve $C^{+}_{\mathrm{HB}}$

Going back to the original 4D system (3), we have that for ϵ small there exists a stable periodic orbit $\mathcal{S}^{+}_{\mathrm{osc}}$ (which belongs to the invariant torus $\mathcal{T}_{\epsilon }$), which disappears at a Hopf bifurcation of the origin in $C^{+}_{\mathrm{HB}}$.

4.1.2 Dynamics of $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$

The fixed point $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ changes from a saddle-focus with a 1-dimensional stable manifold near $C^{-}_{ \mathrm{HB}}$ to a saddle with a 2-dimensional stable manifold for ϵ small and $\bar{\alpha }^{-} > 0$. Moreover, in this case the trace for $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ vanishes. Therefore, if

$$ \beta _{\epsilon 0R} < - C_{\mathrm{det}} + \sqrt{C^{2}_{ \mathrm{det}} + \beta ^{2}_{\epsilon 0I}}, $$

(39)

then $\operatorname{Tr}^{-} = 0$ and $\varDelta ^{-} < 0$ and $\bar{ \mathcal{S}}^{-}_{\mathrm{osc}}$ undergoes a Hopf bifurcation.

So, we will distinguish two cases as follows.

(1) Case $\beta _{\epsilon 0R} < - C_{\mathrm{det}} + \sqrt{C^{2}_{\mathrm{det}} + \beta ^{2}_{\epsilon 0I}}$. For $\bar{\alpha }^{+} \geq 0$, λ fixed and ϵ small, the fixed point $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ is a saddle point with a 1-dimensional unstable manifold (in the Δφ direction) contained in the invariant curve $\bar{\mathcal{T}}_{ \epsilon }$ (region D in Fig. 8). When crossing the curve $\operatorname{Det}^{-}=0$ (region C), the point $\bar{ \mathcal{S}}^{-}_{\mathrm{osc}}$ becomes a stable node. As the coupling ϵ is increased, $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ crosses the curve $\varDelta ^{-}=0$ and $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ becomes a stable focus (region B). When the parameters cross the curve $\operatorname{Tr}^{-} = 0$, $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ undergoes a Hopf bifurcation $\bar{\mathcal{H}}$ in the $d, \varDelta \varphi$ directions and $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ becomes a saddle focus with a 1-dimensional unstable manifold (region A). At this bifurcation there appears or disappears a periodic orbit $\bar{\mathcal{T}}^{-}$ depending whether the Hopf bifurcation is supercritical or subcritical. Finally, the fixed point $\bar{ \mathcal{S}}^{-}_{\mathrm{osc}}$ disappears at a pitchfork bifurcation of the origin in the s-direction occurring at the curve $C^{-}_{ \mathrm{HB}}$.

Figure 8

Bifurcation diagram for $\bar{\mathcal{S}}^{-}_{ \mathrm{osc}}$ in the case $\beta _{\epsilon 0R} > 0$, $C_{ \mathrm{det}} + \beta _{\epsilon 0R} > 0$ and $\beta _{\epsilon 0R} < - C_{\mathrm{det}} + \sqrt{C^{2}_{\mathrm{det}} + \beta ^{2}_{\epsilon 0I}}$. The fixed point $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ appears at a supercritical pitchfork bifurcation of the origin occurring at the curve $C^{-}_{\mathrm{HB}}$, undergoes a Hopf bifurcation $\bar{ \mathcal{H}}$ at the curve $\operatorname{Tr}^{-}=0$ and becomes unstable at the curve $\operatorname{Det}^{-}=0$

Going back to the original full 4D system (3), for ϵ small enough, there exists an unstable periodic orbit $\mathcal{S}^{-}_{\mathrm{osc}}$, belonging to the torus ${\mathcal{T}} _{\epsilon }$, which will become stable at the curve $\operatorname{Det}^{-}=0$. The periodic orbit undergoes a torus bifurcation and $\mathcal{S}^{-}_{\mathrm{osc}}$ becomes unstable at the curve $\operatorname{Tr}^{-}=0$, and a new torus $\mathcal{T}^{-}$ appears or disappears depending whether the torus bifurcation is subcritical or supercritical. Finally, $\mathcal{S}^{-}_{\mathrm{osc}}$ will disappear at a Hopf bifurcation of the origin occurring at $C^{-}_{\mathrm{HB}}$.

(2) Case $\beta _{\epsilon 0R} > - C_{\mathrm{det}} + \sqrt{C^{2}_{\mathrm{det}} + \beta ^{2}_{\epsilon 0I}}$. For $\bar{\alpha }^{+} \geq 0$, λ fixed and ϵ small, the fixed point $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ is a saddle point with a 1-dimensional unstable manifold (in the Δφ direction) contained in the invariant curve $\bar{\mathcal{T}}_{ \epsilon }$ (region C in Fig. 9). As ϵ increases, $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ becomes a saddle with a 2-dimensional unstable manifold at the curve $\operatorname{Det}^{-}=0$ (region B). When further increasing the coupling ϵ, $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ becomes a saddle-focus point at the curve $\varDelta ^{-}=0$ (region A), which disappears at a pitchfork bifurcation of the origin in the s-direction occurring at the curve $C^{-}_{\mathrm{HB}}$.

Figure 9

Bifurcation diagram for $\bar{\mathcal{S}}^{-}_{ \mathrm{osc}}$ in the case $\beta _{\epsilon 0R} > 0$, $C_{ \mathrm{det}} + \beta _{\epsilon 0R} > 0$ and $\beta _{\epsilon 0R} > - C_{\mathrm{det}} + \sqrt{C^{2}_{\mathrm{det}} + \beta ^{2}_{\epsilon 0I}}$. The fixed point $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ appears at a supercritical pitchfork bifurcation of the origin in the s direction occurring at the curve $C^{-}_{\mathrm{HB}}$ and undergoes a bifurcation at the curve $\operatorname{Det}^{-}=0$

Going back to the original full 4D system (3), for ϵ small enough, there exists an unstable periodic orbit $\mathcal{S}^{-}_{\mathrm{osc}}$ belonging to the torus ${\mathcal{T}} _{\epsilon }$. The periodic orbit undergoes a bifurcation at the curve $\operatorname{Det}^{-}=0$ in which a stable manifold becomes unstable. Finally, $\mathcal{S}^{-}_{\mathrm{osc}}$ will disappear at a Hopf bifurcation of the origin occurring at $C^{-}_{\mathrm{HB}}$.

4.1.3 Regions of bistability

Since $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$ is always stable, bistability between fixed points will appear in those regions where $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ is also stable. As in the case $\beta _{\epsilon 0R} > - C_{\mathrm{det}} + \sqrt{C^{2}_{ \mathrm{det}} + \beta ^{2}_{\epsilon 0I}}$, the fixed point $\bar{ \mathcal{S}}^{-}_{\mathrm{osc}}$ is never stable, it is not possible to find bistability regions. By contrast, if $\beta _{\epsilon 0R} < - C _{\mathrm{det}} + \sqrt{C^{2}_{\mathrm{det}} + \beta ^{2}_{\epsilon 0I}}$, there exists a region in the ($\lambda , \epsilon$) parameter space defined as

$$ \operatorname{Tr}^{-}(\lambda , \epsilon ) < 0\quad \text{and}\quad \operatorname{Det}^{-}(\lambda , \epsilon ) > 0, $$

(40)

in which $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ can be either a stable node or a stable focus (see Fig. 8). Thus, the system is bistable in region (40).

Moreover, in the case $\beta _{\epsilon 0R} < - C_{\mathrm{det}} + \sqrt{C ^{2}_{\mathrm{det}} + \beta ^{2}_{\epsilon 0I}}$, the point $\bar{ \mathcal{S}}^{-}_{\mathrm{osc}}$ undergoes a Hopf bifurcation $\bar{\mathcal{H}}$. If the Hopf bifurcation is supercritical, then $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ becomes unstable and a stable limit cycle $\bar{\mathcal{T}}^{-}$ appears, generating bistability between $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$ and $\bar{\mathcal{T}} ^{-}$. The detailed analysis of this situation is beyond the scope of this paper.

Finally, we remark that the same bistable scenarios can be found in the full system (3) replacing the fixed points $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$ by the limit cycles $\mathcal{S}^{\pm }_{\mathrm{osc}}$ and the periodic orbit $\bar{ \mathcal{T}}^{-}$ by the torus $\mathcal{T}^{-}$.

4.2 Case $\beta _{\epsilon 0R} < 0$ and $C_{\mathrm{det}} + \beta _{\epsilon 0R} > 0$ (or $\beta _{\epsilon 0R} > 0$ and $C_{\mathrm{det}} + \beta _{\epsilon 0R} < 0$ )

4.2.1 Dynamics of $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$

In this case the trace for $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$ vanishes ($\operatorname{Tr}^{+} = 0$). Therefore, as

$$ \beta _{\epsilon 0R} < - C_{\mathrm{det}} < - C_{\mathrm{det}} + \sqrt{C ^{2}_{\mathrm{det}} + \beta ^{2}_{\epsilon 0I}}, $$

(41)

then $\operatorname{Tr}^{+} = 0$ and $\varDelta ^{+} < 0$ and $\bar{ \mathcal{S}}^{+}_{\mathrm{osc}}$ will always undergo a Hopf bifurcation $\bar{\mathcal{H}}$.

For $\bar{\alpha }^{+} \geq 0$, λ fixed and ϵ small, the fixed point $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$ is a stable node (region C in Fig. 10) and becomes a stable focus when the parameters cross the curve $\varDelta ^{+}=0$ (region B). For larger values of ϵ, the fixed point $\bar{\mathcal{S}}^{+} _{\mathrm{osc}}$ undergoes a Hopf bifurcation $\bar{\mathcal{H}}$ at the curve $\operatorname{Tr}^{+}=0$ and becomes a saddle-focus point (region A). At this bifurcation there appears or disappears a limit cycle $\bar{\mathcal{T}}^{+}$ depending whether this Hopf bifurcation is subcritical or supercritical. For larger values of ϵ, the fixed point $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$ disappears at a pitchfork bifurcation of the origin in the s-direction at the curve $C^{+}_{\mathrm{HB}}$.

Figure 10

Phase space for $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$ in the case $\beta _{\epsilon 0R} < 0$ and $C_{\mathrm{det}} + \beta _{\epsilon 0R} > 0$. The fixed point $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$ appears at a supercritical pitchfork bifurcation of the origin in the s direction occurring at the curve $C^{+}_{\mathrm{HB}}$ and undergoes a Hopf bifurcation $\bar{\mathcal{H}}$ at the curve $\operatorname{Tr} ^{+}=0$

Going back to the original 4D system (3), for ϵ small enough, there exists a stable periodic orbit $\mathcal{S}^{+}_{\mathrm{osc}}$. This stable periodic orbit will lose its stability across a torus bifurcation occurring at the curve $\operatorname{Tr}^{+} = 0$. At this bifurcation there appears or disappears a torus $\mathcal{T}^{+}$ depending whether the torus bifurcation is subcritical or supercritical. Finally, the unstable limit cycle $\mathcal{S}^{+}_{\mathrm{osc}}$ collapses to the origin at a Hopf bifurcation occurring at the curve $C^{+}_{\mathrm{HB}}$.

4.2.2 Dynamics of $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$

For $\bar{\alpha }^{-} \geq 0$, λ fixed and ϵ small, the fixed point $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ of system (13) is a saddle point with a 1-dimensional unstable manifold in the Δφ direction contained in $\bar{ \mathcal{T}}_{\epsilon }$ (region C in Fig. 11), and as ϵ increases it becomes a stable node when ϵ crosses the curve $\operatorname{Det} ^{-} = 0$ (region B). For larger values of ϵ, the fixed point $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ becomes a stable focus at the curve $\varDelta ^{-} = 0$ (region A) and disappears at a pitchfork bifurcation of the origin in the s direction at the curve $C^{-}_{\mathrm{HB}}$.

Figure 11

Bifurcation diagram for $\bar{\mathcal{S}}^{-}_{ \mathrm{osc}}$ in the case $\beta _{\epsilon 0R} < 0$ and $C_{ \mathrm{det}} + \beta _{\epsilon 0R} > 0$. The fixed point $\bar{ \mathcal{S}}^{-}_{\mathrm{osc}}$ appears at a supercritical pitchfork bifurcation of the origin in the s direction occurring at the curve $C^{-}_{\mathrm{HB}}$ and undergoes a bifurcation at $\operatorname{Det}^{-}=0$

Going back to the original 4D system (3), for ϵ small, there exists an unstable periodic orbit $\mathcal{S}^{-}_{\mathrm{osc}}$. This unstable periodic orbit becomes stable at the curve $\operatorname{Det}^{-} = 0$. Finally, the stable limit cycle $\mathcal{S}^{-}_{\mathrm{osc}}$ collapses to the origin at a Hopf bifurcation occurring at the curve $C^{-}_{\mathrm{HB}}$.

4.2.3 Regions of bistability

There exists a region in the ($\lambda , \epsilon$)-parameter space given by

$$ \operatorname{Tr}^{+}(\lambda , \epsilon ) < 0 \quad\text{and}\quad \operatorname{Det}^{-}(\lambda , \epsilon )> 0, $$

(42)

in which both fixed points $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$ are stable. If the Hopf bifurcation is supercritical, then $\bar{ \mathcal{S}}^{+}_{\mathrm{osc}}$ becomes unstable and a stable limit cycle $\bar{\mathcal{T}}^{+}$ appears, generating bistability between $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ and $\bar{\mathcal{T}}^{+}$. The detailed analysis of this situation is beyond of the scope of this paper.

Finally, we remark that the same bistable scenarios can be found in the full system (3) replacing the fixed points $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$ by the limit cycles $\mathcal{S}^{\pm }_{\mathrm{osc}}$ and the periodic orbit $\bar{ \mathcal{T}}^{+}$ by the torus $\mathcal{T}^{+}$.

4.3 The degenerated case $\beta _{\epsilon 0R} = 0$ and $C_{\mathrm{det}} > 0$ (or $\beta _{\epsilon 0R} = 0$ and $C_{\mathrm{det}} < 0$ )

In this case, the curves $C^{\pm }_{\mathrm{HB}}$ coincide. Moreover, the trace in (31) is identically zero for $( \lambda , \epsilon ) \in C^{\pm }_{\mathrm{HB}}$. To obtain the sign of Tr^±, we compute Tr^± when $\lambda + \epsilon \alpha _{\epsilon 0R} \rightarrow 0^{+}$. We have

$$ \operatorname{Tr}(\lambda , \epsilon ) = (\lambda + \epsilon \alpha _{\epsilon 0R}) \bigl( -2 + \mathcal{O}_{2}(\epsilon ) \bigr). $$

(43)

So, near the $C_{\mathrm{HB}}^{\pm }$ curves, both fixed points $\bar{\mathcal{S}}^{\pm }_{\mathrm{osc}}$ are stable.

4.3.1 Dynamics of $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$

For $\bar{\alpha }^{+} \geq 0$, λ fixed and ϵ small, the fixed point $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$ is a stable node (region B in Fig. 12), and as ϵ increases it becomes a stable focus when the parameters cross the curve $\varDelta ^{+} = 0$ (region A). For larger values of ϵ, the fixed point $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$ disappears at a pitchfork bifurcation of the origin in the s direction at the curve $C^{+}_{\mathrm{HB}}$.

Figure 12

Bifurcation diagram for $\bar{\mathcal{S}}^{+}_{ \mathrm{osc}}$ in the case $\beta _{\epsilon 0R} = 0$ and $C_{ \mathrm{det}} > 0$. The fixed point $\bar{\mathcal{S}}^{+}_{ \mathrm{osc}}$ appears at a supercritical pitchfork bifurcation of the origin in the s direction occurring at the curve $C^{+}_{ \mathrm{HB}}$

Going back to the original 4D system (3), for ϵ small, there exists a stable periodic orbit $\mathcal{S} ^{+}_{\mathrm{osc}}$, which collapses to the origin at a Hopf bifurcation occurring at the curve $C^{+}_{\mathrm{HB}}$.

4.3.2 Dynamics of $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$

For $\bar{\alpha }^{+} \geq 0$, λ fixed and ϵ small, the fixed point $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ is a saddle point with a 1-dimensional unstable manifold (region C in Fig. 13), and as ϵ increases it becomes a stable node when the parameters cross the curve $\operatorname{Det} ^{-} = 0$ (region B). For larger ϵ values, the fixed point $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ becomes a stable focus at $\varDelta ^{-}=0$ (region A) which collapses at a pitchfork bifurcation of the origin in the s direction at the curve $C^{-}_{\mathrm{HB}}$.

Figure 13

Phase space for the $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ fixed point in the case $\beta _{\epsilon 0R} = 0$ and $C_{ \mathrm{det}} > 0$. The fixed point $\bar{\mathcal{S}}^{-}_{ \mathrm{osc}}$ appears at a supercritical pitchfork bifurcation of the origin in the s direction occurring at the curve $C^{-}_{ \mathrm{HB}}$ and undergoes a bifurcation at the curve $\operatorname{Det}^{-}=0$

Going back to the original 4D system (3), for ϵ small, there exists an unstable periodic orbit $\mathcal{S}^{-}_{\mathrm{osc}}$ which changes stability at the curve $\operatorname{Det}^{-}=0$. Finally, the stable periodic orbit $\mathcal{S}^{-}_{\mathrm{osc}}$ collapses to the origin at a Hopf bifurcation at the curve $C^{-}_{\mathrm{HB}}$.

4.3.3 Regions of bistability

In the region in the ($\lambda , \epsilon$)-parameter space given by

$$ \operatorname{Det}^{-}(\lambda , \epsilon ) > 0 $$

(44)

both fixed points $\bar{\mathcal{S}}^{+}_{\mathrm{osc}}$ and $\bar{\mathcal{S}}^{-}_{\mathrm{osc}}$ are stable.

We remark that the same bistability scenarios can be found in the full system (3) replacing the fixed points $\bar{ \mathcal{S}}^{\pm }_{\mathrm{osc}}$ by the limit cycles $\mathcal{S} ^{\pm }_{\mathrm{osc}}$.

5 Wilson–Cowan models for perceptual bistability

Wilson–Cowan oscillators are biophysically motivated neural oscillators providing a population-averaged firing rate description of neural activity, which have been widely used to study cortical dynamics and cortical oscillations [40, 51]. The Wilson–Cowan oscillator (an excitatory (E), inhibitory (I) pair) considered here has dynamics described by

$$ \begin{aligned}& \tau \dot{E}= -E + S(aE - bI), \\ &\tau \dot{I}= -I + S(cE - dI), \end{aligned} $$

(45)

where τ is a time constant and the constants $a,b,c$ and d are the intrinsic E to E, I to E, E to I and I to I coupling weights, respectively. The function S is the sigmoidal response function

$$ S(x) = \frac{1}{1 + e^{-\lambda x + \theta }} - \frac{1}{1 + e^{ \theta }}, $$

(46)

which has threshold θ and slope λ with the convenient property $S(0) = 0$. The function S has the property $S'(0) = \lambda S_{1}$, where $S_{1} = \frac{e^{\theta }}{(1 + e^{\theta })^{2}}$, and λ is treated as a bifurcation parameter playing the equivalent role to λ in the previous sections.

The system generically has a steady state $(E,I)=(0,0)$, which undergoes a Hopf bifurcation at $\lambda _{c}=\frac{2}{(a-d)S_{1}}$. When coupled with a second, identical oscillator the 4-dimensional pair of Wilson–Cowan oscillators (E–I pairs) coupled with strength ϵ are given by

$$ \begin{aligned} &\tau \dot{E}_{1}= -E_{1} + S(aE_{1} - bI_{1}), \\ &\tau \dot{I}_{1}= -I_{1} + S \bigl(cE_{1} - dI_{1} + \epsilon (E_{2} - b _{\mathrm{sp}}I_{2}) \bigr), \\ &\tau \dot{E}_{2}= -E_{2} + S(aE_{2} - bI_{2}), \\ &\tau \dot{I}_{2}= -I_{2} + S \bigl(cE_{2} - dI_{2} + \epsilon (E_{1} - b _{\mathrm{sp}}I_{1}) \bigr), \end{aligned} $$

(47)

whose dynamics will be explored in this section.

For this study, we will consider the following set of parameters:

$$ \mathcal{P} = \{ a=7, b=5.25, c=5, d=0.7, \theta =2, \tau =1 \}, $$

(48)

whereas λ and ϵ will be the bifurcation parameters. By considering $b_{\mathrm{sp}} = -0.03, 0.03, 0.0$, we will study different types of dynamics. For each case we will write system (47) in the normal form (3) by numerically computing its corresponding coefficients (see Appendix B). Next, by using numerical continuation, we will compute bifurcation diagrams for system (47), so we can check the theoretical predictions in Sect. 4 and complete the bifurcation diagrams for large values of λ and ϵ, where the normal form approximation breaks down.

5.1 Case $b_{\mathrm{sp}} < 0$

We consider the case $b_{\mathrm{sp}} = -0.03$. The coefficients of the normal form, which were computed using the techniques described in Appendix B, are given in Table 2 and satisfy the conditions $\beta _{\epsilon 0R} > 0$, $C_{\mathrm{det}} + \beta _{\epsilon 0R} > 0$ and $\beta _{\epsilon 0R} < - C_{\mathrm{det}} + \sqrt{C^{2}_{\mathrm{det}} + \beta ^{2} _{\epsilon 0I}}$. Therefore, this case corresponds to the one considered in Sect. 4.1. Figure 14 shows the bifurcation diagram of system (47) for $b_{\mathrm{sp}} = -0.03$ obtained numerically. The results match the theoretical predictions obtained in Sect. 4.1. More precisely, for a fixed ϵ value and varying the bifurcation parameter λ, we have:

• A stable in-phase (IP) solution corresponding to $\mathcal{S}^{+}_{ \mathrm{osc}}$ will emerge from the Hopf bifurcation at $C^{+}_{ \operatorname{HB}}$. Moreover, when varying the bifurcation parameter, the IP solution will maintain its stability (see Fig. 7).

• An unstable anti-phase (AP) solution corresponding to $\mathcal{S} ^{-}_{\mathrm{osc}}$ will emerge from the Hopf bifurcation at $C^{-}_{\mathrm{HB}}$. For fixed ϵ and varying the bifurcation parameter, AP solution gains stability across a torus bifurcation, but when further increasing the bifurcation parameter, it will lose it again across a pitchfork bifurcation (corresponding respectively to the lines $\operatorname{Tr}^{-}=0$ and $\operatorname{Det}^{-}=0$ in Fig. 8).

Figure 14

Bifurcation diagram with parameters $\mathcal{P}$ and $b_{\mathrm{sp}}=-0.03$ in (47) (corresponding to the case $\beta _{\epsilon 0R} > 0$, $C_{\mathrm{det}} + \beta _{\epsilon 0R} > 0$ and satisfying $\beta _{\epsilon 0R} < - C _{\mathrm{det}} + \sqrt{C^{2}_{\mathrm{det}} + \beta ^{2}_{\epsilon 0I}}$ as described in Section 4.1). (A): Two-parameter bifurcation diagram in the $(\lambda , \epsilon )$-plane. The legend indicates bifurcations of a fixed point (FP) or a limit cycle (LC) giving rise to or involving the $\varDelta \varphi = 0$ in-phase (IP) or $\varDelta \varphi = \pi$ anti-phase (AP) solution branches; PD: period doubling; PF: pitchfork; TR: torus bifurcation. Text labels indicate the solutions that are stable in a given region, e.g. ‘IP+AP’ is a region with coexisting, stable IP and AP solutions. (B): One-parameter bifurcation diagram at $\varepsilon =0.05$ showing the FP branch, IP branch and AP branch; dashed segments are unstable. The IP and AP branches bifurcate from the FP branch in subsequent Hopf bifurcations (bullet) for λ increasing. The IP branch emerges stable and remains stable. For increasing λ, the AP branch is initially unstable, gains stability at a torus bifurcation (star) and loses stability at a pitchfork bifurcation (diamond). (C): Coexisting solutions at $\lambda \approx 3.05$ and $\epsilon =0.05$ in the $(E_{1},E_{2})$-plane. Motion on the diagonal (blue) corresponds to in-phase oscillations. (D): As (C) in the $(E_{1},I_{1})$-plane for one E–I oscillator. (E): As (C) at $\epsilon =0.5$, where a torus bifurcation (star) is on an unstable branch that gains stability at a fold of limit cycle (square)

Table 2 Coefficients of the normal form (3) for the three considered cases, namely $b_{\mathrm{sp}} = -0.03, 0.03$ and 0. These coefficients have been computed using the procedure described in Appendix B

5.2 Case $b_{\mathrm{sp}} > 0$

We consider the case $b_{\mathrm{sp}} = 0.03$. The coefficients of the normal form, which were computed using the techniques described in Appendix B, are given in Table 2 and satisfy the conditions $\beta _{\epsilon 0R} < 0$ and $C_{\mathrm{det}} + \beta _{\epsilon 0R} > 0$. Therefore, this case corresponds to the one considered in Sect. 4.2. Figure 15 shows the bifurcation diagram of system (47) for $b_{\mathrm{sp}} = 0.03$ obtained numerically. The results match the theoretical predictions in Sect. 4.2. More precisely, for a fixed ϵ value and varying the bifurcation parameter λ, we have:

• A stable anti-phase (AP) solution corresponding to $\mathcal{S}^{-} _{\mathrm{osc}}$ will emerge from a Hopf bifurcation at $C^{-}_{ \mathrm{HB}}$, whereas an unstable in-phase (IP) solution corresponding to $\mathcal{S}^{+}_{\mathrm{osc}}$ will emerge from the Hopf bifurcation at $C^{+}_{\mathrm{HB}}$.

• The stability of both solutions is reversed as the bifurcation parameter grows. Moreover, the bifurcations giving rise to these stability changes are of the same type as we predicted: IP solution becomes stable across a torus bifurcation (corresponding to the Hopf bifurcation $\bar{ \mathcal{H}}$ at the $\operatorname{Tr}^{+}=0$ line in Fig. 10), whereas the AP solution loses stability across a pitchfork bifurcation of limit cycles (corresponding to the $\operatorname{Det}^{-}=0$ line in Fig. 11).

Figure 15

Bifurcation diagram with parameters $\mathcal{P}$ and $b_{\mathrm{sp}}=0.03$ in (47) (corresponding to the case $\beta _{\epsilon 0R} < 0$ and $C_{\mathrm{det}} + \beta _{\epsilon 0R} > 0$, as described in Section 4.2). (A): Two-parameter bifurcation diagram in the $(\lambda ,\epsilon )$-plane. Legends and labelling as in Fig. 14; TR: torus bifurcation. (B): One-parameter bifurcation diagram at $\varepsilon =0.05$ showing the FP branch, IP branch and AP branch; dashed segments are unstable. The AP and IP branches bifurcate from the FP branch in subsequent Hopf bifurcations (bullet) for λ increasing. The AP branch loses stability in a pitchfork bifurcation (diamond). The IP branch is initially unstable and gains stability at a torus bifurcation (star)

5.3 Case $b_{\mathrm{sp}} = 0$

We consider the case $b_{\mathrm{sp}} = 0.0$. The coefficients of the normal form, which were computed using the techniques described in Appendix B, are given in Table 2 and satisfy the conditions $\beta _{\epsilon 0R} = 0$ and $C_{\mathrm{det}} > 0$. Therefore, this case corresponds to the “degenerated case” discussed in Sect. 4.3. Figure 16 shows the bifurcation diagram of system (47) for $b_{\mathrm{sp}} = 0$ obtained numerically. Notice that it matches the theoretical predictions, namely:

Figure 16

Bifurcation diagram with parameters $\mathcal{P}$ and $b_{\mathrm{sp}}=0$ in (47) (corresponding to the “degenerated case” in Section 4.3). (A): Two-parameter bifurcation diagram where curves are the locus of bifurcations in the $(\mu ,\epsilon )$-plane. The legend indicates bifurcations of a fixed point (FP) or a limit cycle (LC) giving rise to or involving the $\varDelta \varphi = 0$ in-phase (IP) or $\varDelta \varphi = \pi$ anti-phase (AP) solution branches; PD: period doubling; PF: pitchfork. Text labels indicate the solutions that are stable in a given region, e.g. ‘IP+AP’ is a region with coexisting, stable IP and AP solutions. (B): One-parameter bifurcation diagram for fixed $\epsilon =0.05$ showing the fixed point branch, IP branch and AP branch; dashed segments are unstable. The IP and AP branches bifurcation from the FP branch at a degenerate Hopf bifurcation (bullet). The AP branch loses stability in a pitchfork bifurcation (diamond)

• Both Hopf bifurcation curves $C^{\pm }_{\mathrm{HB}}$ coincide and give rise to a bistable situation. On one side of the double Hopf curve, there exists bistability between the in-phase (IP) solution $\varDelta \varphi = 0$ corresponding to $\mathcal{S}^{+}_{\mathrm{osc}}$ and the anti-phase (AP) $\varDelta \varphi = \pi$ solution corresponding to $\mathcal{S}^{-}_{\mathrm{osc}}$.

• For ϵ fixed and increasing the bifurcation parameter λ, the $\mathcal{S}^{-}_{\mathrm{osc}}$ (AP) solution loses stability across a pitchfork bifurcation of limit cycles that we found for the 3D system as the line having $\operatorname{Det}^{-}=0$ (see Fig. 13).

5.4 Dynamics beyond the weak coupling limit

Our numerical bifurcation analysis has revealed the possibility for richer dynamics, whilst noting a wide range of parameters for which the IP and AP solutions are stable and coexist. Furthermore, a Bautin bifurcation on the AP Hopf branch for $\epsilon _{\mathrm{BT}}\approx 0.4$, as seen in Figs. 14, 15 and 16, gives rise to a region of parameter space for $\lambda \lesssim \lambda _{c}$ where a stable AP solution coexists with a stable fixed point. The bifurcation point $\epsilon _{ \mathrm{BT}}$ separates branches of sub- and supercritical Hopf bifurcations in the parameter space. As we can see, for nearby $\lambda , \epsilon$ parameter values, the system has two limit cycles which collide and disappear via a fold bifurcation of periodic orbits. Although the analysis done in Sects. 3 and 4 is restricted to the weak coupling case, we briefly discuss how the reduced system (13) can provide some insight about this bifurcation.

In the weak coupling regime, the denominator in formula (26) for the $s^{\pm }_{\mathrm{osc}}$ solutions is given by $\alpha _{01R} + \epsilon K_{\mathrm{stb}}^{ \pm }$ and is assumed to be negative. Therefore, $s^{\pm }_{ \mathrm{osc}}$ solutions appear for $\alpha ^{\pm }= \lambda + \epsilon (\alpha _{\epsilon 0R} \pm \beta _{\epsilon 0R}) > 0$ at a supercritical pitchfork bifurcation of the origin (see Fig. 5). Nevertheless, writing the equation for s in (25) in the following way:

$$ \dot{s}=A(\lambda ,\epsilon )s + B(\lambda ,\epsilon )s^{3}, $$

we clearly see that at the curve $A(\lambda ,\epsilon )=0$ the origin undergoes a pitchfork bifurcation that is supercritical or subcritical depending on the sign of $B(\lambda ,\epsilon )$. Consequently, the point ($\lambda , \epsilon$) satisfying $A(\lambda ,\epsilon ) = 0$ and $B(\lambda ,\epsilon ) = 0$ corresponds to a Bautin bifurcation. Thus, using the expression for A and B (which are known up to first order in ϵ and λ), we can estimate that a Bautin bifurcation occurs for

$$ \epsilon _{\mathrm{BT}} \approx -\frac{\alpha _{01R}}{K^{-}_{ \mathrm{stb}}}, $$

(49)

assuming that $K^{-}_{\mathrm{stb}}>0$ and for $\lambda _{\mathrm{BT}}$ such that $(\lambda _{\mathrm{BT}}, \epsilon _{\mathrm{BT}}) \in C_{ \mathrm{HB}}^{-}$. Although an accurate derivation is beyond the scope of this work, this transition from subcritical to supercritical involves the appearance of a curve of saddle-node bifurcations of fixed points for system (13) for nearby values of the parameters. More precisely, if we consider the exact expression of the determinant of the $2\times 2$ block of Jacobian matrix (29) given by

$$ \operatorname{Det} \bigl(s_{\mathrm{osc}}^{-} \bigr) = c_{d}^{d} c_{\varphi }^{ \varphi }- c_{d}^{\varphi }c_{\varphi }^{d}, $$

(50)

where the constants are given by Eqs. (56) in Appendix A with $s = s_{\mathrm{osc}}^{-}$ in (26), one can see that it is singular at $B(\lambda ,\epsilon ) = 0$. Therefore, we consider the curve

$$ B(\lambda ,\epsilon ) \operatorname{Det} \bigl(s_{\mathrm{osc}}^{-} \bigr) = 0, $$

and one can see that the Bautin point ($\lambda _{\mathrm{BT}}, \epsilon _{\mathrm{BT}}$) belongs to it. Moreover, for $\epsilon > \epsilon _{\mathrm{BT}}$ as $B(\lambda , \epsilon ) > 0$, this curve corresponds to the saddle-node bifurcations of the solutions $s_{\mathrm{osc}}^{-}$ outside the $C_{\mathrm{HB}}^{-}$ curve.

Using the numerical values given in Table 2, $K_{\mathrm{stb}}^{-}>0$. Thus, we can estimate from the normal form that the Bautin bifurcation occurs for $\epsilon _{\mathrm{BT}} \approx 0.42, 0.43, 0.42$ for $b_{\mathrm{sp}} = -0.03, 0.03, 0$, respectively, which matches the results obtained numerically (see Figs. 14, 15 and 16). Recall that in the original 4D system (3) the pitchfork and saddle-node bifurcations correspond to Hopf and fold of limit cycles bifurcations, respectively.

Besides this previous behaviour, we also remark that the IP solution undergoes a period-doubling bifurcation for large ϵ and λ leading to richer dynamical behaviour away from the analytically-investigated uncoupling limit.

5.5 Periodically forced coupled Wilson–Cowan equations

With the aim of finding coexisting IP and AP solutions (corresponding to “percept 1” and “percept 2” as described in Sect. 1), we now introduce periodic forcing terms to the coupled WC system given by (47). We consider anti-phase inputs with forcing frequency $f=2.5$ Hz and amplitude A which will be varied as a bifurcation parameter:

$$ \begin{aligned} &\tau \dot{E}_{1}= -E_{1} + S \bigl(aE_{1} - bI_{1}+A\sin ^{2n}(2\pi f t)+(1-h)A \cos ^{2n}(2\pi f t) \bigr), \\ &\tau \dot{I}_{1}= -I_{1} + S \bigl(cE_{1} - dI_{1} + \epsilon (E_{2} - b _{\mathrm{sp}}I_{2}) \bigr), \\ &\tau \dot{E}_{2}= -E_{2} + S \bigl(aE_{2} - bI_{2}+A\cos ^{2n}(2\pi f t)+(1-h)A \sin ^{2n}(2\pi f t) \bigr)), \\ &\tau \dot{I}_{2}= -I_{2} + S \bigl(cE_{2} - dI_{2} + \epsilon (E_{1} - b _{\mathrm{sp}}I_{1}) \bigr), \end{aligned} $$

(51)

where the parameters $\mathcal{P}$ (with the exception of τ) and nonlinearity (46) are as above. The input asymmetry parameter h controls the balance of inputs across the two oscillators; when $h=1$ the oscillators receive exclusive inputs (the case typically considered in competition models [24, 28, 43, 49]), and when $h=0$ the oscillators receive identical inputs (the case considered here). The forcing terms are raised to an even power 2n with $n=5$ to be positive and sharpened such that the anti-phase inputs do not overlap in time. Noting that the isolated Wilson–Cowan oscillator undergoes a supercritical Hopf bifurcation at $\lambda =\frac{2}{(a-d)S _{1}}=3.025$, we set $\lambda =2.6$ before this bifurcation. Further, noting that the bifurcating branch emerges with period

$$ T=\tau \frac{2\pi }{\sqrt{\lambda ^{2}S_{1}^{2}( b c - a d)+ \lambda S_{1} (d-a) + 1 }} $$

(52)

and fixing $T=\frac{1}{2f}$, we can set $\tau =\frac{\sqrt{\lambda ^{2}S_{1}^{2}( b c - a d)+ \lambda S_{1} (d-a) + 1 }}{4f\pi }$ such that the frequencies of oscillations produced at the Hopf match the forcing frequency.

Figure 17 shows a bifurcation diagram for the pair periodically-forced Wilson–Cowan oscillators. Each E–I oscillator receives the same input ($h=0$). Panel (A) shows regions of the $(\epsilon ,A)$ plane in which different types of oscillatory behaviours are stable. For low forcing amplitude, there are only low-amplitude oscillations, effectively modulating the FP solution in the unforced system. As A is increased, pitchfork bifurcations give rise to stable IP and AP branches that coexist (see panel (B)) for small ϵ approaching the uncoupling limit. For large ϵ, the IP solution persists at intermediate values of A. For large A, there is a saturated high-amplitude solution.

Figure 17

Bifurcation diagram with parameters $\mathcal{P}$ whilst setting λ and τ as described in the text. (A): Two-parameter bifurcation diagram in the $(A,\epsilon )$ plane showing locus of bifurcations with legends and labelling as in (16); LA is symmetric ($(E_{1},I_{1})=(E _{2},I_{2})$) low-amplitude limit cycle oscillations (following the periodic input) and HA is a symmetric high-amplitude limit cycle. The IP and AP solutions coexist in the region up to the dashed fold curve to the right. (B): One-parameter bifurcation diagram at fixed $\epsilon =0.5$; dashed curve segments are unstable. Diamonds are pitchfork bifurcations and squares are fold bifurcations. The stable IP branch exists between a pitchfork bifurcation to the left and fold to the right. The AP branch emerges unstable and is stable between a secondary pitchfork bifurcation on the left and a fold bifurcation to the right

The key result here is that the behaviour found in the unforced system is preserved for sufficiently small coupling strength and for weak forcing (IP and AP solutions persist close to the uncoupling limit, IP+AP region in Fig. 17(A)). For larger forcing amplitude, the intrinsic dynamics is overwhelmed and the forcing modulates a symmetrical fixed point (HA region in Fig. 17(A)). This bifurcation analysis demonstrates the possibility for coexisting in-phase and anti-phase responses of the coupled Wilson–Cowan oscillators to encode network states corresponding to “percept 1” (IP) and “percept 2” (AP) as described in Sect. 1. This is possible without strong mutual inhibition (i.e. in the uncoupling limit) between abstract representations of the two possible percepts.

6 Discussion and conclusions

The study of identical coupled oscillators near a Hopf bifurcation is applicable to a wide range of systems where near-identical units undergo oscillatory instability. These systems may in general be represented by very different vector fields. Using the normal form theory in [8], we are able to predict universal aspects of the mathematical behaviour for such systems. The analysis performed in this work for two oscillators reveals that, as is often the case in normal forms, although (3) involves a big number of parameters, in the weak coupling limit just a few of them govern and determine the possible bifurcations of the system.

Because of the symmetries of the system, there are usually two phase-locked oscillating solutions corresponding to in-phase ($\varDelta \varphi = 0$) and anti-phase ($\varDelta \varphi = \pi$). Depending on parameters, we find that all possible combinations between different stabilities of both solutions are possible. Our numerical analysis has shown that away from the coupling limit, richer dynamical behaviour is possible, with secondary bifurcations from the anti-phase branch and regions of coexistence between fixed-point and anti-phase solutions mediated by a fold of cycles. These scenarios can include modulated states that appear at torus bifurcations (see for example Fig. 15). Furthermore, we find that the coexistence of in-phase and anti-phase solutions persists even in the presence of periodic forcing.

6.1 Implications for models of perceptual bistability and neural competition

Models of perceptual bistability are widely based on the assumption of strong mutual inhibition between populations of neurons that encode different perceptual interpretations of ambiguous stimuli. In general, this assumes that populations associated with different percepts are separated in some feature space (e.g. orientation in binocular rivalry) and that these populations enter into competition through mutual inhibition. However, when stimuli are periodic and the two possible perceptual interpretations involve the same features, it is less clear how competition between percepts might arise. For example, for the visual (auditory) stimulus in Fig. 1 both “percept 1” and “percept 2” involve the left spatial location (higher pitch A tone). It is therefore unclear how mutual inhibition between “percept 1” and “percept 2” could be implemented in neural hardware (although see [39] where population pooling inputs from an intermediate feature location were proposed). Another possibility, proposed and demonstrated to be feasible in this study, involves oscillatory neural activity. Indeed, encoding of perceptual interpretations through oscillations allows for complete synchronisation of the network with all incoming inputs (like “percept 1”) or for partial synchronisation of different parts of a network with separate elements (here in anti-phase). Furthermore, such an encoding mechanism does not rely on strong mutual inhibition, widely assumed between the abstracted percept-based neural populations in competition models with little supporting evidence.

6.2 Future perspectives

An obvious extension of the bifurcation analysis would be to the forced symmetry broken case. If there is no assumed symmetry between percepts 1 and 2, this will result in a separation of Hopf bifurcations in the uncoupled limit and presumably mode locking and torus breakup scenarios familiar from the non-symmetric Hopf–Hopf interaction case [14]. Finally, one can consider the periodically forced system. Periodic forcing of the oscillators considered here (e.g. [21] for a single oscillator) will bring us to potentially much more complex bifurcation problems.

The study has demonstrated the potential role of oscillations in encoding different interpretations of periodically modulated ambiguous stimuli. It remains to explore the further role of feature space (say spatial location or tone frequency) and its interaction with oscillatory mechanisms. Additionally, as bistable perception involves spontaneous switching between perceptual interpretations, the mechanisms for these switches in the light of oscillatory stimuli remain to be explored.

Perceptual bistability with periodically modulated stimuli is robust over a range of input rates for the stimulus, whereas the simple network motif studied here has a fixed preferred input rate. So-called gradient networks of coupled oscillators have been proposed as a framework to understand many elements of early auditory processing and for perception of musical rhythm and beat [25, 26]. Such a framework could be extensible to the study of perceptual bistability, relying on the dynamic mechanisms proposed here in the simple case of only two coupled oscillators.

Appendices

Appendix 1: Coupling terms

The coupling terms $f_{r}$ and $f_{\varphi }$ for system (9) are given by

$$ \begin{aligned} &f_{r} (r_{1}, r_{2}, \varDelta \varphi ) \\ &\quad=r^{2}_{1}r_{2} \bigl[( \beta _{\epsilon 1R} + \alpha _{\epsilon 3R})\cos (\varDelta \varphi ) - ( \beta _{\epsilon 1I} - \alpha _{\epsilon 3I})\sin (\varDelta \varphi ) \bigr] \\ &\qquad {}+ r^{2}_{2}r_{1} \bigl[\alpha _{\epsilon 2R} + \beta _{\epsilon 3R} \cos (2\varDelta \varphi ) - \beta _{\epsilon 3I}\sin (2 \varDelta \varphi ) \bigr] + r_{1} \alpha _{\epsilon 0R} + r^{3}_{1} \alpha _{\epsilon 1R} \\ &\qquad{} + r^{3}_{2} \bigl[\beta _{\epsilon 2R}\cos (\varDelta \varphi ) - \beta _{\epsilon 2I}\sin (\varDelta \varphi ) \bigr] + r_{2} \bigl[\beta _{ \epsilon 0R}\cos (\varDelta \varphi ) - \beta _{\epsilon 0I} \sin (\varDelta \varphi ) \bigr], \\ &f_{\varphi }(r_{1}, r_{2}, \varDelta \varphi )\\ &\quad= r^{2}_{1}r_{2} \bigl[ (\beta _{\epsilon 1I} + \alpha _{\epsilon 3I})\cos (\varDelta \varphi ) + ( \beta _{\epsilon 1R} - \alpha _{\epsilon 3R})\sin (\varDelta \varphi ) \bigr] \\ &\qquad{}+ r^{2}_{2}r_{1} \bigl[ \alpha _{\epsilon 2I} + \beta _{\epsilon 3I} \cos (2\varDelta \varphi ) + \beta _{\epsilon 3R}\sin (2 \varDelta \varphi ) \bigr] + r_{1} \alpha _{\epsilon 0I} + r^{3}_{1} \alpha _{\epsilon 1I} \\ &\qquad{} + r^{3}_{2} \bigl[ \beta _{\epsilon 2I}\cos ( \varDelta \varphi ) + \beta _{\epsilon 2R}\sin (\varDelta \varphi ) \bigr] + r_{2} \bigl[ \beta _{\epsilon 0I}\cos (\varDelta \varphi ) + \beta _{\epsilon 0R} \sin (\varDelta \varphi ) \bigr]. \end{aligned} $$

(53)

The coupling term $f_{\varDelta \varphi }$ for system (10) is given by

$$\begin{aligned} &f_{\varDelta \varphi }(r_{1}, r_{2}, \varDelta \varphi ) \\ &\quad= f_{\varphi }(r _{2}, r_{1}, -\varDelta \varphi )/r_{2} - f_{\varphi }(r_{1}, r_{2}, \varDelta \varphi )/r_{1} \\ &\quad =\bigl(r_{1}^{2} - r^{2}_{2} \bigr) \bigl[ \alpha _{\epsilon 2I} - \alpha _{\epsilon 1I} + \beta _{\epsilon 3I}\cos (2 \varDelta \varphi ) \bigr] - 2r_{1}r_{2}( \beta _{\epsilon 1R} - \alpha _{\epsilon 3R})\sin (\varDelta \varphi ) \\ &\qquad{}- \bigl(r_{1}^{2} + r^{2}_{2} \bigr) \beta _{\epsilon 3R}\sin (2 \varDelta \varphi ) + \biggl( \frac{r^{3}_{1}}{r_{2}} - \frac{r^{3}_{2}}{r_{1}} \biggr) \beta _{\epsilon 2I}\cos (\varDelta \varphi ) - \biggl( \frac{r^{3}_{1}}{r _{2}} + \frac{r^{3}_{2}}{r_{1}} \biggr) \beta _{\epsilon 2R}\sin ( \varDelta \varphi ) \\ &\qquad{} + \biggl( \frac{r_{1}}{r_{2}} - \frac{r_{2}}{r _{1}} \biggr)\beta _{\epsilon 0I}\cos (\varDelta \varphi ) - \biggl(\frac{r _{1}}{r_{2}} + \frac{r_{2}}{r_{1}} \biggr)\beta _{\epsilon 0R}\sin ( \varDelta \varphi ). \end{aligned}$$

(54)

The coupling terms $g_{s}$, $g_{d}$ and $g_{\varDelta \varphi }$ for system (13) are given by

$$\begin{aligned} & g_{s}(s, d, \varDelta \varphi ) \\ &\quad = f_{r} \biggl(\frac{s+d}{2}, \frac{s-d}{2}, \varDelta \varphi \biggr) + f_{r} \biggl( \frac{s-d}{2}, \frac{s+d}{2}, -\varDelta \varphi \biggr) \\ &\quad= s \bigl(\cos (\varDelta \varphi )\beta _{\epsilon 0R} + \alpha _{\epsilon 0R} \bigr) + d\sin (\varDelta \varphi )\beta _{\epsilon 0I} + \frac{s}{4} \bigl(s ^{2} + 3d^{2} \bigr) \bigl( \beta _{\epsilon 2R}\cos (\varDelta \varphi ) + \alpha _{\epsilon 1R} \bigr) \\ &\qquad{} + \frac{s}{4} \bigl(s^{2} - d^{2} \bigr) \bigl[ (\beta _{\epsilon 1R} + \alpha _{\epsilon 3R})\cos (\varDelta \varphi ) + \alpha _{\epsilon 2R} + \beta _{\epsilon 3R}\cos (2\varDelta \varphi ) \bigr] \\ &\qquad {}+ \frac{d}{4} \bigl(s^{2} - d^{2} \bigr) \bigl[ \beta _{\epsilon 3I} \sin (2\varDelta \varphi ) - (\beta _{\epsilon 1I} - \alpha _{\epsilon 3I}) \sin (\varDelta \varphi ) \bigr] \\ &\qquad{} + \frac{d}{4} \bigl(3s^{2} + d^{2} \bigr) \beta _{\epsilon 2I}\sin ( \varDelta \varphi ), \\ &g_{d}(s, d, \varDelta \varphi ) \\ &\quad = f_{r} \biggl( \frac{s+d}{2}, \frac{s-d}{2}, \varDelta \varphi \biggr) - f_{r} \biggl(\frac{s-d}{2}, \frac{s+d}{2}, -\varDelta \varphi \biggr) \\ &\quad= {-}d \bigl(\cos (\varDelta \varphi )\beta _{\epsilon 0R} - \alpha _{\epsilon 0R} \bigr) - s\sin (\varDelta \varphi )\beta _{\epsilon 0I} \\ &\qquad{} - \frac{d}{4} \bigl(d ^{2} + 3s^{2} \bigr) \bigl( \beta _{\epsilon 2R}\cos (\varDelta \varphi ) - \alpha _{\epsilon 1R} \bigr) \\ &\qquad{} + \frac{d}{4} \bigl(s^{2} - d^{2} \bigr) \bigl[ (\beta _{\epsilon 1R} + \alpha _{\epsilon 3R})\cos (\varDelta \varphi ) - \alpha _{\epsilon 2R} - \beta _{\epsilon 3R}\cos (2\varDelta \varphi ) \bigr] \\ &\qquad {}- \frac{s}{4} \bigl(s^{2} - d^{2} \bigr) \bigl[ \beta _{\epsilon 3I} \sin (2\varDelta \varphi ) + (\beta _{\epsilon 1I} - \alpha _{\epsilon 3I}) \sin (\varDelta \varphi ) \bigr] \\ &\qquad{} - \frac{s}{4} \bigl(3d^{2} + s^{2} \bigr) \beta _{\epsilon 2I}\sin ( \varDelta \varphi ), \\ &g_{\varDelta \varphi }(s, d, \varDelta \varphi ) \\ &\quad = f_{\varDelta \varphi } \biggl( \frac{s+d}{2}, \frac{s-d}{2}, \varDelta \varphi \biggr) \\ &\quad = \beta _{\epsilon 0I}\cos (\varDelta \phi ) \biggl(\frac{4sd}{s^{2} - d ^{2}} \biggr)- 2\beta _{\epsilon 0R}\sin (\varDelta \phi ) \biggl(\frac{s ^{2} + d^{2}}{s^{2} - d^{2}} \biggr) \\ &\qquad{} - \beta _{\epsilon 2R}\sin (\varDelta \phi ) \biggl(\frac{(s^{2} + d^{2})^{2}}{(s ^{2} - d^{2})} - \frac{(s^{2} - d^{2})}{2} \biggr) + \beta _{\epsilon 2I}\cos (\varDelta \phi ) \frac{2sd(s^{2} + d^{2})}{(s^{2} - d^{2})} \\ &\qquad {}- \beta _{\epsilon 3R}\sin (2\varDelta \phi )\frac{(s^{2} + d^{2})}{2} - (\beta _{\epsilon 1R} - \alpha _{\epsilon 3R})\sin (\varDelta \phi ) \frac{(s ^{2} - d^{2})}{2} \\ &\qquad{}+ \bigl(\alpha _{\epsilon 2I} + \beta _{\epsilon 3I}\cos (2\varDelta \phi ) - \alpha _{\epsilon 1I} \bigr)sd. \end{aligned}$$

(55)

The terms for the Jacobian matrix in (29) are given by

$$ \begin{aligned}& c^{s}_{s}= \lambda + \epsilon (\alpha _{\epsilon 0R} \pm \beta _{\epsilon 0R}) \\ &\phantom{ c^{s}_{s}=}{}+ \frac{3s^{2}}{4} \bigl(\alpha _{01R} + \epsilon \bigl(\alpha _{\epsilon 1R} \pm (\beta _{\epsilon 2R} + \beta _{\epsilon 1R} + \alpha _{\epsilon 3R}) + \alpha _{\epsilon 2R} + \beta _{\epsilon 3R} \bigr) \bigr), \\ &c^{d}_{d}= \lambda + \epsilon (\alpha _{\epsilon 0R} \mp \beta _{\epsilon 0R}) \\ &\phantom{c^{d}_{d}=}{}+ \frac{s^{2}}{4} \bigl(3 \alpha _{01R} + \epsilon \bigl(3(\alpha _{\epsilon 1R} \mp \beta _{\epsilon 2R}) \pm (\beta _{\epsilon 1R} + \alpha _{\epsilon 3R}) - \alpha _{\epsilon 2R} - \beta _{\epsilon 3R} \bigr) \bigr), \\ &c^{d}_{\varDelta \varphi }= \epsilon \biggl(- \frac{s^{3}}{4} \bigl(2\beta _{\epsilon 3I} \pm (\beta _{\epsilon 1I} - \alpha _{\epsilon 3I}) \pm \beta _{\epsilon 2I} \bigr) \mp \beta _{\epsilon 0I}s \biggr), \\ &c^{\varDelta \varphi }_{d}= - \alpha _{01I}s + \epsilon \biggl(s(\alpha _{\epsilon 2I} - \alpha _{\epsilon 1I} + \beta _{\epsilon 3I} \pm 2\beta _{\epsilon 2I}) \pm 4\frac{\beta _{\epsilon 0I}}{s} \biggr), \\ &c^{\varDelta \varphi }_{ \varDelta \varphi }= \epsilon \biggl(\frac{s^{2}}{2} \bigl(\mp ( \beta _{\epsilon 1R} - \alpha _{\epsilon 3R}) - 2\beta _{\epsilon 3R} \mp \beta _{\epsilon 2R} \bigr) \mp 2\beta _{\epsilon 0R} \biggr), \end{aligned} $$

(56)

where $s = s^{\pm }_{\mathrm{osc}}$ and, consistently with the notation used throughout the article, the ± sign corresponds to $\varDelta \varphi = 0, \pi$ respectively.

Appendix 2: Normal form computation

In this appendix we provide a brief description of the numerical procedure used to compute the coefficients of the normal form (3). The procedure is related to normal form techniques in which one constructs a change of variables of the form

$$ z = y + Q_{2}(y, \bar{y}) + Q_{3}(y, \bar{y}), $$

(57)

where $Q_{2}(y, \bar{y})$ and $Q_{3}(y, \bar{y})$ are polynomials or order two and three, respectively, such that system (1) expressed in the variables y and ȳ has the simplest expression possible. That is,

$$ \dot{y} = Ay + f_{2}(y, \bar{y}) + f_{3}(y, \bar{y}) + \mathcal{O} _{4}(y, \bar{y}), $$

(58)

where A is the linearised system (3) around the origin, $f_{2} = 0$ and $f_{3}$ has the same monomials appearing in (3), namely $y_{1}|y_{1}|^{2}$, $y_{1}|y_{2}|^{2}$, $y_{4}y_{1}^{2}$, $y_{2}|y_{1}|^{2}$, $y_{2}|y_{2}|^{2}$ and $y_{3}y_{2}^{2}$.

To that aim we perform the following steps:

• Consider the Taylor expansion of system (1) around the origin:

  $$ \begin{aligned} &\dot{z}= Az + P_{2}(z, \bar{z}) + P_{3}(z, \bar{z}) + \mathcal{O} _{4}(z, \bar{z}), \\ &\dot{\bar{z}}= \bar{A}\bar{z} + \bar{P_{2}}(z, \bar{z}) + \bar{P_{3}}(z, \bar{z}) + \mathcal{O}_{4}(z, \bar{z}), \end{aligned} $$

  (59)

  where $z = (z_{1}, z_{2})$, $\bar{z} = (\bar{z}_{1}, \bar{z}_{2})\in \mathbb{C}^{2}$, $A = \operatorname{diag}(\mu ^{+}, \mu ^{-})$ is a diagonal matrix with $\mu ^{+}, \mu ^{-} \in \mathbb{C}$. $P_{2}$ and $P_{3}$ in (59) correspond to polynomials of degree 2 and 3, respectively. As z̄ is the complex conjugate of z, we will just consider the first equation in (59).

• Compute the $q_{ij}$ coefficients of the polynomial $Q_{2}(y, \bar{y})$ given by

  $$ Q_{2}(y, \bar{y}) = \prod ^{N=4}_{i=1} \prod ^{N=4}_{j=i}q_{ij}y_{i}y_{j}, $$

  (60)

  where $y_{3} = \bar{y}_{1}$ and $y_{4} = \bar{y}_{2}$, by solving the following equation for each monomial:

  $$ AQ_{2}(y, \bar{y}) - D_{y}Q_{2}(y, \bar{y}) Ay - D_{\bar{y}}Q_{2}(y, \bar{y}) \bar{A}\bar{y} = f_{2}(y, \bar{y}) - P_{2}(y, \bar{y}). $$

  (61)

  With this choice, all the monomials in $f_{2}$ in (58) are null.

• Compute $f_{3}(y, \bar{y})$ given by the expression

  $$ f_{3}(y, \bar{y}) = D_{y}P_{2}(y, \bar{y})Q_{2}(y, \bar{y}) + D_{ \bar{y}}P_{2}(y, \bar{y})\bar{Q}_{2}(y, \bar{y}) + P_{3}(y, \bar{y}), $$

  (62)

  thus obtaining the coefficients corresponding to the surviving monomials in (3): $y_{i}|y_{i}|^{2}$, $y^{2}_{i}\bar{y}_{j}$, $y_{j}|y_{i}|^{2}$, $y_{i}|y_{j}|^{2}$, $y^{2}_{j}\bar{y}_{i}$, $y_{j}|y_{j}|^{2}\ (i = 1,2 , j \neq i)$.

• Perform the change of coordinates $y = Cx$ in system (58), where

  $$ C = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & -1 \end{pmatrix} , $$

  (63)

  so that the system is written in the form (3).

Notice that to compute the coefficients of $f_{3}$ in (58) it is enough to compute the change in (57) up to order two. As a final remark, notice that apart from ω and $\alpha _{01}$ all the coefficients $\alpha _{\epsilon i},\beta _{\epsilon i}$ ($i = 0,\ldots,3$) in (3) are multiplied by ϵ. Therefore, to obtain the value of the coefficients, we follow the procedure described above for $\epsilon = 0$, thus obtaining ω and $\alpha _{01}$, and then repeat the same procedure for small $\epsilon \neq 0$, which, using that ω and $\alpha _{01}$ are known, provides the coefficients $\alpha _{\epsilon i}, \beta _{\epsilon i}$.