# Kuramoto model for populations of quadratic integrate-and-fire neurons with
chemical and electrical coupling

Pau Clusella  Department of Experimental and Health Sciences, Universitat
Pompeu Fabra, Barcelona Biomedical Research Park, 08003, Barcelona, Spain
Bastian Pietras  Institute of Mathematics, Technical University Berlin, 10623
Berlin, Germany.  Bernstein Center for Computational Neuroscience Berlin,
10115 Berlin, Germany.  Ernest Montbrió  Neuronal Dynamics Group, Department
of Information and Communication Technologies, Universitat Pompeu Fabra, 08018
Barcelona, Spain.

###### Abstract

We derive the Kuramoto model (KM) corresponding to a population of weakly
coupled, nearly identical quadratic integrate-and-fire (QIF) neurons with both
electrical and chemical coupling. The ratio of chemical to electrical coupling
determines the phase lag of the characteristic sine coupling function of the
KM, and critically determines the synchronization properties of the network.
We apply our results to uncover the presence of chimera states in two coupled
populations of identical QIF neurons. We find that the presence of both
electrical and chemical coupling is a necessary condition for chimera states
to exist. Finally, we numerically demonstrate that chimera states gradually
disappear as coupling strengths cease to be weak.

> The Kuramoto model (KM) is a minimal mathematical model for investigating
> the emergence of collective oscillations in populations of heterogeneous,
> self-sustained oscillators Kuramoto (1975, 1984). Though the KM model was
> not originally intended to describe any specific natural system, an abundant
> body of work applies it to explore large-scale neuronal oscillations; see,
> e.g., Refs. Glomb  _et al._ , 2021; Lynn and Bassett, 2019; Breakspear,
> Heitmann, and Daffertshofer, 2010; Cabral  _et al._ , 2011; Villegas,
> Moretti, and Muñoz, 2014; Ponce-Alvarez  _et al._ , 2015; Sanz-Leon  _et
> al._ , 2015; Schmidt  _et al._ , 2015; Petkoski  _et al._ , 2016; Andrzejak
> _et al._ , 2016; Petkoski, Palva, and Jirsa, 2018; Roberts  _et al._ , 2019;
> Choi and Mihalas, 2019; Ziaeemehr  _et al._ , 2020; Noori  _et al._ , 2020;
> Jung, Eickhoff, and Popovych, 2021; Tewarie  _et al._ , 2021; Ódor, Kelling,
> and Deco, 2021; Pang, Gollo, and Roberts, 2021; Weerasinghe  _et al._ ,
> 2021. Yet, it remains unclear how the parameters of the KM relate to
> parameters—such as chemical or electrical synaptic strengths—critical for
> setting up synchronization in biophysically realistic neuronal models Wang
> (2010); Börgers (2017). Here, we unveil a mathematical relation between a
> popular spiking neuron model, the quadratic integrate-and-fire (QIF)
> Ermentrout (1996); Ermentrout and Kopell (1986); Izhikevich (2007), with a
> well-known variant of the KM Sakaguchi, Shinomoto, and Kuramoto (1988);
> Sakaguchi and Kuramoto (1986). This provides support in favor of the use of
> the KM for modeling studies in computational neuroscience and introduces the
> powerful mathematical framework of the KM Pikovsky, Rosenblum, and Kurths
> (2001); Strogatz (2000); Pikovsky and Rosenblum (2015) for the analysis of
> the dynamics of QIF networks.

##  I Introduction

Large-scale neuronal oscillations emerge due to the synchronous interplay of
ensembles of neurons. These oscillations are successfully replicated by
mathematical models of spiking neurons, which also allow for a mechanistic
understanding of neuronal rhythmogenesis Wang (2010); Börgers (2017).
According to these theories, inhibitory synapses play a central role in
setting up neuronal synchronization either in isolation Wang and Buzsáki
(1996) or due to their interplay with excitatory neurons Wilson and Cowan
(1972). Additionally, inhibitory cells are very often coupled electrically,
and this coupling is usually mediated by so-called gap junctions Nagy, Pereda,
and Rash (2018). Such electrical synapses are well-known to largely favor
synchrony.

Recently, important efforts have been put forward to model the oscillatory
dynamics of so-called whole-brain networks Deco  _et al._ (2015); Lynn and
Bassett (2019); Glomb  _et al._ (2021). To facilitate both the analysis and
the computational work many studies do not use spiking neuron models, but
apply the mathematical framework of the Kuramoto model (KM); see, e.g., Refs.
Glomb  _et al._ , 2021; Lynn and Bassett, 2019; Breakspear, Heitmann, and
Daffertshofer, 2010; Cabral  _et al._ , 2011; Villegas, Moretti, and Muñoz,
2014; Ponce-Alvarez  _et al._ , 2015; Sanz-Leon  _et al._ , 2015; Schmidt  _et
al._ , 2015; Petkoski  _et al._ , 2016; Andrzejak  _et al._ , 2016; Petkoski,
Palva, and Jirsa, 2018; Roberts  _et al._ , 2019; Choi and Mihalas, 2019;
Ziaeemehr  _et al._ , 2020; Noori  _et al._ , 2020; Jung, Eickhoff, and
Popovych, 2021; Tewarie  _et al._ , 2021; Ódor, Kelling, and Deco, 2021; Pang,
Gollo, and Roberts, 2021; Weerasinghe  _et al._ , 2021. Yet, it remains
unclear how to relate the parameters of the KM to bio-physically meaningful
parameters, such as synaptic strengths.

In this paper, we aim to theoretically substantiate the use of the KM for
neuronal modeling, by providing a mathematical link between the quadratic
integrate-and-fire (QIF) model and the KM. We derive a KM for QIF neurons and
subsequently justify its validity in two different ways: First, we compare the
predictions of the KM with those of an exact mean-field model for QIF
neurons—often referred to as firing rate or neural mass model (NMM) Pietras
_et al._ (2019); Montbrió and Pazó (2020); Montbrió, Pazó, and Roxin (2015).
Second, we use two populations of _identical_ Kuramoto oscillators to find so-
called chimera states Montbrió, Kurths, and Blasius (2004); Abrams  _et al._
(2008); Laing (2009); Martens, Bick, and Panaggio (2016). In a chimera state,
one of the two homogeneous populations displays in-phase synchrony, and NMMs
are useless in this case. However, the KM for QIF neurons is perfectly suited
to describe full synchrony, and we exploit this to uncover the existence of
chimera states in two-population networks of QIF neurons.

Our derivation of the KM for QIF neurons mainly builds on a previous work by
Izhikevich Izhikevich (2007) and also on Refs. Pazó and Montbrió, 2014;
Montbrió and Pazó, 2018. In Chap. 10 of Ref. Izhikevich (2007), Izhikevich
applied perturbation methods to derive a simplified model that approximated
the dynamics of two identical QIF neurons, with _either_ chemical or
electrical coupling 111More recently a similar derivation has been obtained
for populations of heterogeneous Winfree oscillators with sinusoidal
infinitesimal Phase Resetting Curves (iPRC) Pazó and Montbrió (2014). When the
oscillators have the iPRC of the QIF neuron, the population of Winfree
oscillators can be well approximated to to a population of QIF neurons with
(weak) chemical synapses Montbrió and Pazó (2018).. Here, we extend the work
of Izhikevich and derive a model that approximates the dynamics of an
_ensemble of heterogeneous_ QIF neurons with _both_ chemical and electrical
coupling. The approximated model turns out to be a well-known version of the
KM Sakaguchi, Shinomoto, and Kuramoto (1988); Sakaguchi and Kuramoto (1986)
and is valid when both heterogeneities and coupling strengths are weak.

This paper is organized as follows: In Section II, we introduce the QIF
population model, and in Section III, we describe the method to reduce the QIF
model to the KM. In Section IV, we analyze the dynamics of the KM and
demonstrate that it correctly describes the collective dynamics of populations
of nearly identical QIF neurons, with weak electrical and chemical synapses.
In Section V, we exploit the KM to uncover the presence of chimera states in
coupled populations of identical QIF neurons. Finally, in Section IV, we
briefly discuss and summarize our results.

##  II Population of QIF neurons with electrical and chemical synapses

We investigate a population of $N$ quadratic integrate-and-fire (QIF) neurons
$i=1,\dots,N$ interacting all-to-all via both electrical and chemical synapses
Kopell and Ermentrout (2004); Laing (2015); Pietras  _et al._ (2019)

$\tau\dot{V}_{i}=V_{i}^{2}+\eta_{i}+\epsilon~{}I_{i,syn}(t),~{}~{}\text{if
}V_{i}>V_{p},\text{ then }V_{i}\leftarrow V_{r},$ (1)

where $V_{i}$ is the membrane potential of neuron $i$, $\tau$ is the membrane
time constant of the neurons, and $\eta_{i}$ represents an external current,
which varies from cell to cell. Due to the quadratic nonlinearity of the QIF
model, the membrane potential blows up in finite time, and a resetting rule is
needed: When the neurons reach the peak value $V_{p}$, they emit a spike and
the voltage is reset to $V_{r}$. We assume symmetric spike resetting,
$V_{p}=-V_{r}$ and $V_{p}\to\infty$, so that the QIF model is equivalent to
the so-called theta-neuron Ermentrout and Kopell (1986); Ermentrout (1996). In
addition, we consider $\eta_{i}>0$ and hence, in the absence of synaptic
inputs ($I_{j,syn}=0$), QIF neurons are self-sustained oscillators. Finally,
synaptic inputs (whose total strength is controlled by the _small_ parameter
$\epsilon\geq 0$) are composed of electrical and chemical synapses,

$I_{i,syn}(t)=g(v(t)-V_{i})+J\tau r(t).$ (2)

Specifically, electrical synapses (of strength $\epsilon g\geq 0$) diffusively
couple each neuron with the mean membrane potential

$v(t)=\frac{1}{N}\sum_{j=1}^{N}V_{j}(t).$ (3)

Electrical synapses mostly connect inhibitory neurons, and hence the chemical
synaptic strength, $J$, is thought of as a negative parameter thereafter.
Finally, chemical synapses (of strength $\epsilon J$) are mediated by the mean
firing rate

$r(t)=\frac{1}{N}\sum_{j=1}^{N}\sum_{k}\delta\left(t-t_{j}^{(k)}\right),$ (4)

where $t_{j}^{(k)}$ is the time of the $k$th spike of the $j$th neuron and
$\delta(t)$ is the Dirac delta function.

##  III Derivation of the Kuramoto model for populations of QIF neurons

In the following we derive the Kuramoto model corresponding to Eq. (1). The
derivation exploits well-known mathematical methods that are reviewed for
example in Refs. Hoppensteadt and Izhikevich, 1997; Izhikevich, 2007; Pietras
and Daffertshofer, 2019.

We perform the derivation of the KM as follows: First, we obtain the phase
resetting curve (PRC) of the QIF model. Second, invoking weak coupling, we
derive the so-called Winfree model corresponding to the QIF model. Finally, we
assume weak heterogeneity and apply the method of averaging to obtain the
Kuramoto model corresponding to Eq. (1).

###  III.1 Phase resetting curve (PRC) of a QIF neuron

We first consider an isolated, regularly spiking QIF neuron; i.e.,
$I_{i,syn}=0$ and $\eta_{i}>0$. The solution of the QIF model immediately
after a spike, $V_{i}(0)=-\infty$, is

$V_{i}(t)=\sqrt{\eta_{i}}\tan\left(t\sqrt{\eta_{i}}/\tau-\pi/2\right).$ (5)

The frequency of the oscillations is

$\Omega_{i}=2\sqrt{\eta_{i}}/\tau,$ (6)

and a phase variable $\theta_{i}$ can be defined in the interval $[0,2\pi)$ as

$\theta_{i}=\Omega_{i}t=2\arctan(V_{i}/\sqrt{\eta_{i}})+\pi.$ (7)

Next, we assume that the neuron is perturbed so that its membrane potential
instantaneously changes from $V_{i}$ to $V_{i}+\delta V$. Then, the new phase
after the perturbation is $\theta_{i,new}=2\arctan((V_{i}+\delta
V)/\sqrt{\eta_{i}})+\pi$.

The PRC measures the phase shift produced by the perturbation, i.e.
PRC$=\theta_{i,new}-\theta_{i}$. Hence, for the QIF model, the PRC is
Izhikevich (2007)

$\text{PRC}(\theta_{i},\delta V)=2\arctan(\delta
V/\sqrt{\eta_{i}}-\cot(\theta_{i}/2))+\pi-\theta_{i}.$ (8)

This function depends on both the strength of the perturbation and the phase
of the neuron at the instant of the perturbation. The PRC Eq. (8) is always
positive, indicating that positive/negative perturbations only produce
positive/negative phase shifts. This characterizes the so-called Class 1
neuronal oscillators Ermentrout (1996).

###  III.2 Weak coupling approximation and the Winfree model

The PRC Eq. (8) exactly characterizes the phase response of the QIF neuron to
a perturbation. Next, we invoke weak coupling, which allows for deriving a new
phase model—called the Winfree model—that approximates the network Eq. (1) for
$\epsilon\ll 1$.

Weak perturbations produce small changes in the membrane potential, $|\delta
V|\ll 1$. Then the PRC scales linearly with the strength of the perturbation
Winfree (1967)

$\text{PRC}(\theta_{i},\delta V)\approx Z(\theta_{i})\delta V,$

where $Z(\theta_{i})$ is called phase sensitivity function or infinitesimal
phase resetting curve (iPRC). For the QIF model, the iPRC is

$Z(\theta_{i})=\left.\frac{\partial\text{PRC}(\theta_{i},\delta
V)}{\partial(\delta V)}\right|_{\delta
V=0}=\frac{1-\cos\theta_{i}}{\sqrt{\eta_{i}}}.$ (9)

When weak perturbations are described by a continuous function $P(t)$ with
$|P(t)|\ll 1$, the infinitesimal change in the phase due to the perturbations
is $d\theta=Z(\theta)P(t)dt$. Accordingly, assuming weak coupling $\epsilon\ll
1$, the population of QIF neurons Eq. (1) is well approximated by the Winfree
model,

$\displaystyle\dot{\theta}_{i}$
$\displaystyle=\Omega_{i}+\frac{\epsilon}{\tau}(1-\cos\theta_{i})\sum_{j=1}^{N}P(\theta_{i},\theta_{j}),$
(10)

where perturbations to neuron $i$ are due to synaptic inputs from neuron $j$
and can be written in terms of the phase variables as

$P(\theta_{i},\theta_{j})=\frac{g}{N}\left(\cot(\theta_{i}/2)-\sqrt{\frac{\eta_{j}}{\eta_{i}}}\cot(\theta_{j}/2)\right)+\frac{2J}{N}\sqrt{\frac{\eta_{j}}{\eta_{i}}}\delta\left(\theta_{j}\right).$
(11)

Recall that $\theta_{j}\in[0,2\pi)$ so that the Dirac delta function in Eq.
(11) has argument zero whenever neuron $j$ fires a spike.

###  III.3 Weak heterogeneity and the averaging approximation

The Winfree model can be further simplified using the method of averaging. We
consider the external currents in Eq. (1) as a common current $\bar{\eta}$
plus a weakly distributed parameter as

$\eta_{i}=\bar{\eta}+\epsilon\chi_{i}.$ (12)

In the derivation of the Winfree model, we already assumed weak coupling,
$\epsilon\ll 1$. Therefore, the smallness of parameter $\epsilon$ implies now
the smallness of both coupling terms and the level of heterogeneity. This
assumption allows for a separation of time scales so that the phases
$\theta_{i}$ can be written as

$\theta_{i}=\Phi+\phi_{i},$ (13)

where $\Phi$ describes the fast, free-running oscillation of period

$T=\tau\pi/\sqrt{\bar{\eta}},$

whereas the phases $\phi_{i}$ describe slow phase drifts produced by weak
heterogeneities and synaptic inputs. Substituting Eq. (13) into the Winfree
model Eqs. (10,11) and collecting terms of order $\epsilon$, we find the
evolution equation for the slow phases

$\displaystyle\dot{\phi}_{i}=\frac{\epsilon\chi_{i}}{\tau\sqrt{\bar{\eta}}}+[1-\cos(\Phi+\phi_{i})]\frac{\epsilon}{\tau
N}\sum_{j=1}^{N}p(\Phi+\phi_{i},\Phi+\phi_{i}+\Delta_{ji}).$ (14)

Here, we defined pairwise phase differences as $\Delta_{ji}=\phi_{j}-\phi_{i}$
and a function describing synaptic perturbations as

$p(x,y)=g\left[\cot(x/2)-\cot(y/2)\right]+2J\delta(y).$

To apply the method of averaging to Eq. (14), we consider that in one period
of the fast oscillation, $T$, the slow phases $\phi_{i}$ can be assumed
constant. Then, Eq. (14) reduces to

$\dot{\phi}_{i}=\frac{\epsilon}{\tau N}\sum_{j=1}^{N}\Gamma(\Delta_{ji}),$
(15)

where the coupling function $\Gamma$ is obtained by averaging the r.h.s. of
Eq. (14) over one period $T$. This involves the evaluation of four integrals
that can be explicitly computed and yields the phase interaction function

$\Gamma(\Delta_{ji})=\frac{\chi_{i}}{\sqrt{\bar{\eta}}}+g\sin\Delta_{ji}+\frac{J}{\pi}(1-\cos\Delta_{ji}).$
(16)

###  III.4 Kuramoto model for populations of QIF neurons

Substituting Eq. (16) into Eq. (15) and expressing the result in terms of the
original phases Eq. (13), we find the Kuramoto model

$\dot{\theta}_{i}=\omega_{i}+\frac{\epsilon}{\tau
N}\sum_{j=1}^{N}\left[g\sin(\theta_{j}-\theta_{i})-\frac{J}{\pi}\cos(\theta_{j}-\theta_{i})\right]+\frac{\epsilon}{\tau}\frac{J}{\pi}$
(17)

with natural frequencies

$\omega_{i}=\frac{2\sqrt{\bar{\eta}}}{\tau}+\epsilon\frac{\chi_{i}}{\tau\sqrt{\bar{\eta}}}.$
(18)

In the absence of electrical synapses $g=0$, Eqs. (17) essentially 222For
$g=0$ the natural frequencies in Refs. Pazó and Montbrió (2014); Montbrió and
Pazó (2018) differ from Eq. (18). The reason for this discrepancy is that
Refs. Pazó and Montbrió (2014); Montbrió and Pazó (2018) consider the Winfree
model with distributed natural frequencies $\Omega_{i}$, while here we study
the QIF model with distributed currents $\eta_{i}$. reduces to the Kuramoto
model with chemical synapses derived in Refs. Pazó and Montbrió, 2014 and
Montbrió and Pazó, 2018. The KM for QIF neurons Eq. (17) generalizes the
results in Refs. Pazó and Montbrió, 2014 and Montbrió and Pazó, 2018 to
networks with both electrical and chemical coupling, and it is our main
result.

Eq. (18) is the linear approximation of Eq. (6) for weak heterogeneity—see
also Eq. (12). The last term of Eq. (17) describes the deviation of the
natural frequencies due to synaptic coupling, which exclusively depends on
chemical coupling. Excitatory coupling ($J>0$) speeds up the frequencies of
the oscillators, and inhibition ($J<0$) slows them down. These frequency
shifts do not qualitatively affect the collective dynamics of Eq. (17), but
they may become relevant if the oscillators are not all-to-all coupled
Sakaguchi, Shinomoto, and Kuramoto (1988); Blasius and Tönjes (2005) or in the
case of interacting excitatory and inhibitory populations Montbrió and Pazó
(2018).

Alternatively, Eq. (17) can be cast in the more transparent form Sakaguchi,
Shinomoto, and Kuramoto (1988); Sakaguchi and Kuramoto (1986)

$\dot{\theta}_{i}=\omega_{i}+\frac{K}{N}\sum_{j=1}^{N}\left[\sin(\theta_{j}-\theta_{i}-\alpha)+\sin\alpha\right],$
(19)

with the coupling constant

$K=\frac{\epsilon}{\tau}\sqrt{(J/\pi)^{2}+g^{2}},$ (20)

and the phase lag parameter

$\alpha=\arctan\left(\frac{J/\pi}{g}\right).$ (21)

The coupling parameters $K$ and $\alpha$ satisfy a simple geometric relation
with the coupling parameters of the QIF model Eq. (1), illustrated in Fig. 1.
Given a particular choice of the QIF coupling parameters, Fig. 1 shows that
electrical coupling and chemical coupling (divided by a factor $\pi$)
contribute equally to the overall coupling strength $K$ of the KM, and that
$K$ is insensitive to the sign of the chemical coupling—in Fig. 1 we consider
an inhibitory network; i.e., $J<0$.

To lighten the notation, we consider $\epsilon=1$ thereafter. Hence, for the
KM Eq. (19) to be a good approximation of Eq. (1), in the following, the
synaptic weights $J$ and $g$ need to be regarded as small quantities.

Figure 1: Geometric relation—determined by Eqs. (20,21)—between the coupling
parameters of the QIF model Eqs. (1) and the Kuramoto model Eqs. (19).

##  IV Analysis of the Kuramoto model for Quadratic Integrate-and-fire neurons

Using Fig. 1—or, equivalently, Eq. (21)—we may infer how chemical and
electrical synapses contribute to synchronization, using well-known results
for the KM. For example, the phase constant $\alpha$ critically determines the
synchronization behavior of Eq. (19) Sakaguchi and Kuramoto (1986). In the
absence of electrical coupling, $g=0$, we find $\alpha=+\pi/2$ for excitatory
coupling and $\alpha=-\pi/2$ for inhibitory coupling. This indicates that
collective synchronization is unreachable—consistent with the well-known fact
that instantaneous chemical coupling is unable to synchronize type 1 neuronal
oscillators Ermentrout (1996); Hansel, Mato, and Meunier (1995); Devalle,
Roxin, and Montbrió (2017). In contrast, in the absence of chemical coupling,
one finds $\alpha=0$, and Eq. (19) reduces to the standard KM, in which
collective synchronization is achieved at a critical degree of heterogeneity
$\Delta=\Delta_{c}(K)$ that depends on the coupling strength Kuramoto (1984).
Between these two extreme cases, that is in networks with both electrical and
chemical synapses, we find the phase lag parameter $|\alpha|\in(0,\pi/2)$, and
synchronization generally depends on both $\alpha$ and the overall shape of
the distribution of natural frequencies Sakaguchi and Kuramoto (1986).

To validate the KM for QIF neurons, in Section IV A, we obtain the mean-field
model corresponding to Eqs. (19) and compare its predictions with those of the
mean-field model derived in Ref.Pietras  _et al._ (2019), which describes the
dynamics of the QIF network Eq. (1) exactly.

###  IV.1 Mean-field model

Figure 2: Synchronization boundary of the QIF network Eq. (1) with the
Lorentzian distribution of currents for various values of the (scaled)
inhibitory coupling strength, $J/\sqrt{\bar{\eta}}$. The solid line
corresponds to the approximated critical width Eq. (24), which is independent
of $J$. Dashed lines correspond to the exact synchronization boundaries,
obtained using Eq. (7) in Ref. Pietras  _et al._ (2019).

In the thermodynamic limit ($N\to\infty$), the dynamics of Eqs. (19) are
greatly simplified assuming $\chi_{i}$ in Eq. (18) to be Lorentzian-
distributed

$G(\chi)=\frac{\Delta/\pi}{\chi^{2}+\Delta^{2}},$

where $\Delta$ is the half-width of the distribution. Then, using the so-
called Ott-Antonsen (OA) ansatz Ott and Antonsen (2008), the KM Eqs. (19) can
be exactly reduced to a mean-field model consisting of two differential
equations for the complex Kuramoto order parameter,

$Z=Re^{i\psi}=\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}},$ (22)

in the limit $N\to\infty$. The mathematical approach to obtain the mean-field
equations corresponding to Eqs. (19) is a standard procedure. Here, we skip
the mathematical details and refer the reader to, for example Ref. Montbrió
and Pazó (2011), where the mean-field model corresponding to Eq. (19) was
derived in detail. Accordingly, using Eqs. (18,20,21), we obtain the mean-
field equations

$\displaystyle\dot{R}$ $\displaystyle=$
$\displaystyle\frac{R}{2\tau}\left(-\frac{2\Delta}{\sqrt{\bar{\eta}}}+g(1-R^{2})\right),$
(23a) $\displaystyle\dot{\psi}$ $\displaystyle=$
$\displaystyle\frac{2\sqrt{\bar{\eta}}}{\tau}+\frac{J}{2\pi\tau}(1-R^{2}),$
(23b)

which approximate the dynamics of the QIF model Eqs. (1) for small $g$ and
$J$. The radial equation Eq. (23a) shows that the incoherent state ($R=0$) is
a stable fixed point above the critical width,

$\Delta_{c}=g\sqrt{\bar{\eta}}/2,$ (24)

which is independent of chemical coupling, $J$. At $\Delta=\Delta_{c}$ a
stable nontrivial solution—corresponding to a partially synchronized
state—bifurcates from incoherence with

$R=\sqrt{(\Delta_{c}-\Delta)/\Delta_{c}},$

and frequency

$\Omega=\frac{2\sqrt{\bar{\eta}}}{\tau}+\frac{\Delta}{\tau\pi\sqrt{\bar{\eta}}}\frac{J}{g}.$
(25)

The solid line in Fig. 2 corresponds to the critical boundary Eq. (24), while
dashed lines correspond to the exact synchronization boundaries of the QIF
network for various degrees of inhibitory coupling—see Eq. (7) in Ref. Pietras
_et al._ , 2019. Note that for weak electrical coupling and/or weak
heterogeneity, all boundaries approach Eq. (24). Furthermore, for weak
heterogeneity, the frequency of the synchronized cluster Eq. (25) agrees with
Eq. (8) in Ref. Pietras  _et al._ (2019), which describes the frequency of the
oscillations near their onset.

In sum, these results confirm the validity of the Kuramoto model Eq. (19) as
an approximation of a population of heterogeneous QIF neurons with electrical
and chemical coupling, Eq. (1).

##  V Chimera states in coupled homogeneous populations of QIF neurons

To further illustrate the appropriateness of the KM to investigate the
dynamics of QIF networks, we investigate the presence of chimera states in
populations of QIF neurons. Our motivation is threefold:

  1. 1.

Chimera states were originally uncovered in a nonlocally coupled network of
identical Kuramoto oscillators Kuramoto and Battogtokh (2002). Given that the
phase dynamics Eqs. (19) are an approximation of Eqs. (1) valid for weak
heterogeneity and weak coupling, we expect QIF networks to display similar
chimera states, at least for weak coupling.

  2. 2.

Several papers have been devoted to investigate chimera states in networks of
spiking neurons; see, e.g., Refs. Sakaguchi, 2006; Olmi, Politi, and Torcini,
2010; Omelchenko  _et al._ , 2013; Vüllings  _et al._ , 2014; Hizanidis  _et
al._ , 2014; Panaggio and Abrams, 2015; Majhi  _et al._ , 2019; Laing, 2019;
Hizanidis  _et al._ , 2016; Majhi  _et al._ , 2019; Gerster  _et al._ , 2020;
Lucchetti, Jensen, and Heltberg, 2021. Some of them provide numerical evidence
that the presence of both chemical and electrical synapses favors the
emergence of chimera states Hizanidis  _et al._ (2016); Majhi  _et al._
(2019). Yet, the relation between chimera states in spiking neuron networks
with the original chimera states uncovered in the KM Kuramoto and Battogtokh
(2002); Montbrió, Kurths, and Blasius (2004); Abrams  _et al._ (2008); Laing
(2009) is lacking Bick  _et al._ (2020); Omel’chenko (2018); Haugland (2021).

  3. 3.

Recently, exact mean-field models for large populations of QIF neurons (often
called neural mass models, NMMs) with electrical and chemical synapses have
been put forward Laing (2015); Pietras  _et al._ (2019); Montbrió and Pazó
(2020); Byrne  _et al._ (2021). However, such NMMs have an important
limitation when neurons are—as in a chimera state—identical and fully
synchronized since both the mean membrane potential and the mean firing rate
diverge at the instant of collective firing 333 NMMs for QIF neurons are
derived after adopting the thermodynamic limit, $N\to\infty$, and under the
assumption $V_{r}\to-\infty$, and $V_{p}\to\infty$, see Montbrió and Pazó
(2020); Montbrió, Pazó, and Roxin (2015). If neurons are fully synchronized,
($V_{i}(t)=V_{j}(t),~{}\forall i,j$) the population of QIF neurons behaves as
a single neuron. Therefore, in QIF-NMM, the mean voltage Eq. (3) diverges when
all neurons fire a spike—in Refs. Laing (2015); Ermentrout (2006) the mean-
field variable $v$ is approximated to avoid this divergence, see also Pietras
_et al._ (2019). . This divergence is avoided using the averaging
approximation, and hence the KM for QIF neurons, Eq. (17), becomes singularly
suited to study collective behavior where neurons are fully synchronized.

Chimera states were originally uncovered in a ring of identical Kuramoto
oscillators with nonlocal coupling when $\alpha\lesssim\pi/2$ Kuramoto and
Battogtokh (2002). Shortly after their discovery, chimera states were also
found in a simpler setup, consisting of two populations of identical Kuramoto
oscillators Montbrió, Kurths, and Blasius (2004); Abrams  _et al._ (2008).
Here, to investigate chimera states in networks of QIF neurons, we adopt the
two-population setup of Refs. Montbrió, Kurths, and Blasius, 2004; Abrams  _et
al._ , 2008.

Specifically, we analyze the dynamics of two identical populations (labeled
$\sigma\in\\{1,2\\}$) of $n=N/2$ identical QIF neurons, interacting all-to-all
via both chemical and electrical synapses

$\tau\dot{V}_{i}^{\sigma}=\left(V_{i}^{\sigma}\right)^{2}+\bar{\eta}+I_{i,syn,s}^{\sigma}+I_{i,syn,c}^{\sigma},$
(26)

with the resetting rule of Eq. (1). Synaptic inputs have a contribution
$I_{i,syn,s}^{\sigma}$ due to self-interactions within each population
$\sigma$ and another contribution $I_{i,syn,c}^{\sigma}$ due to cross-
interactions of population $\sigma=\\{1,2\\}$ with population
$\sigma^{\prime}=\\{2,1\\}$,

$\displaystyle I_{i,syn,s}^{\sigma}$ $\displaystyle=$ $\displaystyle
g_{s}(v^{\sigma}-V_{i}^{\sigma})+J_{s}\tau r^{\sigma},$ $\displaystyle
I_{i,syn,c}^{\sigma}$ $\displaystyle=$ $\displaystyle
g_{c}(v^{\sigma^{\prime}}-V_{i}^{\sigma})+J_{c}\tau r^{\sigma^{\prime}}.$

Here, $v^{\sigma}$ and $r^{\sigma}$ are the mean membrane voltage and mean
firing rate of population $\sigma$, respectively. Using Eqs. (19,20,21), it is
straightforward to write the KM corresponding to Eqs. (26) as

$\displaystyle\dot{\theta}_{i}^{\sigma}=\omega$ $\displaystyle+$
$\displaystyle\frac{K_{s}}{n}\sum_{j=1}^{n}\left[\sin\left(\theta_{j}^{\sigma}-\theta_{i}^{\sigma}-\alpha_{s}\right)+\sin\alpha_{s}\right]$
(27) $\displaystyle+$
$\displaystyle\frac{K_{c}}{n}\sum_{j=1}^{n}\left[\sin\left(\theta_{j}^{\sigma^{\prime}}-\theta_{i}^{\sigma}-\alpha_{c}\right)+\sin\alpha_{c}\right],$

with $\omega=2\sqrt{\bar{\eta}}/\tau$ and

$\displaystyle K_{s,c}$ $\displaystyle=$
$\displaystyle\frac{1}{\tau}\sqrt{(J_{s,c}/\pi)^{2}+g_{s,c}^{2}},$ (28)
$\displaystyle\alpha_{s,c}$ $\displaystyle=$
$\displaystyle\arctan\left(\frac{J_{s,c}/\pi}{g_{s,c}}\right).$ (29)

The KM Eqs. (27) is slightly more general than the model originally
investigated in Refs. Montbrió, Kurths, and Blasius (2004); Abrams  _et al._
(2008); Laing (2009)—which considered $\alpha_{c}=\alpha_{s}$. In the QIF
network, this equality of the phase lag parameters implies that the ratios of
chemical to electrical coupling

$\rho_{s}=\frac{J_{s}/\pi}{g_{s}},\quad\rho_{c}=\frac{J_{c}/\pi}{g_{c}},$ (30)

are identical, $\rho_{s}=\rho_{c}.$ Recent work has also considered the
dynamics of chimera states in populations of Kuramoto oscillators with
distributed phase lags Martens, Bick, and Panaggio (2016); Choe, Ri, and Kim
(2016). Specifically, Martens et al. Martens, Bick, and Panaggio (2016)
investigated chimera states in the two-population model Eqs. (27).

###  V.1 Mean-field model

As we discussed previously, in the thermodynamic limit ($n=N/2\to\infty$) the
KM can be exactly reduced to a low-dimensional mean-field model using the OA
ansatz. In the case of the homogeneous, two-population Kuramoto model Eqs.
(27), the dynamics reduces to six ordinary differential equations using the
Watanabe-Strogatz ansatz Pikovsky and Rosenblum (2008); Watanabe and Strogatz
(1994). Assuming a particular set of initial conditions for the phases, the
system further reduces to four differential equations and it is described by
the OA ansatz Pikovsky and Rosenblum (2008); Abrams  _et al._ (2008). Such
mean-field equations describe the evolution of the complex Kuramoto order
parameters of the two populations,

$Z_{\sigma}=R_{\sigma}e^{i\psi_{\sigma}}=\frac{1}{n}\sum_{j=1}^{n}e^{i\theta^{\sigma}_{j}}.$
(31)

Using the mean-field analysis in Refs. Kawamura  _et al._ , 2010; Martens,
Bick, and Panaggio, 2016 and Eqs. (28,29,30), the mean-field equations for the
complex Kuramoto order parameters can be further reduced (by virtue of the
rotational symmetry of the KM) to the three dimensional system

$\displaystyle\frac{dR_{1}}{d\tilde{t}}$ $\displaystyle=$
$\displaystyle\frac{1-R_{1}^{2}}{2}\left[R_{1}+\frac{g_{c}}{g_{s}}R_{2}\cos\Psi-\rho_{s}\frac{J_{c}}{J_{s}}R_{2}\sin\Psi\right],$
(32a) $\displaystyle\frac{dR_{2}}{d\tilde{t}}$ $\displaystyle=$
$\displaystyle\frac{1-R_{2}^{2}}{2}\left[R_{2}+\frac{g_{c}}{g_{s}}R_{1}\cos\Psi+\rho_{s}\frac{J_{c}}{J_{s}}R_{1}\sin\Psi\right],$
(32b) $\displaystyle\frac{d\Psi}{d\tilde{t}}$ $\displaystyle=$
$\displaystyle\rho_{s}\frac{R_{1}^{2}-R_{2}^{2}}{2R_{1}R_{2}}\left(\frac{J_{c}}{J_{s}}\cos\Psi-
R_{1}R_{2}\right)-$ (32c)
$\displaystyle\frac{g_{c}}{g_{s}}\frac{R_{1}^{2}+R_{2}^{2}+2R_{1}^{2}R_{2}^{2}}{2R_{1}R_{2}}\sin\Psi,$

where the phase difference between the complex order parameters Eq. (31) is
defined as $\Psi=\psi_{1}-\psi_{2}$. In addition, we have rescaled time as
$\tilde{t}=g_{s}t/\tau$ so that the dynamics of Eqs. (27) depends only on
three combinations of parameters: the ratios of cross to self couplings
$g_{c}/g_{s}$ and $J_{c}/J_{s}$ and the ratio of chemical to electrical
coupling $\rho_{s}$ —see Eq. (30). In contrast, the original QIF model Eq.
(26) can, after appropriate rescaling, only be reduced to involve at least
four parameters.

###  V.2 Phase diagram of the mean-field model

Chimera states in two-population Kuramoto networks correspond to symmetry-
broken states where one of the populations is fully synchronized (i.e.
$R_{\sigma}=1$), while the other remains only partially synchronized
($R_{\sigma^{\prime}}<1$). In addition, chimera states in two-population
networks of identical Kuramoto oscillators coexist with the fully synchronized
state, $R_{1}=R_{2}=1$ Montbrió, Kurths, and Blasius (2004); Abrams  _et al._
(2008); Laing (2009).

Figure 3: Phase diagrams of the mean field model Eqs. (32), for two values of
the ratio $\rho_{s}$; see Eqs. (30). Shaded and hatched regions correspond to
regions of steady and unsteady stable chimera states, respectively. Red lines:
Saddle-node (SN) bifurcations. Blue lines: Hopf bifurcations. Green lines:
Homoclinic bifurcations. Filled circles: Takens-Bogdanov points. In panel (a),
the $+$ symbol corresponds to the coordinates used for the numerical
simulations of the QIF network depicted in Fig. 4(f,g,h): $J_{c}/J_{s}=0.15$,
$g_{c}/g_{s}=2.0$. In panel (b), symbols correspond to the coordinates used
for the numerical simulations of the QIF network depicted in Fig. 4(b,c,d) and
Fig. 5: $J_{c}/J_{s}=0.75$, $g_{c}/g_{s}=0.5$ ($+$ symbol) and in Fig. 6:
$J_{c}/J_{s}=0.5$, $g_{c}/g_{s}=0.5$ ($\times$ symbol).

To obtain the phase diagrams depicted in Fig. 3, we set $R_{2}=1$ in Eqs. (32)
and numerically continued Doedel  _et al._ (2007) chimera states using initial
conditions in their basin of attraction—see Fig. 4(a) and Ref. Abrams  _et
al._ , 2008. The diagrams show the regions where steady (shaded) and unsteady
(hatched) chimera states are stable for two different values of the ratio
$\rho_{s}$; see Eqs. (30). These regions lie between a saddle-node (red) and a
Homoclinic (green) bifurcation lines, which—together with a Hopf (blue)
bifurcation line separating steady and unsteady chimera states—meet at two
Takens-Bogdanov (TB) points. For decreasing $|\rho_{s}|$, the region of
chimera states shrinks and eventually disappears when the two TB points
collide.

The phase diagrams in Fig. 3 are qualitatively identical to that of Fig. 4a in
Ref. Martens, Bick, and Panaggio, 2016, but here the regions of chimeras are
represented in the parameter space of the QIF model 444Figures 4(b-f) in Ref.
Martens, Bick, and Panaggio (2016) show a bifurcation scenario with a
transcritical bifurcation that is not observed in Fig. 3. This bifurcation
occurs for $\alpha_{s}>\pi/2$, and these values of the phase lag parameter are
unreachable in the QIF model, where $\alpha_{s}\in[-\pi/2,0)$ for inhibitory
coupling and $\alpha_{s}\in(0,\pi/2]$ for excitatory coupling), see Eq. (29)..
This allows us to determine three necessary conditions for the existence of
chimera states in two-population networks of QIF neurons:

  1. 1.

Chimera states only exist in the presence of _both_ chemical and electrical
coupling.

  2. 2.

Self-chemical coupling needs to be much larger than self-electrical coupling,
$|J_{s}|\gg g_{s}$ —or, equivalently, $|\rho_{s}|\gg 1$. Using Eq. (29), this
implies that $\alpha_{s}$ is close to $\pm\pi/2$ in correspondence with
previous work Martens, Bick, and Panaggio (2016); Montbrió, Kurths, and
Blasius (2004); Abrams  _et al._ (2008).

  3. 3.

The modulus of self-chemical coupling needs to be larger than that of cross-
chemical coupling, $|J_{s}|>|J_{c}|$.

These three conditions are not sufficient conditions to have chimera states in
Eqs. (26) though. Indeed, Eq. (27) and the corresponding mean-field Eqs. (32)
are an approximation of the QIF Eqs. (26) for weak coupling, but it remains to
be seen whether chimera states persist in QIF networks when coupling strengths
become stronger. We numerically explore this issue in the next sections.

Figure 4: Basins of attraction of chimera states in the two-population KM
(dashed lines) and in the two-population QIF model (green-dotted regions).
Panels (a-d): Phase portraits of the (a) KM and (b-d) QIF model, with
$J_{c}/J_{s}=0.75$, $g_{c}/g_{s}=0.5,\rho_{s}=-40/\pi$, see $\times$ symbol in
Fig. 3(b). Panels (e-h): Phase portraits of the (e) KM and (f-h) QIF model,
with $J_{c}/J_{s}=0.15$, $g_{c}/g_{s}=2,\rho_{s}=-120/\pi$, see $+$ symbol in
Fig. 3(a). Symbols and lines in panels (a,e): Solid/open diamond:
Stable/unstable chimera states, respectively; solid dot: In-phase synchronized
state. Solid blue line: Unstable manifold of the saddle point. Dashed black
lines: Stable manifold of the saddle point, basin of attraction of stable
chimeras. The basins of attraction of chimeras have been transformed from the
$(R_{1}\cos\Psi,R_{1}\sin\Psi)$ coordinates (panels a,e) to the
$(r_{1},v_{1})$ coordinates (panels b-d and f-h) using Eq. (33). Green dots
corresponds to initial values leading to chimera states after $t=2500$ time
units in numerical simulations of two populations of $n=200$ QIF neurons.
Parameters: $\tau=1$ and $\eta=1$. Simulations using the Euler scheme with
time step: $dt=10^{-4}$ and symmetric resetting: $V_{p}=-V_{r}=1000$.

###  V.3 Chimera states in populations of QIF neurons

In the following, we numerically investigate the presence of chimera states in
the spiking neuron network model of QIF neurons Eqs. (26). First, we confirm
that, for weak coupling, chimeras are present in QIF networks and they exist
in the parameter range predicted by the phase diagrams of the KM, Fig. 3.
However, then we show that the basin of attraction of chimera states shrinks
as synaptic coupling strengths become stronger.

####  V.3.1 Dynamics of chimera states

Using the OA ansatz, the dynamics of the two-population model Eqs. (27) with
$N\to\infty$ can be exactly reduced to the three-dimensional system Eqs. (32).
In a chimera state we may set $R_{2}=1$ so that Eqs. (32) further reduce to a
planar system with variables $R_{1}$ and $\Psi$. Figures 4(a,e) show the
basins of attraction (dashed lines) of (a) a steady chimera state (solid
diamond symbol) and of (e) an unsteady chimera state (red limit cycle), with
parameters corresponding to $+$ symbols in Fig. 3. As mentioned previously,
chimera states coexist with the stable fully synchronized solution (solid dot
symbols), $Z_{1}=R_{1}=1,~{}\Psi=0$. The basin of attraction of the chimera
state is defined by the stable manifold of a saddle point (open diamond
symbols) Abrams  _et al._ (2008).

To set initial conditions leading to the chimera state of Fig. 4(a) in the
network of QIF neurons, we considered the initial condition $Z_{2}=1$ and
$Z_{1}=R_{1}e^{i\Psi}$ —with $R_{1}$ and $\Psi$ such that the system is in the
basin of attraction of the steady chimera state. Then we used the conformal
map Montbrió, Pazó, and Roxin (2015)

$\pi r_{\sigma}-iv_{\sigma}=\frac{1-Z_{\sigma}}{1+Z_{\sigma}},$ (33)

to transform the mean-field coordinates $Z_{2}$ and $Z_{1}$ into the mean
firing rate $r$ and the mean membrane potential $v$ of the populations of QIF
neurons—for population 2 we find $\pi r_{2}+iv_{2}=0$. Then, we initialized
the membrane voltages of the populations according to the formula
$V_{i}(0)=v+(\pi\tau r)\tan[\pi/2(2i-n-1)/(n+1)],$ for $i=1,\dots,n$ 555 For
$n\to\infty$, this corresponds to a Lorentzian distribution of voltages
$\rho(V)=\pi\tau r/[(V-v)^{2}+(\pi\tau r)^{2}]$; see Ref. Montbrió, Pazó, and
Roxin, 2015, and to uniformly distributed constants of motion in the Watanabe-
Strogatz theory Pikovsky and Rosenblum (2008); Laing (2018); Bick  _et al._
(2020) .

In Fig. 5, we show the results of a numerical simulation of the QIF network
($g_{s}=0.1$). The raster plot in Fig. 5(a) clearly shows the signature of a
chimera state: Neurons in population 1 (blue) are only partially synchronized,
while neurons in population 2 remain fully synchronized. The time evolution of
the firing rate $r_{1}$ and the mean-membrane potential $v_{1}$ for the
incoherent group are displayed in Figs. 5(b) and (c), respectively. These
collective variables indicate a periodic evolution of the incoherent
population, with fluctuations caused by the finite resetting of the QIF
neurons and by finite-size effects. Finally, Fig. 5(d) shows the Kuramoto
order parameter $R_{1}$ (blue) obtained using the time series $r_{1}(t)$ and
$v_{1}(t)$ and the conformal map Eq. (33). In contrast with the steady chimera
state in the mean-field Eqs.(32) (black dotted line), the chimera state in the
network of QIF neurons is not stationary but oscillates periodically in time.
The same unsteady chimeras arise in two-population networks of Winfree
oscillators Pazó and Montbrió (2014) and are the consequence of the lack of
rotational symmetry in the Winfree and QIF models.

In Fig. 6, we also explored how the unsteady chimera states in the Kuramoto
model Eqs. (27) translate to networks of QIF oscillators. To this aim we set
the parameters of the QIF model in the hatched region of the bifurcation
diagram in Fig. 3b ($\times$ symbol), and used the same initial conditions as
in the previous simulation. Here, we find a more complex chimera state that
seems to display macroscopic quasiperiodic dynamics 666We cannot discount that
the collective dynamics of the QIF Eqs.(26) is chaotic, as it also occurs in
coupled populations of Winfree oscillators Pazó and Montbrió (2014).: both the
firing rate and the mean-membrane potential of the incoherent population
oscillate with two characteristic frequencies as can be appreciated in Figs.
6(b,c). Again, the quasiperiodic chimera state in the QIF network corresponds
to a periodic chimera state in the Kuramoto model. In Fig. 5(d) we used Eq.
(33) to represent the Kuramoto order parameter for the QIF network (blue),
which roughly approximates the periodic dynamics of the corresponding Eqs.
(32).

Figure 5: Periodic chimera state in a two-population network of $N=5000$
identical inhibitory QIF neurons ($n=2500$ neurons in each population). (a)
Raster plot of 500 randomly chosen neurons. Neurons in population 1 (blue) are
partially synchronized, and neurons on population 2 are in-phase synchronized
(red). (b) Time series of the mean firing rate $r_{1}$ of population 1,
computed averaging the firing rate at each time step in time windows of
$\delta t=0.05$. (c) Time series of the mean membrane potential $v_{1}$ of
population 1. (d) Time series of the Kuramoto order parameter of population 1,
$R_{1}$, obtained from the mean-field quantities $r_{1},v_{1}$ using the
conformal map Eq. (33) (blue lines) and from direct integration of Eqs. (32)
(black dots). Parameters as in Fig. 4d (see also $+$ symbol in Fig.3b):
$J_{c}/J_{s}=0.75$, $g_{c}/g_{s}=0.5$, $\rho_{s}=-40/\pi$, and $g_{s}=0.1$.
Numerical simulations performed using the Euler scheme with time step:
$dt=10^{-4}$ and symmetric resetting: $V_{p}=-V_{r}=1000$.

Figure 6: Quasiperiodic chimera state in a two-population network of $N=5000$
identical inhibitory QIF neurons. The description of the panels and parameters
is as in Fig. 5 except $J_{c}/J_{s}=0.5$ (see also the $\times$ symbol in
Fig.3b).

####  V.3.2 Chimeras for strong coupling

The derivation of Eqs. (27) from the QIF network Eqs. (26) has been made under
the assumption of weak coupling. Yet, are chimera states in QIF networks
robust against stronger levels of coupling?

To investigate this issue, we used the conformal map Eq. (33) to express the
boundary of the basin of attraction of Fig. 4(a) in terms of $(r_{1},v_{1})$.
This transformed boundary is represented as a dashed line in Figs. 4(b-d).
Then, we performed three sets of numerical simulations of the QIF network for
increasing values of $g_{s}$ while keeping the ratios $g_{c}/g_{s}$,
$J_{c}/J_{s}$, and $\rho_{s}$ constant—note that this implies increasing all
the other coupling parameters.

For very small values of $g_{s}$, we expect the averaging approximation to
hold, and so the stability boundary of chimera states in the QIF model. The
green dotted region in Fig. 4(b) corresponds to the stability boundary of the
steady chimera state in the QIF network for $g_{s}=0.025$. The boundary
approximately agrees with that of the KM (dashed line), although the region is
slightly smaller in the QIF model. Notably, further increases in coupling
strength—see Figs. 4(c,d)—lead to a gradual reduction of the chimera’s basin
of attraction.

To further investigate the reduction of the stability boundary of chimeras in
the QIF model, in Figs. 4(f-h) we computed the basins of attraction of chimera
states in a different parameter regime—corresponding to the $+$ symbol in Fig.
3(a). Here $|\rho_{s}|$ is three times larger than in the previous case, and
hence, chemical couplings are three times stronger than in Figs. 4(b-d).
Correspondingly, the reduction of the basin of attraction of chimera states in
Figs. 4(f-h) is more pronounced than in Figs. 4(b-d). In fact, in Fig 4(h), we
find the complete disappearance of the region of stable chimera states in the
QIF network. This suggests that chimera states in QIF networks are only
observable for weak coupling.

##  VI Conclusions

In this paper we have applied a perturbative approach to simplify a weakly
heterogeneous population of QIF neurons, with weak all-to-all chemical and
electrical coupling Eqs. (1). This approach leads to a classical variant of
the Kuramoto model Eq. (19) Sakaguchi, Shinomoto, and Kuramoto (1988);
Sakaguchi and Kuramoto (1986), whose coupling parameters satisfy a simple
geometric relation with those of the QIF model 777A similar approach has been
applied to the Leaky Integrate-and-Fire model. In this case the approximated
phase model is not the KM, but it contains higher harmonics in the coupling
function Politi and Rosenblum (2015) Fig. 1.

The approximation of the QIF network by Eq. (19) allows one to use the
framework of the KM to investigate the role of chemical and electrical
synapses in setting up synchronization. For example, we find that in the
absence of electrical coupling the phase lag parameter of the KM is
$\alpha=\pm\pi/2$, which prohibits synchronization; see also Ref. Montbrió and
Pazó, 2018. Moreover, for Lorentzian distributions of currents the
synchronization threshold depends only on electrical coupling, Eq. (24),
whereas the oscillation frequency Eq. (25) is determined by the ratio of
chemical to electrical coupling. These results are in consonance with the
exact description provided by so-called neural mass (or firing rate) models
for networks of QIF neurons Pietras  _et al._ (2019); Montbrió and Pazó
(2020).

The framework of the KM allows for uncovering and investigating dynamical
states that are not reachable using neural mass models for QIF neurons
Montbrió, Pazó, and Roxin (2015); Laing (2015); Pietras  _et al._ (2019);
Montbrió and Pazó (2020). Here, we analyzed the case of chimera states in two-
population networks of identical QIF neurons Montbrió, Kurths, and Blasius
(2004); Abrams  _et al._ (2008). Despite the large number of studies devoted
to investigate chimera states in spiking neuron networks—see e.g. Olmi,
Politi, and Torcini (2010); Omelchenko  _et al._ (2013); Vüllings  _et al._
(2014); Hizanidis  _et al._ (2014); Panaggio and Abrams (2015); Majhi  _et
al._ (2019); Laing (2019); Hizanidis  _et al._ (2016); Majhi  _et al._ (2019);
Gerster  _et al._ (2020); Lucchetti, Jensen, and Heltberg (2021)—, the
relation between such states and the original chimera states uncovered in the
KM Kuramoto and Battogtokh (2002); Montbrió, Kurths, and Blasius (2004);
Abrams  _et al._ (2008) is lacking Bick  _et al._ (2020); Haugland (2021). We
showed that chimera states in QIF networks emerge in the presence of both
chemical and electrical couplings but only if chemical coupling is much
stronger than electrical coupling. However, our numerical results suggest that
chimeras in QIF networks are not robust against stronger levels of coupling.

Finally, we introduced a framework for the analysis of QIF networks that can
be readily applied to a variety of extensions of Eqs. (1). In particular, the
derivation of Eq. (19) does not impose constraints on the structure of the
network or the shape of the distribution of heterogeneities. Given that
network structure Sakaguchi, Shinomoto, and Kuramoto (1988); Blasius and
Tönjes (2005) and heterogeneities Pazó (2005); Lafuerza, Colet, and Toral
(2010); Omel’chenko and Wolfrum (2012); Martens  _et al._ (2009); Pazó and
Montbrió (2009); Pietras, Deschle, and Daffertshofer (2018) greatly affect the
dynamics of the KM, it may be interesting to investigate how this translates
to QIF networks.

###### Acknowledgements.

The authors thank Diego Pazó for helpful discussions. PC acknowledges
financial support from the European Union’s Horizon 2020 research and
innovation programme under grant agreement No 101017716 (Neurotwin). EM
acknowledges support by the Agencia Estatal de Investigación under the Project
No. PID2019-109918GB-I00.

## Data availability

Data sharing is not applicable to this article as no new data were created or
analyzed in this study.

## Author Declarations

The authors have no conflicts to disclose.

## References

  * Kuramoto (1975) Y. Kuramoto, “Self-entrainment of a population of coupled non-linear oscillators,” in  _International Symposium on Mathematical Problems in Theoretical Physics_ , Lecture Notes in Physics, Vol. 39, edited by H. Araki (Springer, Berlin, 1975) pp. 420–422.
  * Kuramoto (1984) Y. Kuramoto,  _Chemical Oscillations, Waves, and Turbulence_ (Springer-Verlag, Berlin, 1984).
  * Glomb  _et al._ (2021) K. Glomb, J. Cabral, A. Cattani, A. Mazzoni, A. Raj, and B. Franceschiello, “Computational models in electroencephalography,” Brain Topography , 1–20 (2021).
  * Lynn and Bassett (2019) C. W. Lynn and D. S. Bassett, “The physics of brain network structure, function and control,” Nature Reviews Physics 1, 318–332 (2019).
  * Breakspear, Heitmann, and Daffertshofer (2010) M. Breakspear, S. Heitmann, and A. Daffertshofer, “Generative models of cortical oscillations: neurobiological implications of the kuramoto model,” Frontiers in human neuroscience 4 (2010).
  * Cabral  _et al._ (2011) J. Cabral, E. Hugues, O. Sporns, and G. Deco, “Role of local network oscillations in resting-state functional connectivity,” NeuroImage 57, 130 – 139 (2011).
  * Villegas, Moretti, and Muñoz (2014) P. Villegas, P. Moretti, and M. A. Muñoz, “Frustrated hierarchical synchronization and emergent complexity in the human connectome network,” Sci. Rep. 4 (2014).
  * Ponce-Alvarez  _et al._ (2015) A. Ponce-Alvarez, G. Deco, P. Hagmann, G. L. Romani, D. Mantini, and M. Corbetta, “Resting-state temporal synchronization networks emerge from connectivity topology and heterogeneity,” PLoS Comput. Biol. 11, e1004100 (2015).
  * Sanz-Leon  _et al._ (2015) P. Sanz-Leon, S. A. Knock, A. Spiegler, and V. K. Jirsa, “Mathematical framework for large-scale brain network modeling in the virtual brain,” NeuroImage 111, 385 – 430 (2015).
  * Schmidt  _et al._ (2015) R. Schmidt, K. J. R. LaFleur, M. A. de Reus, L. H. van den Berg, and M. P. van den Heuvel, “Kuramoto model simulation of neural hubs and dynamic synchrony in the human cerebral connectome,” BMC Neuroscience 16, 54 (2015).
  * Petkoski  _et al._ (2016) S. Petkoski, A. Spiegler, T. Proix, P. Aram, J.-J. Temprado, and V. K. Jirsa, “Heterogeneity of time delays determines synchronization of coupled oscillators,” Phys. Rev. E 94, 012209 (2016).
  * Andrzejak  _et al._ (2016) R. G. Andrzejak, C. Rummel, F. Mormann, and K. Schindler, “All together now: Analogies between chimera state collapses and epileptic seizures,” Scientific reports 6, 1–10 (2016).
  * Petkoski, Palva, and Jirsa (2018) S. Petkoski, J. M. Palva, and V. K. Jirsa, “Phase-lags in large scale brain synchronization: Methodological considerations and in-silico analysis,” PLOS Computational Biology 14, 1–30 (2018).
  * Roberts  _et al._ (2019) J. A. Roberts, L. L. Gollo, R. G. Abeysuriya, G. Roberts, P. B. Mitchell, M. W. Woolrich, and M. Breakspear, “Metastable brain waves,” Nature communications 10, 1–17 (2019).
  * Choi and Mihalas (2019) H. Choi and S. Mihalas, “Synchronization dependent on spatial structures of a mesoscopic whole-brain network,” PLOS Computational Biology 15, 1–21 (2019).
  * Ziaeemehr  _et al._ (2020) A. Ziaeemehr, M. Zarei, A. Valizadeh, and C. R. Mirasso, “Frequency-dependent organization of the brain’s functional network through delayed-interactions,” Neural Networks 132, 155–165 (2020).
  * Noori  _et al._ (2020) R. Noori, D. Park, J. D. Griffiths, S. Bells, P. W. Frankland, D. Mabbott, and J. Lefebvre, “Activity-dependent myelination: A glial mechanism of oscillatory self-organization in large-scale brain networks,” Proceedings of the National Academy of Sciences 117, 13227–13237 (2020).
  * Jung, Eickhoff, and Popovych (2021) K. Jung, S. B. Eickhoff, and O. V. Popovych, “Tractography density affects whole-brain structural architecture and resting-state dynamical modeling,” NeuroImage 237, 118176 (2021).
  * Tewarie  _et al._ (2021) P. K. Tewarie, B. Prasse, J. M. Meier, Áine Byrne, M. D. Domenico, C. J. K. Stam, M. J. Brookes, A. Hillebrand, A. Daffertshofer, S. Coombes, and P. V. Mieghem, “Interlayer connectivity reconstruction for multilayer brain networks using phase oscillator models,” New Journal of Physics (2021).
  * Ódor, Kelling, and Deco (2021) G. Ódor, J. Kelling, and G. Deco, “The effect of noise on the synchronization dynamics of the Kuramoto model on a large human connectome graph,” Neurocomputing (2021).
  * Pang, Gollo, and Roberts (2021) J. C. Pang, L. L. Gollo, and J. A. Roberts, “Stochastic synchronization of dynamics on the human connectome,” NeuroImage 229, 117738 (2021).
  * Weerasinghe  _et al._ (2021) G. Weerasinghe, B. Duchet, C. Bick, and R. Bogacz, “Optimal closed-loop deep brain stimulation using multiple independently controlled contacts,” PLOS Computational Biology 17, 1–31 (2021).
  * Wang (2010) X.-J. Wang, “Neurophysiological and computational principles of cortical rhythms in cognition,” Physiological Reviews 90, 1195–1268 (2010).
  * Börgers (2017) C. Börgers,  _An introduction to modeling neuronal dynamics_ , Vol. 66 (Springer, 2017).
  * Ermentrout (1996) B. Ermentrout, “Type i membranes, phase resetting curves, and synchrony,” Neural Comp. 8, 979–1001 (1996).
  * Ermentrout and Kopell (1986) B. Ermentrout and N. Kopell, “Parabolic bursting in an excitable system coupled with a slow oscillation,” SIAM J. Appl. Math. 46, 233–253 (1986).
  * Izhikevich (2007) E. M. Izhikevich,  _Dynamical Systems in Neuroscience_ (The MIT Press, Cambridge, Massachusetts, 2007).
  * Sakaguchi, Shinomoto, and Kuramoto (1988) H. Sakaguchi, S. Shinomoto, and Y. Kuramoto, “Mutual Entrainment in Oscillator Lattices with Nonvariational Type Interaction,” Progress of Theoretical Physics 79, 1069–1079 (1988).
  * Sakaguchi and Kuramoto (1986) H. Sakaguchi and Y. Kuramoto, “A soluble active rotator model showing phase transitions via mutual entrainment,” Prog. Theor. Phys. 76, 576–581 (1986).
  * Pikovsky, Rosenblum, and Kurths (2001) A. S. Pikovsky, M. G. Rosenblum, and J. Kurths,  _Synchronization, a Universal Concept in Nonlinear Sciences_ (Cambridge University Press, Cambridge, 2001).
  * Strogatz (2000) S. H. Strogatz, “From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators,” Physica D 143, 1–20 (2000).
  * Pikovsky and Rosenblum (2015) A. Pikovsky and M. Rosenblum, “Dynamics of globally coupled oscillators: Progress and perspectives,” Chaos 25, 097616 (2015).
  * Wang and Buzsáki (1996) X.-J. Wang and G. Buzsáki, “Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model,” The Journal of Neuroscience 16, 6402–6413 (1996).
  * Wilson and Cowan (1972) H. R. Wilson and J. D. Cowan, “Excitatory and inhibitory interactions in localized populations of model neurons,” Biophys. J. 12, 1–24 (1972).
  * Nagy, Pereda, and Rash (2018) J. I. Nagy, A. E. Pereda, and J. E. Rash, “Electrical synapses in mammalian cns: Past eras, present focus and future directions,” Biochimica et Biophysica Acta (BBA) - Biomembranes 1860, 102 – 123 (2018).
  * Deco  _et al._ (2015) G. Deco, G. Tononi, M. Boly, and M. L. Kringelbach, “Rethinking segregation and integration: contributions of whole-brain modelling,” Nature Reviews Neuroscience 16, 430–439 (2015).
  * Pietras  _et al._ (2019) B. Pietras, F. Devalle, A. Roxin, A. Daffertshofer, and E. Montbrió, “Exact firing rate model reveals the differential effects of chemical versus electrical synapses in spiking networks,” Phys. Rev. E 100, 042412 (2019).
  * Montbrió and Pazó (2020) E. Montbrió and D. Pazó, “Exact mean-field theory explains the dual role of electrical synapses in collective synchronization,” Phys. Rev. Lett. 125, 248101 (2020).
  * Montbrió, Pazó, and Roxin (2015) E. Montbrió, D. Pazó, and A. Roxin, “Macroscopic description for networks of spiking neurons,” Phys. Rev. X 5, 021028 (2015).
  * Montbrió, Kurths, and Blasius (2004) E. Montbrió, J. Kurths, and B. Blasius, “Synchronization of two interacting populations of oscillators,” Phys. Rev. E 70, 056125 (2004).
  * Abrams  _et al._ (2008) D. M. Abrams, R. Mirollo, S. H. Strogatz, and D. A. Wiley, Phys. Rev. Lett. 101, 084103 (2008).
  * Laing (2009) C. R. Laing, “Chimera states in heterogeneous networks,” Chaos: An Interdisciplinary Journal of Nonlinear Science 19, 013113 (2009).
  * Martens, Bick, and Panaggio (2016) E. A. Martens, C. Bick, and M. J. Panaggio, “Chimera states in two populations with heterogeneous phase-lag,” Chaos: An Interdisciplinary Journal of Nonlinear Science 26, 094819 (2016).
  * Pazó and Montbrió (2014) D. Pazó and E. Montbrió, “Low-dimensional dynamics of populations of pulse-coupled oscillators,” Phys. Rev. X 4, 011009 (2014).
  * Montbrió and Pazó (2018) E. Montbrió and D. Pazó, “Kuramoto model for excitation-inhibition-based oscillations,” Phys. Rev. Lett. 120, 244101 (2018).
  * Note (1) More recently a similar derivation has been obtained for populations of heterogeneous Winfree oscillators with sinusoidal infinitesimal Phase Resetting Curves (iPRC) Pazó and Montbrió (2014). When the oscillators have the iPRC of the QIF neuron, the population of Winfree oscillators can be well approximated to to a population of QIF neurons with (weak) chemical synapses Montbrió and Pazó (2018).
  * Kopell and Ermentrout (2004) N. Kopell and B. Ermentrout, “Chemical and electrical synapses perform complementary roles in the synchronization of interneuronal networks,” Proceedings of the National Academy of Sciences 101, 15482–15487 (2004).
  * Laing (2015) C. R. Laing, “Exact neural fields incorporating gap junctions,” SIAM Journal on Applied Dynamical Systems 14, 1899–1929 (2015).
  * Hoppensteadt and Izhikevich (1997) F. C. Hoppensteadt and E. M. Izhikevich,  _Weakly connected neural networks._ (Spinger Verlag, N.Y., 1997).
  * Pietras and Daffertshofer (2019) B. Pietras and A. Daffertshofer, “Network dynamics of coupled oscillators and phase reduction techniques,” Physics Reports 819, 1–105 (2019).
  * Winfree (1967) A. T. Winfree, “Biological rhythms and the behavior of populations of coupled oscillators.” J. Theor. Biol. 16, 15–42 (1967).
  * Note (2) For $g=0$ the natural frequencies in Refs. Pazó and Montbrió (2014); Montbrió and Pazó (2018) differ from Eq. (18). The reason for this discrepancy is that Refs. Pazó and Montbrió (2014); Montbrió and Pazó (2018) consider the Winfree model with distributed natural frequencies $\Omega_{i}$, while here we study the QIF model with distributed currents $\eta_{i}$.
  * Blasius and Tönjes (2005) B. Blasius and R. Tönjes, “Quasiregular concentric waves in heterogeneous lattices of coupled oscillators,” Phys. Rev. Lett. 95, 084101 (2005).
  * Hansel, Mato, and Meunier (1995) D. Hansel, G. Mato, and C. Meunier, “Synchrony in excitatory neural networks,” Neural Comput. 7, 307–337 (1995).
  * Devalle, Roxin, and Montbrió (2017) F. Devalle, A. Roxin, and E. Montbrió, “Firing rate equations require a spike synchrony mechanism to correctly describe fast oscillations in inhibitory networks,” PLoS Computational Biology 13, 1–21 (2017).
  * Ott and Antonsen (2008) E. Ott and T. M. Antonsen, “Low dimensional behavior of large systems of globally coupled oscillators,” Chaos 18, 037113 (2008).
  * Montbrió and Pazó (2011) E. Montbrió and D. Pazó, “Shear diversity prevents collective synchronization,” Phys. Rev. Lett. 106, 254101 (2011).
  * Kuramoto and Battogtokh (2002) Y. Kuramoto and D. Battogtokh, “Coexistence of coherence and incoherence in nonlocally coupled phase oscillators,” Nonlinear Phenom. Complex Syst. 5, 380 (2002).
  * Sakaguchi (2006) H. Sakaguchi, “Instability of synchronized motion in nonlocally coupled neural oscillators,” Phys. Rev. E 73, 031907 (2006).
  * Olmi, Politi, and Torcini (2010) S. Olmi, A. Politi, and A. Torcini, “Collective chaos in pulse-coupled neural networks,” EPL (Europhys. Lett.) 92, 60007 (2010).
  * Omelchenko  _et al._ (2013) I. Omelchenko, O. E. Omel’chenko, P. Hövel, and E. Schöll, “When nonlocal coupling between oscillators becomes stronger: Patched synchrony or multichimera states,” Phys. Rev. Lett. 110, 224101 (2013).
  * Vüllings  _et al._ (2014) A. Vüllings, J. Hizanidis, I. Omelchenko, and P. Hövel, “Clustered chimera states in systems of type-i excitability,” New Journal of Physics 16, 123039 (2014).
  * Hizanidis  _et al._ (2014) J. Hizanidis, V. G. Kanas, A. Bezerianos, and T. Bountis, “Chimera states in networks of nonlocally coupled hindmarsh–rose neuron models,” International Journal of Bifurcation and Chaos 24, 1450030 (2014).
  * Panaggio and Abrams (2015) M. J. Panaggio and D. M. Abrams, “Chimera states: Coexistence of coherence and incoherence in networks of coupled oscillators,” Nonlinearity 28, R67 (2015).
  * Majhi  _et al._ (2019) S. Majhi, B. K. Bera, D. Ghosh, and M. Perc, “Chimera states in neuronal networks: A review,” Physics of Life Reviews 28, 100 – 121 (2019).
  * Laing (2019) C. R. Laing, “Dynamics and stability of chimera states in two coupled populations of oscillators,” Phys. Rev. E 100, 042211 (2019).
  * Hizanidis  _et al._ (2016) J. Hizanidis, N. E. Kouvaris, G. Zamora-López, A. Díaz-Guilera, and C. G. Antonopoulos, “Chimera-like states in modular neural networks,” Scientific reports 6, 19845 (2016).
  * Gerster  _et al._ (2020) M. Gerster, R. Berner, J. Sawicki, A. Zakharova, A. Škoch, J. Hlinka, K. Lehnertz, and E. Schöll, “Fitzhugh–nagumo oscillators on complex networks mimic epileptic-seizure-related synchronization phenomena,” Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 123130 (2020).
  * Lucchetti, Jensen, and Heltberg (2021) A. Lucchetti, M. H. Jensen, and M. L. Heltberg, “Emergence of chimera states in a neuronal model of delayed oscillators,” Phys. Rev. Research 3, 033041 (2021).
  * Bick  _et al._ (2020) C. Bick, M. Goodfellow, C. R. Laing, and E. A. Martens, “Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: a review,” The Journal of Mathematical Neuroscience 10, 1–43 (2020).
  * Omel’chenko (2018) O. E. Omel’chenko, “The mathematics behind chimera states,” Nonlinearity 31, R121 (2018).
  * Haugland (2021) S. W. Haugland, “The changing notion of chimera states, a critical review,” J. Phys. Complex. 2, 032001 (2021).
  * Byrne  _et al._ (2021) A. Byrne, J. Ross, R. Nicks, and S. Coombes, “Mean-field models for eeg/meg: from oscillations to waves,” Brain Topography , 1–18 (2021).
  * Note (3) NMMs for QIF neurons are derived after adopting the thermodynamic limit, $N\to\infty$, and under the assumption $V_{r}\to-\infty$, and $V_{p}\to\infty$, see Montbrió and Pazó (2020); Montbrió, Pazó, and Roxin (2015). If neurons are fully synchronized, ($V_{i}(t)=V_{j}(t),~{}\forall i,j$) the population of QIF neurons behaves as a single neuron. Therefore, in QIF-NMM, the mean voltage Eq. (3) diverges when all neurons fire a spike—in Refs. Laing (2015); Ermentrout (2006) the mean-field variable $v$ is approximated to avoid this divergence, see also Pietras  _et al._ (2019).
  * Choe, Ri, and Kim (2016) C.-U. Choe, J.-S. Ri, and R.-S. Kim, “Incoherent chimera and glassy states in coupled oscillators with frustrated interactions,” Phys. Rev. E 94, 032205 (2016).
  * Pikovsky and Rosenblum (2008) A. Pikovsky and M. Rosenblum, “Partially integrable dynamics of hierarchical populations of coupled oscillators,” Phys. Rev. Lett. 101, 264103 (2008).
  * Watanabe and Strogatz (1994) S. Watanabe and S. H. Strogatz, “Constant of motion for superconducting Josephson arrays,” Physica D 74, 197–253 (1994).
  * Kawamura  _et al._ (2010) Y. Kawamura, H. Nakao, K. Arai, H. Kori, and Y. Kuramoto, “Phase synchronization between collective rhythms of globally coupled oscillator groups: Noiseless nonidentical case,” Chaos: An Interdisciplinary Journal of Nonlinear Science 20, 043110 (2010).
  * Doedel  _et al._ (2007) E. J. Doedel, A. R. Champneys, F. Dercole, T. F. Fairgrieve, Y. A. Kuznetsov, B. Oldeman, R. Paffenroth, B. Sandstede, X. Wang, and C. Zhang, “Auto-07p: Continuation and bifurcation software for ordinary differential equations,” (2007).
  * Note (4) Figures 4(b-f) in Ref. Martens, Bick, and Panaggio (2016) show a bifurcation scenario with a transcritical bifurcation that is not observed in Fig. 3. This bifurcation occurs for $\alpha_{s}>\pi/2$, and these values of the phase lag parameter are unreachable in the QIF model, where $\alpha_{s}\in[-\pi/2,0)$ for inhibitory coupling and $\alpha_{s}\in(0,\pi/2]$ for excitatory coupling), see Eq. (29).
  * Note (5) For $n\to\infty$, this corresponds to a Lorentzian distribution of voltages $\rho(V)=\pi\tau r/[(V-v)^{2}+(\pi\tau r)^{2}]$; see Ref. MPR15, and to uniformly distributed constants of motion in the Watanabe-Strogatz theory Pikovsky and Rosenblum (2008); Laing (2018); Bick  _et al._ (2020).
  * Note (6) We cannot discount that the collective dynamics of the QIF Eqs.(26) is chaotic, as it also occurs in coupled populations of Winfree oscillators Pazó and Montbrió (2014).
  * Note (7) A similar approach has been applied to the Leaky Integrate-and-Fire model. In this case the approximated phase model is not the KM, but it contains higher harmonics in the coupling function Politi and Rosenblum (2015).
  * Pazó (2005) D. Pazó, “Thermodynamic limit of the first order phase transition in the Kuramoto model,” Phys. Rev. E 72, 046211 (2005).
  * Lafuerza, Colet, and Toral (2010) L. F. Lafuerza, P. Colet, and R. Toral, “Nonuniversal results induced by diversity distribution in coupled excitable systems,” Phys. Rev. Lett. 105, 084101 (2010).
  * Omel’chenko and Wolfrum (2012) O. E. Omel’chenko and M. Wolfrum, “Nonuniversal transitions to synchrony in the Sakaguchi-Kuramoto model,” Phys. Rev. Lett. 109, 164101 (2012).
  * Martens  _et al._ (2009) E. A. Martens, E. Barreto, S. H. Strogatz, E. Ott, P. So, and T. M. Antonsen, “Exact results for the kuramoto model with a bimodal frequency distribution,” Phys. Rev. E 79, 026204 (2009).
  * Pazó and Montbrió (2009) D. Pazó and E. Montbrió, “Existence of hysteresis in the kuramoto model with bimodal frequency distributions,” Phys. Rev. E 80, 046215 (2009).
  * Pietras, Deschle, and Daffertshofer (2018) B. Pietras, N. Deschle, and A. Daffertshofer, “First-order phase transitions in the kuramoto model with compact bimodal frequency distributions,” Phys. Rev. E 98, 062219 (2018).
  * Ermentrout (2006) B. Ermentrout, “Gap junctions destroy persistent states in excitatory networks,” Phys. Rev. E 74, 031918 (2006).
  * Laing (2018) C. R. Laing, “The dynamics of networks of identical theta neurons,” The Journal of Mathematical Neuroscience 8, 1–24 (2018).
  * Politi and Rosenblum (2015) A. Politi and M. Rosenblum, “Equivalence of phase-oscillator and integrate-and-fire models,” Phys. Rev. E 91, 042916 (2015).