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

TL;DR

This paper analyzes the dynamics of coupled identical oscillators near a Hopf bifurcation, revealing bistability and synchronization phenomena relevant to perceptual bistability, with theoretical and numerical insights.

Contribution

It introduces a normal form approach for coupled Hopf oscillators and applies it to a perceptual bistability model, highlighting synchronization as a segregation mechanism.

Findings

Bistability between in-phase and anti-phase solutions.

Synchronization can segregate periodic inputs.

Normal form analysis matches numerical bifurcation results.

Abstract

We study the dynamics arising when two identical oscillators are coupled near a Hopf bifurcation where we assume a parameter $\epsilon$ 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 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.

Figures

Tables

Table 1. Table 1: Values for the trace ( T ​ r 𝑇 𝑟 Tr ), the determinant ( D ​ e ​ t 𝐷 𝑒 𝑡 Det ) and the discriminant ( Δ Δ \Delta ) of the linearisation of system ( 13 ) at the fixed points 𝒮 ¯ o ​ s ​ c ± subscript superscript ¯ 𝒮 plus-or-minus 𝑜 𝑠 𝑐 \bar{\mathcal{S}}^{\pm}_{osc} near the curves C H ​ B ± subscript superscript 𝐶 plus-or-minus 𝐻 𝐵 C^{\pm}_{HB} ( α ¯ ± = 0 superscript ¯ 𝛼 plus-or-minus 0 \bar{\alpha}^{\pm}=0 ) and near to the uncoupled case ( α ¯ ± ≥ 0 superscript ¯ 𝛼 plus-or-minus 0 \bar{\alpha}^{\pm}\geq 0 and ϵ italic-ϵ \epsilon small).

	${\overline{𝒮}}_{osc}^{+}$		${\overline{𝒮}}_{osc}^{-}$
	${\overline{\alpha }}^{+}\to 0^{+}$	${\overline{\alpha }}^{+}\ge 0$, $ϵ\to 0^{+}$	${\overline{\alpha }}^{-}\to 0^{+}$	${\overline{\alpha }}^{-}\ge 0$, $ϵ\to 0^{+}$
$Tr$	$-4ϵ{\beta }_{ϵ0R}$	$-2\lambda$	$4ϵ{\beta }_{ϵ0R}$	$-2\lambda$
$Det$	$4{ϵ}^2\left({\beta }_{ϵ0I}^2+{\beta }_{ϵ0R}^2\right)$	$4ϵ\lambda \left(C_{det}+{\beta }_{ϵ0R}\right)$	$4{ϵ}^2\left({\beta }_{ϵ0I}^2+{\beta }_{ϵ0R}^2\right)$	$-4ϵ\lambda \left(C_{det}+{\beta }_{ϵ0R}\right)$
$\Delta$	$-4{ϵ}^2{\beta }_{ϵ0I}^2$	${\lambda }^2$	$-4{ϵ}^2{\beta }_{ϵ0I}^2$	${\lambda }^2$

Table 2. Table 2: Coefficients of the normal form ( 3 ) for the three considered cases, namely b s ​ p = − 0.03 , 0.03 subscript 𝑏 𝑠 𝑝 0.03 0.03 b_{sp}=-0.03,0.03 and 0 0 . These coefficients have been computed using the procedure described in Appendix B .

	$b_{sp}$
	-0.03	0.03	0
${\alpha }_{01R}$	-21.94	-21.94	-21.94
${\alpha }_{01I}$	-20.94	-20.94	-20.94
${\alpha }_{ϵ0R}$	0	0	0
${\alpha }_{ϵ0I}$	0	0	0
${\alpha }_{ϵ1R}$	0	0	0
${\alpha }_{ϵ1I}$	0	0	0
${\alpha }_{ϵ2R}$	8.4	9.02	8.72
${\alpha }_{ϵ2I}$	6.34	6.8	6.57
${\alpha }_{ϵ3R}$	-24.02	-22.3	-23.2
${\alpha }_{ϵ3I}$	-46.36	-44.92	-45.46

Equations

$$\begin{array}{c}\frac{d{x}_1}{dt}=H_{\lambda }(x_1)+ϵh_{\lambda ,ϵ}(x_1;x_2),\\ \frac{d{x}_2}{dt}=H_{\lambda }(x_2)+ϵh_{\lambda ,ϵ}(x_2;x_1),\end{array}{x}_1,x_2\in {\mathbb{R}}^2ϵ,\lambda \in \mathbb{R}$$

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

$$\begin{array}{c}\frac{d{z}_1}{dt}=U_{\lambda }(z_1)+ϵF_N(z_1,z_2,ϵ)+{\mathcal{O}}_{N+1}(z_1,z_2),\\ \frac{d{z}_2}{dt}=U_{\lambda }(z_2)+ϵF_N(z_2,z_1,ϵ)+{\mathcal{O}}_{N+1}(z_2,z_1),\end{array}{z}_1,z_2\in \mathbb{C}$$

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

$$\begin{split}\frac{dz_{1}}{dt}&=z_{1}\Big{(}\lambda+i\omega+\alpha_{01}|z_{1}|^{2}\Big{)}+\epsilon\Big{[}z_{1}\Big{(}\alpha_{\epsilon 0}+\alpha_{\epsilon 1}|z_{1}|^{2}+\alpha_{\epsilon 2}|z_{2}|^{2}+\alpha_{\epsilon 3}\bar{z}_{2}z_{1}\Big{)}\\ &+z_{2}\Big{(}\beta_{\epsilon 0}+\beta_{\epsilon 1}|z_{1}|^{2}+\beta_{\epsilon 2}|z_{2}|^{2}+\beta_{\epsilon 3}\bar{z}_{1}z_{2}\Big{)}\Big{]},\\ \frac{dz_{2}}{dt}&=z_{2}\Big{(}\lambda+i\omega+\alpha_{01}|z_{2}|^{2}\Big{)}+\epsilon\Big{[}z_{2}\Big{(}\alpha_{\epsilon 0}+\alpha_{\epsilon 1}|z_{2}|^{2}+\alpha_{\epsilon 2}|z_{1}|^{2}+\alpha_{\epsilon 3}\bar{z}_{1}z_{2}\Big{)}\\ &+z_{1}\Big{(}\beta_{\epsilon 0}+\beta_{\epsilon 1}|z_{2}|^{2}+\beta_{\epsilon 2}|z_{1}|^{2}+\beta_{\epsilon 3}\bar{z}_{2}z_{1}\Big{)}\Big{]},\end{split}$$

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

$$\left(\begin{array}[]{cccc}\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{array}\right),$$

$$\begin{array}[]{rcl}\mu_{+}=\lambda+i\omega+\epsilon(\alpha_{\epsilon 0}+\beta_{\epsilon 0}),\quad\quad\mu_{-}=\lambda+i\omega+\epsilon(\alpha_{\epsilon 0}-\beta_{\epsilon 0}),\end{array}$$

$$C_{HB}^{+}$$

$$C_{HB}^{-}$$

$$\text{if}(\lambda ,ϵ)\in C_{HB}^{+}\to$$

$$\text{if}(\lambda ,ϵ)\in C_{HB}^{-}\to$$

$$\begin{split}\dot{r}_{1}&=r_{1}\Big{(}\lambda+\alpha_{01R}r^{2}_{1}\Big{)}+\epsilon f_{r}(r_{1},r_{2},\Delta\varphi),\\ \dot{r}_{2}&=r_{2}\Big{(}\lambda+\alpha_{01R}r^{2}_{2}\Big{)}+\epsilon f_{r}(r_{2},r_{1},-\Delta\varphi),\\ r_{1}\dot{\varphi}_{1}&=r_{1}\Big{(}\omega+\alpha_{01I}r^{2}_{1}\Big{)}+\epsilon f_{\varphi}(r_{1},r_{2},\Delta\varphi),\\ r_{2}\dot{\varphi}_{2}&=r_{2}\Big{(}\omega+\alpha_{01I}r^{2}_{2}\Big{)}+\epsilon f_{\varphi}(r_{2},r_{1},-\Delta\varphi),\end{split}$$

$${\dot{r}}_1$$

$${\dot{r}}_2$$

$$\dot{\Delta \phi }$$

$${\dot{\phi }}_1$$

$$K:(r_1,r_2,\Delta \phi )\to (r_2,r_1,-\Delta \phi )\text{and}{K}^2=Id.$$

$$\tilde{K}(s,d,\Delta \phi )\to (s,-d,-\Delta \phi ).$$

$$\begin{split}\dot{s}&=s\Big{(}\lambda+\frac{\alpha_{01R}}{4}(s^{2}+3d^{2})\Big{)}+\epsilon g_{s}(s,d,\Delta\varphi),\\ \dot{d}&=d\Big{(}\lambda+\frac{\alpha_{01R}}{4}(d^{2}+3s^{2})\Big{)}+\epsilon g_{d}(s,d,\Delta\varphi),\\ \dot{\Delta\varphi}&=-\alpha_{01I}sd+\epsilon g_{\Delta\varphi}(s,d,\Delta\varphi),\end{split}$$

$$\begin{split}\dot{s}&=s\Big{(}\lambda+\frac{\alpha_{01R}}{4}(s^{2}+3d^{2})\Big{)},\\ \dot{d}&=d\Big{(}\lambda+\frac{\alpha_{01R}}{4}(d^{2}+3s^{2})\Big{)},\\ \dot{\Delta\varphi}&=-\alpha_{01I}sd.\end{split}$$

$$\displaystyle\Big{(}0,0\Big{)},\quad\Big{(}\sqrt{\frac{-4\lambda}{\alpha_{01R}}},0\Big{)},\quad\Big{(}+\sqrt{\frac{-\lambda}{\alpha_{01R}}},-\sqrt{\frac{-\lambda}{\alpha_{01R}}}\Big{)},\quad\Big{(}+\sqrt{\frac{-\lambda}{\alpha_{01R}}},+\sqrt{\frac{-\lambda}{\alpha_{01R}}}\Big{)}.$$

$$\left(\begin{array}[]{cc}\lambda+\frac{3\alpha_{01R}}{4}(s^{2}+d^{2})&\frac{\alpha_{01R}}{4}6ds\\ \frac{\alpha_{01R}}{4}6ds&\lambda+\frac{3\alpha_{01R}}{4}(s^{2}+d^{2})\end{array}\right),$$

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

$$\bar{\mathcal{S}}_{1}(\Delta\varphi_{0})=\Big{\{}s=\sqrt{\frac{-4\lambda}{\alpha_{01R}}},\quad\quad d=0,\quad\quad\Delta\varphi=\Delta\varphi_{0}\Big{\}},$$

$$\bar{\mathcal{T}}_{0}=\Big{\{}s=\sqrt{\frac{-4\lambda}{\alpha_{01R}}},\quad\quad d=0,\quad\quad\Delta\varphi\in\mathbb{T}\Big{\}},$$

$$\displaystyle\mathcal{S}_{1}(\varphi_{2}^{0})=\Big{\{}z_{1}=\sqrt{\frac{-\lambda}{\alpha_{01R}}}e^{i\varphi_{1}(t)},$$

$$\Delta {\phi }_0={\phi }_2^0-{\phi }_1^0,$$

$$\mathcal{T}_{0}=\bigcup_{\varphi_{2}^{0}\in\mathbb{T}}\mathcal{S}_{1}(\varphi_{2}^{0})=\Big{\{}|z_{1}|=|z_{2}|=\sqrt{\frac{-\lambda}{\alpha_{01R}}},\quad\quad\varphi_{1},\varphi_{2}\in\mathbb{T}\Big{\}},$$

$$\displaystyle\bar{\mathcal{S}}^{2}=\Big{\{}s=d=\sqrt{\frac{-\lambda}{\alpha_{01R}}},$$

$$\displaystyle\bar{\mathcal{S}}^{3}=\Big{\{}s=-d=\sqrt{\frac{-\lambda}{\alpha_{01R}}},$$

$$\begin{array}[]{rcl}\mathcal{S}^{2}&=&\Big{\{}z_{1}=\sqrt{\frac{-\lambda}{\alpha_{01R}}}e^{i\varphi_{1}(t)},\quad z_{2}=0,\quad\varphi_{1}(t)=\varphi^{0}_{1}+\big{(}\omega-\lambda\frac{\alpha_{01I}}{\alpha_{01R}}\big{)}t,\quad t\in\mathbb{R}\Big{\}},\\ \mathcal{S}^{3}&=&\Big{\{}z_{1}=0,\quad z_{2}=\sqrt{\frac{-\lambda}{\alpha_{01R}}}e^{i\varphi_{2}(t)},\quad\varphi_{2}(t)=\varphi^{0}_{2}+\big{(}\omega-\lambda\frac{\alpha_{01I}}{\alpha_{01R}}\big{)}t,\quad t\in\mathbb{R}\Big{\}},\end{array}$$

$$\Xi^{+}=\Big{\{}(s,d,\Delta\varphi)=(s,0,0)\Big{\}},\quad\quad\quad\Xi^{-}=\Big{\{}(s,d,\Delta\varphi)=(s,0,\pi)\Big{\}},$$

$$\begin{split}\dot{s}=&\lambda s+\frac{s^{3}\alpha_{01R}}{4}+\epsilon\Big{[}(\alpha_{\epsilon 0R}\pm\beta_{\epsilon 0R})s+\\ &\frac{s^{3}}{4}(\overbrace{\alpha_{\epsilon 2R}+\alpha_{\epsilon 1R}+\beta_{\epsilon 3R}\pm(\beta_{\epsilon 2R}+\beta_{\epsilon 1R}+\alpha_{\epsilon 3R})}^{K_{stb}^{\pm}})\Big{]},\\ \dot{d}=&0,\\ \dot{\Delta\varphi}=&0,\end{split}$$

Full text

1. 1 Introduction

2. 1.1 Coupled oscillatory instabilities

3. 1.2 Oscillatory models of perceptual bistability

4. 1.3 Outline

5. 2 Two identical Hopf bifurcations with an uncoupling limit

6. 2.1 Truncated Normal Form in Complex Coordinates

7. 3 Dynamical analysis of the truncated normal form

8. 3.1 Hopf bifurcations of the origin

9. 3.2 Truncated Normal Form in Polar Coordinates

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

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

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

13. 4.1 Case $\beta_{\epsilon 0R}>0$ and $C_{det}+\beta_{\epsilon 0R}>0$ ( or $\beta_{\epsilon 0R}<0$ and $C_{det}+\beta_{\epsilon 0R}<0$ )

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

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

16. 4.1.3 Regions of bistability

17. 4.2 Case $\beta_{\epsilon 0R}<0$ and $C_{det}+\beta_{\epsilon 0R}>0$ (or $\beta_{\epsilon 0R}>0$ and $C_{det}+\beta_{\epsilon 0R}<0$ )

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

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

20. 4.2.3 Regions of bistability

21. 4.3 The degenerated case $\beta_{\epsilon 0R}=0$ and $C_{det}>0$ (or $\beta_{\epsilon 0R}=0$ and $C_{det}<0$ )

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

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

24. 4.3.3 Regions of bistability

25. 5 Wilson-Cowan models for perceptual bistability

26. 5.1 Case $b_{sp}<0$

27. 5.2 Case $b_{sp}>0$

28. 5.3 Case $b_{sp}=0$

29. 5.4 Dynamics beyond the weak coupling limit

30. 5.5 Periodically Forced Coupled Wilson-Cowan Equations

31. 6 Discussion and conclusions

32. 6.1 Implications for models of perceptual bistability and neural competition

33. 6.2 Future perspectives

34. A Appendix: Coupling terms

35. B Normal Form Computation

**List of abbreviations

**

[TABLE]

1 Introduction

The Hopf bifurcation is a generic and well-characterized 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 an 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 AshComNic2016 .

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 literature on coupled oscillators, viewed from a wide variety of theoretical view points, and from the point of view of applications, e.g. Pikovsky01 . 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 Acebron2005 , suitable for answering a lot of questions about synchronization 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 Borisyuk1995 .

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 necker:32 ; rubin:21 , binocular rivalry levelt1968binocular ; blake:89 ; blake:01 , random-dot rotating spheres wallach-oconnell:53 , motion plaids hupe2003dynamics and multistable barber-pole illusion meso2016relative . 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 fries1997synchronization , this mechanism has been widely overlooked.

Further background and motivation for the study of coupled oscillatory instabilities close to the uncoupling limit is given in Section 1.1, whilst further background and motivation for the study of oscillatory dynamics in the context of perceptual competition is given in Section 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 golubitsky2003symmetry 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 golubitsky2003symmetry .

A recent paper ashwin2016hopf 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, ashwin2016hopf considers $N$ identical and identically interacting smooth ($C^{\infty}$) vector fields on $x_{i}\in\mathbb{R}^{d}$ $(d\geq 2)$ and present 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 Aronson1990 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 emphasize 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 gavrilov1980some (see also guckenheimer2013nonlinear ); 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 kolers1964illusion ; anstis1980perception . Figure 1A shows two frames of such an apparent motion display111More complex example than the one we’re interested in: https://open-mind.net/videomaterials/kohler-motion-quartet.mp4/view, 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. 1B, this can be perceived as a single square moving from side to side (“percept 1” in Fig. 1C). However, another interpretation is possible, of distinct squares blinking on and off either side of the fixation point (“percept 2” in Fig. 1C). Watching such a display, perception switches between percept 1 and percept 2 every few seconds; see anstis1985adaptation-b ; ramachandran1983perceptual , references within and more recently gilroy2004multiplicative ; muckli2002apparent . Perceptual bistability also occurs for the so-called auditory streaming paradigm van1975temporal ; anstis1985adaptation ; Pressnitzer2006 222https://auditoryneuroscience.com/scene-analysis/streaming-alternating-tones. The stimulus consists of interleaved sequences of tones A and B, separated by a difference in tone frequency $\Delta f$, and repeating in an “ABABAB…” pattern (Fig. 1D). This can be perceived as one stream, integrated into an alternating rhythm (“percept 1” in Fig. 1E) or as two segregated streams (“percept 2” in Fig. 1E); see recent reviews moore2002factors ; snyder2017recent . There are commonalities between these visual and auditory paradigms, in percept 1 (Fig. 1C and E) the stimulus elements are linked into a single percept. 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 anstis1985adaptation-b ; in the range 5–10 Hz for the auditory stimulus van1975temporal ), 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 Laing2002 ; Wilson2003 ; wilson2007minimal , the dependence of these dynamics on stimulus parameters Laing2002 ; moreno2010alternation ; seely2011role ; rankin2015neuromechanistic , mechanisms for attention li2017attention , entrainment to slowly varying stimuli kim2006stochastic and the effects of stimulus perturbations rankin2017stimulus . Generally models are based on competition between abstract, percept-based units Wilson2003 ; Shpiro2009 ; huguet2014noise ; li2017attention , but more recently models with a feature-based representation of competition have been developed Laing2002 ; kilpatrick2013short ; rankin2014bifurcation ; rankin2015neuromechanistic . Some percept-based models have explored how rapidly alternating inputs ($>$2 Hz) can still give rise to stable perception over several seconds wilson2007minimal ; vattikuti2016canonical ; li2017attention . 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 wang2008oscillatory 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 Section 2 we use recent theoretical results in ashwin2016hopf to write the normal form of a system of two weakly coupled identical oscillators near a Hopf bifurcation. In Section 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 Section 4 we identify different dynamical regimes depending on specific coefficients of the normal form and study the bifurcation diagrams. In Section 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 Section 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:

[TABLE]

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 $\lambda$ crosses zero.

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

[TABLE]

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 ashwin2016hopf , 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

[TABLE]

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

[TABLE]

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

[TABLE]

where the constants $\alpha_{01},\alpha_{\epsilon i},\beta_{\epsilon i}\in\mathbb{C}$ with the restriction $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

[TABLE]

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

[TABLE]

and their eigenvalues are given by

[TABLE]

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 $Re(\mu_{+})=0$ and $Re(\mu_{-})=0$. These conditions define the following Hopf bifurcation curves $\mathcal{C}^{\pm}_{HB}$ in the ($\lambda$, $\epsilon$)-parameter space

[TABLE]

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

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}_{HB}$ defined in (7). Thus,

[TABLE]

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

• •

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

• •

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

• •

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

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

3.2 Truncated Normal Form in Polar Coordinates

To perform the analysis of the oscillatory solutions $\mathcal{S}_{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}$:

[TABLE]

where $\Delta\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 $\Delta\varphi$:

[TABLE]

where the expression for the function $f_{\Delta\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 gavrilov1980some (see also guckenheimer2013nonlinear ). 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 $\Delta\varphi=\varphi_{2}-\varphi_{1}$ leading to different generic dynamics than the one in gavrilov1980some , 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 phase space as

[TABLE]

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

[TABLE]

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

[TABLE]

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

The 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 Section 3.2.2, working in the variables $s,d,\Delta\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

[TABLE]

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

[TABLE]

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

[TABLE]

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

Thus, 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)=\Big{(}\sqrt{\frac{-\lambda}{\alpha_{01R}}},0\Big{)} plus two unstable corresponding to (s,d)=\Big{(}\sqrt{\frac{-\lambda}{\alpha_{01R}}},\pm\sqrt{\frac{-\lambda}{\alpha_{01R}}}\Big{)}.

Now let us study the solutions of system (14) obtained from the fixed points (15) when considering the variable $\Delta\varphi$. The (singular) solution

[TABLE]

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

For any value $\Delta\varphi_{0}$, the solution

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

Solutions $\mathcal{S}^{2}$ and $\mathcal{S}^{3}$ are hyperbolic periodic orbits for $\lambda>0$ and $\epsilon=0$. Therefore, for $\lambda>0$ fixed and $\epsilon$ 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 Fenichel71 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\leq\epsilon\leq\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 $\epsilon$ 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 ashwin2016hopf highlighted that continuation of the torus $\mathcal{T}_{\epsilon}$ with $\epsilon$ 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 afraimovich1991invariant .

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

3.2.2 The oscillating solutions $\mathcal{S}_{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 $\tilde{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

[TABLE]

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

[TABLE]

where the $\pm$ sign corresponds to $\Delta\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}_{osc}$ given by

[TABLE]

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}_{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

[TABLE]

which are exactly the conditions defining the curves $C^{\pm}_{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}_{HB}$ we can define, respectively, the following fixed points of system (13)

[TABLE]

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}_{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}_{osc}$ is block diagonal

[TABLE]

where the terms $c^{s}_{s}$, $c^{d}_{d}$, $c^{d}_{\Delta\varphi}$, $c^{\Delta\varphi}_{d}$ and $c^{\Delta\varphi}_{\Delta\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 1x1 block. Thus, the eigenvalue $\bar{\mu}^{\pm}_{1}$ takes the form

[TABLE]

and therefore, the solutions $\bar{\mathcal{S}}^{\pm}_{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).

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

[TABLE]

where

[TABLE]

So, computing the discriminant

[TABLE]

we find that the eigenvalues of the 2x2 block of the Jacobian matrix (29) write as,

[TABLE]

where

[TABLE]

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

[TABLE]

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}_{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}_{HB}$, the eigenvalues of the Jacobian matrix (29) at the fixed points $\bar{\mathcal{S}}^{\pm}_{osc}$ are given by

[TABLE]

Therefore, when the parameters $(\lambda,\epsilon$) cross the curves $C_{HB}^{\pm}$ from left to right, if $\beta_{\epsilon 0R}>0$, $\bar{\mathcal{S}}^{+}_{osc}$ is a stable focus-node whereas $\bar{\mathcal{S}}^{-}_{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 Section 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 $\beta_{\epsilon 0R}$.

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

[TABLE]

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_{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 $\epsilon$ 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}_{osc}$. Consequently, if $\beta_{\epsilon 0R}+C_{det}>0$, $\bar{\mathcal{T}}_{\epsilon}$ consists of the union of the saddle point $\bar{\mathcal{S}}^{-}_{osc}$, its unstable 1-dimensional manifold and the stable node $\bar{\mathcal{S}}^{+}_{osc}$ (and vice versa if $\beta_{\epsilon 0R}+C_{det}<0$) (see Fig. 6). In conclusion, for $\lambda>0$ fixed and $\epsilon$ small enough, the invariant torus $\mathcal{T}_{\epsilon}$ of the system (3) contains the periodic orbits, $\mathcal{S}^{+}_{osc}$ and $\mathcal{S}^{-}_{osc}$ with $\Delta\varphi=0$ and $\Delta\varphi=\pi$, respectively, whose stability depends on the sign of $\beta_{\epsilon 0R}+C_{det}$.

Remark 2

The existence of the invariant torus $\mathcal{T}_{\epsilon}$ is only guaranteed for $\lambda>0$ fixed and $\epsilon$ 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 $\epsilon$ increases is beyond the scope of this paper.

However, in Section 4, using system (13), we study the evolution and bifurcations of the periodic orbits $\mathcal{S}^{\pm}_{osc}$ (corresponding to fixed points $\bar{\mathcal{S}}^{\pm}_{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 afraimovich1991invariant .

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

In the previous Sections we have shown that when $\epsilon$ is small and $\bar{\alpha}^{\pm}\geq 0$ there exist two critical points $\bar{\mathcal{S}}^{\pm}_{osc}$ of system (13) belonging to the curve $\bar{\mathcal{T}}_{\epsilon}$ which disappear at two independent curves: $C^{\pm}_{HB}$. Therefore, the points $\bar{\mathcal{S}}^{\pm}_{osc}$ undergo several bifurcations in the domain $\mathcal{A}^{\pm}$ defined in (36). Table 1 shows the values of the trace $Tr^{\pm}$ in (31), the determinant $Det^{\pm}$ in (32) and the discriminant $\Delta^{\pm}$ in (34) of the Jacobian matrix of system (13) at $\bar{\mathcal{S}}^{\pm}_{osc}$ near the curves $C^{\pm}_{HB}$ (given by the condition $\bar{\alpha}^{\pm}=0$) and for $\bar{\alpha}^{\pm}\geq 0$ and $\epsilon$ small. Notice that the sign of the constants $\beta_{\epsilon 0R}$ and $C_{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}_{osc}$ can have a null trace. For $\beta_{\epsilon 0R}>0$, is $\bar{\mathcal{S}}^{+}_{osc}$, whereas for $\beta_{\epsilon 0R}<0$ is $\bar{\mathcal{S}}^{-}_{osc}$.

• •

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

• •

Moreover, as we increase $\epsilon$ 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}_{osc}$ change from a stable node and a saddle point to a stable focus and a saddle-focus.

Depending on the sign of $\beta_{\epsilon 0R}$ and $C_{det}+\beta_{\epsilon 0R}$ we consider three different cases: (1) $\beta_{\epsilon 0R}>0$, $C_{det}+\beta_{\epsilon 0R}>0$, (2) $\beta_{\epsilon 0R}<0$, $C_{det}+\beta_{\epsilon 0R}>0$, and (3) $\beta_{\epsilon 0R}=0$, $C_{det}>0$. The cases (i) $\beta_{\epsilon 0R}<0$, $C_{det}+\beta_{\epsilon 0R}<0$, (ii) $\beta_{\epsilon 0R}>0$, $C_{det}+\beta_{\epsilon 0R}<0$, and (iii) $\beta_{\epsilon 0R}=0$, $C_{det}<0$ are analogous to (1), (2) and (3), respectively, just replacing $\bar{\mathcal{S}}^{\pm}_{osc}$ by $\bar{\mathcal{S}}^{\mp}_{osc}$. For each case, we study in detail the different bifurcations of the solutions $\bar{\mathcal{S}}^{\pm}_{osc}$ in the $(\lambda,\epsilon)$ parameter space, we link 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_{det}+\beta_{\epsilon 0R}>0$

( or $\beta_{\epsilon 0R}<0$ and $C_{det}+\beta_{\epsilon 0R}<0$ )

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

For $\bar{\alpha}^{+}\geq 0$, $\lambda$ fixed and $\epsilon$ small, the fixed point $\bar{\mathcal{S}}^{+}_{osc}$ for system (13) is a stable node contained in the invariant curve $\bar{\mathcal{T}}_{\epsilon}$ (region B in Fig. 7), and as $\epsilon$ increases it becomes a stable focus at the curve $\Delta^{+}=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 $\epsilon$ small there exists a stable periodic orbit $\mathcal{S}^{+}_{osc}$ (which belongs to the invariant torus $\mathcal{T}_{\epsilon}$), which disappears at a Hopf bifurcation of the origin in $C^{+}_{HB}$.

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

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

[TABLE]

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

So, we will distinguish two cases:

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

For $\bar{\alpha}^{+}\geq 0$, $\lambda$ fixed and $\epsilon$ small, the fixed point $\bar{\mathcal{S}}^{-}_{osc}$ is a saddle point with a $1$-dimensional unstable manifold (in the $\Delta\varphi$ direction) contained in the invariant curve $\bar{\mathcal{T}}_{\epsilon}$ (region D in Fig. 7). When crossing the curve $Det^{-}=0$ (region C), the point $\bar{\mathcal{S}}^{-}_{osc}$ becomes a stable node. As the coupling $\epsilon$ is increased, $\bar{\mathcal{S}}^{-}_{osc}$ crosses the curve $\Delta^{-}=0$ and $\bar{\mathcal{S}}^{-}_{osc}$ becomes a stable focus (region B). When the parameters cross the curve $Tr^{-}=0$, $\bar{\mathcal{S}}^{-}_{osc}$ undergoes a Hopf bifurcation $\bar{\mathcal{H}}$ in the $d,\Delta\varphi$ directions and $\bar{\mathcal{S}}^{-}_{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}}^{-}_{osc}$ 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 $\epsilon$ small enough, there exists an unstable periodic orbit $\mathcal{S}^{-}_{osc}$, belonging to the torus ${\mathcal{T}}_{\epsilon}$, which will become stable at the curve $Det^{-}=0$. The periodic orbit undergoes a Torus bifurcation and $\mathcal{S}^{-}_{osc}$ becomes unstable at the curve $Tr^{-}=0$ and a new torus $\mathcal{T}^{-}$ appears or disappears depending whether the Torus bifurcation is subcritical or supercritical. Finally, $\mathcal{S}^{-}_{osc}$ will disappear at a Hopf bifurcation of the origin occurring at $C^{-}_{HB}$.

2) Case $\beta_{\epsilon 0R}>-C_{det}+\sqrt{C^{2}_{det}+\beta^{2}_{\epsilon 0I}}$

For $\bar{\alpha}^{+}\geq 0$, $\lambda$ fixed and $\epsilon$ small, the fixed point $\bar{\mathcal{S}}^{-}_{osc}$ is a saddle point with a $1$-dimensional unstable manifold (in the $\Delta\varphi$ direction) contained in the invariant curve $\bar{\mathcal{T}}_{\epsilon}$ (region C in Fig. 9). As $\epsilon$ increases, $\bar{\mathcal{S}}^{-}_{osc}$ becomes a saddle with a $2$-dimensional unstable manifold at the curve $Det^{-}=0$ (region B). When further increasing the coupling $\epsilon$, $\bar{\mathcal{S}}^{-}_{osc}$ becomes a saddle-focus point at the curve $\Delta^{-}=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 $\epsilon$ small enough, there exists an unstable periodic orbit $\mathcal{S}^{-}_{osc}$ belonging to the torus ${\mathcal{T}}_{\epsilon}$. The periodic orbit undergoes a bifurcation at the curve $Det^{-}=0$ in which a stable manifold becomes unstable. Finally, $\mathcal{S}^{-}_{osc}$ will disappear at a Hopf bifurcation of the origin occurring at $C^{-}_{HB}$.

4.1.3 Regions of bistability

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

[TABLE]

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

Moreover, in the case $\beta_{\epsilon 0R}<-C_{det}+\sqrt{C^{2}_{det}+\beta^{2}_{\epsilon 0I}}$, the point $\bar{\mathcal{S}}^{-}_{osc}$ undergoes a Hopf bifurcation $\bar{\mathcal{H}}$. If the Hopf bifurcation is supercritical, then $\bar{\mathcal{S}}^{-}_{osc}$ becomes unstable and a stable limit cycle $\bar{\mathcal{T}}^{-}$ appears, generating bistability between $\bar{\mathcal{S}}^{+}_{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}_{osc}$ by the limit cycles $\mathcal{S}^{\pm}_{osc}$ and the periodic orbit $\bar{\mathcal{T}}^{-}$ by the torus $\mathcal{T}^{-}$.

4.2 Case $\beta_{\epsilon 0R}<0$ and $C_{det}+\beta_{\epsilon 0R}>0$

(or $\beta_{\epsilon 0R}>0$ and $C_{det}+\beta_{\epsilon 0R}<0$ )

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

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

[TABLE]

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

For $\bar{\alpha}^{+}\geq 0$, $\lambda$ fixed and $\epsilon$ small, the fixed point $\bar{\mathcal{S}}^{+}_{osc}$ is a stable node (region C in Fig. 10) and becomes a stable focus when the parameters cross the curve $\Delta^{+}=0$ (region B). For larger values of $\epsilon$, the fixed point $\bar{\mathcal{S}}^{+}_{osc}$ undergoes a Hopf bifurcation $\bar{\mathcal{H}}$ at the curve $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 $\epsilon$, the fixed point $\bar{\mathcal{S}}^{+}_{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 $\epsilon$ small enough there exists a stable periodic orbit $\mathcal{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 $\mathcal{T}^{+}$ depending whether the torus bifurcation is subcritical or supercritical. Finally the unstable limit cycle $\mathcal{S}^{+}_{osc}$ collapses to the origin at a Hopf bifurcation occurring at the curve $C^{+}_{HB}$.

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

For $\bar{\alpha}^{-}\geq 0$, $\lambda$ fixed and $\epsilon$ small, the fixed point $\bar{\mathcal{S}}^{-}_{osc}$ of system (13) is a saddle point with a $1$-dimensional unstable manifold in the $\Delta\varphi$ direction contained in $\bar{\mathcal{T}}_{\epsilon}$ (region C in Fig. 11), and as $\epsilon$ increases it becomes a stable node when $\epsilon$ crosses the curve $Det^{-}=0$ (region B). For larger values of $\epsilon$, the fixed point $\bar{\mathcal{S}}^{-}_{osc}$ becomes a stable focus at the curve $\Delta^{-}=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 $\epsilon$ small there exists an unstable periodic orbit $\mathcal{S}^{-}_{osc}$. This unstable periodic orbit becomes stable at the curve $Det^{-}=0$. Finally, the stable limit cycle $\mathcal{S}^{-}_{osc}$ collapses to the origin at a Hopf bifurcation occurring at the curve $C^{-}_{HB}$.

4.2.3 Regions of bistability

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

[TABLE]

in which both fixed points $\bar{\mathcal{S}}^{\pm}_{osc}$ are stable. If the Hopf bifurcation is supercritical, then $\bar{\mathcal{S}}^{+}_{osc}$ becomes unstable and a stable limit cycle $\bar{\mathcal{T}}^{+}$ appears, generating bistability between $\bar{\mathcal{S}}^{-}_{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}_{osc}$ by the limit cycles $\mathcal{S}^{\pm}_{osc}$ and the periodic orbit $\bar{\mathcal{T}}^{+}$ by the torus $\mathcal{T}^{+}$.

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

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

[TABLE]

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

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

For $\bar{\alpha}^{+}\geq 0$, $\lambda$ fixed and $\epsilon$ small, the fixed point $\bar{\mathcal{S}}^{+}_{osc}$ is a stable node (region B in Fig. 12), and as $\epsilon$ increases it becomes a stable focus when the parameters cross the curve $\Delta^{+}=0$ (region A). For larger values of $\epsilon$, the fixed point $\bar{\mathcal{S}}^{+}_{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 $\epsilon$ small there exists a stable periodic orbit $\mathcal{S}^{+}_{osc}$, which collapses to the origin at a Hopf bifurcation occurring at the curve $C^{+}_{HB}$.

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

For $\bar{\alpha}^{+}\geq 0$, $\lambda$ fixed and $\epsilon$ small, the fixed point $\bar{\mathcal{S}}^{-}_{osc}$ is a saddle point with a 1-dimensional unstable manifold (region C in Fig. 13), and as $\epsilon$ increases it becomes a stable node when the parameters cross the curve $Det^{-}=0$ (region B). For larger $\epsilon$ values the fixed point $\bar{\mathcal{S}}^{-}_{osc}$ becomes a stable focus at $\Delta^{-}=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 $\epsilon$ small there exists an unstable periodic orbit $\mathcal{S}^{-}_{osc}$ which changes stability at the curve $Det^{-}=0$. Finally, the stable periodic orbit $\mathcal{S}^{-}_{osc}$ collapses to the origin at a Hopf bifurcation at the curve $C^{-}_{HB}$.

4.3.3 Regions of bistability

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

[TABLE]

both fixed points $\bar{\mathcal{S}}^{+}_{osc}$ and $\bar{\mathcal{S}}^{-}_{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}_{osc}$ by the limit cycles $\mathcal{S}^{\pm}_{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 Wilson1972 ; rinzel1998analysis . The Wilson-Cowan oscillator (an excitatory ($E$), inhibitory ($I$) pair) considered here has dynamics described by

[TABLE]

where $\tau$ 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

[TABLE]

which has threshold $\theta$ and slope $\lambda$ with the convenient property $S(0)=0$. The function $S$ has the property $S^{\prime}(0)=\lambda S_{1}$, where $S_{1}=\frac{e^{\theta}}{(1+e^{\theta})^{2}}$, and $\lambda$ is treated as a bifurcation parameter playing the equivalent role to $\lambda$ in 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 $\epsilon$ are given by

[TABLE]

whose dynamics will be explored in this Section.

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

[TABLE]

whereas $\lambda$ and $\epsilon$ 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 Section 4 and complete the bifurcation diagrams for large values of $\lambda$ and $\epsilon$, where the normal form approximation breaks down.

5.1 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 $\beta_{\epsilon 0R}>0$, $C_{det}+\beta_{\epsilon 0R}>0$ and $\beta_{\epsilon 0R}<-C_{det}+\sqrt{C^{2}_{det}+\beta^{2}_{\epsilon 0I}}$. Therefore, this case corresponds to the one considered in Section 4.1. Fig. 14 shows the bifurcation diagram of system (47) for $b_{sp}=-0.03$ obtained numerically. The results match the theoretical predictions obtained in Section 4.1. More precisely, for a fixed $\epsilon$ value and varying the bifurcation parameter $\lambda$ we have:

• •

A stable in-phase (IP) solution corresponding to $\mathcal{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 $\mathcal{S}^{-}_{osc}$ will emerge from the Hopf bifurcation at $C^{-}_{HB}$. For a fixed $\epsilon$ and varying the bifurcation parameter AP solution gains stability across a Torus bifurcation, but when further increasing the bifurcation parameter it will loose it again across a pitchfork bifurcation (corresponding respectively to the lines $Tr^{-}=0$ and $Det^{-}=0$ in Fig. 8).

5.2 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 $\beta_{\epsilon 0R}<0$ and $C_{det}+\beta_{\epsilon 0R}>0$. Therefore, this case corresponds to the one considered in Section 4.2. Fig. 15 shows the bifurcation diagram of system (47) for $b_{sp}=0.03$ obtained numerically. The results match the theoretical predictions in Section 4.2. More precisely, for a fixed $\epsilon$ value and varying the bifurcation parameter $\lambda$ we have:

• •

A stable anti-phase (AP) solution corresponding to $\mathcal{S}^{-}_{osc}$ will emerge from a Hopf bifurcation at $C^{-}_{HB}$ whereas an unstable in-phase (IP) solution corresponding to $\mathcal{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 bifurcation $\bar{\mathcal{H}}$ at the $Tr^{+}=0$ line in Fig. 10) whereas the AP solution looses stability across a pitchfork bifurcation of limit cycles (corresponding to the $Det^{-}=0$ line in Fig. 11).

5.3 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 $\beta_{\epsilon 0R}=0$ and $C_{det}>0$. Therefore, this case corresponds to the “degenerated case” discussed in Section 4.3. Fig. 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^{\pm}_{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 $\Delta\varphi=0$ corresponding to $\mathcal{S}^{+}_{osc}$ and the anti-phase (AP) $\Delta\varphi=\pi$ solution corresponding to $\mathcal{S}^{-}_{osc}$.

• •

For $\epsilon$ fixed and increasing the bifurcation parameter $\lambda$, the $\mathcal{S}^{-}_{osc}$ (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).

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_{BT}\approx 0.4$ as seen in Figures 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_{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 Sections 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 the formula (26) for the $s^{\pm}_{osc}$ solutions, is given by $\alpha_{01R}+\epsilon K_{stb}^{\pm}$ and is assumed to be negative. Therefore, $s^{\pm}_{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

[TABLE]

we clearly see that at the curve $A(\lambda,\epsilon)=0$, the origin undergoes a pitchfork bifurcation that it 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 $\epsilon$ and $\lambda$), we can estimate that a Bautin bifurcation occurs for

[TABLE]

assuming that $K^{-}_{stb}>0$ and for $\lambda_{BT}$ such that $(\lambda_{BT},\epsilon_{BT})\in 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 the system (13) for nearby values of the parameters. More precisely, if we consider the exact expression of the determinant of the 2x2 block of Jacobian Matrix (29) given by:

[TABLE]

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

[TABLE]

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

Using the numerical values given in Table 2, $K_{stb}^{-}>0$. Thus, we can estimate from the normal form that the Bautin bifurcation occurs for $\epsilon_{BT}\approx 0.42,0.43,0.42$ for $b_{sp}=-0.03,0.03,0$, respectively, which matches the results obtained numerically (see Figures 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 $\epsilon$ and $\lambda$ 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 section 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:

[TABLE]

where the parameters $\mathcal{P}$ (with the exception of $\tau$) 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 Laing2002 ; Wilson2003 ; Shpiro2009 ; li2017attention ) 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

[TABLE]

and fixing $T=\frac{1}{2f}$ we can set $\tau=\frac{\sqrt{\lambda^{2}S_{1}^{2}(bc-ad)+\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 $\epsilon$ approaching the uncoupling limit. For large $\epsilon$ 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 close to the uncoupling limit, IP+AP region in Figure 17A). For larger forcing amplitude the intrinsic dynamics is overwhelmed and the forcing modulates a symmetrical fixed point (HA region in Figure 17A). 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 Section 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 of in ashwin2016hopf , 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 ($\Delta\varphi=0$) and anti-phase ($\Delta\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 Figure 15). Furthermore, we find 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 rankin2015neuromechanistic where a population pooling inputs from an intermediate feature location was 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 gavrilov1980some . Finally, one can consider the periodically forced system. Periodic forcing of the oscillators considered here (e.g. kim2015signal 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 remains 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 large2010canonical ; large2015neural . Such a framework could be extensible to the study of perceptual bistability, relying in the dynamic mechanisms proposed here in the simple case of only two coupled oscillators.

Declarations

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of data and material

Not applicable

Acknowledgements

This work has been partially funded by the grants MINECO-FEDER MTM2015-65715-P, MDM-2014-0445, PGC2018-098676-B-100 AEI/FEDER/UE, the Catalan grant 2017SGR1049, (GH, AP, TS), the MINECO-FEDER-UE MTM-2015-71509-C2-2-R (GH), and the Russian Scientific Foundation Grant 14-41-00044 (TS). GH acknowledges the RyC project RYC-2014-15866. TS is supported by the Catalan Institution for research and advanced studies via an ICREA academia price 2018. AP acknowledges the FPI Grant from project MINECO-FEDER-UE MTM2012-31714. We thank T. Lázaro for providing us valuable references to compute the normal form coefficients. We also acknowledge the use of the UPC Dynamical Systems group’s cluster for research computing.333See: https://dynamicalsystems.upc.edu/en/computing. PA and JR acknowledge the financial support of the EPSRC Centre for Predictive Modelling in Healthcare, via grant EP/N014391/1. JR acknowledges support from an EPSRC New Investigator Award (EP/R03124X/1).

Competing Interests

The authors declare they have no competing interests.

Authors’ Contributions

PA, JR, AP formulated the problem. All authors were involved in the theoretical analysis, discussion of the results and writing the manuscript. AP and JR performed numerical simulations. All authors have read and approved the final version.

Appendix A Appendix: Coupling terms

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

Appendix B 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

[TABLE]

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

[TABLE]

where $A$ is the linearized sytem (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:

[TABLE]

where $z=(z_{1},z_{2})$, $\bar{z}=(\bar{z}_{1},\bar{z}_{2})$ $\in$ $\mathbb{C}^{2}$, $A=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 $\bar{z}$ is the complex conjugate of $z$, we will just consider the first equation in (59). 2. 2.

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

[TABLE]

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

[TABLE]

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

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

[TABLE]

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}\quad(i=1,2\quad j\neq i)$. 4. 4.

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

[TABLE]

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 $\omega$ and $\alpha_{01}$ all the coefficients $\alpha_{\epsilon i},\beta_{\epsilon i}$ ($i=0,...,3$) in (3) are multiplied by $\epsilon$. Therefore, to obtain the value of the coefficients we follow the procedure described above for $\epsilon=0$, thus obtaining $\omega$ abd $\alpha_{01}$, and then repeat the same procedure for a small $\epsilon\neq 0$, which, using that $\omega$ and $\alpha_{01}$ are known, provides the coefficients $\alpha_{\epsilon i},\beta_{\epsilon i}$.