Phase-locked states in oscillating neural networks and their role in neural communication

TL;DR

This paper models neural oscillations to understand phase locking and its role in neural communication, revealing optimal timing for effective signaling and mechanisms for switching communication modes without structural changes.

Contribution

It introduces a population rate model analyzing phase locking and bistability, providing insights into neural communication mechanisms proposed by the CTC theory.

Findings

Optimal phase locking naturally emerges in the model.

Bistable regions enable switching between communication regimes.

Input gain depends critically on timing relative to inhibitory activity.

Abstract

The theory of communication through coherence (CTC) proposes that brain oscillations reflect changes in the excitability of neurons, and therefore the successful communication between two oscillating neural populations depends not only on the strength of the signal emitted but also on the relative phases between them. More precisely, effective communication occurs when the emitting and receiving populations are properly phase locked so the inputs sent by the emitting population arrive at the phases of maximal excitability of the receiving population. To study this setting, we consider a population rate model consisting of excitatory and inhibitory cells modelling the receiving population, and we perturb it with a time-dependent periodic function modelling the input from the emitting population. We consider the stroboscopic map for this system and compute numerically the fixed and periodic points of this map and their bifurcations as the amplitude and the frequency of the perturbation are varied. From the bifurcation diagram, we identify the phase-locked states as well as different regions of bistability. We explore carefully the dynamics emphasizing its implications for the CTC theory. In particular, we study how the input gain depends on the timing between the input and the inhibitory action of the receiving population. Our results show that naturally an optimal phase locking for CTC emerges, and provide a mechanism by which the receiving population can implement selective communication. Moreover, the presence of bistable regions, suggests a mechanism by which different communication regimes between brain areas can be established without changing the structure of the network

Equations

$$\begin{array}{rl}\dot{r_e}& =-r_e+S_e(c_1{r}_e-c_2{r}_i+P+Ap(t)),\\ \dot{r_i}& =-r_i+S_i(c_3{r}_e-c_4{r}_i+Q),\end{array}$$

$$S_k(x)=\frac{1}{1+e^{-a_k(x-{\theta }_k)}},\text{for }k=e,i,$$

$$\begin{array}{r}\mathcal{P}=\{c_1=13,c_2=12,a_e=1.3,{\theta }_e=4,c_3=6,c_4=3,a_i=2,{\theta }_i=1.5\},\end{array}$$

$$p(t)=1+\cos \left(\frac{2\pi t}{T^{\prime }}\right).$$

$$F_A:{\mathbb{R}}^2$$

$$x$$

$$[F_A(x){]}^q={\varphi }_A(t_0+q{T}^{\prime };t_0,x)={\varphi }_A(t_0+pT;t_0,x)=x.$$

$$\Delta \theta =\frac{t_{inh}-t_{pert}}{T^{\prime }},\Delta \theta \in [-0.5,0.5),$$

$$\Delta \alpha =\frac{{\alpha }_A}{{\alpha }_0}.$$

$$\Delta {\theta }_1=\frac{t_{inh}^{(1)}-t_{pert}^{(1)}}{T^{\prime }},\Delta {\alpha }_1=\frac{{\alpha }_A^{(1)}}{{\alpha }_0},$$

$$\Delta {\theta }_2=\frac{t_{inh}^{(2)}-t_{pert}^{(2)}}{T^{\prime }},\Delta {\alpha }_2=\frac{{\alpha }_A^{(2)}}{{\alpha }_0},$$

$${\tilde{F}}_A:$$

$$G(x,T^{\prime })=\{\begin{array}{l}{\tilde{F}}^x(x,T^{\prime })-x=0\\ {\Phi }_{BIF}(x,T^{\prime })=0.\end{array}$$

$${\Phi }_{SN}(x,T^{\prime })$$

$${\Phi }_{PD}(x,T^{\prime })$$

$${\Phi }_{NS}(x,T^{\prime })$$

$$\stackrel{DG(z)}{\overbracket{\left[\begin{array}{ccc}{\partial }_{x_1}{\tilde{F}}^{x_1}(z)-1& {\partial }_{x_2}{\tilde{F}}^{x_1}(z)& {\partial }_{T^{\prime }}{\tilde{F}}^{x_1}(z)\\ {\partial }_{x_1}{\tilde{F}}^{x_2}(z)& {\partial }_{x_2}{\tilde{F}}^{x_2}(z)-1& {\partial }_{T^{\prime }}{\tilde{F}}^{x_2}(z)\\ {\partial }_{x_1}{\Phi }_{BIF}(z)& {\partial }_{x_2}{\Phi }_{BIF}(z)& {\partial }_{T^{\prime }}{\Phi }_{BIF}(z)\end{array}\right]}}\stackrel{\Delta z}{\overbracket{\left[\begin{array}{c}\Delta x_1\\ \Delta x_2\\ \Delta T^{\prime }\end{array}\right]}}=-\stackrel{E}{\overbracket{\left[\begin{array}{c}{\tilde{F}}^{x_1}(z)-x_1\\ {\tilde{F}}^{x_2}(z)-x_2\\ {\Phi }_{BIF}(z)\end{array}\right]}}$$

$${\partial }_{x_j}{\tilde{F}}^{x_k}(z)$$

$${\partial }_{T^{\prime }}{\tilde{F}}^{x_k}(z)$$

$$=X^{x_k}(\tilde{\varphi }(T^{\prime };z))+\frac{\partial {\tilde{\varphi }}^{x_k}}{\partial T^{\prime }}(t;x,T^{\prime }){|}_{t=T^{\prime }}$$

$$\frac{d}{dt}{D}_z\tilde{\varphi }(t;z)=A(t)D_z\tilde{\varphi }(t;z)$$

$$A(t)={\left[\begin{array}{ccc}{\partial }_{x_1}{X}^{x_1}& {\partial }_{x_2}{X}^{x_1}& {\partial }_{T^{\prime }}{X}^{x_1}\\ {\partial }_{x_1}{X}^{x_2}& {\partial }_{x_2}{X}^{x_2}& {\partial }_{T^{\prime }}{X}^{x_2}\\ 0& 0& 0\end{array}\right]}_{|\tilde{\varphi }(t,z)}$$

$${\Phi }_{SN}=\det (D_x{\tilde{F}}^x-I{d}_2)=({\tilde{\varphi }}_{x_1}^{x_1}-1)({\tilde{\varphi }}_{x_2}^{x_2}-1)-{\tilde{\varphi }}_{x_1}^{x_2}{\tilde{\varphi }}_{x_2}^{x_1},$$

$${\partial }_{x_1}{\Phi }_{SN}=$$

$${\partial }_{x_2}{\Phi }_{SN}=$$

$${\partial }_{T^{\prime }}{\Phi }_{SN}=$$

$$\frac{d}{dt}{\tilde{\varphi }}_{z_i{z}_j}^{z_k}(t;z)=\sum _{p,q=1}^3{\partial }_{z_p{z}_q}{X}^{z_k}(\tilde{\varphi }(t;z))\frac{\partial {\tilde{\varphi }}^{z_p}(t;z)}{\partial z_i}\frac{\partial {\tilde{\varphi }}^{z_p}(t;z)}{\partial z_j}+\sum _{p=1}^3{\partial }_{z_p}{X}^{z_k}(\tilde{\varphi }(t;z))\frac{{\partial }^2{\tilde{\varphi }}^{z_p}(t;z)}{\partial z_i\partial z_j},$$

$${\dot{r}}_{e_1}$$

$${\dot{r}}_{i_1}$$

$${\dot{r}}_{e_2}$$

$${\dot{r}}_{i_2}$$

Full text

Phase-locked states in oscillating neural networks and their role in neural communication

Alberto Pérez-Cervera1, Tere M. Seara1 and Gemma Huguet1

1

Departament de Matemàtiques, Universitat Politècnica de Catalunya, Avda. Diagonal 647, 08028 Barcelona. BGSMATH

Corresponding author: Alberto Pérez-Cervera, alberto.perez@upc.edu

Keywords: Oscillatory dynamics, Wilson-Cowan equations, communication through coherence, phase-locking.

MSC2010 codes: 92B25, 65P30, 37N25

Abstract

The theory of communication through coherence (CTC) proposes that brain oscillations reflect changes in the excitability of neurons, and therefore the successful communication between two oscillating neural populations depends not only on the strength of the signal emitted but also on the relative phases between them. More precisely, effective communication occurs when the emitting and receiving populations are properly phase locked so the inputs sent by the emitting population arrive at the phases of maximal excitability of the receiving population. To study this setting, we consider a population rate model consisting of excitatory and inhibitory cells modelling the receiving population, and we perturb it with a time-dependent periodic function modelling the input from the emitting population. We consider the stroboscopic map for this system and compute numerically the fixed and periodic points of this map and their bifurcations as the amplitude and the frequency of the perturbation are varied. From the bifurcation diagram, we identify the phase-locked states as well as different regions of bistability. We explore carefully the dynamics of particular phase-locking regimes emphasizing its implications for the CTC theory. In particular, we study how the input gain depends on the timing between the input and the inhibitory action of the receiving population. Our results show that naturally an optimal phase locking for CTC emerges, and provide a mechanism by which the receiving population can implement selective communication. Moreover, the presence of bistable regions, suggests a mechanism by which different communication regimes between brain areas can be established without changing the structure of the network.

**List of abbreviations

**

[TABLE]

1 Introduction

Neural oscillations are ubiquitous in the brain. Since they were first observed in 1929 by Hans Berger [1], they have been profusely studied to unveil their link with brain function. Nowadays, they are classified in the following bands: delta (1-4 Hz), theta (4-8 Hz), alpha (8-13 Hz), beta (13-30 Hz) and gamma (30-70 Hz). Although some of these frequency bands have been associated to specific tasks or behaviours, their functional role is still not completely understood [2].

Fast brain oscillations in the gamma frequency band have been hypothesized to occur in local neural networks composed by excitatory pyramidal neurons and inhibitory interneurons (E-I networks) [3, 4]. There is an increasing number of studies which link oscillations in the gamma band frequencies with cognitive processes and communication between brain areas [5, 6, 7]. In this context, the Communication Through Coherence (CTC) Theory [8] conjectures that oscillations can account for a flexible mechanism of communication between neural populations. More precisely, oscillations generated across the interaction of excitatory and inhibitory cells cause that the excitability of the excitatory population is not the same for all the phases of the cycle due to the inhibitory action [9, 10]. Indeed, when the excitatory population receives an external input at the phase in which the inhibition is not present, the excitatory cells can respond effectively, thus promoting communication, while if the inhibition is present, the input might be ignored, thus preventing communication (see Fig. 1). Therefore, according to the CTC theory, two neuronal populations with underlying oscillatory activity communicate much effectively when they are coherent, that is, they are properly phase locked so that the output sent by the pre-synaptic (emitting) population reaches the post-synaptic (receiving) population in its peaks of excitability.

Different predictions following the CTC theory have been experimentally tested [11]. On one hand, different studies link the phase of the inhibitory receiving population with the modulation of the input gain [12]. On the other hand, different studies support that selective communication, that is, the ability of the post-synaptic group to respond to a given input and ignore the others, is implemented through selective coherence [13, 14].

Furthermore, besides the experimental studies, the CTC framework has also been studied by means of computational models, which focus on the link between gamma oscillations and stimulus selection [15, 16, 17]. Conclusions agree with the CTC hypothesis that the phase relationship which is established between a rhythmic input and the post-synaptic group turns to be optimal for the CTC scheme [18]. Most of these computational studies are based on E-I networks of spiking neurons. Nevertheless, mean field approaches are also useful to complement these results and to gain insight into the mechanisms underlying the CTC hypotheses. Indeed, they allow for a more manageable analytical treatment, while keeping the essential processes involved [19].

In this work we propose a theoretical approach to the CTC theory by means of a phenomenological description of the population activity in terms of the mean firing rate. More precisely, in order to deepen in the mechanism underlying phase-locking in a neuronal cycle having different excitability phases, we will consider the effect of an external periodic input onto a network model consisting of a single population of excitatory neurons and a single population of inhibitory neurons (E-I network). For such setting we consider the simplest canonical model describing the mean firing rates of an E-I network: the Wilson-Cowan equations [20]. The parameters of this model will be chosen so that the system shows oscillations [21]. In particular, we focus on oscillations of the Wilson-Cowan model arising from a Hopf bifurcation -we also provide a preliminary exploration of the case close to a Saddle-Node on an Invariant Circle (SNIC) bifurcation. The goal is to study the different phase-locking patterns that emerge between the oscillatory E-I network and the external periodic input for different input parameters. We remark that our setting explores the dynamics emerging from unidirectional communication. Indeed, some studies have conjectured that a given brain area has neurons receiving inputs and different neurons sending outputs [22]. Nevertheless, we point out that schemes based on bidirectional communication have also been proposed to play a role in the context of CTC theory [23, 24].

To determine and compute the phase-locked states, we consider the stroboscopic map for this system and compute numerically its fixed and periodic points and their bifurcations, as the amplitude and the frequency of the perturbation are varied. The techniques that we use to do the bifurcation analysis have no restriction neither on the amplitude nor on the frequency of the perturbation, or how close the limit cycle of the Wilson-Cowan equations is from a bifurcation. From the bifurcation diagram, we can identify the phase-locked states as well as different regions of bistability between different invariant objects. We explore carefully the dynamics on these invariant objects and we discuss the implications of these results for the CTC theory, paying attention to the phase-locking and amplitude of the response of the oscillatory neuronal population to the external input. Notably, our results provide a mechanism by which the receiving population can implement selective communication, as well as a mechanism by which different communication regimes between areas can be established (communication can be turned on and off) without changing the connectivity of the network.

The structure of the paper is as follows. In Section 2, we introduce the mathematical model by which we explore the theoretical basis of CTC. Section 3 contains the mathematical analysis of the model. More precisely, in Section 3.1, we introduce the stroboscopic map and compute the bifurcation diagram of its fixed points as the frequency and the amplitude of the perturbation are varied. In Section 3.2 we provide a complete dynamical analysis of the different phase-locking regions including the bistability regions identified in the bifurcation diagram. In Section 4 we discuss the implications for CTC theory of the different dynamical scenarios found in Section 3.2 and we finish with a discussion in Section 5. The Appendix A contains details of the numerical algorithms used to compute the bifurcation diagram, Appendix B contains a preliminary exploration of the Wilson-Cowan equations in the oscillatory regime close to a SNIC bifurcation and Appendix C includes an exploratory result for non-sinusoidal type of inputs.

2 Mathematical setting for CTC

In this Section we present the theoretical setting based on a canonical population firing rate model that implements mathematically the CTC framework. We consider the classical Wilson-Cowan model, which describes the behaviour of a coupled network of excitatory and inhibitory neurons [20], and we perturb it with a time-periodic function $p(t)$ which models the input from an external oscillating source. The perturbed Wilson-Cowan equations have the form

[TABLE]

where the variables $r_{e}$ and $r_{i}$ are the firing rate activity of the excitatory and inhibitory populations, respectively, and

[TABLE]

is the input-output function.

As we will use the Wilson-Cowan equations to model the receiving population, we need to choose parameters in (1) such that for $A=0$ they show a limit cycle. The conditions for the Wilson-Cowan equations to display oscillations have been studied in classical papers [20, 21]. Typically, the external currents $P$ and $Q$ are set as the bifurcation parameters. The reason is because they translate the nullclines of system (1) and thus determine the position and number of the critical points. For this problem we will use the following set of parameters

[TABLE]

for which, as the bifurcation diagram in Fig. 2 shows, the system (1) for $A=0$ displays a limit cycle denoted by $\Gamma_{0}$ for some $(P,Q)$ values. In particular, in this paper we choose $(P,Q)=(2.5,0)$, so the unperturbed limit cycle is near a Hopf bifurcation.

Besides the Wilson-Cowan equations, we model the external input $p(t)$ to the excitatory population by means of a positive $T^{\prime}$-periodic function. In this paper we have chosen:

[TABLE]

3 Dynamical Analysis

In this Section we will define the stroboscopic map of system (1) and study the bifurcations of its fixed points as the amplitude and the frequency of the perturbation $p(t)$ in (4) are varied. In particular, we will focus on the study of the 1:1 and 1:2 phase-locked states, as we will see in Section 4, they can serve to interpret several aspects of the CTC theory.

3.1 The stroboscopic map

To study the $T^{\prime}$-periodic system (1) we use the stroboscopic map defined by

[TABLE]

where $\phi_{A}(t;t_{0},x)$ is the solution of (1) such that $\phi_{A}(t_{0};t_{0},x)=x$. Calculations in this paper will always assume that $t_{0}$ = 0.

As it is well known, periodic orbits of system (1) are given by the fixed and periodic points of the stroboscopic map (3.1), whereas the quasi-periodic solutions correspond to its invariant curves [25]. More precisely, if $\gamma(t)=\phi_{A}(t;t_{0},x)$ is a solution of system (1) and $[F_{A}(x)]^{q}=x$, then $\phi_{A}(t_{0}+qT^{\prime};t_{0},x)=x$ and therefore $\gamma(t)$ is a periodic orbit of system (1) with period $qT^{\prime}$. Analogously, if $\gamma(t)=\phi_{A}(t;t_{0},x)$ is a periodic orbit of period $T$ of (1) with $\frac{T^{\prime}}{T}=\frac{p}{q}$, $p,q\in\mathbb{N}$, then

[TABLE]

Otherwise, if $T^{\prime}/T\in\mathbb{R}\setminus\mathbb{Q}$, then the iterates of $F_{A}$ fill densely an invariant curve denoted by $\Gamma_{A}$.

In the perturbed Wilson-Cowan model (1), the relationship (6) indicates a p:q phase locked state between the population and the perturbation. This means that the neuronal population variables $r_{e}$ and $r_{i}$ have completed p revolutions in the same time that the perturbation $p(t)$ has completed q revolutions.

3.2 Dynamics of the stroboscopic map $F_{A}$

Computing bifurcations of the fixed points of the stroboscopic map becomes relevant to identify different synchronous regimes as well as asynchronous ones. Using the techniques described in Appendix A we can compute the bifurcation diagram for the fixed points of the stroboscopic map (3.1) of system (1) as the amplitude and the frequency of the perturbation $p(t)$ in (4) are varied. As Fig. 3 shows, the fixed points of the map $F_{A}$ undergo different bifurcations, namely, saddle-node, Neimark-Sacker and period doubling bifurcations, which bound the 1:1 and 1:2 phase locking areas, and allow for a natural identification of the synchronous regimes of interest. The yellow and pink regions correspond to 1:1 and 1:2 phase locked states of system (1), respectively. We recall that they correspond to fixed points (1:1) and 2-periodic points (1:2) of the map $F_{A}$. The white regions may contain other p:q phase-locked states, as well as asynchronous states. The orange regions contain more than one stable invariant object for the map $F_{A}$.

Next, we study in detail the dynamics predicted by the bifurcation diagram in Fig. 3. In particular, we will consider different $T^{\prime}/T$ intervals and study in detail the dynamics as the amplitude $A$ of the perturbation $p(t)$ is increased. We recall that by choosing the set of parameters $\mathcal{P}$ given in (3), ($P$, $Q$) = (2.5, 0) and $A=0$, the phase space for system (1) shows the limit cycle $\Gamma_{0}$ and an unstable focus $P_{1}$ (see Fig. 2 bottom left). As both objects are normally hyperbolic we expect them to persist for weak enough amplitudes as an invariant curve $\Gamma_{A}$ and a fixed point $P_{1}$, respectively, for the corresponding stroboscopic map $F_{A}$ in (3.1). The analysis that we perform focuses on 1:1 and 1:2 phase-locked states because they occupy the largest regions of the parameter space. Furthermore, as we will see in Section 4, the dynamics emerging in the 1:1 and 1:2 phase-locking regions can be interpreted in terms of the CTC theory which motivates this study. In particular, the 1:1 phase-locking pattern allows for the study of the modulation of the input gain whereas the 1:2 phase-locking pattern accounts for selective communication.

Dynamics close to the saddle-node bifurcation at the 1:1 phase-locking region

The dynamics for values of $T^{\prime}$ such that 0.9388 $<\frac{T^{\prime}}{T}<$ 1.04, was studied in [26]. For the sake of completeness we recall here the main results, which are shown in Fig. 4. For $A$ small, an attracting invariant curve $\Gamma_{A}$ and an unstable focus $P_{1}$ inside it exist (regions $A_{1}$ and $B$). In region $A_{1}$ the invariant curve $\Gamma_{A}$ has no fixed points and once the saddle-node bifurcation (solid blue curve) is crossed (region B), there appear two fixed points on the invariant curve: a stable node $P_{2}$ and a saddle $P_{3}$, thus a SNIC bifurcation occurs, so $\Gamma_{A}$ consists of the union of the saddle $P_{3}$, its unstable invariant manifolds, and the stable node $P_{2}$. If the amplitude is increased (region C), $P_{1}$ becomes an unstable node (dashed gray curve). Furthermore, if the amplitude is increased further, $P_{1}$ coalesces with $P_{3}$ in an unstable saddle-node bifurcation (dashed blue curve), causing the disappearance of the invariant curve $\Gamma_{A}$ and the stable node $P_{2}$ remains as the unique fixed point (region D). Observe that it is possible to pass from region $A_{1}$ to region C, without passing through region B. In the region $A_{2}$ the unstable focus $P_{1}$ becomes an unstable node before crossing the saddle-node bifurcation curve (solid blue curve).

Dynamics close to the Neimark-Sacker bifurcation at the 1:1 phase-locking region

The dynamics for values of $T^{\prime}$ such that 0.51 $<\frac{T^{\prime}}{T}<$ 0.9388, was studied in [26]. For the sake of completeness we recall here the main results, which are shown in Fig. 5. For $A$ small, the attracting invariant curve $\Gamma_{A}$ has no fixed points of $F_{A}$, and an unstable focus $P_{1}$ exists inside $\Gamma_{A}$ (region A). If the amplitude $A$ is further increased, a Neimark-Sacker bifurcation occurs (green curve). At this point, the curve $\Gamma_{A}$ collapses to $P_{1}$ and disappears, while $P_{1}$ becomes a stable focus (region B).

Dynamics on the left hand side of the 1:2 phase-locking region

For values of $T^{\prime}$ such that 0.32 $<\frac{T^{\prime}}{T}<$ 0.42, the phase portrait for the map $F_{A}$ in different regions of the parameter space is shown in Fig. 6. For $A$ small, there exists a stable invariant curve and an unstable focus $P_{1}$ (region A). When the amplitude increases, $P_{1}$ becomes an unstable node when it crosses the dashed grey curve (region B). If the amplitude increases more, one finds, depending on the $T^{\prime}$ value considered, different bifurcation curves where there appear unstable fixed points for the map $F^{2}_{A}$. Next, we describe the dynamics in the zoomed region in Fig. 6 containing a sketch of such bifurcations.

For 0.355 $<\frac{T^{\prime}}{T}<$ 0.4, as $A$ increases an unstable saddle-node bifurcation for $F^{2}_{A}$ is crossed (dashed blue curve) and two saddles ($P_{3}$ and $P_{5}$) and two unstable nodes ($P_{2}$ and $P_{4}$) appear as fixed points for the map $F^{2}_{A}$ (region C). By contrast for 0.32 $<\frac{T^{\prime}}{T}<$ 0.355, one finds a period doubling bifurcation (dashed purple curve), at which there appear two unstable nodes ($P_{2}$ and $P_{4}$) and the unstable node $P_{1}$ becomes a saddle (region G). In both cases, a slight increase of the amplitude $A$ causes the unstable nodes $P_{2}$ and $P_{4}$ of $F^{2}_{A}$ to become unstable focuses at the dashed grey line (regions D and H). These focuses change their stability at a subcritical Neimark-Sacker bifurcation (green curve). Thus, an unstable invariant curve appears surrounding each of the stable focuses (regions E and I, respectively), generating a situation of bistability between the invariant curve $\Gamma_{A}$ and the fixed points $P_{2}$ and $P_{4}$. These unstable invariant curves undergo a homoclinic bifurcation (not shown) generating a unique unstable invariant curve which collides with the stable invariant curve $\Gamma_{A}$ at a saddle-node bifurcation of invariant curves (brown dashed curve). Therefore, the stable focuses $P_{2}$ and $P_{4}$ remain as the unique attractors (regions F and J). Additionally, the saddles $P_{3}$ and $P_{5}$ in region F disappear at a period doubling bifurcation (dashed purple curve) and transitioning from region F to J.

So, as the amplitude $A$ increases, by means of different bifurcation routes, that depend on $T^{\prime}/T$, the map $F^{2}_{A}$ shows the phase portrait depicted in region J, consisting of a saddle $P_{1}$ and two stable focuses $P_{2}$ and $P_{4}$. Finally, for large enough amplitudes, the stable focuses $P_{2}$ and $P_{4}$ for the map $F^{2}_{A}$ become stable nodes when crossing the dashed grey line (region K). Increasing the amplitude further, both points collapse at a period doubling bifurcation (solid purple line) where the saddle $P_{1}$ becomes a stable node (region L).

Dynamics on the right hand side of the 1:2 phase-locking region

For values of $T^{\prime}$ such that 0.42 $<\frac{T^{\prime}}{T}<$ 0.51, the phase portrait for the map (3.1) in different regions of the parameter space is shown in Fig. 7. For $A$ small, the attracting invariant curve $\Gamma_{A}$ generated from the unperturbed limit cycle $\Gamma_{0}$ has no fixed points, and an unstable focus $P_{1}$ exists inside $\Gamma_{A}$ (region A). When the amplitude is increased, a saddle node bifurcation curve of $F^{2}_{A}$ is crossed (blue curve), and there appear four fixed points on the invariant curve for the map $F_{A}^{2}$: two stable nodes ($P_{2}$ and $P_{4}$), and two saddles ($P_{3}$ and $P_{5}$) (region B). The invariant curve consists of the union of both saddles and their unstable invariant manifolds with the fixed points $P_{2}$ and $P_{4}$.

For values of $T^{\prime}$ such that 0.48 $<\frac{T^{\prime}}{T}<$ 0.51, as the amplitude is increased, $P_{2}$, $P_{3}$, $P_{4}$ and $P_{5}$ pair-collide again on a saddle-node bifurcation (blue curve) and disappear leaving an attracting invariant curve without periodic points and the unstable fixed point $P_{1}$ (region A). The amplitude of this invariant curve decreases as the amplitude $A$ increases until it reaches a supercritical Neimark-Sacker bifurcation (green curve) for the map $F_{A}$ leaving just a stable focus as the unique fixed point $P_{1}$ (Region F).

For values of $T^{\prime}$ such that 0.42 $<\frac{T^{\prime}}{T}<$ 0.48, as the amplitude $A$ increases, very close to the saddle node bifurcation curve, $P_{2}$ and $P_{4}$ become stable focuses at the grey dashed curve (region C). In this region the only stable objects are the focuses $P_{2}$ and $P_{4}$. As the amplitude is increased further, this situation is maintained until $P_{2}$ and $P_{4}$ cross again the grey dashed curve and become stable nodes, and the invariant curve is the union of the saddle points $P_{3}$ and $P_{5}$ and their unstable invariant manifolds with the fixed points $P_{2}$ and $P_{4}$ (region B’).

As the amplitude is increased further, a homoclinic bifurcation is crossed (red curve) and an invariant curve appears (region D). Therefore, we have found a region where our system presents bistability between an attracting invariant curve and the fixed points $P_{2}$ and $P_{4}$ for $F^{2}_{A}$. If the amplitude is increased further, a Neimark-Sacker bifurcation is crossed (green curve), thus the invariant curve disappears and $P_{1}$ changes stability (region E). This situation of bistability between a 2-periodic orbit and a fixed point $P_{1}$ of the map $F_{A}$ persists as the amplitude increases further until it reaches a saddle-node bifurcation (blue curve) for $F^{2}_{A}$ when the stable focus $P_{1}$ remains as the only fixed point (region F). See Appendix A for the description of the procedure followed to compute the homoclinic bifurcation curve.

Dynamics on the bottom right of the 1:1 phase-locking region

For values of $T^{\prime}$ such that 1.04 $<\frac{T^{\prime}}{T}<$ 1.125, the phase portrait for the map (3.1) in different regions of the parameter space can be seen in Fig. 8. The invariant curve $\Gamma_{A}$ (region A) evolves as the amplitude increases until crossing a saddle-node bifurcation (blue curve). At this bifurcation a stable node $P_{2}$ and a saddle $P_{3}$ are born and the invariant curve consists of the union of the saddle $P_{3}$ and its unstable invariant manifolds with the stable node $P_{2}$ (region B). Increasing the amplitude, a homoclinic bifurcation is crossed (red curve) and there appears a stable invariant curve without fixed points generating bistability between the invariant curve itself and the fixed point $P_{2}$ (region C). This invariant curve collapses at a Neimark-Sacker bifurcation (green curve), where the focus becomes stable, generating bistability between the fixed points $P_{1}$ and $P_{2}$ (region D). This situation persists until $P_{1}$ becomes a stable node at the grey dashed line (region E) which coalesces with the saddle $P_{3}$ at a saddle node bifurcation (blue curve) and disappears leaving $P_{2}$ as the unique (stable) fixed point (region F). See Appendix A for the description of the procedure followed to compute the homoclinic bifurcation curve.

Dynamics on the top right of the 1:1 phase-locking region

For values of $T^{\prime}$ such that 1.125 $<\frac{T^{\prime}}{T}<$ 1.255, the phase portrait for the map (3.1) in different regions of the parameter space can be seen in Fig. 9. As the amplitude is increased, the invariant curve $\Gamma_{A}$ (region A) collapses at a Neimark-Sacker bifurcation (green curve), where the stability of the focus $P_{1}$ changes (region B). If the amplitude is increased further, a stable node $P_{2}$ and a saddle $P_{3}$ appear at a saddle-node bifurcation (blue curve), generating a situation of bistability between the focus $P_{1}$ and the node $P_{2}$ (region C). As the amplitude $A$ increases further, $P_{1}$ becomes a stable node at the grey dashed line (region D) and coalesces with $P_{3}$ at a saddle-node bifurcation, leaving node $P_{2}$ as the unique (stable) fixed point (region E).

4 Implications for CTC Theory

In this section, we interpret the results obtained in Section 3 in terms of the CTC framework. More precisely, we explore the implications of the 1:1 phase-locked states in the modulation of the input gain, the 1:2 phase-locked states in selective communication and the bistability regions in regulating communication.

According to the CTC theory, phase locking between the emitting and receiving populations is required to establish an effective communication. Nevertheless, as an effective communication is characterized by a noticeable increase of the response of the excitatory receiving population, it turns out that the timing between the input and the inhibitory response might modulate the response of the excitatory receiving population. Indeed, inputs preceding inhibition may participate effectively in the response of the receiving population, thus increasing it. By contrast, inputs following the inhibitory action may be partially or totally silenced and, thus, have almost no noticeable effect in the receiving population. Next, in order to explore the relationship between the timing of the inhibition and the increase of the response of the excitatory population to the input we define and compute two magnitudes, $\Delta\theta$ and $\Delta\alpha$, on the phase-locking areas of interest.

In particular, we define $\Delta\theta$ to compute the phase difference between the maximum of the inhibitory population and the maximum of the perturbation, that is,

[TABLE]

where $t_{inh}$ and $t_{pert}$ are the times at which $I(t)$ and $p(t)$ of system (1) achieve a maximum inside a cycle. Notice that when $\Delta\theta$ is positive, the perturbation precedes the activation of the inhibitory population, so we expect that the excitatory receiving population is sensitive to the input. On the contrary, when $\Delta\theta$ is negative, the perturbation follows the activation of the inhibitory population, so we expect that the excitatory receiving population is less sensitive to the input due to the presence of inhibition (see Fig. 10).

In addition, $\Delta\alpha$ computes the maximum of the activity of the excitatory population $\alpha_{A}$, normalized by the maximum of the activity of the unperturbed excitatory population $\alpha_{0}$, that is,

[TABLE]

Notice that when the rate $\Delta\alpha$ is greater than one, the perturbation increases the amplitude of $r_{e}(t)$. Therefore, the larger $\Delta\alpha$ the more effective the input.

Next, we compute both magnitudes $\Delta\theta$ and $\Delta\alpha$ for each of the two phase locked regions considered. Notice that results for both regions can be interpreted differently. In the 1:1 region, as there is just one input per period, the input can precede or follow the inhibitory action, whereas in the 1:2 case, as there are two inputs per period, we expect one to precede and the other to follow the inhibitory action.

4.1 Modulation of the Input Gain (1:1 phase-locking region)

By looking at the bifurcation diagram in Fig. 3, we observe that there is a large region of 1:1 entrainment, which corresponds to the yellow region. To investigate the features of this entrainment we compute the quantities $\Delta\theta$ and $\Delta\alpha$ described above (see Fig. 11).

We observe that, in general, $0\leq\Delta\theta\leq 0.5$, indicating that inhibition typically follows the input. Predominance of positive values of $\Delta\theta$ (blue in the left panel of Fig. 11) seems to indicate that the perturbation will have a positive effect onto the activity of the excitatory population. Indeed, we observe that the activity of the excitatory population increases since $\Delta\alpha>1$ (see the right panel for Fig. 11). Nevertheless, this increase is not the same for all the points in the 1:1 phase-locking region. Notice that, as it is expected, the response of the receiving population is larger as the amplitude of the input increases. Nevertheless, for a fixed forcing amplitude $A$, the factor $\Delta\alpha$ is lower near the borders of the 1:1 phase-locking region (white and red regions), where the inhibition action precedes the input ($\Delta\theta<0$) and it can suppress totally or partially the input effect.

In conclusion, the 1:1 phase-locking pattern naturally produces a stable phase relationship that is optimal for CTC in the sense that it promotes an increase in the firing rate activity of the receiving population. Interestingly, near the boundaries of the 1:1 region this situation is reversed ($\Delta\theta<0$) so the perturbation follows the inhibitory action.

4.2 Selective communication (1:2 phase-locking region)

In the previous Section, we have shown that for the 1:1 phase-locking region, the input typically precedes the inhibitory response. This is especially interesting when studying the forcing with higher frequencies, as it is the case of the 1:2 phase-locking region. In this region, the input undergoes two cycles for one cycle of activity of the receiving population. Because of this, we expect that one of the input cycles precedes the inhibitory action whereas the other one follows it. Indeed, we can interpret the input in the 1:2 phase-locking region as two identical inputs, $I_{1}(t)$ and $I_{2}(t)$, from two different emitting neural populations, which arrive to the receiving population separated by a half-period (see Fig. 4.2), and study competition between inputs. That is, we explore whether one input phase-locks at an optimal phase so that it increases the post-synaptic response, while the other one is ignored [27, 28]. In the CTC context, this situation is known as selective communication (only one pre-synaptic population communicates effectively) [8]. We refer to input $I_{1}(t)$ as the one that precedes the main inhibitory response and $I_{2}(t)$ as the one that follows it. Thus, we expect that input $I_{1}$ produces an increase in the activity of the receiving population whereas the other one $I_{2}$, is ignored.

Similarly as the procedure followed in the 1:1 phase-locking case, we will study selective communication by computing the timing between the input and the inhibitory response $\Delta\theta$ and the rate change in the response of the excitatory receiving population $\Delta\alpha$. Since in this case the input is interpreted as the sum of two inputs (Fig. 4.2), we will compute $\Delta\theta$ and $\Delta\alpha$ for each input, provided that each input generates a response of the population. For the first input $I_{1}(t)$, as we assume that it will always precede the inhibitory response, we can always compute the following magnitudes

[TABLE]

where $t^{(1)}_{inh}$ is the time at which $r_{i}(t)$ achieves a maximum inside a cycle and $t^{(1)}_{pert}$ is the time at which $I_{1}(t)$ achieves a maximum inside the interval $0<t<T^{\prime}$. Moreover, we denote by $\alpha^{(1)}_{A}$ the value of the excitatory activity at the main maximum.

Nevertheless, for the second input $I_{2}(t)$ the situation is not so straightforward. We will consider that this input elicits a response from the receiving population if the activity of the excitatory/inhibitory population shows, apart from the main maximum, a second peak (see Fig. 12). In that case, we will also compute

[TABLE]

where $t^{(2)}_{inh}$ is the time at which $r_{i}(t)$ achieves a second local maximum inside a cycle and $t^{(2)}_{pert}$ is the time at which $I_{2}(t)$ achieves a maximum inside a cycle, so $T^{\prime}<t^{(2)}_{pert}<2T^{\prime}$. Moreover, we denote by $\alpha^{(2)}_{A}$ the value of the excitatory activity at the second local maximum.

We remark that, differently from the 1:1 case in which the perturbation can either follow or precede the inhibition, so $-0.5<\Delta\theta<0.5$, in this case both magnitudes $\Delta\theta_{1}$ and $\Delta\theta_{2}$ are defined in such a way that they always precede an inhibitory response, so $\Delta\theta_{1},\Delta\theta_{2}\in[0,1]$.

Fig 13 shows the magnitudes $\Delta\theta_{1}$ and $\Delta\alpha_{1}$ defined in (9), for the first input $I_{1}$. Observe that as $I_{1}$ was defined as always preceding the main inhibitory response, then $0<\Delta\theta_{1}<1$ (see Fig. 13 left panel), so similarly to the 1:1 phase-locking case, the effect of this input is to increase the activity of the excitatory population ($\Delta\alpha_{1}>1$) (see right panel in Fig. 13). By contrast, the second input only elicits a response of the receiving population for large values of the amplitude (see coloured region in Fig. 14 and Fig 15), and this is smaller than the one produced by the first input. Indeed, both inputs only elicit the same response just at the upper boundary of the 1:2 phase-locking region.

In conclusion, the 1:2 phase-locking pattern naturally establishes a stable phase relationship so that one of the inputs enhances the response of the excitatory neurons while preventing the second one to elicit a response (except at the upper boundary of the 1:2 phase-locking region). Notice that because of the symmetry of the problem, no input is preferred, so phase-shifts can change the input selected for effective communication.

4.3 Bistable regions

The analysis in Section 3.2 revealed the existence of different bistable regions which can be interpreted in terms of the CTC framework. Bistability suggests that, for a given input, the population may operate in different regimes depending on the initial conditions (which in fact correspond to the initial phase difference between oscillators). More interestingly, the bistability regions that we have found can generate situations in which two different synchronous regimes or co-existence of synchronous and asynchronous regimes are possible. Namely,

• •

Bistability between a 2-periodic orbit and an invariant curve without fixed or 2-periodic points for the map $F_{A}$ (panel D in Fig. 7 and panels E and I in Fig. 6).

• •

Bistability between a 2-periodic orbit and a fixed point of $F_{A}$ (panel E in Fig. 7).

• •

Bistability between a fixed point and an invariant curve without fixed points of the map $F_{A}$ (panel C in Fig. 8).

• •

Bistability between two stable fixed points of the map $F_{A}$ (panels D and E in Fig. 8 and panels C and D in Fig. 9).

Figures 16 and 17 show the main bistable regions. In particular, bistable situations between fixed points of the stroboscopic map imply a defined phase locking relationship, suggesting that there can exist different encodings of the input by the receiving population depending on the initial phase difference, as Fig 16 left illustrates. In this case, one of the solutions shows a larger variation in the activity of the E cells. We also observe bistability between 1:1 and 1:2 entrainment, see Fig.17 left, where the receiving population can either select only one input or respond to both at the price of reducing its effects.

By contrast, bistable situations between fixed points and attracting invariant curves, as illustrated in Fig 16 right and Fig. 17 right, suggest that there might exist or not coherence between the emitting and receiving neural groups depending on the initial conditions. The absence of coherence (asynchronous regimes) prevents the communication between them.

5 Discussion

In this paper we have introduced a mathematical framework based on a phenomenological description of the population activity to study some aspects of the CTC theory, namely coherence or selective communication. Our approach considers an oscillating population of excitatory and inhibitory neurons described by the Wilson-Cowan equations, which models the receiving population, submitted to an external time-periodic input, which models the effect of the emitting population. By varying the amplitude and the frequency of the external forcing, we studied the phase-locking regions between the forcing and the system, and interpreted the dynamics in terms of the CTC theory.

To do so, we considered the stroboscopic map $F_{A}$ and computed the bifurcation diagram of the fixed points in terms of the amplitude and the frequency. We have focused in the regions corresponding to 1:1 and 1:2 phase locking since they are the largest ones and have stronger implications for CTC theory. In general, our analysis revealed the existence in the parameter space of only one attracting object, either a fixed point for the map $F_{A}$ (or $F^{2}_{A}$) or an invariant curve without fixed points on it, which correspond to stable synchronous or asynchronous regimes, respectively. By performing a detailed analysis of the boundaries of these regions, we have found rich dynamics, like bistability between invariant objects.

Once we have identified the stable phase-locked states of the system, we have analysed those aspects of the dynamics that have important implications for the CTC theory. In particular, the phase relationship between the inhibition and the input and the increase in the activity of the target excitatory population due to the external input. Indeed, for an effective communication (a positive effect of the input onto the activity of the excitatory target population), the theory requires the input to precede activation of the inhibition, since the inhibition may partially or totally silence the input effect. In general, we have found that the entrainment of the postsynaptic population to the rhythmic input from the presynaptic population naturally sets up a phase relationship that is optimal for CTC, in the sense that the input precedes the inhibition, leading to an effective communication. Interestingly, we have found that near the borders of the phase-locking region, where the transition from synchronous to asynchronous dynamics occurs, the general tendency is reversed and the inhibition precedes the input. This result suggests a relationship between the loss of effective communication and the loss of phase-locking.

By repeating the analysis for the 1:2 phase-locking region, we explored a different aspect of the CTC theory, which is selective communication, that is, for a population that receives two identical inputs from different sources, how it can respond to one while ignoring the other [11]. To do so, for this region, we interpreted the input to the receiving population as the sum of two identical inputs arriving in anti-phase (see Fig 4.2), and observed that the phase automatically sets up so that one of the inputs precedes and the other one follows the inhibitory action. Our results confirm the hypothesis that the input following the inhibitory action has almost no effect onto the target population.

Moreover, we have found regions with bistability. Bistability suggests that depending on the initial conditions, the population may operate in different regimes without changing the structure or the connections of the network. Thus, bistability between synchronous and asynchronous solutions suggests that communication between brain areas cannot be predicted by the actual network but depends on the current state of the network (initial conditions), indicating that communication between neuronal populations can be switched on and off by means of possibly top-down influences [11].

Thus, we stress that thanks to the low-dimensionality of our system, we have been able to obtain a bifurcation diagram which provides a broad picture of the dynamics in a large parameter space. Thus, our results confirm previous computational results on spiking networks, while they provide new results that suggest further research. More precisely, our study corroborates the general result for network approaches [11, 18], in which the phase established between the input and the receiving network is optimal for CTC. Moreover, we have also observed that selective communication occurs and depends on the amplitude and frequency (see also [17, 16]). Moreover, because of the wide range of parameters explored, we have found new interesting results. Namely, we have been able to detect the break up of the phase-locking patterns and relate them to the phase relationship between excitation and inhibition, as well as to detect the possibility of having bistability between synchronous regimes or synchronous and asynchronous regimes, which can motivate future research.

The CTC theory proposes several hypothesis and we have not included all of them here. Other papers in the neuroscience literature regarding the mathematical implementation of the CTC theory, focus on different aspects than the ones considered here. As the only requirement of our mathematical analysis is the periodicity of the input, it can be applied to explore other hypothesis and propositions of the CTC theory. For instance, some studies related to our problem have explored the coherence (in the sense of width) of the input in selective communication [16, 17]. Instead of a sinusoidal input, we could have considered other types of functions for the input, like Gaussian-shaped with a parameter that controls the width of the input coherence (see, for instance, [16, 18]) to study this feature of CTC by means of our setting. We refer the reader to the Appendix C for preliminary results on other types of inputs. Moreover, we could have also explored the implications of inputs to the inhibitory cells (see, for instance, [29]), as they have been also proposed to be a mechanisms of phase-shifting in models of cortical networks [15]. Indeed, regarding the role of inhibition in the CTC, the Wilson-Cowan equations that we used describe oscillations across the PING (Pyramidal Interneuron Network Gamma) mechanism [30]. As the ING (Interneuron Network Gamma) mechanism can also account for gamma E-I oscillations [31, 32], we can use our setting to explore the implications of ING mechanism in the CTC [7]. Furthermore the extension of our setting to networks of two populations, can lead to the exploration of bi-directional communication [23, 24, 33, 34] instead of the unidirectional communication that we have explored in this paper.

Finally, the model can be adjusted in order to provide a more biologically-based mathematical framework. As we stated in the introduction, the framework that we introduce in this paper is an alternative to large-dimensional network approaches to the CTC. A macroscopic observable measuring the mean rate of such neuronal networks has been described by means of PDE equations such as the Fokker-Planck equation [35] or the refractory density equation [36] which have been used to describe emerging oscillations in such neuronal networks. Interestingly, recently several groups have developed exact firing rate models solving the Fokker-Planck PDE equation, thus leading to new neural mass models described by ODEs [37, 38], to which we can apply the methodology described in this paper. As the results from these models can be automatically checked on the spiking network, these new models appear as a feasible and more realistic alternative to heuristic equations such as the Wilson-Cowan equations that we use. Furthermore, independently of the rate model which is being used, the addition of synaptic equations can account for the description of realistic gamma band scenarios in which our scheme can also be applied [39].

Besides the neuroscience considerations, we would like also to emphasize several mathematical implications of our work. Perturbations of non-linear oscillators have been extensively studied in the mathematical literature [40, 41]. In particular, this paper explores the dynamics of a perturbed oscillator close to a Hopf bifurcation. Our results are based on a rigorous numerical study of a particular system (the Wilson-Cowan equations), which does not have any restrictions regarding the size of the amplitude or the forcing frequency [42]. The results obtained numerically match the theoretical predictions for perturbations of a Hopf bifurcation obtained in [43] using the normal form. We highlight here that our method does not require the system to be close to the Hopf bifurcation. Indeed, we also have applied it to the case close to the SNIC bifurcation (see Appendix B). Interestingly, although the regions of strong resonances show a different shape, we obtain the same bifurcations of fixed points. A detailed study of the implications for CTC theory of this scenario is an interesting topic for future research.

In conclusion, we have shown that a simplified setting that extracts the essence of the CTC theory allows a basic understanding of the processes involved in the generation of communication through rhythms according to the CTC theory. We expect that it would shed light into the field and open the door for future studies.

Appendix A Numerical computation of the bifurcation Diagram

In this Section we provide a brief description of the numerical procedure used to compute the bifurcations of the fixed points of the map $F_{A}$ defined in (3.1). We highlight that we have developed our own numerical software in Python to compute the bifurcation diagram instead of relying on existing software packages, thus providing more control on the calculations performed. We consider three types of bifurcations of fixed points: Saddle-Node (SN), Period Doubling (PD) and Neimark-Sacker (NS) (see [44]).

We first note that the parameters that we will vary for the map $F_{A}$ in (3.1) will be the amplitude $A$ and the period $T^{\prime}$ of the perturbation. The computation of the bifurcation curves is based on the numerical methods provided in [45]. Here we summarize the main steps of the procedure.

We first assume that $A$ is fixed and we look for a point $x\in\mathbb{R}^{2}$ and a period $T^{\prime}\in\mathbb{R}$ that satisfy that $x$ is a fixed point of the map $F_{A}$ together with a bifurcation condition $\Phi_{BIF}(x)=0$.

To set the mathematical formalism, we consider the extended system consisting of system (1) (that we will denote generically $\dot{x}=X(x,t,T^{\prime})$ where $x=(x_{1},x_{2})\in\mathbb{R}^{2}$) with the extra equation $\dot{T}^{\prime}=0$. Let us denote $\tilde{\phi}_{A}(t;t_{0},x,T^{\prime})$ the flow of the extended system and let us introduce the map

[TABLE]

Notice that we have set $t_{0}=0$ and abusing of notation we have removed the dependence on $t_{0}$ from the expression of the flow $\tilde{\phi}_{A}$. From now on we will also remove the subscript $A$ to avoid stodgy notation. The superscript $\tilde{F}^{w}$ refers to the $w$-component of the map $\tilde{F}$, where $w=x_{1},x_{2},T^{\prime}$.

Remark A.1**.**

We consider the extended system with the trivial equation $\dot{T^{\prime}}=0$ because, as we will see later, we need to know how the solutions of system (1) vary with respect to the parameter $T^{\prime}$.

Thus, to detect the bifurcation values we need to look for zeroes of the equation

[TABLE]

The conditions which must be satisfied at the bifurcation values for a Saddle-Node (SN), Period Doubling (PD) and Neimark-Sacker (NS) bifurcations are, respectively,

[TABLE]

where we denote by $D_{x}\tilde{F}^{x}$ the Jacobian matrix of the map $\tilde{F}^{x}$ restricted to the first two components.

We have implemented a Newton method to solve system (12). Suppose we have an approximate solution $z=(x,T^{\prime})$ of system (12) and we look for an improved solution $z^{*}=z+\Delta z$. Straightforward calculations show that $\Delta z$ has to satisfy:

[TABLE]

Notice that for the first two rows of the matrix $DG$ we have

[TABLE]

for $j,k=1,2$. In order to obtain the derivatives of the flow $\tilde{\phi}$ with respect to initial conditions $(x_{1},x_{2},T^{\prime})$ at time $T^{\prime}$, we need to integrate the first variational equations given by

[TABLE]

with initial condition $D_{z}\tilde{\phi}(0;z)=Id_{3}$, where

[TABLE]

Recall that these equations need to be integrated together with the flow $\tilde{\phi}(t,z)$ up to time $T^{\prime}$.

Finally, one needs to compute the terms in the third row of $DG$ in (14). The exact expression depends on the bifurcation for which we look for (see Eq. (13)). Next, we will derive the expression for the $SN$ case to illustrate the method, and the other cases are analogous. The determinant of $D_{x}\tilde{F}^{x}-Id_{2}$ writes as

[TABLE]

where $\tilde{\phi}^{x_{k}}_{x_{j}}=\partial_{x_{j}}\tilde{\phi}^{x_{k}}(T^{\prime};x,T^{\prime})$ and therefore

[TABLE]

where $\tilde{\phi}^{z_{k}}_{z_{i}z_{j}}=\frac{\partial^{2}}{\partial_{z_{j}}\partial_{z_{i}}}\tilde{\phi}^{z_{k}}(T^{\prime};z)$ for $k=1,2$ and $i,j=1,2,3$. Notice that $z_{3}=T^{\prime}$.

The computation of the terms in equations (18) requires to solve the second order variational equations, which are given in components as

[TABLE]

with initial condition $\tilde{\phi}^{z_{k}}_{z_{i}z_{j}}(0;z)=0$ and $k=1,2$, $i,j=1,2,3$.

Once a point $z^{*}$ satisfies (12) with the established tolerance (in our case $10^{-8}$), the procedure to continue the bifurcation curve consists in selecting a new value for the amplitude $A_{new}=A+\Delta_{A}$, use the computed point $z^{*}$ as initial seed, and repeat the above procedure for the new value $A_{new}$. This procedure –which is the one used in this manuscript– requires to change the sign of $\Delta_{A}$ by hand when the derivative of the bifurcation curve with respect to $A$ is zero. Nevertheless, alternative strategies as the Lagrange multipliers or the pseudo arclength method (see [45, 46], respectively) can be used so the tangent vector $v$ to the bifurcation curve evaluated at $z^{*}$ is obtained and can be used to provide $z^{*}+\epsilon v$ as initial seed for $\epsilon$ small enough.

**Homoclinic bifurcation

**

In this Section, we explain the method used to compute the homoclinic bifurcation curve in Section 3.2. The crossing of an homoclinic bifurcation implies the appearance/disappearance of an attracting invariant curve which surrounds the unstable focus $P_{1}$. Therefore, the crossing of this bifurcation implies a qualitative change in the asymptotic solutions of system (3.1) when using an initial condition near the unstable focus $P_{1}$. More precisely, in the case corresponding to “Dynamics on the right hand side of the 1:2 phase-locking region” (see Fig. 7), if an attracting invariant curve exists, an initial condition near the unstable focus $P_{1}$ will tend to the invariant curve for a large enough time of integration (see panel D in Fig. 7). Otherwise, the 2-periodic orbit $(P_{2},P_{4})$ will be the asymptotic solution (see panel B’ in Fig. 7). Similarly, in the case “Dynamics on the bottom right of the 1:1 phase-locking region” (see Fig. 8), if an attracting invariant curve exists, an initial condition near the unstable focus $P_{1}$ will tend to the invariant curve for a large enough time of integration (see panel C in Fig. 8). Otherwise, the fixed point $P_{2}$ will be the asymptotic solution (see panel B in Fig. 8). In conclusion, in both cases, the homoclinic bifurcation curve is delimited by considering different values of the fraction $\frac{T^{\prime}}{T}$, and slightly varying the amplitude while checking whether there is a qualitative change in the asymptotic solution of points near $P_{1}$.

Appendix B Bifurcation Diagram for oscillations close to a SNIC bifurcation

The unperturbed Wilson-Cowan equations (1) with the set of parameters $\mathcal{P}$ given in (3), can also have oscillations which are born from a Saddle-Node on Invariant Curve (SNIC) bifurcation (see Fig. 2). Thus, when we pick the values $(P,Q)=(1.4,-0.75)$, identically as to the Hopf case, the phase space for system (1) shows a limit cycle and an unstable focus (see Fig. 18). We have carried out a numerical exploration of the bifurcations that occur when we perturb this limit cycle close to a SNIC bifurcation. Preliminary results are shown in Figure 18 bottom. The regions with p:q phase-locking, start out of the point $(p/q,0)$ but, as the amplitude increases, they tilt towards lower values of $p/q$ as predicted by the Phase Response Curves [47]. Indeed, the PRC for limit cycles close to a SNIC are mainly positive indicating that the phase can only be advanced by an external excitatory perturbation. Therefore, they can only synchronize with external inputs with higher frequency.

Moreover, when we compare with the Hopf diagram, in both cases, the phase-locking regions are bounded by saddle-node bifurcation curves (for small values of $A$) and Neimark-Sacker and period-doubling bifurcation curves for other values of the amplitude. So, our preliminary study suggests a qualitatively similar dynamics as in the Hopf case. A thorough study is left for future work.

Appendix C Non-sinusoidal inputs

Our analysis has been carried out for sinusoidal inputs while the oscillations generated by the E-I network are not. Thus, in order to consider a more realistic approach, in this Section, we present preliminary results on the dynamical effects of unidirectional coupling between two Wilson-Cowan populations. Thus, we use the value of $r_{e}$ of one population as the input for the second one, which is periodic but non-sinusoidal. Mathematically, we consider the following system of equations:

[TABLE]

where the parameters are the same for both populations and correspond to the set $\mathcal{P}$ in (3).

When we simulate this system, the network establishes a 1:1 phase-locking and a phase relationship $\Delta\theta>0$ that increases the amplitude of oscillations $\Delta\alpha>1$ (see Fig. 19), in agreement with the results of the paper presented in Section 4.1. A complete investigation varying the amplitude and frequency of the input, or even the coherence (in the sense of sharpness) of the input (see [18]) is left for future work. Indeed, regarding the coherence of the input, preliminary results with inputs of the form $p(t)=\cos^{n}(2\pi t/T^{\prime})$ with $n=4,8$, which have higher coherence than the input (4) considered in this paper, show qualitatively the same bifurcation diagram as in Fig 3 (results not shown).

Acknowledgements

This work has been partially funded by the Spanish grants MTM2015-65715-P, MDM-2014-0445, PGC2018-098676-B-100, 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.

Competing Interests

The authors declare they have no competing interests.