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

We study the dynamics arising when two identical oscillators are coupled near a Hopf bifurcation where we assume a parameter ϵ uncouples the system at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\epsilon =0$\end{document}ϵ=0. Using a normal form for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$N=2$\end{document}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.


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 spatiotemporal 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).

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 ∞ ) vector fields on x i ∈ R d (d ≥ 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 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) 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.

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, a 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]. b 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).
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.

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.

Two identical Hopf bifurcations with an uncoupling limit
We will study systems consisting of two identically coupled oscillators of the form: having S 2 permutation symmetry. We assume that when system (1) is uncoupled ( = 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 ∈ R 2 given by has a stable focus at x = 0 for λ < 0 that undergoes a supercritical Hopf bifurcation for λ = 0, which gives rise to a small amplitude stable limit cycle for λ > 0. For simplicity we assume that the eigenvalues of DH λ (0) are λ ± iω with ω = 0. Moreover, without loss of generality, we assume that (x 1 , x 2 ) = (0, 0) is an equilibrium point for (λ, ) in some neighbourhood of (0, 0) for system (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 λ = 0, can be written in the following normal form: where F N is an N -degree polynomial function that is equivariant under the rotational symmetries F N z 1 e iφ , z 2 e iφ , = e iφ F N (z 1 , z 2 , ).
If we consider the normal form up to order three and ignore the O 4 (z) terms, we obtain the truncated normal form where the constants α 01 , α i , β i ∈ C with the restriction Re(α 01 ) < 0 because the Hopf bifurcation is supercritical.
At each curve C ± HB , a limit cycle, that will be denoted by S ± osc , is born. To study the stability of the origin of system (3), we analyse the sign of the real part of its eigenvalues μ + and μgiven in (6) at the Hopf bifurcation curves C ± HB defined in (7). Thus, Therefore, we conclude that (see Fig. 2): -If β 0R > 0, for (λ, ) ∈ C + HB , the solution S 0 changes from a stable focus to a saddle-focus and a stable limit cycle S + osc emerges from C + HB . Moreover, when (λ, ) ∈ C -HB , the solution S 0 changes from a saddle-focus to an unstable focus and a saddle limit cycle Sosc appears. -If β 0R < 0, for (λ, ) ∈ C -HB , the solution S 0 changes from a stable focus to a saddle-focus and a stable limit cycle Sosc emerges from C -HB . Moreover, when (λ, ) ∈ C + HB , the solution S 0 changes from a saddle-focus to an unstable focus and a saddle limit cycle S + osc appears. -If β 0R = 0, for (λ, ) ∈ C -HB = C + HB , the solution S 0 changes from a stable focus to an unstable focus and two stable limit cycles S + osc and Sosc appear.  (2). If β 0R > 0, a stable limit cycle emerges from C + HB , whereas a saddle limit cycle emerges from C -HB . The case β 0R < 0 is analogous just reversing ± by ∓. For the special case β 0R = 0, two stable limit cycles emerge at the coincident curves C + HB and C -HB . For these plots, we assume β 0R > α 0R > 0 In the next section we analyse the oscillatory solutions S ± osc that arise from the bifurcation curves C ± HB of system (3).

Truncated normal form in polar coordinates
To perform the analysis of the oscillatory solutions S ± osc , we express the normal form in (3) in polar coordinates, that is, we write z n = r n e iϕ n with r n > 0 and ϕ n ∈ T: r 1 = r 1 λ + α 01R r 2 1 + f r (r 1 , r 2 , Δϕ), r 2 = r 2 λ + α 01R r 2 2 + f r (r 2 , r 1 , -Δϕ), r 1φ1 = r 1 ω + α 01I r 2 1 + f ϕ (r 1 , r 2 , Δϕ), where Δϕ = ϕ 2ϕ 1 and the subscript X = R, I in α 01 refers to its real and imaginary parts, respectively. The expression for the functions f r and f ϕ can be found in Eq. (53) in the Appendix. System (9) can be also written using the variable Δϕ: where the expression for the function f Δϕ 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 ϕ 1 , ϕ 2 . However, in our case (see system (10)) the equations for the amplitudes r 1 , r 2 depend on Δϕ = ϕ 2ϕ 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 , Δϕ) → (r 2 , r 1 , -Δϕ) and K 2 = Id. (11) This action can be diagonalised using sum and difference variables s = r 1 + r 2 , d = r 1r 2 , with s, d ∈ R + × R: in this casẽ Thus, expressing the first three equations of system (10) in the variables (s, d, Δϕ), we havė where the expressions for functions g s , g d and g Δϕ 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, Δϕ has the advantage that the linearised system about the solutions of interest becomes block diagonal.

Dynamical analysis of the reduced system in the uncoupled case ( = 0)
The general picture of the uncoupled case can be obtained straightforwardly from the original system (3) for = 0. Indeed, as we consider two identical systems having a supercritical Hopf bifurcation at λ = 0, the solutions of system (3) for λ > 0 will correspond to the Cartesian product of solutions of each 2-dimensional system. In this section we show how the solutions for = 0 are seen in the reduced system (13) so that we can explore how they evolve for = 0. System (13) for = 0 writeṡ 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 R + × R, the fixed points of the Then, as the Jacobian matrix for the two first equations of system (14) is given by it is straightforward to see that the eigenvalues of (16) for (s, d) = (0, 0) are λ (double), for (s, d) = ( -4λ α 01R , 0) are -2λ (double) and for (s, d) = ( -λ α 01R , ± -λ α 01R ) are λ and -2λ. Thus, as Fig. 3 shows, when λ = 0 the origin undergoes a bifurcation and changes from stable to unstable, while three new fixed points appear: one stable corresponding to (s, d) = ( -λ α 01R , 0) plus two unstable corresponding to (s, d) = ( -λ α 01R , ± -λ α 01R ). Now let us study the solutions of system (14) obtained from the fixed points (15) when considering the variable Δϕ. The (singular) solution corresponds to the origin S 0 in (4) of system (3), which is a focus with eigenvalues λ ± iω (double) (see Sect. 3.1).
For any value Δϕ 0 , the solution is a fixed point of system (14) with eigenvalues -2λ (double) and 0. These fixed points fill up the invariant curvē whose characteristic exponents are -2λ (double). The fixed pointsS 1 (Δϕ 0 ) and the invariant curveT 0 correspond in the original system (3) for = 0 to the periodic orbits and the 2-dimensional invariant torus T 0 respectively. Notice that the periodic orbits S 1 (ϕ 0 2 ) fill the torus T 0 . The characteristic exponents of T 0 are the eigenvalues of the fixed point (s, d) = ( -4λ α 01R , 0) of the first two equations of system (14) which are -2λ (double).
The invariant 2-torus T 0 is the product of two periodic orbits with the same period in the uncoupled case = 0. Note that 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λt , whereas the internal dynamics is just a rotation. Therefore the torus T 0 is normally hyperbolic.
The last two fixed points in (15) give rise to the following periodic orbits of system (14): whose characteristic exponents are λ and -2λ, so they are of saddle type. These solutions correspond to the periodic solutions of the original system (3) for = 0 which have characteristic exponents -2λ, λ ± iω. Therefore, they are hyperbolic periodic orbits of saddle type for λ > 0.
In conclusion, for = 0, the 4D solutions S 0 , T 0 and S 2 and 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,S 0 andT 0 , and two saddle periodic orbits,S 2 andS 3 , respectively (see Fig. 4).

Figure 4
Phase space for the unperturbed system (14) for λ > 0. There are two invariant curves,S 0 (which is unstable) andT 0 (which is stable), filled with fixed points. Moreover, there exist two saddle periodic orbitsS 2 andS 3 Solutions S 2 and S 3 are hyperbolic periodic orbits for λ > 0 and = 0. Therefore, for λ > 0 fixed and small enough, there exist periodic orbits S 2 and S 3 that are C 1 -close to the unperturbed ones.
To ensure the persistence of the torus 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
The analytic continuation when increases of the periodic orbits S 2 , S 3 and the invariant torus T provided by Lemma 1 is beyond the scope of this paper. We note that the periodic orbits S 2 and S 3 are limited only by hyperbolicity. Moreover, previous work [8] highlighted that continuation of the torus T with in Lemma 1 is only possible for = o(λ). 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 (λ, ) small, of the periodic solutions S ± osc that are born at the bifurcation curves C ± HB (see (7)). In Remark 2 we relate these periodic orbits S ± osc with the invariant torus T 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 S ± osc .

3.2.2
The oscillating solutions S ± osc in the coupled case ( > 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 mapK in (12). Notice that by denoting r 1 = r 2 = r * , the curves (r * , r * , 0) and (r * , r * , π) are invariant for system (13). Then, if we write these curves in the (s, d) coordinates the dynamics for system (13) when restricted to Ξ ± reduces tȯ where the ± sign corresponds to Δϕ = 0, π , respectively. It is straightforward to check that the equation for s in (25) has three steady solutions, namely s = 0 (which corresponds to the solutionS 0 studied before) and s ± osc given by Notice that since s ∈ R + we have discarded the negative solutions for the square root. Taking into account that α 01R < 0, solutions s ± osc in (26) are only admissible whenᾱ ± = λ + (α 0R ± β 0R ) > 0. This restriction defines the following conditions for the bifurcation: which are exactly the conditions defining the curves C ± HB in (7) corresponding to the Hopf bifurcations of the origin.
Therefore, for (λ, )-values on the right-hand side of curves C ± HB , we can define, respectively, the following fixed points of system (13): 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 S ± 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 pointsS ± osc is block diagonal where the terms c s s , c d d , c d Δϕ , c Δϕ d and c Δϕ Δϕ 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 × 1 block. Thus, the eigenvaluē μ ± 1 takes the form and therefore, the solutionsS ± osc are always stable in the s direction as they appear forᾱ ± = (α 0R ± β 0R ) + λ > 0. Therefore, the pitchfork bifurcations of the origin are supercritical (see Fig. 5).
As the solutionsS ± osc are always stable in the s direction, one has to consider the eigenvalues of the 2 × 2 block, corresponding to the transverse directions in order to study possible bifurcations of the symmetric solutionsS ± osc . The trace (Tr ± ) and the determinant (Det ± ) of the 2 × 2 block of (29) atS ± osc are given up to order two in λ, by Det ± (λ, ) = ±4 λ + (α 0R ± β 0R ) (C det + β 0R ) + 4 2 β 2 0I + β 2 0R , where So, computing the discriminant we find that the eigenvalues of the 2 × 2 block of the Jacobian matrix (29) write as where Next, we study the stability of the solutionsS ± osc given in (28) when the parameters λ, lie in the domain whereᾱ ± are defined in (27). Notice that the domain A ± corresponds to the region on the right-hand side of curves C ± 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ᾱ ± = 0, that is (λ, ) ∈ C ± HB , the eigenvalues of the Jacobian matrix (29) at the fixed pointsS ± osc are given bȳ Therefore, when the parameters (λ, ) cross the curves C ± HB from left to right, if β 0R > 0,S + osc is a stable focus-node whereasSosc is a saddle-focus with a 1-dimensional stable manifold (corresponding to the s direction which is always stable) and vice versa if β 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 + HB and C -HB given in (7), and the stability of these periodic orbits depends on the sign of β 0R .
For small andᾱ ± ≥ 0, the eigenvalues of the Jacobian matrix (29) at the fixed points S ± osc are given bȳ which are O( )-close to the ones of the uncoupled case, -2λ (double) and 0. In particular, depending on the sign of (β 0R + C det ), one fixed point is a stable node whereas the other is a saddle with a 1-dimensional unstable manifold. We remark that, for λ > 0 fixed and small enough, we know that there exists an invariant curveT corresponding to the invariant torus T obtained in Lemma 1. Since this invariant curve is provided by Fenichel theory, it will contain the invariant pointsS ± osc . Consequently, if β 0R + C det > 0,T consists of the union of the saddle pointSosc , its unstable 1-dimensional manifold and the stable nodeS + osc (and vice versa if β 0R + C det < 0) (see Fig. 6). In conclusion, for λ > 0 fixed and small enough, the invariant torus T of system (3) contains the periodic orbits S + osc and Sosc with Δϕ = 0 and Δϕ = π , respectively, whose stability depends on the sign of β 0R + C det .
Remark 2 The existence of the invariant torus T is only guaranteed for λ > 0 fixed and small enough by Lemma 1. The evolution and eventual breakdown of this torus T (or, equivalently, the invariant curveT ) 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 S ± osc (corresponding to fixed pointsS ± osc ) for (λ, ) small and no assumption on = o(λ). 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].  (13) for β 0R + C det > 0, λ > 0, and 0 < < 0 (λ) (in particular λ + (α 0R ± β 0R ) > 0). There exist two fixed pointsS ± osc , a stable node and a saddle point, respectively, which together with the unstable invariant manifold of the saddle point form the invariant curveT . Due to the coupling term, there are only two fixed points onT , whereas we had an infinite number in the unperturbed case. Notice that the dynamics on the s direction is always attracting Table 1 Values for the trace (Tr), the determinant (Det) and the discriminant (Δ) of the linearisation of system (13) at the fixed pointsS ± osc near the curves C ± HB (ᾱ ± = 0) and near to the uncoupled case (ᾱ ± ≥ 0 and small)

Bifurcation diagrams of the oscillating S ± osc solutions
In the previous sections we have shown that when is small andᾱ ± ≥ 0 there exist two critical pointsS ± osc of system (13) belonging to the curveT which disappear at two independent curves C ± HB . Therefore, the pointsS ± osc undergo several bifurcations in the domain A ± defined in (36). Table 1 shows the values of the trace Tr ± in (31), the determinant Det ± in (32) and the discriminant Δ ± in (34) of the Jacobian matrix of system (13) atS ± osc near the curves C ± HB (given by the conditionᾱ ± = 0) and forᾱ ± ≥ 0 and small. Notice that the sign of the constants β 0R and C det + β 0R is relevant to determine the local dynamics around the fixed points. In particular, β 0R determines which of the two solutionsS ± osc can have a null trace. For β 0R > 0, it isS + osc , whereas for β 0R < 0, it isSosc . -The sign of C det + β 0R determines which of the two solutionsS ± osc can have a null determinant. For C det + β 0R > 0, it isSosc , whereas for C det + β 0R < 0, it isS + 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 pointsS ± osc change from a stable node and a saddle point to a stable focus and a saddle-focus. Depending on the sign of β 0R and C det + β 0R , we consider three different cases: (1) β 0R > 0, C det + β 0R > 0, (2) β 0R < 0, C det + β 0R > 0, and (3) β 0R = 0, C det > 0. The cases (i) β 0R < 0, C det + β 0R < 0, (ii) β 0R > 0, C det + β 0R < 0, and (iii) β 0R = 0, C det < 0 are analogous to (1), (2) and (3), respectively, just replacingS ± osc byS ∓ osc . For each case, we study in detail the different bifurcations of the solutionsS ± osc in the (λ, ) 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.

Dynamics ofS + osc
Forᾱ + ≥ 0, λ fixed and small, the fixed pointS + osc for system (13) is a stable node contained in the invariant curveT (region B in Fig. 7), and as increases it becomes a stable focus at the curve Δ + = 0 (region A in Fig. 7). It disappears at a pitchfork bifurcation of the origin in the s-direction at C + HB . Going back to the original 4D system (3), we have that for small there exists a stable periodic orbit S + osc (which belongs to the invariant torus T ), which disappears at a Hopf bifurcation of the origin in C + HB .

Dynamics ofSosc
The fixed pointSosc changes from a saddle-focus with a 1-dimensional stable manifold near C -HB to a saddle with a 2-dimensional stable manifold for small andᾱ -> 0. Moreover, in this case the trace forSosc vanishes. Therefore, if then Tr -= 0 and Δ -< 0 andSosc undergoes a Hopf bifurcation. So, we will distinguish two cases as follows.
(1) Case β 0R < -C det + C 2 det + β 2 0I . Forᾱ + ≥ 0, λ fixed and small, the fixed pointSosc is a saddle point with a 1-dimensional unstable manifold (in the Δϕ direction) contained in the invariant curveT (region D in Fig. 8). When crossing the curve Det -= 0 (region C), the pointSosc becomes a stable node. As the coupling is increased,Sosc crosses the curve Δ -= 0 andSosc becomes a stable focus (region B). When the parameters cross the curve Tr -= 0,Sosc undergoes a Hopf bifurcationH in the d, Δϕ directions andSosc becomes a saddle focus with a 1-dimensional unstable manifold (region A). At this bifurcation there appears or disappears a periodic orbitTdepending whether the Hopf bifurcation is supercritical or subcritical. Finally, the fixed pointSosc disappears at a pitchfork bifurcation of the origin in the s-direction occurring at the curve C -HB . Going back to the original full 4D system (3), for small enough, there exists an unstable periodic orbit Sosc , belonging to the torus T , which will become stable at the curve Det -= 0. The periodic orbit undergoes a torus bifurcation and Sosc becomes unstable at the curve Tr -= 0, and a new torus Tappears or disappears depending whether the torus Figure 8 Bifurcation diagram forSosc in the case β 0R > 0, C det + β 0R > 0 and β 0R < -C det + C 2 det + β 2 0I . The fixed pointSosc appears at a supercritical pitchfork bifurcation of the origin occurring at the curve C -HB , undergoes a Hopf bifurcationH at the curve Tr -= 0 and becomes unstable at the curve Det -= 0 bifurcation is subcritical or supercritical. Finally, Sosc will disappear at a Hopf bifurcation of the origin occurring at C -HB . (2) Case β 0R > -C det + C 2 det + β 2 0I . Forᾱ + ≥ 0, λ fixed and small, the fixed pointSosc is a saddle point with a 1-dimensional unstable manifold (in the Δϕ direction) contained in the invariant curveT (region C in Fig. 9). As increases,Sosc becomes a saddle with a 2dimensional unstable manifold at the curve Det -= 0 (region B). When further increasing the coupling ,Sosc becomes a saddle-focus point at the curve Δ -= 0 (region A), which disappears at a pitchfork bifurcation of the origin in the s-direction occurring at the curve C -HB . Going back to the original full 4D system (3), for small enough, there exists an unstable periodic orbit Sosc belonging to the torus T . The periodic orbit undergoes a bifurcation at the curve Det -= 0 in which a stable manifold becomes unstable. Finally, Sosc will disappear at a Hopf bifurcation of the origin occurring at C -HB .

Regions of bistability
SinceS + osc is always stable, bistability between fixed points will appear in those regions whereSosc is also stable. As in the case β 0R > -C det + C 2 det + β 2 0I , the fixed pointSosc is never stable, it is not possible to find bistability regions. By contrast, if β 0R < -C det + C 2 det + β 2 0I , there exists a region in the (λ, ) parameter space defined as Tr -(λ, ) < 0 and Det -(λ, ) > 0, in whichSosc 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 β 0R < -C det + C 2 det + β 2 0I , the pointSosc undergoes a Hopf bifur-cationH. If the Hopf bifurcation is supercritical, thenSosc becomes unstable and a stable limit cycleTappears, generating bistability betweenS + osc andT -. 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 pointsS ± osc by the limit cycles S ± osc and the periodic orbitTby the torus T -.

Figure 9
Bifurcation diagram forSosc in the case β 0R > 0, C det + β 0R > 0 and β 0R > -C det + C 2 det + β 2 0I . The fixed pointSosc appears at a supercritical pitchfork bifurcation of the origin in the s direction occurring at the curve C -HB and undergoes a bifurcation at the curve Det -= 0

Figure 10
Phase space forS + osc in the case β 0R < 0 and C det + β 0R > 0. The fixed pointS + osc appears at a supercritical pitchfork bifurcation of the origin in the s direction occurring at the curve C + HB and undergoes a Hopf bifurcationH at the curve Tr + = 0 4.2 Case β 0R < 0 and C det + β 0R > 0 (or β 0R > 0 and C det + β 0R < 0 )

Dynamics ofS + osc
In this case the trace forS + osc vanishes (Tr + = 0). Therefore, as then Tr + = 0 and Δ + < 0 andS + osc will always undergo a Hopf bifurcationH. Forᾱ + ≥ 0, λ fixed and small, the fixed pointS + osc is a stable node (region C in Fig. 10) and becomes a stable focus when the parameters cross the curve Δ + = 0 (region B). For larger values of , the fixed pointS + osc undergoes a Hopf bifurcationH at the curve Tr + = 0 and becomes a saddle-focus point (region A). At this bifurcation there appears or disappears a limit cycleT + depending whether this Hopf bifurcation is subcritical or supercritical. For larger values of , the fixed pointS + osc disappears at a pitchfork bifurcation of the origin in the s-direction at the curve C + HB . Going back to the original 4D system (3), for small enough, there exists a stable periodic orbit S + osc . This stable periodic orbit will lose its stability across a torus bifurcation occurring at the curve Tr + = 0. At this bifurcation there appears or disappears a torus T + depending whether the torus bifurcation is subcritical or supercritical. Finally, the unsta- Figure 11 Bifurcation diagram forSosc in the case β 0R < 0 and C det + β 0R > 0. The fixed pointSosc appears at a supercritical pitchfork bifurcation of the origin in the s direction occurring at the curve C -HB and undergoes a bifurcation at Det -= 0 ble limit cycle S + osc collapses to the origin at a Hopf bifurcation occurring at the curve C + HB .

Dynamics ofSosc
Forᾱ -≥ 0, λ fixed and small, the fixed pointSosc of system (13) is a saddle point with a 1dimensional unstable manifold in the Δϕ direction contained inT (region C in Fig. 11), and as increases it becomes a stable node when crosses the curve Det -= 0 (region B). For larger values of , the fixed pointSosc becomes a stable focus at the curve Δ -= 0 (region A) and disappears at a pitchfork bifurcation of the origin in the s direction at the curve C -HB . Going back to the original 4D system (3), for small, there exists an unstable periodic orbit Sosc . This unstable periodic orbit becomes stable at the curve Det -= 0. Finally, the stable limit cycle Sosc collapses to the origin at a Hopf bifurcation occurring at the curve C -HB .

Regions of bistability
There exists a region in the (λ, )-parameter space given by Tr + (λ, ) < 0 and Det -(λ, ) > 0, (42) in which both fixed pointsS ± osc are stable. If the Hopf bifurcation is supercritical, thenS + osc becomes unstable and a stable limit cycleT + appears, generating bistability betweenSosc andT + . 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 pointsS ± osc by the limit cycles S ± osc and the periodic orbitT + by the torus T + .

4.3
The degenerated case β 0R = 0 and C det > 0 (or β 0R = 0 and C det < 0 ) In this case, the curves C ± HB coincide. Moreover, the trace in (31) is identically zero for (λ, ) ∈ C ± HB . To obtain the sign of Tr ± , we compute Tr ± when λ + α 0R → 0 + . We have So, near the C ± HB curves, both fixed pointsS ± osc are stable. Forᾱ + ≥ 0, λ fixed and small, the fixed pointS + osc is a stable node (region B in Fig. 12), and as increases it becomes a stable focus when the parameters cross the curve Δ + = 0 (region A). For larger values of , the fixed pointS + osc disappears at a pitchfork bifurcation of the origin in the s direction at the curve C + HB . Going back to the original 4D system (3), for small, there exists a stable periodic orbit S + osc , which collapses to the origin at a Hopf bifurcation occurring at the curve C + HB .

Dynamics ofSosc
Forᾱ + ≥ 0, λ fixed and small, the fixed pointSosc 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 Det -= 0 (region B). For larger values, the fixed pointSosc becomes a stable focus at Δ -= 0 (region A) which collapses at a pitchfork bifurcation of the origin in the s direction at the curve C -HB . Going back to the original 4D system (3), for small, there exists an unstable periodic orbit Sosc which changes stability at the curve Det -= 0. Finally, the stable periodic orbit Sosc collapses to the origin at a Hopf bifurcation at the curve C -HB .

Regions of bistability
In the region in the (λ, )-parameter space given by both fixed pointsS + osc andSosc are stable. We remark that the same bistability scenarios can be found in the full system (3) replacing the fixed pointsS ± osc by the limit cycles S ± osc .

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 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) = 1 1 + e -λx+θ - which has threshold θ and slope λ with the convenient property S(0) = 0. The function S has the property S (0) = λS 1 , where S 1 = e θ (1+e θ ) 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 λ c = 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 τĖ 1 = -E 1 + S(aE 1 -bI 1 ), whose dynamics will be explored in this section.
For this study, we will consider the following set of parameters: whereas λ and will be the bifurcation parameters. By considering b 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.

Case b sp < 0
We consider the case b 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 β 0R > 0, C det + β 0R > 0 and β 0R < -C det + C 2 det + β 2 0I . Therefore, this case corresponds to the one considered in Sect. 4.1. Figure 14 shows the bifurcation diagram of system (47) for b 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 S + osc will emerge from the Hopf bifurcation at C + 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 Sosc will emerge from the Hopf bifurcation at C -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 Tr -= 0 and Det -= 0 in Fig. 8).

Case b sp > 0
We consider the case b 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 β 0R < 0 and C det + β 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 sp = 0.03 Coexisting solutions at λ ≈ 3.05 and = 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 = 0.5, where a torus bifurcation (star) is on an unstable branch that gains stability at a fold of limit cycle (square) 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 Sosc will emerge from a Hopf bifurcation at C -HB , whereas an unstable in-phase (IP) solution corresponding to S + osc will emerge from the Hopf bifurcation at C + 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 bifurcationH at the Tr + = 0 line in Fig. 10), whereas the AP solution loses stability across a pitchfork bifurcation of limit cycles (corresponding to the Det -= 0 line in Fig. 11).

Case b sp = 0
We consider the case b 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 β 0R = 0 and C 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 sp = 0 obtained numerically. Notice that it matches the theoretical predictions, namely: -Both Hopf bifurcation curves C ± 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 Δϕ = 0 corresponding to S + osc and the anti-phase (AP) Δϕ = π solution corresponding to Sosc .
-For fixed and increasing the bifurcation parameter λ, the Sosc (AP) solution loses stability across a pitchfork bifurcation of limit cycles that we found for the 3D system as the line having Det -= 0 (see Fig. 13).

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 BT ≈ 0.4, as seen in Figs. 14, 15 and 16, gives rise to a region of parameter space for λ λ c where a stable AP solution coexists with a stable fixed point. The bifurcation point BT separates branches of suband supercritical Hopf bifurcations in the parameter space. As we can see, for nearby λ, 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 ± osc solutions is given by α 01R + K ± stb and is assumed to be negative. Therefore, s ± osc solutions appear for α ± = λ + (α 0R ± β 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: we clearly see that at the curve A(λ, ) = 0 the origin undergoes a pitchfork bifurcation that is supercritical or subcritical depending on the sign of B(λ, ). Consequently, the point (λ, ) satisfying A(λ, ) = 0 and B(λ, ) = 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 assuming that Kstb > 0 and for λ BT such that (λ BT , BT ) ∈ C -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 × 2 block of Jacobian matrix (29) given by where the constants are given by Eqs. Using the numerical values given in Table 2, Kstb > 0. Thus, we can estimate from the normal form that the Bautin bifurcation occurs for BT ≈ 0.42, 0.43, 0.42 for b 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 perioddoubling bifurcation for large and λ leading to richer dynamical behaviour away from the analytically-investigated uncoupling limit.

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: where the parameters 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 λ = 2 (a-d)S 1 = 3.025, we set λ = 2.6 before this bifurcation. Further, noting that the bifurcating branch emerges with period and fixing T = 1 2f , we can set τ = √ λ 2 S 2 1 (bc-ad)+λS 1 (d-a)+1 4f π 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 ( , 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.
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  (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 = 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 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.

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 (Δϕ = 0) and anti-phase (Δϕ = π ). 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 antiphase 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.

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.

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.
With this choice, all the monomials in f 2 in (58) are null. -Compute f 3 (y,ȳ) given by the expression f 3 (y,ȳ) = D y P 2 (y,ȳ)Q 2 (y,ȳ) + DȳP 2 (y,ȳ)Q 2 (y,ȳ) + P 3 (y,ȳ), thus obtaining the coefficients corresponding to the surviving monomials in (3): y i |y i | 2 , y 2 iȳ j , y j |y i | 2 , y i |y j | 2 , y 2 jȳ i , y j |y j | 2 (i = 1, 2, j = i). -Perform the change of coordinates y = Cx in system (58), where 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 α 01 all the coefficients α i , β i (i = 0, . . . , 3) in (3) are multiplied by . Therefore, to obtain the value of the coefficients, we follow the procedure described above for = 0, thus obtaining ω and α 01 , and then repeat the same procedure for small = 0, which, using that ω and α 01 are known, provides the coefficients α i , β i .