]Equal contribution.]Equal contribution.

]Correponding authors.

]Correponding authors.

# Multi-band oscillations emerge from a simple spiking network

Tianyi Wu  [  School of Mathematical Sciences, Peking University, Beijing
100871, China;  Center for Quantitative Biology, Peking University, Beijing,
100871, China;  Yuhang Cai  [  Department of Mathematics, University of
California, Berkeley, CA, USA 94720;  Ruilin Zhang  Center for Quantitative
Biology, Peking University, Beijing, 100871, China;  Yuanpei College, Peking
University, Beijing, 100871, China;  Zhongyi Wang  School of Mathematical
Sciences, Peking University, Beijing 100871, China;  Center for Quantitative
Biology, Peking University, Beijing, 100871, China;  Louis Tao  [  Center for
Quantitative Biology, Peking University, Beijing, 100871, China;  Center for
Quantitative Biology, Peking University, Beijing, 100871, China;
<EMAIL_ADDRESS>Zhuo-Cheng Xiao  [  Courant Institute of
Mathematical Sciences, New York University, NY, USA 10003<EMAIL_ADDRESS>
###### Abstract

In the brain, coherent neuronal activities often appear simultaneously in
multiple frequency bands, e.g., as combinations of alpha (8-12 Hz), beta
(12.5-30 Hz), gamma (30-120 Hz) oscillations, among others. These population
rhythms are believed to underlie information processing and cognitive
functions and have been subjected to intense experimental and theoretical
scrutiny. Computational modeling has provided a framework for the emergence of
network-level oscillatory behavior from the interaction of spiking neurons.
However, because of the strong nonlinear interactions between highly recurrent
spiking populations, the dynamic interplay between neuronal rhythms in
multiple frequency bands has rarely been theoretically investigated. Many
studies invoke multiple physiological timescales (e.g., various ion channels
or multiple types of inhibitory neurons) or oscillatory inputs to produce
rhythms in multi-bands. Here we demonstrate the emergence of multi-band
oscillations in a simple network consisting of a single excitatory and a
single inhibitory neuronal population driven by constant input. First, we
construct a data-driven, Poincaré section theory for robust numerical
observations of single-frequency oscillations bifurcating into multiple bands.
Then we develop model reductions of the stochastic, nonlinear, high-
dimensional neuronal network to capture the appearance of multi-band dynamics
and the underlying bifurcations theoretically. Furthermore, when viewed within
the dimensionally-reduced state space, our analysis reveals conserved
geometrical features of the bifurcations on low-dimensional dynamical
manifolds. These results suggest a simple geometric mechanism behind the
emergence of multi-band oscillations without appealing to oscillatory inputs
or multiple synaptic or neuronal timescales. Thus our work points to
unexplored regimes of stochastic competition between excitation and inhibition
behind the generation of dynamic, patterned neuronal activities.

††preprint: Submit to CHAOS

> Coherent neural activities in the brain are well studied because of the
> belief that they are the dynamical substrate underlying information
> processing and cognitive functions. Frequently, these neural dynamical
> patterns appear as synchronous, stochastic, population-level spiking
> activities. Because of the strong, rapid and recurrent interactions between
> excitatory and inhibitory neurons, theoretical modeling of these coherent
> spiking patterns has faced enormous challenges. Here we study a dynamical
> regime where these synchronous activities appear as collective oscillations
> in multiple frequencies, even in the presence of temporally constant inputs.
> Making use of the appearance of nearly periodic dynamical spiking patterns,
> we combine classical Poincare theory with a data-driven approach to
> construct a stochastic bifurcation theory to explain the prominence of
> multi-band oscillations. We develop model reductions to deconstruct the
> underlying bifurcations. Furthermore, applying our recently developed
> dimensional reduction methods, these bifurcations can be viewed as conserved
> geometrical features on low-dimensional manifolds. In particular, we show
> that multi-band oscillatory activities can be described by a geometry of
> “leaves and roots,” where each “leaf” corresponds to a particular beat, and
> the “root” encapsulates the complex initial conditions that trigger each
> stochastic synchronous event.

##  I Introduction

Coherent neuronal activities are found in many brain regions, ranging from
visual and auditory cortices Gray  _et al._ (1989); Azouz and Gray (2000);
Logothetis  _et al._ (2001); Henrie and Shapley (2005); Brosch, Budinger, and
Scheich (2002), to frontal and parietal cortices Buschman and Miller (2007);
Gregoriou  _et al._ (2009); Siegel, Warden, and Miller (2009); Sohal  _et al._
(2009); Canolty  _et al._ (2010); Sigurdsson  _et al._ (2010); van Wingerden
_et al._ (2010); Pesaran  _et al._ (2002); Medendorp  _et al._ (2007), to the
hippocampus Bragin  _et al._ (1995); Csicsvari  _et al._ (2003); Colgin  _et
al._ (2009); Colgin (2016), the amygdala Popescu, Popa, and Paré (2009), and
the striatum van der Meer and Redish (2009), and are often manifested as
neuronal oscillations. Dynamically, these rhythms are likely to be a
reflection of the many possible interactions between excitation and
inhibition, between disparate synaptic and neuronal timescales, and between
local and long-range recurrent connectivities. It is also commonly believed
that they form a neurophysiological basis of sensory processing and cognitive
functions, such as attention Fries  _et al._ (2001, 2008), learning Bauer,
Paz, and Paré (2007), and memory Pesaran  _et al._ (2002). Therefore, neuronal
oscillations have been subjected to intense experimental and theoretical
examinations, ranging from the specific functions of these brain rhythms to
inferences of their roles in information processing and brain computations.

Much experimental evidence has correlated various population oscillations to
enhanced cognition or behavioral tasks. For instance, the enhancement of gamma
oscillations has been shown to sharpen visual feature detection in V1 Azouz
and Gray (2000, 2003); Frien  _et al._ (2000); Womelsdorf  _et al._ (2012) and
direction selectivity in MT Liu and Newsome (2006); Khayat, Niebergall, and
Martinez-Trujillo (2010). Attentional states are often accompanied by more
significant gamma-band synchronization Fries  _et al._ (2001); Engell and
McCarthy (2010); Benchenane, Tiesinga, and Battaglia (2011); Poch, Campo, and
Barnes (2014). Theta-band oscillations have been found to be correlated with
visual spatial memory Kaplan  _et al._ (2014); Herweg, Solomon, and Kahana
(2020); Goyal  _et al._ (2020); Vivekananda  _et al._ (2021), while
physiologically and pharmacologically-enhanced theta rhythms improve learning
and memory Sweet  _et al._ (2014); Stoiljkovic  _et al._ (2015); Roberts  _et
al._ (2018); Herweg, Solomon, and Kahana (2020). Alpha oscillations have been
linked to perceptual learning Sigala  _et al._ (2014); Poch, Campo, and Barnes
(2014); Bays  _et al._ (2015); Brickwedde, Krüger, and Dinse (2019); He  _et
al._ (2021), and recent work has suggested a role in coordinating activity
between brain regions Hindriks  _et al._ (2015); Herweg  _et al._ (2016);
Zhang  _et al._ (2018). Cortical beta rhythms have been correlated with
working memory, decision making, and sensorimotor activities Haegens  _et al._
(2011); Kilavik  _et al._ (2013); Tan, Wade, and Brown (2016); Schmidt  _et
al._ (2019).

Many studies, experimental, theoretical, and modeling, have focused on the
emergence of collective rhythms of a single frequency band. Numerical
simulations of large-scale models of the visual cortex have shown dynamics
consistent with the experimentally observed gamma rhythms Chariker, Shapley,
and Young (2018). Studies demonstrated that coherent neuronal oscillations can
be driven by external inputs or may emerge as internal states of the network
Rubin and Terman (2004); Kumar, Rotter, and Aertsen (2008). Theoretical
investigations have pointed to the importance of recurrent couplings between
interneurons and excitatory pyramidal cells and the role of synaptic
timescales Wang and Buzsáki (1996); Whittington  _et al._ (2000); Borgers and
Kopell (2003); Chariker and Young (2015); Bezaire  _et al._ (2016); Chariker,
Shapley, and Young (2018); Keeley  _et al._ (2019). Substantial theoretical
progress has been made by mapping the so-called theta neurons to systems of
phase oscillators, which, in the weakly coupled limit, are amenable to mean-
field reductions demonstrating how collective dynamics can lead to phase
synchronizationAshwin, Coombes, and Nicks (2016); Bick  _et al._ (2020).

Compared to single-frequency band neuronal rhythms, there has been much less
work on the mechanisms underlying multi-band oscillations. Theoretically,
studies have built upon single frequency synchronization by introducing
multiple synaptic or neuronal timescales or by the interaction between
multiple neuronal pools Bezaire  _et al._ (2016); Ter Wal and Tiesinga (2021);
Zonca and Holcman (2021). Researchers have also pointed out that different
frequency input forcing near the Hopf bifurcation of a Hodgkin-Huxley type
neuron can result in nested oscillations, similar to the theta-gamma nested
rhythms observed in the hippocampus and during rapid eye movements Jensen and
Colgin (2007); Segneri  _et al._ (2020).

Here we demonstrate that multi-band rhythms can emerge from a small neuronal
network consisting of a single excitatory population recurrently coupled with
a single inhibitory population, both driven by constant input. Recently, we
investigated how stochastic gamma oscillations form a low-dimensional manifold
in a dimension-reduced state space Cai  _et al._ (2021). Here we show that a
more detailed exploration of the dynamics reveals regions in parameter space
where single-frequency, stochastic oscillations break up into different
frequency bands, even in the presence of temporally stationary inputs.

In the next section, we introduce our leaky integrate-and-fire (LIF) neuronal
network model and explore the regime of repetitive, nearly periodic multiple-
firing events (MFEs) as a function of model parameters. Focusing first on the
repetition, we formulate and examine a data-driven Poincaré section theory.
The Poincaré sections allow us to map the iterative, nearly synchronous
coherent spiking events in a reduced, low-dimensional space. Our analysis
reveals that the multi-band oscillations emerge from different modes of
excitation-inhibition (E-I) competitions between neural populations, which are
sufficiently captured by membrane potential profiles.

##  II Results

Experimentally, neuronal oscillations are typically studied by
electrophysiological methods, and rhythmic oscillatory behavior is usually
observed in the recordings of local field potentialsRay and Maunsell (2015);
Henrie and Shapley (2005). Mathematically, neuronal oscillations were usually
modeled by temporally reoccurring activity patterns of a neuronal network.
Well-known theories include the interneuronal network gamma (ING Wang and
Buzsáki (1996); Livingstone (1996); Geisler, Brunel, and Wang (2005)), the
pyramidal interneuronal network gamma (PING Whittington  _et al._ (2000);
Borgers and Kopell (2003)), multiple-firing events (MFEs, Chariker and Young
(2015); Chariker, Shapley, and Young (2018)), among others. Although these
models include different levels of biological details and describe various
types of oscillations in different brain areas, they generally share two
common views: First, neuronal oscillations are reflected by reoccurring
synchronized spiking patterns generated by neural ensembles; Second, the
reoccurring spiking patterns are induced by the rapid, recurrent interactions
between neurons. Therefore, the informative aspects of neuronal oscillations
(including frequencies, amplitudes, and phases) can be represented by the
statistics of the spiking patterns. Furthermore, any study of the mechanism
rests on the understanding of the neuronal interactions within the ensemble.

This paper investigates multi-band neuronal oscillations emerging from the
dynamics of a simple spiking neuronal network containing only two neuronal
populations: excitatory (E) and inhibitory (I), which are recurrently coupled.
Oscillatory or any other temporally inhomogeneous inputs are absent. Following
previous work Zhang  _et al._ (2014), we focus on the so-called multiple-
firing events (MFEs). We first present a few different temporally repetitive,
synchronized spiking patterns found during our parameter space explorations.
Appearing with different amplitudes and phases, these reoccurring MFEs form
different beats robustly in many parts of parameter space and some have as
many as three oscillatory bands: alpha (8-12 Hz), beta (12.5-30 Hz), and gamma
(30-90 Hz).

To understand the emergence of different beats, we need to explain how the
repetitions and alternations of MFEs are shaped by the tight competition
between E/I populations. Here, we develop a novel data-driven Poincaré section
theory for the network dynamics. Most remarkably, our Poincaré section reveals
that the E-I competition is essentially captured by the membrane potential
profiles just before the occurrence of each MFE. Furthermore, by considering
only the most informative features of the membrane potential distributions, a
stochastic bifurcation analysis illustrates that

  1. 1.

Different beats are shaped by different balances of the E-I competition
arising from the bifurcations;

  2. 2.

The multi-band oscillations are prevalent in the high-dimensional parameter
space.

In addition, when presented in an appropriate dimensionally-reduced state
space, these MFE beats collectively form notable geometrical structures on the
principal manifolds. This novel representation of spiking network dynamics
reveals global features emerging from local bifurcations.

The above physics & data-informed model reduction provides a top-down view of
the emergence of multi-band oscillations from the detailed LIF model. To the
opposite, we also develop an autonomous ODE model from a bottom-up perspective
with abbreviations of most of the biophysical details, only preserving the
most crucial features for oscillations. Notably, it qualitatively reproduces
the MFE beat patterns and bifurcations by using only the low-order moments of
the membrane potential distributions. In all, the model reduction and ODE
model meet in the half way and exhibit surprisingly similar geometrical
features in the reduced phase space, confirming the tight relation between
multi-band oscillations and delicated competitions between E and I
populations.

Taken together, these results provide a theoretical framework to analyze the
complicated oscillatory dynamics produced by the rapid competition between the
excitatory and inhibitory neuronal populations.

Figure 1: A network model producing different oscillatory dynamics. A. Network
structure of the Leaky-Integrate-Fire network. B. Raster plots exhibiting
three different classical temporal repetitive patterns of MFEs: periods 1, 2,
and 3. Excitatory/inhibitory spikes are indicated by red/blue dots. C.
Temporally averaged frequency spectrums of the three MFE patterns in B. The
error bars indicate the standard error.

###  II.1 Synchronized firing in a small integrate-and-fire network

We numerically simulate a small, 400-neuron network that has been shown to
produce stochastic, oscillatory dynamics in the gamma band Chariker and Young
(2015); Li, Chariker, and Young (2019); Li and Xu (2019); Cai  _et al._
(2021). The parameter values were taken from our previous studies Cai  _et
al._ (2021) of an input layer (Layer 4C$\alpha$) of the macaque primary visual
cortex (V1). We will demonstrate that significantly different temporal
repetitions of spike patterns can be produced by varying the neuronal and
network parameters systematically near the parameter values that well-model
many V1 electrophysiological data.

####  II.1.1 A simple network model

The network contains two populations of spiking neurons: 300 excitatory (E)
cells and 100 inhibitory (I) cells. The network structure is assumed to be
spatially homogeneous, while the synaptic couplings and single neuron dynamics
are intentionally fixed near previous models Li, Chariker, and Young (2019);
Li and Xu (2019); Cai  _et al._ (2021). All neurons are driven by independent,
identically distributed Poisson spike trains. Here the strength of this
feedforward input maintains each neuron at a physiologically realistic firing
rate.

The network dynamics is determined by the membrane potential ($v$) of each
neuron, which is modeled by a linear, leaky integrate-and-fire (LIF) equation:

$\frac{{\rm d}{v}}{{\rm d}t}=-g^{L}(v-V_{r})-g^{E}(v-V^{E})-g^{I}(v-V^{I}).$
(1)

The membrane potential $v$ is a dimensionless variable ranging between the
reversal potential $V^{I}=-2/3$ and the spiking threshold $V_{th}=1$. $v$ is
driven towards $V_{th}$ by the excitatory current $g^{E}(v-V^{E})$, and away
from it through the leaky and inhibitory currents ($g^{L}(v-V_{r})$ and
$g^{I}(v-V^{I})$). In LIF dynamics, when the membrane potential $v$ arrives at
$V_{th}$, a spike is emitted by this neuron to all its postsynaptic cells,
concurrently $v$ is reset to the rest potential $V_{r}=0$, and then held there
for an absolute refractory period $\tau^{\textrm{ref}}$. The E and I
conductance terms are sums of Green’s functions of the spiking series received
by this neuron. The key features of the computation of the LIF neurons are
depicted in Fig. 1A and detailed in Methods.

####  II.1.2 Different MFE beats

The network described above has been shown to produce stochastic gamma-band
oscillations in many previous studies Li, Chariker, and Young (2019); Li and
Xu (2019); Cai  _et al._ (2021). Namely, spikes of E and I neurons cluster in
time as MFEs, which reoccur within timescales of roughly $\sim$40 Hz (see,
e.g., the Fig. 1B, upper panel). Despite the stochastic external drive,
previous work demonstrated that the robust recurrence of MFEs is due to the
synaptic timescales during the dynamical competition between the E and I
populations Chariker, Shapley, and Young (2018); Keeley  _et al._ (2019); Li,
Chariker, and Young (2019); Cai  _et al._ (2021). Specifically, an MFE is
triggered by the stochastic firing of a couple of E cells, leading to the
rapid recruitment of both E and I neuron firings. However, the I-current peaks
later than the E-current because of the longer timescale of the inhibitory
synapse (see Methods). Therefore, in typical network dynamics involving gamma
oscillations, both E and I populations initially undergo highly correlated
firing, then followed by a dominance of inhibition, resulting in the
termination of the MFE. After the inhibition wears off, the network again
enters an excitable state, allowing for the next occurrence of MFE. We refer
the readers to Chariker and Young (2015) for more details of this recurrent
excitation-inhibition mechanism of stochastic gamma oscillations.

However, the recurrence of MFE does not necessarily result in consecutively
similar spiking patterns. Instead, we find that MFEs may recur in different
beats and rhythms. During our parameter space investigation, the amplitudes of
MFEs (i.e., number of neurons involved) repetitively alter in different
fashions. We describe the three most prevalent repetitive patterns of MFEs
below.

The first pattern depicts nearly regular recurrence of MFEs in the gamma band,
where consecutive MFE amplitudes do not vary significantly (Fig. 1B, upper).
This classical 1-beat repetitive pattern has been well studied in previous
literature associated with gamma oscillations Chariker and Young (2015); Li,
Chariker, and Young (2019); Li and Xu (2019); Cai  _et al._ (2021). To the
opposite, consecutive MFEs can alternate with large and small amplitudes and
form multiple beats. For example, in a 2-beat rhythm, a strong MFE is usually
followed by a weak one and vice versa (Fig. 1B, middle). On the other hand, a
weak MFE may occurs after every two consecutive strong MFEs (Fig. 1B, bottom),
forming a 3-beat repetitive pattern. We remark that the three distinct MFE
rhythms can take place when varying only one physiological parameter (in Fig.
1B, the cortical coupling I-to-E weights $S^{EI}$: 1-beat, $S^{EI}=2.01$;
3-beat, $S^{EI}=2.07$; 2-beat, $S^{EI}=2.16$).

These repetitive patterns appear as different numbers of peaks in the temporal
frequency spectrum (Fig. 1C). Compared to the single peak $\sim$45 Hz for the
1-beat rhythm, the alternation of MFE amplitudes in the 2-beat rhythm
introduces an additional peak $\sim$25 Hz (high beta band), and the 3-beat
rhythm results in one more $\sim$15 Hz (low beta band). At the same time, the
gamma-band peak is shifted towards higher frequencies as the network recovers
more rapidly from the weaker inhibition induced by weaker MFEs.

Needless to say, the rhythms listed above was only a small subset of all
possible network dynamics. On the other hand, the patterns were the most
prominent ones during our parameter space exploration (via 1-dimensional
sweeps of a few critical parameters; see Sect II.2.2). Thus, in the remainder
of this paper, we focus on these highly robust and persistent MFE repetitive
patterns and leave more thorough investigations of the high-dimensional
parameter space to future work.

###  II.2 Multi-bands arise from iterations of MFEs

The regularity of MFE recurrence as different beats suggests the importance of
dynamical paths between consecutive MFEs. Specifically, how is the next MFE
induced by the E-I competition determined by the current MFE? We demonstrate
that the dynamical paths can be captured via an iterative scheme. The idea is
intuitively simple: If one could predict the next MFE given all the
information of the current MFE, say, from a mapping $F$, then the long-term
network dynamics are revealed by $F^{n}$, i.e., all subsequent MFEs are
predicted iteratively.

To address the mechanism behind different MFE beats, this section proposes a
data-driven Poincaré section theory regarding MFE iterations. We will perform
a top-down analysis: Starting from the high-dimensional network dynamics, we
deduce and simplify the MFE mapping $F$ by combining first principles
arguments and data-informed numerical simulations.

Figure 2: A theoretical framework to analyze the repetitive patterns of MFEs.
A. The raster plots are divided into MFEs (the accumulated dots) and inter-MFE
intervals (relatively quiescence periods) by a MFE-detection algorithm. The
initiation/termination of an MFE is labeled by a solid/dash green line. Two
different repetitive patterns are depicted in this panel to demonstrate the
effectiveness of the algorithm. B. The iteration from one MFE to another
equivalents to a stochastic, high-dimensional discrete dynamical system. C.
The two E/I population potential profile distributions assumed to be truncated
Gaussian before the initiation of MFEs. The distribution with a higher average
has a support whose upper bound intersects the firing threshold ($V_{th}=1$).
The difference between the averages is indicated ($m$).

####  II.2.1 Reduction of a high-dimensional MFE mapping

We start by dissecting the network dynamics onto a sequence of time intervals
$\\{I_{n}=[t_{n},t_{n+1})\\}$, where $t_{n}$ is the starting time of the
$n$-th MFE (see definition in Methods). Without justification, we first assume
that there exists an MFE mapping $F$ capturing the dynamical paths from one
MFE to another. Thus, the first natural question is: How $F$ predicts the
dynamical path on $I_{n+1}$ given the network dynamics on interval $I_{n}$?

Despite its complexity, the dynamical path on $I_{n}$ is determined by two
components:

  * •

The network state on time section $t_{n}$ (denoted as $\mathbb{M}_{n}$).
$\mathbb{M}_{n}$ is the initial condition of the dynamical path.

  * •

The realization of all stochasticity $\xi$ within interval $I_{n}$, including
noises in external drives and internal randomness in synaptic couplings.

Consequently, the mapping $F$ between dynamical paths is equivalent to a
stochastic Poincaré section mapping $F_{1}$ concerning the network states

$\mathbb{M}_{n+1}=F_{1}(\mathbb{M}_{n},\xi),$ (2)

The network state $\mathbb{M}_{n}$ includes membrane potentials and E/I
synaptic conductances $(v,g^{E},g^{I})$ of all neurons. Even for our small
400-neuron network, $\mathbb{M}_{n}$ contains at least $3\times 400=1200$
variables. Therefore, it would be hard to study $F_{1}$ without a drastic
dimension reduction. To make it worse, MFEs are produced by rapid, transient
competitions between E/I populations Chariker and Young (2015); Li, Chariker,
and Young (2019); Li and Xu (2019); Cai  _et al._ (2021). Therefore, each MFE
can be very sensitive to random events as small as a couple of synaptic
failures. In all, a first-principles-based analytic description of MFE
mappings is unfeasible to the best of our knowledge.

Instead, we propose to investigate $F_{1}$ in a data-informed manner: We first
demonstrate a two-level model reduction of the original mapping $F_{1}$, then
reconstruct the reduced version of $F_{1}$ from surrogate data.

_Step 1: Coarse-Graining._ Due to the homogeneous network structure, we
coarse-grain the spiking network by assuming the interchangeability between
neurons with similar potential profiles. Instead of concerning the triplet
$(v,g^{E},g^{I})$ for every neuron, the coarse-grained network state is
completely determined by the population statistics of the membrane potentials
($\rho^{E}(v)$, $\rho^{I}(v)$) and the total conductances $g^{EE}$, $g^{EI}$,
$g^{IE}$, $g^{II}$. Here, $g^{PQ}$ represents the sum of $Q$-conductances of
all type-$P$ neurons ($P,Q\in{E,I}$). Therefore, the mapping $F_{1}$ is
simplified to

$M_{n+1}=F_{2}(M_{n},\xi),$ (3)

where $M_{n}=(\rho^{E}(v),\rho^{I}(v),g^{EE},g^{EI},g^{IE},g^{II})$ at time
section $t_{n}$. Fig. 2B illustrates the schematics of the iteration in this
coarse-grained state space. We remark that, through Eq. 2 and 3, the E-I
competition during MFEs are effectively represented by the distributions of
E/I membrane potential profiles.

_Step 2: Gaussianity of $\rho^{E,I}(v)$._ The coarse-grained network state
$M_{n}$ still includes four variables (conductances) and two distributions (of
membrane potential profiles). To reduce further the complexity of MFE
mappings, we impose three additional assumptions to network states on time
section $t_{n}$.

  1. 1.

$g^{PQ}\approx 0$ at the beginning of an MFE ($P,Q\in\\{E,I\\}$) since the
durations of inter-MFE intervals are generally 3-6 times of the longest
synaptic timescale (see details in Methods);

  2. 2.

The membrane potential distributions $\rho^{E}(v)$ and $\rho^{I}(v)$ are well
approximated by truncated Gaussian distributions whose variances are dominated
by the noise of external input (Fig. 2C);

  3. 3.

There is at least one neuron whose membrane potential is near the firing
threshold $V_{\textrm{th}}$.

Assumptions 1&2 come from the lack of recurrent drives during the inter-MFE
intervals (Fig. 2A). As observed in the network simulations, the number of
spikes produced outside MFEs is small compared to spikes during MFEs, and
therefore, the effects of the network spikes produced during inter-MFE
intervals on potential distributions and conductances can be neglected.
Assumption 3 comes from the definition of $t_{n}$: it would take more time
after $t_{n}$ to trigger the first spike of MFE if no neuron had a potential
profile at $V_{\textrm{th}}$. According to the truncated Gaussian assumption
(assumption 2), at $t_{n}$, the upper bound of at least one support of
$\rho^{E,I}(v)$ should be $V_{\textrm{th}}$ (Fig. 2C).

Therefore, for the coarse-grained network state $M_{n}$, the information
provided by $\rho^{E,I}(v)$ equivalents to the positions of the two truncated
Gaussians. In other words, the only relevant statistics is the difference
between the averages of the two membrane potentials profiles:
$m_{n}=\mathbb{E}[v^{E}(t_{n})]-\mathbb{E}[v^{I}(t_{n})]$. In all, $F_{2}$ can
be further reduced to a 1-dimensional mapping $f$, such that

$m_{n+1}=f(m_{n},\xi),$ (4)

Eq. 4 is a drastic simplification compared to Eq. 3. We hope $m_{n}$ alone can
capture sufficient information to reproduce the observed MFE beats, whose
dynamics are dictated by the strong E-I competition during and after MFEs.

_Step 3: Reconstruct the reduced MFE mapping._ We now turn to reconstruct $f$
from surrogate data. (We denote the reconstructed version as $\hat{f}$.) For
any $m_{0}$ in the relevant domain, we perform a series of Monte Carlo
simulations of the transient network dynamics between the initiations of two
MFEs. Each simulation contains one independent realization of random paths
$\xi$ on $I_{0}=[t_{0},t_{1})$.

Specifically, for each Monte Carlo simulation, we generate membrane potentials
for each E/I neuron as the initial condition, then simulate the network
dynamics until the next MFE is detected. Since it is important to ensure that
an MFE is triggered at $t_{0}$, we sample membrane potentials randomly from
two truncated Gaussian distributions $\rho^{E}(v,t_{0})$ and
$\rho^{I}(v,t_{0})$ respectively, where
$m_{0}=\mathbb{E}[v^{E}(t_{0})]-\mathbb{E}[v^{I}(t_{0})]$, and the maximum of
membrane potentials is fixed as $V_{\textrm{th}}$. At the start of the second
MFE ($t_{1}$), we collect the empirical distributions of E/I population
potential profiles ($\hat{\rho}^{E}(v,t_{1})$ and $\hat{\rho}^{I}(v,t_{1})$).
Finally, the output of $\hat{f}$ from the Monte Carlo simulation is

$\displaystyle\hat{m}_{1}=\mathbb{E}[\hat{v}^{E}(t_{1})]-\mathbb{E}[\hat{v}^{I}(t_{1})],$
(5)

where $\hat{v}^{E,I}\sim\hat{\rho}^{E,I}(v,t_{1})$.

Taking together the outcomes of all Monte Carlo simulations, the output of the
stochastic MFE mapping $\hat{f}$ is an empirical distribution of
$\hat{m}_{1}$. For more details, see Methods.

Here we must point out that, without the reduction of Step 2, a much larger
ensemble of surrogate data would be necessary to produce an effective, low-
dimensional representation of the MFE mapping. Specifically, if the
construction goal was $F_{2}$, the corresponding Monte Carlo simulations
should cover the high-dimensional state space of $M_{0}$ to estimate the two
potential profile distributions and four total conductances. The surrogate
data set would be unfeasibly large.

We now reconstruct $\hat{f}$ aiming for the three MFE repetitive patterns
depicted in Fig. 1B. For each specific $m_{n}$, 24 Monte Carlo simulations are
performed, and the corresponding $\hat{m}_{n+1}$ values are indicated by blue
dots in the left plots in Fig. 3A-C. We find that the invariant sets of
mapping $\hat{f}$ (denoted $\chi_{\hat{f}}$) are clustered as 1, 2 or 3
subgroups (pink, green, and cyan circles, Fig. 3A-C left panels). Remarkably,
the number of subgroups predicts the beat number of MFE repetitive patterns
(Fig. 3A-C, upper right raster plots).

Moreover, the invariant sets $\chi_{\hat{f}}$ agree in detail with the network
dynamics. More precisely, for each parameter set, $\chi_{\hat{f}}$ fits the
membrane potentials E/I population at the initiation of each MFE. (Fig. 3 A-C,
lower right panels. Colored circles correspond to subgroups of
$\chi_{\hat{f}}$ indicated by the same colors). These results confirm that, in
the selected parameter regimes, the 1-dim mapping $\hat{f}$ successfully
captures the E-I competition and the underlying mechanism of MFE beats by only
considering the averages of potential profiles.

Furthermore, different MFE beats are well predicted by the geometrical
features of $\hat{f}$. We note that the lower branch of $\hat{f}$ (when
$m_{n}>0$) undergoes downward spatial translation when $S^{EI}$ increases
(Fig. 3 left panels. $S^{EI}$ values: $A<C<B$). Namely, when $m_{n}$ is fixed,
a stronger $S^{EI}$ coupling means inhibition produced by the MFEs has a more
significant impact on the E neurons. Therefore, the average of E-population
potential profiles $\mathbb{E}[\hat{v}^{E}]$ will be lower when the next MFE
is initiated, leading to a lower $m_{n+1}$. Collectively, the downward
translation of $\hat{f}$ explains the different MFE repetitive patterns: From
A to C, the only cluster (pink) loses stability and introduces two other
clusters (cyan and green) into $\chi_{\hat{f}}$ as its predecessor and
successor in the MFE iteration, giving rise to 3-beat in the oscillation; From
C to B, the further lowered branch exhibits no intersections with the diagonal
$m_{n+1}=m_{n}$ (red line), and the remaining two clusters form the 2-beat
rhythm oscillation.

Figure 3: Different repetitive patterns are explained by reconstructed MFE
mapping $\hat{f}(m_{n},\xi_{n})$. Panels A-C correspond to the three
repetitive patterns illustrated in Fig. 1B, generated by different choices of
$S_{EI}$. A. $S_{EI}=2.01$. Left: Stochastic mapping $\hat{f}(m_{n},\xi_{n})$
represented by Monte Carlo simulations (blue dots). Line $x=y$ (red diagonal
line) and the only subgroup of the invariant set $\chi_{\hat{f}}$ (pink
circle) also indicated. Right Up: Raster plot exhibiting MFEs. Right Bottom:
Time trace of $m(t)=\mathbb{E}[v^{I}(t)]-\mathbb{E}[v^{E}(t)]$. The first few
$m_{n}$s are indicated by pink circles at the initiation of each MFE. B.
$S_{EI}=2.16$. Same as A, but the $\chi_{\hat{f}}$ includes two subgroups
(green and pink circles). C. $S_{EI}=2.07$. $\chi_{\hat{f}}$ includes three
points (green, cyan, and pink circles) and indicates the 3-beat rhythm
repetitive pattern.

####  II.2.2 Multi-band oscillations caused by bifurcations

Figure 4: MFE beat bifurcation maps indicated by the distribution of
$\rm\Delta m$. A: Bifurcation map of $S^{{EI}}$. $S^{{EI}}$ values of the
three patterns in Fig. 3 are indicated by arrows. B-D: Bifurcation maps of
$S^{\rm{ext}}$, $1/p$ and $\tau^{\rm{ref}}$, where $S^{{EI}}$ is fixed as
indicated by the red line in panel A.

The connection between $\hat{f}$ and MFE beats suggests that multi-band neural
oscillations are induced when the MFE mapping undergoes critical bifurcations
in the parameter space. Here, we argue and examine this point using the
reconstructed $\hat{f}$ to investigate 1D sections of the parameter
landscapes. This subsection aims to demonstrate the prominence of multi-band
oscillations in the parameter space of the simple network.

Here, we focus on a handful of parameters selected from the following groups,
since a thorough investigation of high-dimensional parameter landscape is
computationally unfeasible.

  1. 1.

Cortical synaptic couplings. We choose $S^{EI}$, the I-to-E synaptic strength.

  2. 2.

External drive. We choose $S^{\textrm{ext}}$, the strength of each external
kick. The mean of external stimuli $S^{\textrm{ext}}\cdot\lambda^{\rm{ext}}$
is fixed during the variation of strength.

  3. 3.

Network architecture. In a homogeneous network, this is specified by $p$, he
connection probability $P$ between each pair of neurons. Likely, $P\cdot
S^{Q^{\prime}Q}$ is fixed for $Q,Q^{\prime}\in\\{E,I\\}$ to control the total
postsynaptic current induced by each spike.

  4. 4.

Neuronal parameters that govern single-cell dynamics. We choose
$\tau^{\textrm{ref}}$, the refractory period of a cell.

The four parameters listed above
$\\{S^{EI},S^{\textrm{ext}},p,\tau^{\textrm{ref}}\\}$ are varied individually
to investigate the 1-dimensional bifurcations of $\hat{f}$ . For each
parameter set involved, we reconstruct $\hat{f}$ and expect the invariant set
$\chi_{\hat{f}}$ to reflect the MFE beats shown in Sect II.2.1. In addition,
we propose a differential visualization of $\chi_{\hat{f}}$ by concerning

$\Delta m_{n}=m_{n+1}-m_{n}=(\hat{f}-\mathbbm{1})(m_{n})$

instead of $m_{n}$, where $\mathbbm{1}$ stands for the identity function. The
reason is that some subgroups of $\chi_{\hat{f}}$ may occupy similar locations
when displayed on the 1D space of $m_{n}$. In the differential visualization,
we use the subgroup number of $\chi_{\hat{f}-\mathbbm{1}}$ to represent the
subgroup number of $\chi_{\hat{f}}$ (Fig. 4). In all, the MFE beats
bifurcations are depicted by the heat maps of empirical distributions of
$\chi_{\hat{f}-\mathbbm{1}}$.

Fig. 4A illustrates a systematic investigation of MFE beats bifurcation for
different $S^{EI}$, while all other parameters are fixed (same as Fig. 3).
This confirms our previous finding of the “1-3-2" beat bifurcation when
$S^{EI}$ increases: $\chi_{\hat{f}-\mathbbm{1}}$ exhibits only one subgroup
for $S^{EI}<2.0$, two subgroups for $S^{EI}>2.1$, and three subgroups in the
transition period between them. Notably, the unstable subgroup of the 3-beat
rhythm (pink, Fig. 3C) is affected by the stochasticity nature of $\hat{f}$.
Therefore, instead of standing alone, the 3-beat rhythms is always mixed with
1-beat or 2-beat rhythms (see e.g., Fig. 3C raster plot, 2650-2800 ms). This
leads to an uneven distribution of mass within the three subgroups of
$\chi_{\hat{f}-\mathbbm{1}}$. For example, when $S^{EI}\approx 2.0$, though
two other subgroups already emerge (blue and green circles, Fig. 3C), the
subgroup around $(0,0)$ still carries most of the mass (pink circle), i.e.,
$\Delta m_{n}=m_{n+1}-m_{n}$ is mostly 0 for any consecutive two MFEs. As a
result, the network dynamics are dominated by 1-beat rhythms, with 3-beat
rhythms taking place occasionally. On the opposite, when $S^{EI}\approx 2.1$,
as indicated by the branches of $\Delta m_{n}$, the subgroup around $(0,0)$
almost vanishes, and the network dynamics is dominated by 2-beat rhythm
patterns instead.

The bifurcation map of $S^{EI}$ suggests that 1- and 2-beat rhythms are more
robust, whereas the 3-beat is a more transient phenomenon. Therefore, we next
focus on how other parameters may affect the bifurcations for 3-beat rhythms.
Three 1-dim parameter scan are carried out ($S^{\textrm{ext}}$, $p$, and
$\tau^{\textrm{ref}}$) by fixing $S^{EI}=2.08$ where 3-beat rhythm patterns
are prominent.

We first examine the bifurcations induced by the noise of external and
internal stimulus ($S^{\textrm{ext}}$ and $p$): The strength of external
stimuli $S^{\textrm{ext}}$ represents the size of external noise, whereas
projection probability between neurons $p$ measures the internal variability
of network dynamics. Remarkably, the network dynamics are dominated by 1-beat
when external and internal noises are large (high $S^{\textrm{ext}}$ and
$1/p$), and dominated by 3-beat rhythm patterns when noises are small (low
$S^{\textrm{ext}}$ and $1/p$, Fig. 4BC). These beat bifurcations cohere with
our intuition that 3-beat rhythm are disrupted by larger randomness.

The bifurcations induced by $\tau^{\textrm{ref}}$ are much more intriguing
(Fig. 4D). While dominated by 3-beat rhythm patterns when
$\tau^{\textrm{ref}}\sim 0$, the cluster $\Delta m\approx 0$ carries most mass
of $\chi_{\hat{f}-\mathbbm{1}}$ for $0.2<\tau^{\textrm{ref}}<0.8$ ms. Hence
MFEs reoccur in 1-beat rhythms due to $m_{n+1}\approx m_{n}$. When the
$\tau^{\textrm{ref}}$ becomes longer, the dynamics is dominated by 2-beat
rhythms for $0.8<\tau^{\textrm{ref}}<2.5$. More interestingly, 3-beat rhythms
come back robustly when the refractory becomes sufficiently long. Though it is
so far difficult to provide even a qualitative explanation of this
sophisticated bifurcation phenomena, we point out that the bifurcation points
of $\tau^{\textrm{ref}}$ locate at the same magnitudes of the synaptic
timescales (1-10 ms, see Koch (2004)). This intuition suggests the refractory
periods plays an essential role in the interplay of E-I competition, leading
to different MFE beats. In all, Fig. 4D illustrates a crucial dependence of
oscillation features on the interplay between physiological timescales, as we
have investigated in our previous study Cai  _et al._ (2021).

####  II.2.3 Geometrical features of the bifurcations

Figure 5: Multiband oscillations produced by network simulation are presented
in a reduced 3D space classified by the coarse-grained model in Cai  _et al._
(2021). Trajectories of different beats are colored according to values
$\Delta m$ at the initiations of each MFE. Weak beats are specially labeled
with deep blue color. A. Four points in the 1D parameter space of $S^{EI}$ are
selected. The corresponding network dynamics form a "root-leaf" structure, and
the number of leaf branches indicates the MFE beat number demonstrated by the
bifurcation map in Fig. 4A. B-D. Similar as A, but the varying parameters
correspond to Fig. 4B-D.

The MFE beat bifurcation lies at the center of our multi-band oscillation
theory and is worth more meticulous investigations. Here, we present a novel
geometrical representation of the bifurcations by projecting network MFE
dynamics on low-dimensional manifolds. In our previous study, we developed a
first-principle-based coarse-grain model consisting of the total E/I
conductances ($g^{E,I}$) and the numbers of subthreshold neurons near the
firing threshold $V_{th}$ (E and I “gate" neurons). Though drastically reduced
from the high-dimensional spiking network model, this coarse-grain model
successfully captured gamma oscillations. In the reduced state space, gamma
oscillations are restricted on a nearly 2D manifold except for the initiation
of each MFE.

To describe the geometrical features, we use an analogy of “root-leaf"
structure: “root" for the MFE initiations and “leaf" for trajectories of MFEs
and post-MFE dynamics represented on the 2D manifold. Specifically, at the end
of each inter-MFE interval, the inertial drives produced by recurrent synaptic
couplings mostly dissipate. Therefore, the network states at each MFE
initiation are more governed by the randomness from external stimuli,
composing a common area with large stochasticity and dimensionality (“roots").
The states of roots were found to dictate the following dynamics of MFEs. On
the other hand, despite the nonlinear and highly random evolution, MFEs and
post-MFE dynamics exhibited sufficiently low dimensional structures that the
trajectories are restricted as loops on a 2D manifold (“leaves"). We refer
readers to Cai  _et al._ (2021) for the rest of the details.

The 2D manifold depicting Gamma oscillations also captures multi-band
oscillations by its “root-leaf" structure. Notably, the MFE beats form a one-
to-one mapping to the different branches of leaves: When multiband
oscillations are illustrated in the same state space, we find that the number
of “leaf" clusters equals to the MFE beat number. Furthermore, the emergence
of multiple leaves is dictated by the MFE beat bifurcations (Fig. 5). To
visualize different branches of leaves, we dissect the trajectories of network
dynamics on $I_{n}$ (as shown in Sect. II.2.1), then color each section
according to $m_{n}$ (weak MFEs are specially colored in dark blue). On the
MFE manifold, similarly colored trajectories cluster together on the same
leaves, signaling an overall bifurcation structure.

For each MFE beat bifurcation map in Fig. 4, we select four points in the
scanned 1D parameter space and demonstrate network dynamics in the reduced
space. Taking $S^{EI}$ for example (Fig. 5A), recall the “1-3-2" pattern of
the bifurcation map. Correspondingly, for the 1-beat rhythm, the trajectories
of $S^{EI}=0.0196$ only form one leaf with light green ($m\approx 0$, Fig. 5A
first plot). In the second and third plots, the three-beat rhythms are
represented by the three separate leaves (dark blue, orange, and light green)
in the reduced space. For $S^{EI}=0.0214$, the 2-beat rhythm is echoed by the
merging of the two large leaves. Similar observations of $S^{\rm ext}$ and $p$
(Fig. 5BC) confirm the relation between network randomness and “leaf"
branches. Surprisingly, for different refractory periods, the 2D manifold
capturing network dynamics undergoes substantial changes (Fig. 5D). The
reduced space reveals that the single-neuron dynamics, especially the
parameters governing physiological timescales, play a crucial role in the
emergence of MFE beats and multi-band oscillations.

###  II.3 MFE beats captured by reduced ODE models

Figure 6: Comparison between the dynamics generated by the spiking neuronal
network and the reduced ODE model. Row 1-3 displays the MFE beats produced by
both models, where the parameter sets are chosen the same as Fig. 3. Column 1:
Raster plots of the spiking neuronal network. Column 2&3: Firing rates of both
models. Column 4&5 Differences between the averages of E/I potential profiles
($m$) of both models. Figure 7: The bifurcation map of MFE beats produced by
the ODE model. Only cortical synaptic coupling strength $S^{EI}$ is varied,
and this map is comparable to Fig. 4A.

We have been studying the MFE beats by a combination of first principles and
data-driven model reductions. Although the reduced mapping $f$ provides a
solid explanation for the emergence of multi-band oscillatory dynamics,
implicitly it also shares the common weaknesses of data-driven models.
Specifically, the soundness of the reduced mapping $f$ is dictated by two
factors: First, the validity of the model reduction hypothesis; and second,
the size and quality of the surrogate data from network simulations. The
former leads to a drastic reduction of the complicated E-I competition to a
scalar ($m_{n}$), requiring more careful investigation to make sure it is not
a “blind trust". For the latter, the reliability of $f$ is weakened by the
large fluctuation in some surrogate data. On the other hand, a large size of
high-quality surrogate data from network simulation can be computationally
expensive, and even call into question the meaning of the reduced model
itself.

To overcome the weakness mentioned above, we develop an autonomous ODE model
based only on the averages and variances of E/I potential profiles. The ODE
model is governed by the mean-field representations of E-I competitions
combined with population fluctuations. Opposed to the previous data-driven
model, this ODE model is completely based on first principles and free of any
surrogate data, hence independent from the numerical simulations of network
dynamics.

The key idea of the ODE model is to approximate the distributions of E/I
potential profiles with multiple (instead of one) Gaussian peaks to echo the
observation that E/I populations are only partially involved in the weak MFEs.
During the evolution of the model, the potential profiles will undergo
splitting and merging, which are caused by the threshold-crossing during MFE
and reactivation after MFE. Specifically, the potential profile of a neural
population would be split if a fraction of neurons fire during an MFE but the
others remain subthreshold. The potential profiles of the spiked neuron would
reset to $V_{r}$ and be modeled as an additional truncated Gaussian in the
probability density function. On the other hand, we assume two peaks would
merge into one if they get close enough. Each peak is described as a truncated
Gaussian due to the same reason as above (both E/I-neurons are primarily
driven by external Poissonian inputs in the absence of recurrent stimuli). In
all, the autonomous ODE model describe the E-I competition by the dynamics of
the splitting and merging of Gaussian peaks. We refer readers to Methods for
more details.

Crude as it is, the autonomous ODE model provides a remarkably good
approximation to MFE beats, especially when the variance of the each Gaussian
peak is small. We find that ODE model produces similar crucial dynamical
features as the MFE beats produced by the network (Fig. 6), including firing
rates, firing synchronicity, and oscillation frequency. Most importantly, the
ODE model also exhibits qualitatively similar dynamics of potential profile
averages ($\Delta m$). We selectively demonstrate the performances of the ODE
model for the three parameter sets displayed in Fig. 3. Furthermore, compared
to Fig. 4A, the ODE model generates similar bifurcation map of MFE beats when
the cortical synaptic coupling weights varies (Fig. 7).

Notably, the ODE model reveals a similar “root-leaf" structure and leaf
branching during MFE bifurcations (Fig. 8) when illustrated in the same
reduced space concerning E and I “gate" neurons and conductances. Since the
ODE model evolves the distributions of E/I potential profiles, one can deduce
the numbers of “gate" neurons from the numbers and locations of the Gaussian
peaks. Remarkably, not only the number of “leaves" matches the number of MFE
beats, but the ODE model trajectories shares similar geometrical features with
the representation of the network simulation (Fig. 5). This suggests that the
complex network dynamics can be successfully captured by considering only the
most important features of potential profiles, even though the splitting and
merging of Gaussian peaks in the ODE model are very coarse descriptions of the
former.

Taking together, the autonomous ODE model supports our main conclusions that

  1. 1.

Multi-band oscillation, i.e., the various MFE beats, emerge from the E-I
competition in different fashions;

  2. 2.

The E-I competition in network dynamics can be captured by the E/I potential
profiles, and further sufficiently represented by their important features.

Figure 8: Similar as Fig. 5A, multi-band oscillations produced by the
autonomous ODE model are presented in the same reduced 3D space.

##  III Discussion

The emerging picture of brain cognition is one in which coherent neural
activities form the physiological basis of many cognitive functions. Prominent
amongst the patterned population activities are the various oscillatory
rhythms which have been implicated in sensory, perceptions, attention, and
learning and memory. While the experimental investigations of coherent neural
activities have formed the theoretical underpinnings of information coding in
the brain Salinas and Sejnowski (2001); Colgin  _et al._ (2009); Quiroga and
Panzeri (2009); Quian Quiroga and Panzeri (2013); Sigala  _et al._ (2014);
Herweg  _et al._ (2016), the understanding of how neural circuits can generate
patterned activity remain incomplete. Previous theoretical studies have
largely focused on single-band oscillatory activities. In addition to the
theories for Gamma oscillations referred to earlier in this manuscript,
various theoretical studies the mechanism behind theta Varga  _et al._ (2008);
Hangya  _et al._ (2009); Bezaire  _et al._ (2016); Colgin (2016); Zhang  _et
al._ (2018); Segneri  _et al._ (2020), alpha Zhang  _et al._ (2018); Zonca and
Holcman (2021), beta oscillations Pavlides, Hogan, and Bogacz (2015); Tan,
Wade, and Brown (2016); Yu and Wang (2019), and others. By combining these
ideas with oscillatory inputs, studies have shown that, in general, multi-band
oscillations can emerge near Hopf bifurcations of the driven population Jensen
and Colgin (2007); Segneri  _et al._ (2020).

Here we take a different tack. Recent work have shown that a wide variety of
coherent activity can be generated by homogeneously structured, fluctuation-
driven neuronal networks, even under temporally constant drive. In particular,
experimental and simulational studies have demonstrated the ubiquity of self-
organized, nearly synchronous population spiking activity of variable size and
duration. Theoretically, the main obstacle to a full analysis is a clear
delineation of the strong stochastic competition between excitatory and
inhibitory sub-populations. Already, the mathematical description of the so-
called multiple-firing events needed the development of a specific ensemble
average, i.e., a ‘partitioned ensemble average’ that treated MFEs differently
from nearly homogeneous spiking patterns Zhang  _et al._ (2014); Zhang and
Rangan (2015) (which can be faithfully captured by classical coarse-graining
techniques, e.g, Fokker-Planck approaches Shao, Zhang, and Tao (2020)). In
recent work, we also uncovered a nearly two-dimensional manifold underlying
stochastic gamma oscillations by projecting onto a suitable, low-dimensional
state space Cai  _et al._ (2021).

Here we focus on a dynamical regime between single-band, stochastic
oscillatory rhythms and general heterogeneous patterned activities. During our
exploration of the dynamical possibilities of homogeneous networks that
generate MFEs, we found that many parameters led to the emergence of nearly
periodic patterns consisting of more than one single frequency. Because of the
prevalent periodicity, we found it natural to combine classical Poincaré
sectional theories with modern data-driven techniques. After a few steps of
model reductions, the complex MFE dynamics were captured by the neuronal
potential profiles on the time sections of MFE initiations, which were further
computed from a series of Monte Carlo simulations of transient network
dynamics. The Poincaré sectional theory revealed the crucial role of the
delicate E-I competition in MFE bifurcations, leading to multiple frequency
bands in network oscillations.

While a data-driven Poincaré return map was instrumental in identifying the
bifurcations behind the emergence of multi-band patterns, a more detailed
understanding of the roles played by various dynamical variables emerged in
our dimensional reduced modeling. Viewed in our framework, a simple geometry
of “leaves and roots" described well the multi-band oscillatory regime. When
displayed in a suitably reduced state space, the stochastic neuronal
oscillations coherently formed a nearly 2D manifold, on which a collection of
“leaves" (describing each MFEs and the post-MFE dynamics) and “roots"
(indicating the complex initial conditions triggering each MFE) delineate the
manifold of multi-band oscillations. Furthermore, based on these observations,
and through well-motivated assumptions, we managed to develop an ODE model
description that captured the essentials of the population membrane potential
dynamics. The multi-band oscillations emerging from spiking networks and the
ODE model both exhibits similar geometrical representations of MFE cycles,
with the number of MFE beats corresponds to the branches of “leaves" in the
state space.

The modern view of cognition as neuronal population dynamics has also
considered these neural network computations as trajectories on
manifoldsChurchland  _et al._ (2010); Mante  _et al._ (2013); Janak and Tye
(2015); Gallego  _et al._ (2017); Saxena and Cunningham (2019); Cueva  _et
al._ (2020). Recent dynamical theories of stochastic neuronal network dynamics
have progressed beyond f-I curves, static input-output relations, thresholded-
linear filters, and weakly nonlinear expansions Brunel and Hakim (1999); Cai
_et al._ (2006); Keeley  _et al._ (2019). Amongst the various neuronal
dynamical activities, periodic synchronization has received the most
attention, both experimentally and theoretically. We view our work here as
complementary to previous studies that focused on identifying mechanisms of
synchronization in heterogeneous networks Acebrón  _et al._ (2005); Montbrio,
Pazo, and Roxin (2015); Ashwin, Coombes, and Nicks (2016); Bick  _et al._
(2020). By identifying the population membrane potential distribution as the
major player in the generation of multi-frequency oscillations, we have
uncovered a novel dynamical route to multi-band oscillations in a simple,
homogeneous neuronal network, without appealing to oscillatory inputs or more
cell types and synaptic time-scales. The central role of the membrane
potential distribution points out the possibility of maintaining and
manipulating these multi-frequencied, synchronizing populations, which we hope
to investigate more fully in ongoing and future work.

## Methods

### Integrate-and-Fire Network

We consider an integrate-and-fire (I&F) neuronal network composed by $N_{E}$
excitatory neurons ($E$) and $N_{I}$ inhibitory neurons ($I$). For neuron $i$,
its membrane potential ($v_{i}$) is driven by excitatory and inhibitory (E/I)
synaptic currents:

$\begin{split}\frac{{\rm d}{v_{i}}}{{\rm
d}t}&=\left(g_{i}^{\textrm{ext}}+g_{i}^{E}\right)\cdot\left(V^{E}-v_{i}\right)+g_{i}^{I}\left(V^{I}-v_{i}\right),\\\
g_{i}^{\textrm{ext}}&=S_{i}^{\textrm{ext}}\sum_{\mu_{i}^{\textrm{ext}}}G^{E}(t-t_{\mu_{i}^{\textrm{ext}}}),\\\
g_{i}^{E}&=\sum_{\begin{subarray}{c}j\in E\\\ j\neq
i\end{subarray}}S_{ij}^{E}\sum_{\mu_{j}^{E}}G^{E}(t-t_{\mu_{j}^{E}}),\\\
g_{i}^{I}&=\sum_{\begin{subarray}{c}j\in I\\\ j\neq
i\end{subarray}}S_{ij}^{I}\sum_{\mu_{j}^{I}}G^{I}(t-t_{\mu_{j}^{I}}).\end{split}$
(6)

In Eq. 6, $g_{i}^{\\{\textrm{ext},E,I\\}}$ indicate the external, recurrent
excitatory and inhibitory conductances of neuron $i$, which are induced by
external stimulus ($\mu_{i}^{\textrm{ext}}$) and E/I-spikes produced by other
neurons in the network ($\mu_{j}^{\\{E,I\\}}$) respectively. The corresponding
strength of synaptic couplings are denoted by $S_{i}^{\textrm{ext}}$ and
$S_{ij}^{\\{E,I\\}}$. The postsynaptic conductances are considered as temporal
convolution of the spiking series via a Green’s function,

$\begin{split}G^{E}(t)&=\frac{1}{\tau^{E}}e^{-\nicefrac{{t}}{{\tau^{E}}}}h(t),\\\
G^{I}(t)&=\frac{1}{\tau^{I}}e^{-\nicefrac{{t}}{{\tau^{I}}}}h(t),\end{split}$
(7)

where $h(t)$ is the Heaviside function. The time constants,
$\tau^{\\{E,I\\}}$, model the synaptic timescales of the excitatory and
inhibitory synapses (such as AMPA and GABA Gerstner  _et al._ (2014)),
respectively.

Neuron $i$ fires when its membrane potential $v_{i}$ reaches the firing
threshold $V^{th}$, after which neuron $i$ immediately resets the rest
potential ($V^{r}$) and remains in a refractory state for a fixed time of
$\tau^{\textrm{ref}}$. It is conventional in many previous studies to use
reduced-dimensional units and choose $V^{th}=1$ and $V^{r}=0$ Koch (2004); Cai
_et al._ (2006). Accordingly, $V^{E}=\nicefrac{{14}}{{3}}$ and
$V^{I}=-\nicefrac{{2}}{{3}}$ are the excitatory and inhibitory reversal
potentials.

While Eq. (6) can model a network with arbitrary connectivity structure, in
this paper, we focus on a spatially homogeneous network architecture. That is,
whether a certain spike released by a neuron of type $Q$ is received by
another neuron of type $Q^{\prime}$ is only determined by an independent coin
flip with a probability $P^{Q^{\prime}Q}$, where $Q,Q^{\prime}\in\\{E,I\\}$.
Furthermore, $S_{ij}^{\\{E,I\\}}$ only depends on the subtypes of the
post-/pre-synaptic neuronal pairs, i.e., $S_{ij}=S^{Q^{\prime}Q}$ if $i$ is a
type-$Q^{\prime}$ and $j$ is a type-$Q$ neuron. We point out that this makes
any neurons sharing the same potential profile ($v$) interchangeable in the
evolution of transient network dynamics.

The choice of network parameters are summarized below:

  * •

Frequencies of external input: $\lambda^{E}=\lambda^{I}=7000$ Hz; Strength of
external stimuli: $S_{i}^{\textrm{ext}}=0.001$.

  * •

Strengths of cortical synaptic couplings: $S^{EE}=S^{EI}=S^{II}=0.02$, and
$S^{IE}=0.008$;

  * •

Probability of spike projections: $P^{EE}=0.15$, $P^{IE}=P^{EI}=0.5$, and
$P^{II}=0.4$;

  * •

Synaptic timescales: $\tau^{I}=4.5$ ms. $\tau^{E}=1.4$ ms to reflect the fact
that AMPA is faster than GABA.

### Power spectrum of neuronal oscillations.

The power spectrum density (PSD) measures of the variance in a signal as a
function of frequency. In this study, the PSD is computed as follows:

A time interval $(0,T)$ is divide into time bins $B_{n}=[(n-1)\Delta t,n\Delta
t],n=1,2,...$, the spike density $\mu_{n}$ per neuron in $B_{n}$ is given by
$\mu_{n}=\nicefrac{{m_{n}}}{{N\Delta t}}$ where $m_{n}$ is the total number of
spikes fired in bin $B_{n}$. Hence, the discrete Fourier transform of
$\\{\mu_{n}\\}$ on $(0,T)$ is given as:

$\displaystyle\hat{\mu}(k)=\frac{1}{\sqrt{T}}\sum_{n=1}^{\nicefrac{{T}}{{\Delta
t}}}\mu_{n}\Delta te^{-k\cdot(2\pi i)\cdot(n\Delta t)}.$ (8)

Finally, as a function of $k$, PSD is the “power" concentrated at frequency
$k$, i.e., $|\hat{\mu}(k)|^{2}$.

### Detection of MFE

MFE lacks a rigorous definition from previous studies. On the other hand, it
is always visually categorized as highly correlated spiking patterns on the
temporal rasterplot. In this study, we design an generic algorithm based on
the key idea that each MFE is always triggered by recurrent E-E spikes.

  1. 1.

When more than two E-E spikes shows up in a 2 ms time interval, we set the end
of the interval as the initiation of a MFE.

  2. 2.

After a MFE initiation, when less than two E-E spikes are detected in a 2 ms
time interval, we set the begining of the interval as the termination of the
MFE.

The two steps above temporally dissect the network dynamics into intervals for
MFEs and inter-MFE periods. To finalize the MFE detection, we further combine
two MFE intervals into one and discard the inter-MFE-interval between them, if
the latter is shorter than 1 ms.

It is known that the lengths of MFE and inter-MFE-intervals are related to the
net gain of synaptic current from external stimuli Cai  _et al._ (2021). On
the other hand, for all parameter regimes we investigated, the majority of
inter-MFE-intervals are longer than 10 ms. Although we do not claim it
necessarily coherent with all previous studies, we find this simple algorithm
captures the great majority of the MFEs in this study and is consistently
reliable during our primary investigations to the network parameter space.

### MFE iteration plots

To reconstruct the reduced MFE mapping $m_{1}=f(m_{0},\xi)$, we perform the
following Monte Carlo type simulations of transient network dynamics. For a
fixed value of $m_{0}$, we randomly draw membrane potentials of every E/I
neuron in the network, from two distributions
$D_{E}\sim\mathcal{N}(m_{E},\sigma^{2}_{E})$ and
$D_{I}\sim\mathcal{N}(m_{E}-m_{0},\sigma^{2}_{I})$. The variances of both
distributions are collected from all times sections of MFE initiations during
a 100s simulation of network dynamics. These potential profiles are treated as
initial conditions in the subsequent simulations of the network dynamics until
the detection of the second MFE, at which the E/I distributions of potential
profiles are collected. After that, we compute the difference between the
averages of E/I potential profiles ($\hat{m}_{1}$) as one output of the MFE
mapping $f$. 50 simulations are carried out for each $m_{0}$, contributing to
to the data points in the iteration map on the plane of $(m_{0},\hat{m}_{1})$.
The probabilistic output of $f(m_{0},\xi)$ is therefore estimated as unimodal
or multimodal Gaussians from the empirical distributions of $\hat{m}_{1}$,
producing a data-driven reconstruction of the MFE mapping as $\hat{f}$.

The invariant measures of mapping $\hat{f}$ and $\hat{f}-\mathbbm{1}$
($\mathbbm{1}$ stands for identity function) are denoted as $\chi_{\hat{f}}$
and $\chi_{\hat{f}-\mathbbm{1}}$, and computed by the average of last 500
iterations from a 1000-time iteration of $\hat{f}$.

### ODE model

From the LIF neuronal network model, our ODE model is a further coarse-
graining type reduction. First of all, due to the interchangeability between
neurons with the same potential profiles, the distributions of E/I neuronal
membrane potential profiles ($\rho^{E,I}$) are sufficient for the evolution of
network dynamics (rather than tracing $v_{i}$ for each neuron). Then, we
approximate the distribution of neuronal membrane potential profiles by sum of
a collection of truncated Gaussian peaks, and the ODE model is governed the
splitting and merging of the peaks. According to classical results
[CaiTaoetal06], the E/I potential profiles are well approximated by Gaussian
functions whose variances are dictated by the noise introduced by the external
and recurrent input during the absences of MFEs. Readers should note that the
following assumptions of the ODE model are designed to reproduce the MFE
dynamics in a qualitative and drastically simplified fashion, instead of
mimicking a rigorous evolution of a Fokker-Planck type equation concerning the
distributions of potential profiles.

_Splitting._ We first assume $\rho^{E,I}$ are both initially Gaussian which
are truncated at 3$\sigma$. Being stimulated by external drive, the membrain
potentials increases and their distributions $\rho^{E,I}$ may cross the
threshold the firing threshold $V^{th}$, hence be divided into two parts. We
further assume that the probability mass passing threshold are released from
$V^{r}$ as a new truncated Gaussian peak after the refractory period (and
neurons passing threshold forms the MFE). On the other hand, the remaining
probability mass can be pushed back from the threshold due to the inhibition
produced by the I-spikes involved in the MFE. For simplicity, the shape of the
distribution of the remaining probability mass is fixed during the
suppression, i.e., it remains as a truncated Gaussian peak minus the right
tails that already crossed the threshold. In all, $\rho^{E}$ and $\rho^{I}$
(which are originally unimodal) are divided into two parts after splitting: a
Gaussian peak without tail on the high side and another Gaussian peak released
from $V^{r}$. Splitting may occur for a number of times before the merging of
any two peaks (defined momentarily), producing multiple truncated Gaussian
peaks (with or without tails on the high side).

_Merging._ Two Gaussian peaks are merged manually when they share common
support. The original two Gaussian peaks are replaced by a new Gaussian peak
whose mean and variance are averaged from the two merged peaks.

Therefore, the evolution of the distributions of potential profiles is now
simplified to the averages and variances of (multiple) Gaussian peaks. For
E-neuron profiles, we use $D^{E}_{1},D^{E}_{2},\dots,D^{E}_{k^{E}(t)}$ to
denote those normal distributions. Here, $k^{E}(t)$ stands for the number of
peaks contained by the distribution and
$D^{E}_{i}\sim\mathcal{N}(m_{E,i},\sigma^{2}_{E,i})$ for
$i=1,2,\dots,k^{E}(t)$. The peaks are ranked in a descending order of their
averages. Note that the first peak $D^{E}_{1}$ may lack its right tail due to
threshold-crossing during MFEs. The probability mass of each peak are
$p^{E}_{1},p^{E}_{2},\dots,p^{E}_{k^{E}(t)}$. Similarly, we use
$D^{I}_{1},D^{I}_{2},\dots,D^{I}_{k^{I}(t)}$ and
$D^{I}_{i}\sim\mathcal{N}(m_{I,i},\sigma^{2}_{I,i})$ to denote the Gaussian
peaks of I-neuron profiles and $p^{I}_{1},p^{I}_{2},\dots,p^{I}_{k^{I}(t)}$ to
denote their probability mass.

During the absence of the spliting and merging, each peak is independently
driven by external and recurrent stimulus. Therefore, the evolution of
averages and variances of each peak are dipicted by the following ODEs:

$\displaystyle\left\\{\begin{split}\frac{{\rm d}{\sigma^{2}_{E,i}}}{{\rm
d}t}&=\lambda^{E}(S^{\textrm{ext}})^{2}-2\cdot\frac{\sigma^{2}_{E,i}}{V^{th}-V^{I}}\frac{S^{EI}}{\tau^{I}}g^{EI},\\\
\frac{{\rm d}{m_{E,i}}}{{\rm
d}t}&=\lambda^{E}S^{\textrm{ext}}+\frac{S^{EE}}{\tau^{E}}g^{EE}-\frac{S^{EI}}{\tau^{I}}\frac{m_{E,i}-V^{I}}{V^{th}-V^{I}}g^{EI},\\\
\end{split}\right.$ (9) $\displaystyle\left\\{\begin{split}\frac{{\rm
d}{\sigma^{2}_{I,i}}}{{\rm
d}t}&=\lambda^{I}(S^{\textrm{ext}})^{2}-2\cdot\frac{\sigma^{2}_{I,i}}{V^{th}-V^{I}}\frac{S^{II}}{\tau^{I}}g^{II},\\\
\frac{{\rm d}{m_{I,i}}}{{\rm
d}t}&=\lambda^{I}S^{\textrm{ext}}+\frac{S^{IE}}{\tau^{E}}g^{IE}-\frac{S^{II}}{\tau^{I}}\frac{m_{I,i}-V^{I}}{V^{th}-V^{I}}\cdot
g^{II},\end{split}\right.$

In Eq. 9, we assume that the variances of each peak are dominated by the noise
from external and recurrent inhibition. The term
$\frac{m_{Q,i}-V^{I}}{V^{th}-V^{I}}$ normalize the contribution of I current
driving of $m_{Q,i}$ due to the conductance-based setup in Eq. 6, where
$Q\in\\{E,I\\}$. The corresponding normalization term for E current is saved
for simplicity since $\frac{m_{Q,i}-V^{E}}{V^{th}-V^{E}}$ are closer to 1.

## Acknowledgments

This work was partially supported by the Natural Science Foundation of China
through grants 31771147 (T.W., R.Z., Z.W., L.T.) and 91232715 (L.T.). Z.X. is
supported by the Courant Institute of Mathematical Sciences through Courant
Instructorship.

We thank Lai-Sang Young and David McLaughlin (New York University) and Yao Li
(The University of Massachusetts Amherst) for their useful comments.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Data Availability Statement

Pending data and code available upon the acceptance of the manuscript.

## References

  * Gray  _et al._ (1989) C. Gray, P. König, A. Engel, and W. Singer, “Oscillatory responses in cat visual cortex exhibit inter-columnar synchronization which reflects global stimulus properties,” Nature 338, 334–337 (1989).
  * Azouz and Gray (2000) R. Azouz and C. Gray, “Dynamic spike threshold reveals a mechanism for synaptic coincidence detection in cortical neurons in vivo,” Proc. Natl. Acad. Sci. USA 97, 8110–8115 (2000).
  * Logothetis  _et al._ (2001) N. K. Logothetis, J. Pauls, M. Augath, T. Trinath, and A. Oeltermann, “Neurophysiological investigation of the basis of the fMRI signal,” Nature 412, 150–157 (2001).
  * Henrie and Shapley (2005) J. Henrie and R. Shapley, “LFP power spectra in V1 cortex: The graded effect of stimulus contrast,” J. Neurophysiol. 94, 479–490 (2005).
  * Brosch, Budinger, and Scheich (2002) M. Brosch, E. Budinger, and H. Scheich, “Stimulus-related gamma oscillations in primate auditory cortex,” J. Neurophysiol. 87, 2715–2725 (2002).
  * Buschman and Miller (2007) T. Buschman and E. Miller, “Top-down versus bottom-up control of attention in the prefrontal and posterior parietal cortices,” Science 315, 1860–1862 (2007).
  * Gregoriou  _et al._ (2009) G. Gregoriou, S. Gotts, H. Zhou, and R. Desimone, “High-frequency, long-range coupling between prefrontal and visual cortex during attention,” Science 324, 1207–1210 (2009).
  * Siegel, Warden, and Miller (2009) M. Siegel, M. Warden, and E. Miller, “Phase-dependent neuronal coding of objects in short-term memory,” Proc. Natl. Acad. Sci. USA 106, 21341–21346 (2009).
  * Sohal  _et al._ (2009) V. Sohal, F. Zhang, O. Yizhar, and K. Deisseroth, “Parvalbumin neurons and gamma rhythms enhance cortical circuit performance,” Nature 459, 698–702 (2009).
  * Canolty  _et al._ (2010) R. Canolty, K. Ganguly, S. Kennerley, C. Cadieu, K. Koepsell, J. Wallis, and J. Carmena, “Oscillatory phase coupling coordinates anatomically dispersed functional cell assemblies,” Proc. Natl. Acad. Sci. USA 107, 17356–17361 (2010).
  * Sigurdsson  _et al._ (2010) T. Sigurdsson, K. Stark, M. Karayiorgou, J. Gogos, and J. Gordon, “Impaired hippocampal-prefrontal synchrony in a genetic mouse model of schizophrenia,” Nature 464, 763–767 (2010).
  * van Wingerden  _et al._ (2010) M. van Wingerden, M. Vinck, J. V. Lankelma, and C. M. Pennartz, “Learning-associated gamma-band phase-locking of action–outcome selective neurons in orbitofrontal cortex,” Journal of Neuroscience 30, 10025–10038 (2010).
  * Pesaran  _et al._ (2002) B. Pesaran, J. Pezaris, M. Sahani, P. Mitra, and R. Andersen, “Temporal structure in neuronal activity during working memory in macaque parietal cortex,” Nat. Neurosci. 5, 805–811 (2002).
  * Medendorp  _et al._ (2007) W. Medendorp, G. Kramer, O. Jensen, R. Oostenveld, J. Schoffelen, and P. Fries, “Oscillatory activity in human parietal and occipital cortex shows hemispheric lateralization and memory effects in a delayed double-step saccade task,” Cereb. Cortex 17, 2364–2374 (2007).
  * Bragin  _et al._ (1995) A. Bragin, G. Jandó, Z. Nádasdy, J. Hetke, K. Wise, and G. Buzsáki, “Gamma ($40-100$ Hz) oscillation in the hippocampus of the behaving rat,” J. Neurosci. 15, 47–60 (1995).
  * Csicsvari  _et al._ (2003) J. Csicsvari, B. Jamieson, K. Wise, and G. Buzsáki, “Mechanisms of gamma oscillations in the hippocampus of the behaving rat,” Neuron 37, 311–322 (2003).
  * Colgin  _et al._ (2009) L. Colgin, T. Denninger, M. Fyhn, T. Hafting, T. Bonnevie, O. Jensen, M. Moser, and E. Moser, “Frequency of gamma oscillations routes flow of information in the hippocampus,” Nature 462, 75–78 (2009).
  * Colgin (2016) L. L. Colgin, “Rhythms of the hippocampal network,” Nat Rev Neurosci 17, 239–249 (2016).
  * Popescu, Popa, and Paré (2009) A. Popescu, D. Popa, and D. Paré, “Coherent gamma oscillations couple the amygdala and striatum during learning,” Nat. Neurosci. 12, 801–807 (2009).
  * van der Meer and Redish (2009) M. van der Meer and A. Redish, “Low and high gamma oscillations in rat ventral striatum have distinct relationships to behavior, reward, and spiking activity on a learned spatial decision task,” Front. Integr. Neurosci. 3, 9 (2009).
  * Fries  _et al._ (2001) P. Fries, J. Reynolds, A. Rorie, and R. Desimone, “Modulation of oscillatory neuronal synchronization by selective visual attention,” Science 291, 1560–1563 (2001).
  * Fries  _et al._ (2008) P. Fries, T. Womelsdorf, R. Oostenveld, and R. Desimone, “The effects of visual stimulation and selective visual attention on rhythmic neuronal synchronization in macaque area V4,” J. Neurosci. 28, 4823–4835 (2008).
  * Bauer, Paz, and Paré (2007) E. P. Bauer, R. Paz, and D. Paré, “Gamma oscillations coordinate amygdalo-rhinal interactions during learning,” Journal of Neuroscience 27, 9369–9379 (2007).
  * Azouz and Gray (2003) R. Azouz and C. Gray, “Adaptive coincidence detection and dynamic gain control in visual cortical neurons in vivo,” Neuron 37, 513–523 (2003).
  * Frien  _et al._ (2000) A. Frien, R. Eckhorn, R. Bauer, T. Woelbern, and A. Gabriel, “Fast oscillations display sharper orientation tuning than slower components of the same recordings in striate cortex of the awake monkey,” Eur. J. Neurosci. 12, 1453–1465 (2000).
  * Womelsdorf  _et al._ (2012) T. Womelsdorf, J. Schoffelen, R. Oostenveld, W. Singer, R. Desimone, A. Engel, and P. Fries, “Orientation selectivity and noise correlation in awake monkey area V1 are modulated by the gamma cycle,” Proc. Natl. Acad. Sci. USA 109, 4302–4307 (2012).
  * Liu and Newsome (2006) J. Liu and W. Newsome, “Local field potential in cortical area MT: Stimulus tuning and behavioral correlations,” J. Neurosci. 26, 7779–7790 (2006).
  * Khayat, Niebergall, and Martinez-Trujillo (2010) P. S. Khayat, R. Niebergall, and J. C. Martinez-Trujillo, “Frequency-dependent attentional modulation of local field potential signals in macaque area mt,” Journal of neuroscience 30, 7037–7048 (2010).
  * Engell and McCarthy (2010) A. D. Engell and G. McCarthy, “Selective attention modulates face-specific induced gamma oscillations recorded from ventral occipitotemporal cortex,” Journal of Neuroscience 30, 8780–8786 (2010).
  * Benchenane, Tiesinga, and Battaglia (2011) K. Benchenane, P. H. Tiesinga, and F. P. Battaglia, “Oscillations in the prefrontal cortex: a gateway to memory and attention,” Current opinion in neurobiology 21, 475–485 (2011).
  * Poch, Campo, and Barnes (2014) C. Poch, P. Campo, and G. R. Barnes, “Modulation of alpha and gamma oscillations related to retrospectively orienting attention within working memory,” European Journal of Neuroscience 40, 2399–2405 (2014).
  * Kaplan  _et al._ (2014) R. Kaplan, D. Bush, M. Bonnefond, P. A. Bandettini, G. R. Barnes, C. F. Doeller, and N. Burgess, “Medial prefrontal theta phase coupling during spatial memory retrieval,” Hippocampus 24, 656–665 (2014).
  * Herweg, Solomon, and Kahana (2020) N. A. Herweg, E. A. Solomon, and M. J. Kahana, “Theta oscillations in human memory,” Trends in cognitive sciences 24, 208–227 (2020).
  * Goyal  _et al._ (2020) A. Goyal, J. Miller, S. E. Qasim, A. J. Watrous, H. Zhang, J. M. Stein, C. S. Inman, R. E. Gross, J. T. Willie, B. Lega, _et al._ , “Functionally distinct high and low theta oscillations in the human hippocampus,” Nature communications 11, 1–10 (2020).
  * Vivekananda  _et al._ (2021) U. Vivekananda, D. Bush, J. A. Bisby, S. Baxendale, R. Rodionov, B. Diehl, F. A. Chowdhury, A. W. McEvoy, A. Miserocchi, M. C. Walker, _et al._ , “Theta power and theta-gamma coupling support long-term spatial memory retrieval,” Hippocampus 31, 213–220 (2021).
  * Sweet  _et al._ (2014) J. A. Sweet, K. C. Eakin, C. N. Munyon, and J. P. Miller, “Improved learning and memory with theta-burst stimulation of the fornix in rat model of traumatic brain injury,” Hippocampus 24, 1592–1600 (2014).
  * Stoiljkovic  _et al._ (2015) M. Stoiljkovic, C. Kelley, D. Nagy, and M. Hajós, “Modulation of hippocampal neuronal network oscillations by $\alpha$7 nach receptors,” Biochemical Pharmacology 97, 445–453 (2015).
  * Roberts  _et al._ (2018) B. M. Roberts, A. Clarke, R. J. Addante, and C. Ranganath, “Entrainment enhances theta oscillations and improves episodic memory,” Cognitive neuroscience 9, 181–193 (2018).
  * Sigala  _et al._ (2014) R. Sigala, S. Haufe, D. Roy, H. R. Dinse, and P. Ritter, “The role of alpha-rhythm states in perceptual learning: insights from experiments and computational models,” Frontiers in computational neuroscience 8, 36 (2014).
  * Bays  _et al._ (2015) B. C. Bays, K. M. Visscher, C. C. Le Dantec, and A. R. Seitz, “Alpha-band eeg activity in perceptual learning,” Journal of vision 15, 7–7 (2015).
  * Brickwedde, Krüger, and Dinse (2019) M. Brickwedde, M. C. Krüger, and H. R. Dinse, “Somatosensory alpha oscillations gate perceptual learning efficiency,” Nature communications 10, 1–9 (2019).
  * He  _et al._ (2021) Q. He, B. Gong, K. Bi, and F. Fang, “The causal role of transcranial alternating current stimulation at alpha frequency in boosting visual perceptual learning,” bioRxiv (2021).
  * Hindriks  _et al._ (2015) R. Hindriks, M. Woolrich, H. Luckhoo, M. Joensson, H. Mohseni, M. L. Kringelbach, and G. Deco, “Role of white-matter pathways in coordinating alpha oscillations in resting visual cortex,” Neuroimage 106, 328–339 (2015).
  * Herweg  _et al._ (2016) N. A. Herweg, T. Apitz, G. Leicht, C. Mulert, L. Fuentemilla, and N. Bunzeck, “Theta-alpha oscillations bind the hippocampus, prefrontal cortex, and striatum during recollection: evidence from simultaneous eeg–fmri,” Journal of Neuroscience 36, 3579–3587 (2016).
  * Zhang  _et al._ (2018) H. Zhang, A. J. Watrous, A. Patel, and J. Jacobs, “Theta and alpha oscillations are traveling waves in the human neocortex,” Neuron 98, 1269–1281 (2018).
  * Haegens  _et al._ (2011) S. Haegens, V. Nácher, A. Hernández, R. Luna, O. Jensen, and R. Romo, “Beta oscillations in the monkey sensorimotor network reflect somatosensory decision making,” Proceedings of the National Academy of Sciences 108, 10708–10713 (2011).
  * Kilavik  _et al._ (2013) B. E. Kilavik, M. Zaepffel, A. Brovelli, W. A. MacKay, and A. Riehle, “The ups and downs of beta oscillations in sensorimotor cortex,” Experimental neurology 245, 15–26 (2013).
  * Tan, Wade, and Brown (2016) H. Tan, C. Wade, and P. Brown, “Post-movement beta activity in sensorimotor cortex indexes confidence in the estimations from internal models,” Journal of Neuroscience 36, 1516–1528 (2016).
  * Schmidt  _et al._ (2019) R. Schmidt, M. H. Ruiz, B. E. Kilavik, M. Lundqvist, P. A. Starr, and A. R. Aron, “Beta oscillations in working memory, executive control of movement and thought, and sensorimotor function,” Journal of Neuroscience 39, 8231–8238 (2019).
  * Chariker, Shapley, and Young (2018) L. Chariker, R. Shapley, and L.-S. Young, “Rhythm and synchrony in a cortical network model,” Journal of Neuroscience 38, 8621–8634 (2018).
  * Rubin and Terman (2004) J. Rubin and D. Terman, “High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model,” J. Comp. Neurosci. 16, 211–235 (2004).
  * Kumar, Rotter, and Aertsen (2008) A. Kumar, S. Rotter, and A. Aertsen, “Conditions for propagating synchronous spiking and asynchronous firing rates in a cortical network model,” J. Neurosci. 28, 5268–5280 (2008).
  * Wang and Buzsáki (1996) X.-J. Wang and G. Buzsáki, “Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model,” Journal of neuroscience 16, 6402–6413 (1996).
  * Whittington  _et al._ (2000) M. A. Whittington, R. Traub, N. Kopell, B. Ermentrout, and E. Buhl, “Inhibition-based rhythms: experimental and mathematical observations on network dynamics,” International journal of psychophysiology 38, 315–336 (2000).
  * Borgers and Kopell (2003) C. Borgers and N. Kopell, “Synchronization in networks of excitatory and inhibitory neurons with sparse, random connectivity,” Neural Comput 15, 509–538 (2003).
  * Chariker and Young (2015) L. Chariker and L.-S. Young, “Emergent spike patterns in neuronal populations,” Journal of computational neuroscience 38, 203–220 (2015).
  * Bezaire  _et al._ (2016) M. J. Bezaire, I. Raikov, K. Burk, D. Vyas, and I. Soltesz, “Interneuronal mechanisms of hippocampal theta oscillations in a full-scale model of the rodent ca1 circuit,” Elife 5, e18566 (2016).
  * Keeley  _et al._ (2019) S. Keeley, Á. Byrne, A. Fenton, and J. Rinzel, “Firing rate models for gamma oscillations,” J Neurophysiol 121, 2181–2190 (2019).
  * Ashwin, Coombes, and Nicks (2016) P. Ashwin, S. Coombes, and R. Nicks, “Mathematical Frameworks for Oscillatory Network Dynamics in Neuroscience,” J Math Neurosci 6, 2 (2016).
  * 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,” J Math Neurosci 10, 9 (2020).
  * Ter Wal and Tiesinga (2021) M. Ter Wal and P. H. Tiesinga, “Comprehensive characterization of oscillatory signatures in a model circuit with pv-and som-expressing interneurons,” Biological Cybernetics 115, 487–517 (2021).
  * Zonca and Holcman (2021) L. Zonca and D. Holcman, “Emergence and fragmentation of the alpha-band driven by neuronal network dynamics,” PLoS Computational Biology 17, e1009639 (2021).
  * Jensen and Colgin (2007) O. Jensen and L. L. Colgin, “Cross-frequency coupling between neuronal oscillations,” Trends in cognitive sciences 11, 267–269 (2007).
  * Segneri  _et al._ (2020) M. Segneri, H. Bi, S. Olmi, and A. Torcini, “Theta-nested gamma oscillations in next generation neural mass models,” Frontiers in computational neuroscience 14, 47 (2020).
  * Cai  _et al._ (2021) Y. Cai, T. Wu, L. Tao, and Z. C. Xiao, “Model Reduction Captures Stochastic Gamma Oscillations on Low-Dimensional Manifolds,” Front Comput Neurosci 15, 678688 (2021).
  * Ray and Maunsell (2015) S. Ray and J. H. Maunsell, “Do gamma oscillations play a role in cerebral cortex?” Trends Cogn Sci 19, 78–85 (2015).
  * Livingstone (1996) M. Livingstone, “Oscillatory firing and interneuronal correlations in squirrel monkey striate cortex,” J. Neurophysiol. 66, 2467–2485 (1996).
  * Geisler, Brunel, and Wang (2005) C. Geisler, N. Brunel, and X.-J. Wang, “Contributions of intrinsic membrane dynamics to fast network oscillations with irregular neuronal discharges,” Journal of neurophysiology 94, 4344–4361 (2005).
  * Zhang  _et al._ (2014) J. Zhang, D. Zhou, D. Cai, and A. V. Rangan, “A coarse-grained framework for spiking neuronal networks: between homogeneity and synchrony,” Journal of Computational Neuroscience 37, 81–104 (2014).
  * Li, Chariker, and Young (2019) Y. Li, L. Chariker, and L.-S. Young, “How well do reduced models capture the dynamics in models of interacting neurons?” Journal of mathematical biology 78, 83–115 (2019).
  * Li and Xu (2019) Y. Li and H. Xu, “Stochastic neural field model: multiple firing events and correlations,” Journal of mathematical biology 79, 1169–1204 (2019).
  * Koch (2004) C. Koch,  _Biophysics of computation: information processing in single neurons_ (Oxford university press, 2004).
  * Salinas and Sejnowski (2001) E. Salinas and T. Sejnowski, “Correlated neuronal activity and the flow of neural information,” Nat. Rev. Neurosci. 2, 539–550 (2001).
  * Quiroga and Panzeri (2009) R. Q. Quiroga and S. Panzeri, “Extracting information from neuronal populations: information theory and decoding approaches,” Nature Reviews Neuroscience 10, 173–185 (2009).
  * Quian Quiroga and Panzeri (2013) R. Quian Quiroga and S. Panzeri, eds.,  _Principles of neural coding_ (CRC Press, London, 2013).
  * Varga  _et al._ (2008) V. Varga, B. Hangya, K. Kránitz, A. Ludányi, R. Zemankovics, I. Katona, R. Shigemoto, T. F. Freund, and Z. Borhegyi, “The presence of pacemaker hcn channels identifies theta rhythmic gabaergic neurons in the medial septum,” The Journal of physiology 586, 3893–3915 (2008).
  * Hangya  _et al._ (2009) B. Hangya, Z. Borhegyi, N. Szilágyi, T. F. Freund, and V. Varga, “Gabaergic neurons of the medial septum lead the hippocampal network during theta activity,” Journal of Neuroscience 29, 8094–8102 (2009).
  * Pavlides, Hogan, and Bogacz (2015) A. Pavlides, S. J. Hogan, and R. Bogacz, “Computational models describing possible mechanisms for generation of excessive beta oscillations in parkinson’s disease,” PLoS computational biology 11, e1004609 (2015).
  * Yu and Wang (2019) Y. Yu and Q. Wang, “Oscillation dynamics in an extended model of thalamic-basal ganglia,” Nonlinear Dynamics 98, 1065–1080 (2019).
  * Zhang and Rangan (2015) J. W. Zhang and A. V. Rangan, “A reduction for spiking integrate-and-fire network dynamics ranging from homogeneity to synchrony,” Journal of Computational Neuroscience 38, 355–404 (2015).
  * Shao, Zhang, and Tao (2020) Y. Shao, J. Zhang, and L. Tao, “Dimensional reduction of emergent spatiotemporal cortical dynamics via a maximum entropy moment closure,” PLoS computational biology 16, e1007265 (2020).
  * Churchland  _et al._ (2010) M. M. Churchland, M. Y. Byron, J. P. Cunningham, L. P. Sugrue, M. R. Cohen, G. S. Corrado, W. T. Newsome, A. M. Clark, P. Hosseini, B. B. Scott, _et al._ , “Stimulus onset quenches neural variability: a widespread cortical phenomenon,” Nature neuroscience 13, 369–378 (2010).
  * Mante  _et al._ (2013) V. Mante, D. Sussillo, K. V. Shenoy, and W. T. Newsome, “Context-dependent computation by recurrent dynamics in prefrontal cortex,” Nature 503, 78–84 (2013).
  * Janak and Tye (2015) P. H. Janak and K. M. Tye, “From circuits to behaviour in the amygdala,” Nature 517, 284–292 (2015).
  * Gallego  _et al._ (2017) J. A. Gallego, M. G. Perich, L. E. Miller, and S. A. Solla, “Neural Manifolds for the Control of Movement,” Neuron 94, 978–984 (2017).
  * Saxena and Cunningham (2019) S. Saxena and J. P. Cunningham, “Towards the neural population doctrine,” Curr Opin Neurobiol 55, 103–111 (2019).
  * Cueva  _et al._ (2020) C. J. Cueva, A. Saez, E. Marcos, A. Genovesio, M. Jazayeri, R. Romo, C. D. Salzman, M. N. Shadlen, and S. Fusi, “Low-dimensional dynamics for working memory and time encoding,” Proc Natl Acad Sci U S A 117, 23021–23032 (2020).
  * Brunel and Hakim (1999) N. Brunel and V. Hakim, “Fast global oscillations in networks of integrate-and-fire neurons with low firing rates,” Neural Comput. 11, 1621–1671 (1999).
  * Cai  _et al._ (2006) D. Cai, L. Tao, A. Rangan, and D. McLaughlin, “Kinetic theory for neuronal network dynamics,” Comm. Math. Sci. 4, 97–127 (2006).
  * Acebrón  _et al._ (2005) J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler, “The Kuramoto model: A simple paradigm for synchronization phenomena,” Rev. Mod. Phys. 77, 137–185 (2005).
  * Montbrio, Pazo, and Roxin (2015) E. Montbrio, D. Pazo, and A. Roxin, “Macroscopic description for networks of spiking neurons,” Phys. Rev. X 5, 021028 (2015).
  * Gerstner  _et al._ (2014) W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski,  _Neuronal dynamics: From single neurons to networks and models of cognition_ (Cambridge University Press, 2014).