Symmetry breaking in two interacting populations of quadratic integrate-and-fire neurons

TL;DR

This paper investigates the complex dynamics of two coupled populations of quadratic integrate-and-fire neurons, revealing various non-symmetric states and bifurcations through a novel reduction method validated by simulations.

Contribution

It introduces an exact reduction method for analyzing coupled neuron populations with heterogeneous neurons, uncovering diverse dynamical patterns.

Findings

Discovery of splay, antiphase, and chimera-like states.

Identification of chaotic oscillations and bistabilities.

Validation of the reduction method with numerical simulations.

Abstract

We analyze the dynamics of two coupled identical populations of quadratic integrate-and-fire neurons, which represent the canonical model for class I neurons near the spiking threshold. The populations are heterogeneous; they include both inherently spiking and excitable neurons. The coupling within and between the populations is global via synapses that take into account the finite width of synaptic pulses. Using a recently developed reduction method based on the Lorentzian ansatz, we derive a closed system of equations for the neuron's firing rates and the mean membrane potentials in both populations. The reduced equations are exact in the infinite-size limit. The bifurcation analysis of the equations reveals a rich variety of non-symmetric patterns, including a splay state, antiphase periodic oscillations, chimera-like states, also chaotic oscillations as well as bistabilities between various states. The validity of the reduced equations is confirmed by direct numerical simulations of the finite-size networks.

Figures

Equations

$${\dot{V}}_j=V_j^2+{\eta }_j+I_j^{syn}.$$

$$I_j^{syn}=-K(V_j-V_s)\frac{1}{N}\sum _{l=1}^N{s}_l.$$

$$I^{syn}=J\frac{V_{th}}{N}\sum _{l=1}^{N}H(V_l-V_{th}).$$

$${\dot{V}}_{j,k}=V_{j,k}^2+{\eta }_j+{\mathcal{I}}_k.$$

$${\mathcal{I}}_k=\left(J_{in}{S}_k+J_{ex}{S}_{1-k}\right)V_{th},$$

$$S_k=\frac{1}{N}\sum _{l=1}^{N}H(V_{l,k}-V_{th})$$

$$\frac{\partial }{\partial t}{\rho }_k=-\frac{\partial }{\partial V}\left[{\rho }_k\left\{V^2+\eta +{\mathcal{I}}_k\right\}\right],$$

$$S_k=\underset{-\infty }{\overset{+\infty }{\int }}g(\eta )\underset{-\infty }{\overset{+\infty }{\int }}{\rho }_k(V|\eta ,t)H(V-V_{th})\text{d}V\text{d}\eta .$$

$${\rho }_k(V|\eta ,t)=\frac{1}{\pi }\frac{x_k(\eta ,t)}{[V-y_k(\eta ,t){]}^2+x_k(\eta ,t{)}^2},$$

$$r_k(t)$$

$$v_k(t)$$

$${\dot{x}}_k(\eta ,t)$$

$${\dot{y}}_k(\eta ,t)$$

$$\dot{w_k}(\eta ,t)=i\left[\eta -w_k^2(\eta ,t)+{\mathcal{I}}_k\right].$$

$$g(\eta )=\frac{1}{\pi }\frac{\Delta }{(\eta -\bar{\eta }{)}^2+{\Delta }^2}$$

$$\pi r_k(t)+i{v}_k(t)=w_k(\bar{\eta }-i\Delta ,t)$$

$$S_k(t)=\frac{1}{\pi }\left[\frac{\pi }{2}-\arctan \left(\frac{V_{th}-v_k(t)}{\pi r_k(t)}\right)\right].$$

$${\dot{r}}_k$$

$${\dot{v}}_k$$

$$R$$

$$P$$

$$Q$$

$$M$$

$$\dot{Q}$$

$$\dot{M}$$

$$\left(\begin{array}[]{c}\delta\dot{R}\\ \delta\dot{P}\end{array}\right)=\boldsymbol{A}\left(\begin{array}[]{c}\delta R\\ \delta P\end{array}\right)$$

$$\boldsymbol{A}=2\left(\begin{array}[]{cc}M&Q\\ -\pi^{2}Q+S_{M}(M-V_{th})&M-S_{M}Q\end{array}\right)$$

$$S_M=\frac{V_{th}(J_{ex}-J_{in})}{2\left[{\pi }^2{Q}^2+(M-V_{th}{)}^2\right]}.$$

$$V_{j,k}=\tan ({\theta }_{j,k}/2)$$

$${\dot{\theta }}_{j,k}=(1-\cos {\theta }_{j,k})+(1+\cos {\theta }_{j,k})({\eta }_j+{\mathcal{I}}_k).$$

$${\mathcal{I}}_k=J_{in}{S}_k{V}_{th}+J_{ex}{v}_{1-k},$$

$$v_k=\frac{1}{N}\sum _{l=1}^N{V}_{l,k}$$

Full text

I Introduction

Studies of collective motion in networks of nonlocally or globally coupled oscillators or excitable elements are the focus of current research in diverse fields from physics to neuroscience. Starting from pioneering works of Winfree Winfree (1967) and Kuramoto Kuramoto (1975) various striking effects such as synchronization, collective chaos, and chimera states have been detected Pikovsky and Rosenblum (2015).

The existence of chimera states was first reported by Kuramoto and Battogtokh Kuramoto and Battogtokh (2002). They considered a ring of identical nonlocally coupled oscillators and showed that they could spontaneously split into synchronized and desynchronized subpopulations. Though the oscillators were coupled symmetrically and the system possessed a translational symmetry, for the chimera states this symmetry was broken. By now chimera states were detected and analyzed in diverse systems with different types of topology, various types of oscillators and different types of couplings. For a recent review of the subject we refer to Ref. Panaggio and Abrams (2015).

A major breakthrough in an analytical treatment of chimera states was achieved by Abrams et al Abrams et al. (2008). They considered the simplest setup that supports chimera states: a pair of oscillator populations in which each oscillator is coupled equally to all the others in its group, and with different strength to those in the other group. Such a system is symmetric with respect to exchange of populations and here chimera states represent the symmetry broken solutions with the coherent and incoherent behavior of oscillators in different populations. The authors solved the problem in the infinite-size (thermodynamic) limit by applying Ott and Antonsen ansatz Ott and Antonsen (2008). They derived a simple system of two ordinary differential equations that characterize the macroscopic dynamics of the network and obtained the exact results about the stability, dynamics, and bifurcations of chimera states. Inspired by the work Abrams et al. (2008), different authors employed the setup with two globally coupled populations to analyze chimera states in a large variety of models Laing (2009); Pazó and Montbrió (2014); González-Avella et al. (2014); Panaggio et al. (2016); Olmi et al. (2010); Bolotov et al. (2016); Shepelev et al. (2017).

Chimera states are of particular interest in neural models Olmi et al. (2010); Bolotov et al. (2016); Shepelev et al. (2017); Sakaguchi (2006); Omelchenko et al. (2013); Hizanidis et al. (2014); Omelchenko et al. (2015); Bera et al. (2016a, b); Glaze et al. (2016); Semenova et al. (2016); Isele et al. (2016). Many creatures like birds, reptiles and sea mammals sleep with only half their brain at a time Rattenborg et al. (2000). In such a unihemispheric sleep the awake side of the brain shows desynchronized electrical activity, whereas the sleeping side is highly synchronized. The authors of Ref. Abrams et al. (2008) suggested that chimera states in two coupled neural populations might serve as a model of unihemispheric sleep. An asymmetric brain activity has been also observed in human sleep apnea patients Rial et al. (2013). Thus the study of neural chimera states may have clinical relevance as well.

In this paper, we analyze two globally coupled populations of quadratic integrate-and-fire (QIF) neurons. The isolated QIF neuron is the canonical model for the class I neurons near the spiking threshold. The spiking instability in such neurons appears through a saddle-node bifurcation on an invariant curve (SNIC), in which a pair of fixed points on a closed curve coalesce to disappear, converting the curve to a periodic orbit. The peculiarities of our model are as follows. Unlike to typical models, which consider chimera states in systems of identical oscillators, here we analyze two populations of heterogeneous neurons. Each population contains both excitable and spiking neurons. The interactions between neurons are provided either by synaptic coupling or mean field potential. Our model admits an analytical treatment via a recently developed Lorentzian ansatz (LA) method Montbrió et al. (2015); Ratas and Pyragas (2016). In the thermodynamic limit, we derive a simple system of ordinary differential equations, which describe the macroscopic dynamics of the firing rates and mean membrane potentials in both neural populations. This macroscopic model enables us to perform a thorough bifurcation analysis of the system. As a result, we detect two types of chimera-like states. In one of them, the majority of neurons in one population are quenched, while in another population they spike synchronously. In the second type, the majority of neurons produce spikes in both populations, but with a different synchronization level.

The paper is organized as follows. The microscopic model of two synaptically coupled populations of QIF neurons and derivation of macroscopic model equations in the thermodynamic limit is described in Sec. II. Section III is devoted to the stability analyzes of symmetric solutions of the macroscopic model. The analysis of non-symmetric solutions and their bifurcations is presented in Sec. IV. In Sec. V we consider the case when the coupling between populations is defined by a mean field potential rather than the synaptic interaction. The paper is concluded by discussion in Sec. VI.

II The Model

Our consideration of two neural populations is based on the heterogeneous model of all-to-all synaptically coupled quadratic integrate-and-fire neurons, which are the canonical representatives for a class I neurons near the spiking threshold. One such population has been thoroughly studied in Ref. Ratas and Pyragas (2016). The membrane potential of each neuron $V_{j}$, (here ${1\leq j\leq N}$, $N$ is the size of the population) is described by the following equation Bard Ermentrout and Terman (2010):

[TABLE]

Here the constants $\eta_{j}$ specify the behavior of individual neurons. For $\eta_{j}<0$ the neuron is in an excitable regime and for $\eta_{j}>0$ it is in the spiking regime. We assume that the values of the parameters $\eta_{j}$ are distributed according to some defined density function $g(\eta)$ and that system (1) contains both excitable ($\eta_{j}<0$) and spiking neurons ($\eta_{j}>0$). Whenever the membrane potential $V_{j}$ reaches the peak value $V_{peak}$ its voltage is reset to the value $V_{reset}$. In order to treat the system (1) analytically, the peak and reset voltages are set to infinity ${V_{peak}=-V_{reset}\to\infty}$.

The term $I^{syn}_{j}$ stands for the synaptic current. We assume all-to-all homogeneous neural coupling and write this term in the form:

[TABLE]

Here $K$ is the maximal conductance of postsynaptic receptors and $V_{s}$ is the reversal potential of synapse. In the case of fast synaptic processes, the fraction of open ion channels in the neuron membrane is described by sigmoid function, $s_{l}=\{1+\exp[-\sigma(V_{l}-V_{th})]\}^{-1}$, with steepness parameter $\sigma$ and threshold potential $V_{th}$. This oft-used coupling form is called fast threshold modulation Somers and Kopell (1993). We consider the limiting version $\sigma\to 0$, when the sigmoid function transforms into the Heaviside step function, $s_{l}=H(V_{l}-V_{th})$. Moreover, to reduce the number of parameters, we consider the limit $V_{s}\to\infty$ and $K\to 0$ with the product $KV_{s}$ remaining finite. Then by defining the new parameter $J=KV_{s}/V_{th}$, we obtain the simplified expression for the synaptic current that does not depend on the index $j$, $I_{j}^{syn}=I^{syn}$,

[TABLE]

This expression is a good approximation for small excitatory synapses on a large compartment Roth and van Rossum (2009). In that case, the depolarization of the membrane is small and the difference $V_{j}-V_{s}$ is little changed during the excitatory postsynaptic potential. Note that the latter approximation is not necessary for the analytical treatment of the model. The reduced system of macroscopic equations for the neuron’s firing rate and the mean membrane potential can be derived in the thermodynamic limit without recourse to this approximation Ratas and Pyragas (2016). However, below we will use this approximation, since it simplifies the bifurcation analysis of the macroscopic equations. Moreover, as shown in Ref. Ratas and Pyragas (2016), the bifurcation diagrams obtained with the synaptic currents (2) and (3) are qualitatively similar. In both cases the model demonstrates a large variety of dynamical regimes, including a single steady state solution, bistability between two different steady states, macroscopic self-oscillations as well as bistability between the steady state and self-oscillations.

In this paper we consider two interacting populations of neurons of the described type

[TABLE]

We assume that the populations are identical. The index $j=1,\ldots,N$ labels the neurons inside of each population, while the index $k=0,1$ marks the populations. The term $\mathcal{I}_{k}$ describes both the interaction between the neurons inside of each population and between the populations. In this paper we mainly focus on the analysis of the situation when the internal and external interactions between neurons are modeled by the synaptic current in the form (3). The expression for term $\mathcal{I}_{k}$ in this case is

[TABLE]

where

[TABLE]

and the parameters $J_{in}$ and $J_{ex}$ define the coupling strengths within and between the populations, respectively. In Sec. V, we briefly describe another situation, when the interaction between the populations is provided by a mean field coupling rather than the synaptic coupling.

II.1 Thermodynamic limit

In the thermodynamic limit $N\to\infty$, the system (4) of two interacting neural populations can be reduced to a system of only four ordinary differential equations, which defines the dynamics of firing rates and mean membrane potentials of individual populations. Such macroscopic equations can be derived by a recently developed reduction method based on the Lorentzian ansatz. The idea of this method has been first proposed for a single network of QIF neurons interacting via instantaneous pulses Montbrió et al. (2015) and then extended to a more realistic model of synaptic interaction that takes into account the finite width of synaptic pulses Ratas and Pyragas (2016). The technique of the derivation of the macroscopic equations for the system (4) is similar to that described in Refs. Montbrió et al. (2015) and Ratas and Pyragas (2016), and thus we present this derivation in abbreviated form.

In the infinite-$N$ limit, the macroscopic state of each population in the system (4) can be described by continuous density functions $\rho_{k}(V|\eta,t)$, $k=0,1$. The product $\rho_{k}(V|\eta,t)\text{d}V$ defines the fraction of neurons in the $k$th population with the membrane potential between $V$ and $V+\text{d}V$ and parameter $\eta$ at time $t$. These density functions satisfy continuity equations

[TABLE]

where $\mathcal{I}_{k}$ is defined in Eq. (5). In the continuous limit, the sum $S_{k}$ defined in Eq. (6) becomes a double integral

[TABLE]

The main assumption of the LA ansatz is that the solutions of Eqs. (7), for any initial conditions, converge to a Lorentzian-shaped function (see Ref. Montbrió et al. (2015) for the relation between the LA ansatz and Ott-Antonsen ansatz Ott and Antonsen (2008))

[TABLE]

where time-dependent parameters $x_{k}(\eta,t)$ and $y_{k}(\eta,t)$ define the half-width and the centre of the distribution. The parameters $x_{k}(\eta,t)$ and $y_{k}(\eta,t)$ characterize all relevant dynamics of the system in a reduced subspace. They are related to the total firing rate $r_{k}(t)$ and the mean membrane potential $v_{k}(t)$ via integrals

[TABLE]

Substituting the LA (9) into the continuity Eqs. (7), one can derive a system of differential equations for $x_{k}(\eta,t)$ and $y_{k}(\eta,t)$

[TABLE]

which for the complex variable $w_{k}(\eta,t)\equiv x_{k}(\eta,t)+iy_{k}(\eta,t)$ can be written as

[TABLE]

A simplification can be gained by choosing the density distribution of the $\eta$ parameter in the Lorentzian function form

[TABLE]

with the width $\Delta$ and the center at $\bar{\eta}$. In this case the integrals (8) and (10) can be solved by extending $\eta$ to the complex plane and computing a contour integral over an infinitely large semicircle in the lower half-plane Montbrió et al. (2015). The values of these integrals are defined by the pole ${\eta=\bar{\eta}-i\Delta}$ of $g(\eta)$ function. This enables us to relate the complex variable $w_{k}$ with the firing rate and the mean membrane potential

[TABLE]

as well as obtain the explicit expression for the synaptic function

[TABLE]

Taking into account Eqs. (12) and (14), the firing rates and the mean membrane potentials satisfy the system of four differential equations

[TABLE]

where $k=0,1$. These equations together with Eqs. (5) and (15) form the closed macroscopic model for the network consisting of two synaptically coupled populations of QIF neurons.

III Symmetric solutions and their stability

The macroscopic Eqs. (16) possess permutational symmetry: they are invariant under the change of variables $(r_{0},v_{0},r_{1},v_{1})\to(r_{1},v_{1},r_{0},v_{0})$. This symmetry admits the existence of the symmetric solutions $(r_{1},v_{1})=(r_{0},v_{0})$. To analyze the stability of such solutions, it is convenient to introduce new variables

[TABLE]

where $(R,P)$ and $(Q,M)$ are the transverse and longitudinal coordinates, respectively. In the new coordinates $(R,P,Q,M)$, the trajectories of the symmetric solutions are placed in the invariant subspace $(0,0,Q,M)$, where the variables $Q$ and $M$ satisfy differential equations

[TABLE]

These equations are identical to the equations that describe the dynamics of a single population of QIF neurons with a modified coupling strength $J=J_{in}+J_{ex}$. The solutions of this system have previously been analyzed in Ref. Ratas and Pyragas (2016). It has been shown that the system (18) has two types of asymptotically stable solutions: fixed points (steady states) and limit cycles. For some values of the parameters, these solutions can coexist in different combinations giving rise to the bistability. The stable solutions of the system  (18) constitute longitudinally stable solutions of the system (16) in the invariant subspace $(0,0,Q,M)$, while their transverse stability requires a special analysis.

The transverse stability of the symmetric solutions is defined by the variational equations of the $(R,P)$ variables

[TABLE]

with the matrix

[TABLE]

and the parameter

[TABLE]

The variables $Q$ and $M$ in Eqs. (20) and (21) satisfy Eqs. (18). Equation (19) governs the dynamics of transverse deviations $(\delta R,\delta P)$ from the invariant subspace. When these deviations decay in time than the corresponding symmetric solution defined by Eq. (18) is transverse stable, otherwise it is unstable.

The analysis of transverse stability of fixed points and limit cycles is different. For the fixed points, the matrix $\boldsymbol{A}$ is constant and its eigenvalues $\lambda_{1,2}$ define the stability. If $\mathrm{Re}\lambda_{1,2}<0$, the symmetric fixed point solution of the system (16) is transverse stable. For the limit cycles, the matrix $\boldsymbol{A}$ depends on time periodically, $\boldsymbol{A}(t)=\boldsymbol{A}(t+T)$, where $T$ is the period of the limit cycle. In this case we recourse to the Floquet theory. We solve differential equations for the fundamental matrix $\boldsymbol{\Phi}(t)$: $\dot{\boldsymbol{\Phi}}(t)=\boldsymbol{A}(t)\boldsymbol{\Phi}(t)$ with the initial condition equal to the identity matrix, $\boldsymbol{\Phi}(0)=\mathbb{I}$. Then we compute the eigenvalues $\mu_{1,2}$ of the monodromy matrix $\boldsymbol{\Phi}(T)$. If $|\mu_{1,2}|<1$, the symmetric limit cycle solution of the system (16) is transverse stable. Note that the transverse stability of the symmetric solutions depends on the difference $J_{ex}-J_{in}$, while the solutions themselves are defined by the sum $J_{ex}+J_{in}$.

In Fig. 1, we show the results of the above stability analysis in the plain of parameters $(J_{ex},J_{in})$. The values of the parameters $\bar{\eta}=0$, $V_{th}=50$ and $\Delta=1$ are chosen so that the system  (18) is monostable for any values of $J_{ex}+J_{in}$. The red diagonal line $J_{ex}+J_{in}=J_{H}$ in Fig. 1 indicates the Hopf bifurcation of the system  (18). For $J_{ex}+J_{in}<J_{H}\approx 14.7$, the only attractor of the system is the fixed point and for $J_{ex}+J_{in}>J_{H}$ — a limit cycle. The white color in the figure corresponds to the regions where the above symmetric solutions are transverse stable. In the blue regions, these solutions are transverse unstable and thus here the system (16) has no symmetric attractors. In the next section, we analyze symmetry-broken solutions that appear in the blue regions of Fig. 1.

IV Non-symmetric solutions and bifurcations

The behavior of the interacting neural populations is different for inhibitory ($J_{ex}<0$) and excitatory ($J_{ex}>0$) couplings. Below we present the analysis of non-symmetric solutions for these two different cases in separate sections.

IV.1 Inhibitory coupling

The bifurcation diagram of non-symmetric solutions in the region of parameters $(J_{ex},J_{in})$ relevant to the inhibitory coupling between populations is shown in Fig. 2(a). Here as well as in Fig. 1 the white area corresponds to monostable symmetric states. In the colored region, there are non-symmetric attractors, while in the striped area the system exhibits coupling-induced bistability. Here, depending on the initial conditions, the system can approach either symmetric or non-symmetric state.

The variation of the coupling strengths within and between the populations leads to a rich variety of bifurcations, including limit point (LP), Hopf (H), Neimark-Sacker (NS), limit point of cycles (LPC), and period doubling (PD) bifurcations. In Fig. 2(a) these bifurcations are presented by colored curves marked with corresponding acronyms. A more detailed visualization of the bifurcations is given in the right-hand side of Fig. 2. Here the local maxima $r_{0}^{max}$ and $r_{1}^{max}$ of spiking rates of different populations are shown as functions of a smoothly varying coupling strength $J_{ex}$ for three different fixed values of $J_{in}=10,16$ and $20$. These values are presented in the panel (a) by dotted horizontal lines.

For low coupling strength $J_{in}=10$, the non-interacting populations ($J_{ex}=0$) have a symmetric fixed point attractor with equal time-independent spiking rates $r_{0}=r_{1}=\textit{const}$. They are marked in panels (b) and (e) by red points. The increase of the inhibitory coupling between the populations results in symmetry breaking through a limit point bifurcation. The symmetric fixed point attractor (red dots) exists in the interval $J_{ex}\in(-3.31,0)$. The non-symmetric fixed point attractor (black dots) with different stationary spiking rates in different populations, $r_{0}\neq r_{1}$, is in the interval $J_{ex}\in(-6,-2.35)$. Following Ref. Olmi et al. (2010) we refer to this solution as a splay state. The splay state is characterized by the absence of any collective dynamics, since both the spiking rates and mean fields of populations are constant. In the interval $J_{ex}\in(-3.31,-2.35)$, the system exhibits bistability with coexisting symmetric and non-symmetric fixed point attractors.

More interesting symmetry breaking scenarios are observed at higher coupling strength $J_{in}>J_{H}\approx 14.7$, when the synchronization between neurons within the isolated populations ($J_{ex}=0$) causes macroscopic limit cycle oscillations of their spiking rates and mean fields. In panels (c) and (f), we fix $J_{in}=16$ and continue the symmetric limit cycle solution to the region $J_{ex}<0$ (green dots). With the increase of the inhibitory coupling between the populations, the symmetry of the limit cycle is conserved up to the value $J_{ex}=-0.76$. At this point the system undergoes the Neimark-Sacker bifurcation. Further increase of the coupling strength $|J_{ex}|$ results in quasiperiodic oscillations, which then transform into periodic oscillations and finally become the symmetric fixed point attractor (red dots) following the Hopf bifurcation at $J_{ex}=-3.15$. This symmetric fixed point attractor disappears via limit point bifurcation at $J_{ex}=-5.28$. In addition to the above solutions, the system has a non-symmetric limit cycle attractor presented by blue dots. The attractor appears via a limit point cycle bifurcation at $J_{ex}=-2.41$ and exists for any $J_{ex}<-2.41$. This solution is most interesting since it represents a chimera-like state. Here the spiking rate in one of the populations oscillates at high amplitude, while the spiking rate of other population has a low amplitude of oscillations. Below we discuss such solutions in more details. Note that the system is bistable in the interval $J_{ex}\in(-5.28,-2.41)$. In this interval, the chimera-like state may coexist either with a symmetric fixed point or with a non-symmetric limit cycle or with quasiperiodic oscillations.

Further increase of the coupling strength $J_{in}$ within the populations revokes the bistability. In panels (d) and (g) this is demonstrated for the fixed value $J_{in}=20$. The continuation of the symmetric limit cycle solution into the region $J_{ex}<0$ shows that the system has no coexisting attractors for any $J_{ex}<0$. Unlike in the previous case, here the symmetry breaking appears through a period doubling bifurcation, which occurs at $J_{ex}=-1.46$. Latter this bifurcation leads to a chaotic regime in which high and low activities of neurons change irregularly between the populations. The microscopic dynamics of the system in this regime will be presented below. The chaotic regime is replaced by a periodic solution via a limit point cycle bifurcation at $J_{ex}=-3.25$. The non-symmetric periodic solution that exists for $J_{ex}<-3.25$ represents the above-mentioned chimera-like state.

To better visualize the chimera-like state, in Fig. 3 we fix $J_{ex}=-4$ and show the projections of the trajectory to the $(v_{0},r_{0})$ and $(v_{1},r_{1})$ plains as well as the dynamics of the spiking rates $r_{0}(t)$ and $r_{1}(t)$ in different populations. We see that the spiking rates in different populations oscillate with considerably different amplitudes. In one of the populations the spiking rate is close to zero (blue curve in the figure) and its oscillations are almost unremarkable on the scale of variation of the spiking rate of other population (red curve).

IV.1.1 Simulation of the microscopic model

Now we consider the problem of two interacting QIF neuron populations on the microscopic level. Numerical simulations of the microscopic model Eqs. (4) are interesting for two reasons. First, we can verify whether the macroscopic Eqs. (16) derived in the thermodynamic limit $N\to\infty$ predict well the dynamical regimes of a finite-size network. Second, such simulations allow us to observe the behavior of individual neurons for any particular macroscopic regime.

Direct numerical integration of Eqs. (4) is problematic because the membrane potential of the QIF neuron tends to infinity at the moment when the neuron fires. This problem can be avoided by the change of variables

[TABLE]

that transforms the QIF neurons into theta neurons, where $\theta_{j,k}$ is the phase of the $j$th neuron in the $k$th population. Then the model (4) in the theta representation reads as follows:

[TABLE]

When the QIF neuron fires, its membrane potential approaches infinity, $V_{j,k}\to\infty$, and then its value is reset to minus infinity, $V_{j,k}\to-\infty$. In the theta representation, this process is smooth: the phase $\theta_{j,k}$ simply crosses the value $\pi$.

The parameters $\eta_{j}$, satisfying the Lorentzian distribution Eq. (13), were generated deterministically by using formula $\eta_{j}=\bar{\eta}+\Delta\tan\left[(\pi/2)(2j-N-1)/(N+1)\right]$, where $j=1,\ldots,N$ and $\Delta=1$. Such a numeration of neurons means that the isolated neurons with the index $j<j_{c}=(N+1)/2-(2N+1)\arctan(\bar{\eta})/\pi$ are excitable and the neurons with the index $j>j_{c}$ are spiking. At each step of integration of Eqs. (23), the synaptic variables (6) were estimated as $S_{k}(t)=dN_{k}^{S}/N$, where $dN_{k}^{S}$ is the number of neurons in the $k$th population whose phases are in the interval $\theta_{j,k}\in[2\arctan(V_{th}),\pi]$. Similarly, the firing rates were estimated as $r_{k}=dN_{k}^{r}/(Ndt)$, where $dN^{r}_{k}$ is the number of neurons in the $k$th population whose phases are in the interval $\theta_{j,k}\in(\pi-2dt,\pi)$. Such estimations are based on the assumption that the time step $dt$ is small and thus the phase speed of neurons close to the firing phase $\theta=\pi$ can be approximated as $\dot{\theta}_{j,k}\approx 2$. Because of the finite number of neurons, the quantities $dN_{k}^{r}$ fluctuate in time and the firing rates vary nonsmoothly. For better visualization, we smoothed these quantities by using a moving average with a time window of the size $\delta t=5\cdot 10^{-2}$.

The three columns in Fig. 4 show the dynamics of the microscopic model for three different dynamical regimes defined by different choices of the coupling strengths $(J_{ex},J_{in})$. In order from left to right, they are $(-4,10),(-4,20)$ and $(-2.5,20)$. These values were marked in Fig. 2(a) by crosses. The upper (a)-(c) and bottom (d)-(f) panels in Fig. 4 show, respectively the dynamic of the spiking rates and phases of individual neurons in both populations.

The panels (a) and (d) correspond to the splay state. According to the macroscopic model, this state is characterized by stationary spiking rates, which are different in different populations. For the given values of the parameters, the macroscopic model predicts the values of the spiking rates $r_{0}=0.09$ and $r_{1}=0.98$. In panel (a) they are shown by horizontal lines. The spiking rates obtained from the microscopic model fluctuates around these predicted values. The fluctuations occur due to the finite size of the network. In the population with a small spiking rate, almost all neurons are quenched, while in the population with a larger spiking rate almost all neurons spike, however, their spikes are incoherent and they do not produce any macroscopic oscillations. Note that for the given value $\bar{\eta}=0$, the half of the isolated ($J_{ex}=J_{in}=0$) neurons in each population are spiking and another half are excitable.

The chimera-like state reproduced by the microscopic model is shown in panels (b) and (e). The values of the parameters are the same as in Fig. 3, where this state is demonstrated via the macroscopic model. Comparing Figs. 3(b) and 4(b), we see that the dynamics of the spiking rates derived from the microscopic and macroscopic models are in good agreement. In the chimera-like state, the majority of neurons in one of the populations are quenched and their spiking rate is close to zero, while in another population the majority of neurons spike synchronously and produce large amplitude oscillations of the spiking rate.

Finally, in panels (c) and (f) we demonstrate the solutions of microscopic equations for a chaotic chimera-like state. Here synchronous spiking erratically jumps from one to another population. For the given values of the parameters, the macroscopic model exhibits similar behavior (not shown).

IV.2 Excitatory coupling

We turn now to the situation when couplings within and between the populations are both excitatory, $J_{in}>0$ and $J_{ex}>0$. In this case the bifurcation scenarios are less diverse than for the inhibitory coupling considered above. Now the qualitative change of solutions is mainly defined by only two bifurcations, namely, the Neimark-Sacker and branch point of cycles (BPC) bifurcations.

In Fig. 5(a), we present the bifurcation diagram of the macroscopic model (16) in the region of parameters $(J_{ex},J_{in})$ relevant to the excitatory coupling between populations. Here as well as in Fig. 2 the white area corresponds to monostable symmetric states and the colored region represents non-symmetric attractors. At the border of the colored region (red curve) the system undergoes BPC bifurcation. At this bifurcation the symmetric limit cycle looses stability and there appears a symmetrical pair of stable non-symmetric limit cycles. The two green curves in the colored region indicate Neimark-Sacker bifurcation.

In panels (b) and (c), we fix $J_{in}=20$ and continue the solutions of Eqs.  (16) into the region $J_{ex}>0$. Unlike to Figs. 2(b)-(g), here we plot not only the local maxima $r_{0}^{max}$ and $r_{1}^{max}$ of the spiking rates of different populations, but also the values $r_{1}$ and $r_{0}$, respectively, at the moments when the neighbor population has reached its local maximum. This representation allows us to distinguish antiphase periodic oscillations $r_{0}(t)=r_{1}(t+T/2)$, where $T$ is the period of the limit cycle. Such oscillations are common for the excitatory coupling. For fixed $J_{in}=20$, the antiphase periodic oscillations take place in the region $0<J_{ex}<4.5$. In the region $4.5<J_{ex}<8.3$, between two Neimark-Sacker bifurcations, the system exhibits complex dynamics, including high-periodic cycles, quasiperiodicity and chaos. Then for $8.3<J_{ex}<11.9$, the chimera-like state appears. This state is destroyed via a branch point of cycles bifurcation at $J_{ex}=11.9$, and for $J_{ex}>11.9$ the symmetric limit cycle oscillations are established.

The microscopic dynamics of the system for the periodic chimera-like state at $(J_{ex},J_{in})=(10,20)$ [in Fig. 5(a) these values are marked by a cross] is demonstrated in Fig. 6. In panel (a) the projections of the microscopic dynamics to the plains $(r_{0},v_{0})$ and $(r_{1},v_{1})$ are compared with the corresponding projections obtained from the macroscopic model. Again, we see that the macroscopic model predicts well the dynamics of the finite-size network. Panels (b) and (c) show the time-dependence of the spiking rates and phases of neurons, respectively, for both populations. Comparing Figs. 6(c) and 4(e), we reveal that the microscopic dynamics of the system in chimera-like state are different for the inhibitory and excitatory coupling between populations. For the inhibitory coupling, almost all neurons in one population spike in synchrony and in another population the majority of neurons are quenched. For the excitatory coupling, in both populations almost all neurons spike, but in one population, they spike synchronously and in another – asynchronously. These two different microscopic dynamics lead to the similar macroscopic result: for one population, the spiking rate exhibits large amplitude oscillations, while for another population it is close to zero.

V Mean field coupling

In this section, we additionally consider another type of coupling between the populations. As well as above, we assume that the interaction within the populations is defined by the synaptic current (3), however the coupling between the populations is provided by mean field of membrane potential. Experimentally, such a situation can be imagined as a control problem. Assume we have two non-interacting populations of neurons and can separately measure their mean membrane potentials. Then we stimulate the first population by a signal proportional to the mean field measured from the second population and vice versa. In that case the total current $\mathcal{I}_{k}$ that defines all interactions between and within the populations in Eq. (4) takes the form:

[TABLE]

where

[TABLE]

is the mean membrane potential of the $k$th population. The macroscopic dynamics of this system are again described by Eqs. (16), but the total current is now defined by Eq. (24).

An example of chimera-like state obtained with an inhibitory mean field coupling between the populations is demonstrated in Fig. 7. In numerical simulations of the microscopic model, the mean membrane potential (25) were estimated by ignoring the contribution from neurons with extremely large values of membrane potential, when $V_{l,k}>200$. This allowed us to avoid the divergence of the sum (25) due to the finite-size effect. The dynamics of the system presented in Fig. 7(b)-(c) is similar to that shown in Fig. 4(b)-(e), which has been obtained with the inhibitory synaptic coupling.

VI Discussion

In this paper, we analyzed the dynamics of a neural network consisting of two identical populations of quadratic integrate-and-fire neurons. Both populations are heterogeneous; they include inherently spiking and excitable neurons. The main part of the paper is devoted to the synaptic coupling within and between the populations. The coupling is global and takes into account the finite width of synaptic pulses. We assumed that the interactions within the populations are excitatory, and considered two cases, with the inhibitory and excitatory coupling between the populations.

Using a recently developed reduction technique Montbrió et al. (2015) based on the Lorentzian ansatz, we derived a macroscopic model that describes the neuron spiking rates and mean membrane potentials in different populations. The model is defined by only four ordinary differential equations, which are exact in the limit of the infinite-size network. Such a simplification allowed us to perform a thorough bifurcation analysis of the system. In the parameter plane defined by the coupling strengths within and between the populations, we identified the areas where the symmetric solutions (with the identical dynamics in both populations) lose their stability and non-symmetrical solutions are established. Our analysis showed that the competition of neural interactions within and between the populations may lead to a rich variety of non-symmetric patterns, including a splay state, antiphase periodic oscillations, chimera-like states, also chaotic oscillations as well as bistabilities between various states. The most interesting non-symmetric pattern is the chimera-like state. Here the neurons in one population behave synchronously and produce high amplitude oscillations of the spiking rate, while the neurons in another population are quenched or desynchronized and their spiking rate is close to zero. The chimera-like state exists for both the inhibitory and excitatory synaptic coupling between the populations. In addition, we showed that the chimera-like state may appear when the synaptic coupling between the populations is replaced with a mean field coupling.

To verify the validity of the macroscopic model we performed numerical simulations of the microscopic model equations for a finite-size network. As a result, we were convinced that the macroscopic model predicts well the behavior of a finite-size network consisting of only thousand neurons in each population. Thus the macroscopic models of the type considered here are natural candidates for use in future large-scale brain simulations. Such models can be considered as an alternative to neural mass models Destexhe and Sejnowski (2009), which are especially useful for understanding brain rhythms. The neural mass models also employ only several differential equations to describe the coarse-grained activity of large-scale neural networks, however, they are phenomenological in nature. In contrast, the approach discussed in this paper provides an exact macroscopic description of an underlying microscopic spiking neurodynamics.