# Desynchronization of large-scale neural networks by stabilizing unknown
unstable incoherent equilibrium states

Tatjana Pyragienė and Kestutis Pyragas  Center for Physical Sciences and
Technology, Sauletekio al. 3, LT-10257 Vilnius, Lithuania

###### Abstract

In large-scale neural networks, coherent limit cycle oscillations usually
coexist with unstable incoherent equilibrium states, which are not observed
experimentally. We implement a first-order dynamic controller to stabilize
unknown equilibrium states and suppress coherent oscillations. The
stabilization of incoherent equilibria associated with unstable focus and
saddle is considered. The algorithm is demonstrated for networks composed of
quadratic integrate-and-fire (QIF) neurons and Hindmarsh-Rose neurons. The
microscopic equations of an infinitely large QIF neural network can be reduced
to an exact low-dimensional system of mean-field equations, which makes it
possible to study the control problem analytically.

###### keywords:

Neural network; Mean-field equations; Synchronization control; Quadratic
integrate-and-fire neurons; Hindmarsh-Rose neurons

##  1 Introduction

Synchronization studies in large populations of coupled oscillatory or
excitable elements are relevant in fields ranging from physics to neuroscience
[1, 2, 3, 4]. The role of synchronization in neural systems can be twofold. In
a healthy state, it is responsible for learning and cognition [5, 6], however,
excessive synchronization can cause a variety of neurological conditions such
as Parkinson’s disease [7], epilepsy [8, 9], tinnitus [10], and others. High-
frequency (HF) deep brain stimulation (DBS) is a standard procedure for the
treatment of neurological disorders [11, 12]. The mechanisms of DBS are not
yet well understood [13, 14]. Simple models show that the HF DBS effect can be
explained either as the result of stabilizing the resting state of individual
neurons [15] or as suppressing synchronized oscillations without forcing
individual neurons into silence [16]. HF DBS may cause side effects and its
therapeutic effect may decrease over time, so there is a significant clinical
need for less invasive and more effective stimulation methods [17]. In open
loop control systems such as HF DBS, adverse effects on neural tissue can be
reduced by optimizing the waveform of the stimulus signal [18, 19].

However, a number of theoretical works show that the desynchronization of
coherent oscillations is especially effective with the help of closed-loop
(feedback) control algorithms. Various control strategies based on linear [20,
21, 22, 23, 24] and nonlinear [25, 26, 27] time-delayed feedback, linear
feedback bandpass filters [28, 29, 30], proportional-integro-differential
feedback with a separate stimulation-registration setup [31], act-and-wait
time-delayed feedback [32, 33] and others [34, 35, 36] were considered.

Recent advances in the theory of nonlinear dynamical systems have provided the
neuroscience community with simple, low-dimensional models of neural networks
referred to as next-generation neural mass models [37]. Such models are useful
objects for developing, testing, and understanding various synchronization
control algorithms. Here we show that these models can naturally explain the
desynchronization mechanism of our feedback control algorithm in terms of
stabilizing unknown unstable incoherent states. The next-generation models are
derived directly from the microscopic dynamics of individual neurons and are
accurate in the thermodynamic limit of infinite network size. These models
represent a closed system of mean-field equations for biophysically relevant
parameters such as mean membrane potential and firing rate. Low-dimensional
dynamics in a large population of coupled oscillatory elements was first
discovered by Ott and Antonsen [38] in the Kuramoto model [2]. Later, this
discovery was successfully applied to derive a low-dimensional system of mean-
field equations for a certain class of networks consisting of all-to-all
pulse-coupled QIF neurons [39], which are canonical models of class I neurons
[40].

In recent years, next-generation models have been obtained for a large number
of different modifications of QIF neural networks [41, 42, 43, 44, 45, 46, 47,
48, 49, 50, 51]. These models make it possible to carry out their detailed
bifurcation analysis and reveal synchronization mechanisms. It has been shown
that synchronized limit cycle oscillations can arise from various
bifurcations, such as the Hopf bifurcation in Refs. [16, 42] or the homoclinic
bifurcation in Ref. [42]. However, stable limit cycles are always accompanied
by unstable fixed points, which correspond to unstable incoherent equilibrium
states of the network. These unstable states are not observed experimentally.
Here, we show that a priori unknown unstable incoherent states can be
stabilized using the control algorithm proposed in Refs. [52, 53]. Initially,
this algorithm was developed and tested to stabilize unknown unstable
equilibrium states of low-dimensional dynamical systems, and recently it has
been implemented to stabilize unstable pedestrian flows in the collective
behavior of large crowds of people [54]. Here, we implement this algorithm to
stabilize unstable incoherent states in large-scale neural networks consisting
of QIF and Hinmarsh-Rose [55] neurons. We demonstrate effective control of two
types of equilibrium states associated with an unstable focus and a saddle
point. As far as we know, the control of saddle equilibrium states in neural
networks has not been considered in the literature.

The paper is organized as follows. Section 2 describes the control algorithm.
In Sec. 3, we apply this algorithm to a population of synaptically coupled
excitatory QIF neurons. Here we stabilize incoherent states associated with an
unstable focus and a saddle fixed point. The latter is stabilized by an
unstable controller. Section 4 is devoted to the control of two interacting
populations of excitatory and inhibitory QIF neurons. In Sec. 5, we apply our
algorithm to a population of chaotically spiking Hindmarsh-Rose neurons, whose
microscopic model equations cannot be reduced to a low-dimensional system. The
conclusions are presented in Sec. 6.

##  2 Control algorithm

We consider a large network of coupled neurons generating collective coherent
oscillations. We assume that, along with the synchronous mode of coherent
oscillations, the network has an unstable equilibrium state characterized by
incoherent oscillations of individual neurons. Our goal is to stabilize the
incoherent state and transition the network from synchronous to incoherent
mode. To achieve this goal, we turn to the algorithm for stabilizing unknown
unstable equilibrium points of low-dimensional dynamical systems, developed in
Refs. [52, 53]. The algorithm uses a simple first-order dynamic controller
based on a low-pass filter (LPF). The block diagram of this algorithm, adapted
for neural networks, is shown in Fig. 1.

Figure 1: Block diagram of stabilization of unknown incoherent states in
neural networks. The mean membrane potential $v(t)$ represents the output of
the network. The network is stimulated by the input current $I_{c}(t)$. In a
feedback loop, LPF stands for low-pass filter.

We assume that the mean membrane potential $v(t)$ of the entire or some part
of the neural population can be measured at the output of the network. In
addition, we assume that all or part of the population of neurons can be
stimulated by the input current $I_{c}(t)$. In general, the measured and
stimulated subpopulations may differ. The input and output of the network are
connected by a feedback loop described by the following equations:

$\displaystyle\dot{w}$ $\displaystyle=$ $\displaystyle\omega_{c}(v-w),$ (1a)
$\displaystyle I_{c}$ $\displaystyle=$ $\displaystyle k(w-v),$ (1b)

where $w$ is a dynamic variable of the controller (LPF). The control algorithm
has two adjustable parameters: the cutoff frequency $\omega_{c}$ of the LPF
and the control gain $k$. Let us denote the average membrane potential of the
free network in a state of unstable equilibrium as $v=v^{*}$, which in the
thermodynamic limit should be a constant, $v^{*}=const$. We assume that this
value is a priori unknown. The control algorithm is designed in such a way
that the equilibrium value of $v^{*}$ remains unchanged in the stationary
state of the closed loop system. Indeed, at $\dot{w}=0$ the control variable
coincides with the mean membrane potential $w=w^{*}=v^{*}$, and the feedback
perturbation vanishes, $I_{c}=0$. However, feedback perturbation affects the
stability of the incoherent state. The examples below show that this state can
be stabilized by adjusting control parameters $\omega_{c}$ and $k$
accordingly.

This algorithm has a number of advantages. Firstly, it is weakly invasive.
Below, we will show that the feedback perturbation $I_{c}$ decreases according
to a power law with increasing network size and vanishes as the network size
tends to infinity. Secondly, this algorithm does not require knowledge of the
mean membrane potential $v^{*}$ of an unstable equilibrium state and, thirdly,
the algorithm provides tracking of the equilibrium state in the case of slowly
varying system parameters [53].

Note that the control algorithm with an ordinary LPF ($\omega_{c}>0$) has a
limitation. It works well for unstable equilibrium points like focuses but
doesn’t work for saddles. More precisely, Ref. [52] gives a theorem that a
stable controller cannot stabilize unstable equilibrium points with an odd
number of real positive eigenvalues. This limitation can be avoided by using
an unstable controller in the same way as it is done in the delayed feedback
control algorithm [56] when stabilizing a certain type of unstable periodic
orbits [57]. Here, to stabilize an unstable incoherent state of the saddle
type, we will use an unstable LPF with the parameter $\omega_{c}<0$. An
unstable LPF can be implemented using an RC circuit with a negative resistor.

In the following sections, we will demonstrate the performance of this
algorithm for three examples of neural networks. The first two examples deal
with large populations of synaptically coupled QIF neurons. In the limit of
infinite size, microscopic models of these networks can be reduced to exact
low-dimensional systems of mean-field equations. In the first example, one
population of excitatory neurons is considered, and in the second example, two
interacting populations of excitatory and inhibitory neurons are analyzed. The
third example is devoted to electrically coupled chaotic Hindmarsh-Rose
neurons.

##  3 Controlling a population of synaptically coupled excitatory QIF neurons

First, we apply the algorithm described above to a heterogeneous population of
QIF excitatory neurons interacting via finite-width synaptic pulses [42]. The
microscopic state of the population is defined by the set of $N$ neurons’
membrane potentials $\\{V_{j}\\}_{j=1,\ldots,N}$. They satisfy the following
set of equations [40]:

$\displaystyle\dot{V}_{j}$ $\displaystyle=$ $\displaystyle
V_{j}^{2}+\eta_{j}+Js(t)+I_{c}(t),\;\;$ (2)
$\displaystyle\;\;\text{if}\;V_{j}\geq V_{p}\;\text{then}\;V_{j}\leftarrow
V_{r}.$

Here, $\eta_{j}$ is a heterogeneous excitability parameter that specifies the
behavior of individual neurons and the term $Js(t)$ stands for the synaptic
coupling, where $J$ is the synaptic weight and $s(t)$ is the normalized mean
synaptic current emitted by spiking neurons. The term $I_{c}(t)$ describes an
external current, which we interpret as a control variable. In this model, the
membrane time constant of QIF neurons is assumed to be unity. This means that
time here is measured in units of the membrane time constant.

For $J=0$ and $I_{c}=0$, the neurons with the parameter $\eta_{j}<0$ are at
rest, and the neurons with the parameter $\eta_{j}>0$ generate spikes. When
the potential $V_{j}$ reaches the threshold value $V_{p}$, it is instantly
reset to the value $V_{r}$. We choose thresholds in the form
$V_{p}=-V_{r}=\infty$, which allows us to transform QIF neurons into theta
neurons and obtain an accurate system of reduced mean-field equations [39]. We
consider the case when the heterogeneous parameter $\eta$ is distributed
according to the Lorentzian density function

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

where $\Delta$ is the half-widths and $\bar{\eta}$ is the center of the
distribution. For the Lorentzian heterogeneity, the reduction of microscopic
equations is the most efficient. Note that other distributions of the
heterogeneous parameter have been considered in recent publications. [49, 50].

Here we use the model of global coupling in which neurons emit synaptic pulses
of finite width with the mean synaptic current defined as [42]

$s(t)=\frac{V_{th}}{N}\sum_{i=1}^{N}H(V_{i}(t)-V_{th}),$ (4)

where $H(\cdot)$ is the Heaviside step function and $V_{th}$ is a threshold
potential that determines the height and width of synaptic pulses.

In the limit $N\to\infty$, the above microscopic model reduces to an exact
system of two ordinary differential equations (ODEs) [42]

$\displaystyle\dot{r}$ $\displaystyle=$ $\displaystyle\Delta/\pi+2rv,$ (5a)
$\displaystyle\dot{v}$ $\displaystyle=$
$\displaystyle\bar{\eta}+v^{2}-\pi^{2}r^{2}+Js(t)+I_{c}(t)$ (5b)

for two biophysically relevant parameters, the mean spiking rate $r(t)$ and
the mean membrane potential $v(t)$. In the infinite size limit, the mean
synaptic current (4) is expressed in terms of the parameters $r(t)$ and $v(t)$
as [42]

$s(t)=\frac{V_{th}}{\pi}\left[\frac{\pi}{2}-\arctan\left(\frac{V_{th}-v(t)}{\pi
r(t)}\right)\right].$ (6)

This expression closes the system of mean-field Eqs. (5). The bifurcation
analysis of these equations without control $I_{c}(t)=0$ was carried out in
Ref. [42]. This analysis showed that synchronous limit cycle oscillations can
occur through two types of bifurcations: the Hopf bifurcation and the
homoclinic bifurcation. In the first case, the system (5) has a stable focus
before the bifurcation. On a microscopic level, this corresponds to a stable
equilibrium state of the network with incoherent dynamics of individual
neurons. After the bifurcation, the incoherent equilibrium state becomes an
unstable focus, and neurons exhibit coherent limit cycle oscillations. Our
goal here is to bring back the incoherent dynamics by stabilizing the unstable
equilibrium state. In the case of a homoclinic bifurcation, the limit cycle
touches the saddle point and becomes a homoclinic orbit. Near this
bifurcation, we will suppress coherent oscillations by using an unstable
controller to stabilize the incoherent state of the saddle equilibrium.

We begin the application of our control algorithm from the case of limit cycle
oscillations arising from the Hopf bifurcation. We use typical system
parameters corresponding to this mode [42]: $\Delta=1$, $V_{th}=50$,
$\bar{\eta}=2$, and $J=20$. For these parameters, the only attractor in the
two-dimensional phase space ($r,v$) of the free ($I_{c}=0$) system (5) is the
limit cycle. Inside this cycle there is an unstable focus with the coordinates
$(r^{*},v^{*})\approx(2.1081,-0.0755)$ and two complex-conjugate eigenvalues
$\lambda_{1,2}\approx 0.2621\pm 9.6190$. Let us now estimate how the local
properties of this fixed point change in the presence of a control defined by
the Eq. (1). Due to the additional variable $w$, the phase space of the closed
loop system is expanded to three dimensions: $(r,v,w)$. The coordinates of the
fixed point in the three-dimensional phase space are $(r^{*},v^{*},v^{*})$,
i.e. its projection onto the original two-dimensional phase space remains
unchanged. However, the stability properties of this fixed point now depend on
the controller parameters $\omega_{c}$ and $k$ and are determined by the
eigenvalue problem

$\det(A-\lambda I)=0$ (7)

of the linearized system of Eqs. (5) and (1). Here

$A=\begin{pmatrix}a_{11}&a_{12}&0\\\ a_{21}&a_{22}-k&k\\\
0&\omega_{c}&-\omega_{c}\end{pmatrix}$ (8)

is the Jacobian matrix of this system, $a_{ij}$ are the coefficients of the
Jacobian matrix of of the system (5) without control evaluated at the fixed
point $(r^{*},v^{*})$. Specifically, $a_{11}=2v^{*}$, $a_{12}=2r^{*}$,
$a_{21}=-2\pi^{2}r^{*}+JV_{th}(\pi r^{*})^{-2}(V_{th}-v^{*})c^{-1}$ and
$a_{22}=2v^{*}+JV_{th}\pi^{-2}(cr^{*})^{-1}$, where
$c=1+\left[(V_{th}-v^{*})/(\pi r^{*})\right]^{2}$. Finally, $I$ is the
identity matrix, and $\lambda$ is the eigenvalue.

For a given fixed point, the dependence of the solutions of the Eq. (7) on the
parameters $\omega_{c}$ and $k$ is shown in Fig. 2.

Figure 2: The performance of the control algorithm depending on the control
parameters $\omega_{c}$ and $k$. The results for an unstable focus in a
population of synaptically coupled QIF neurons are presented. The contour
lines and colors indicate the maximum real part of the eigenvalues
$\max[\mathrm{Re}(\lambda)]$ obtained from the Eq. (7). The thick red contour
line corresponds to $\max[\mathrm{Re}(\lambda)]=0$. It separates stable and
unstable regions. The originally unstable focus is stabilized in the region
$\max[\mathrm{Re}(\lambda)]<0$. Network parameters: $\Delta=1$, $V_{th}=50$,
$\bar{\eta}=2$, and $J=20$.

The colors encode the values of $\max[\mathrm{Re}(\lambda)]$. The thick red
contour line corresponds to $\max[\mathrm{Re}(\lambda)]=0$. It separates
regions of a stable and unstable fixed point. We see that the control
algorithm is robust to the choice of control parameters $\omega_{c}$ and $k$.
The algorithm provides stabilization of the unstable focus for any
$\omega_{c}>0$ and $k\gtrapprox 0.55$.

Figure 3 shows the performance of the control algorithm for fixed values of
$\omega_{c}=1$ and $k=2$.

Figure 3: Suppression of coherent oscillations by stabilization of an unstable
focus in a population of synaptically coupled QIF neurons. For $t<5$, there is
no control and the network generates collective coherent oscillations. For
$t>5$, the control is turned on and the system goes into a previously unstable
incoherent state. The dynamics of (a) mean membrane potential, (b) mean
spiking rate, and (d) control perturbation derived from the mean-field Eqs.
(5) are shown as thick gray curves. The thin red curves show the same results
derived from the microscopic model (10). (c) Raster plot of $200$ randomly
selected neurons. The spike moments for each neuron are shown by dots. The
neuron numbers are shown on the vertical axis. The parameters of the network
are the same as in Fig. 2. Controller parameters: $\omega_{c}=1$ and $k=2$.
The microscopic model was simulated using $N=10^{4}$ neurons.

The thick gray curves show the dynamics of the free and controlled neuronal
population obtained from the mean-field equations (5). During the time $t<5$
the control is switched off and the system is in the mode of limit cycle
oscillations. The mean membrane potential [Fig. 3(a)] and the mean spiking
rate [Fig. 3(b)] show periodic oscillations. At $t>5$ the control is activated
and the oscillations are damped. The system approaches a stabilized
equilibrium state. The control perturbation [Fig. 3(d)] experiences transient
damped oscillations and vanishes asymptotically.

As a next step, we tested the performance of our algorithm for networks of
finite size, described by the microscopic Eqs. (2). Unlike the low-dimensional
mean-field Eqs. (5), the microscopic model is defined by a huge number of
differential equations. The typical population sizes we model here are $N\sim
10^{4}$ neurons. There is of course no a priori guarantee whether the control
algorithm will work for such high-dimensional systems. Numerical simulation is
more convenient after changing variables

$V_{j}=\tan(\theta_{j}/2),$ (9)

which transforms QIF neurons into theta neurons. The advantage of theta
neurons is that they avoid the discontinuty problem. When the membrane
potential $V_{j}$ of the QIF neuron rises to $+\infty$ and falls to $-\infty$,
the theta neuron simply crosses the phase $\theta_{j}=\pi$. For theta neurons,
the Eqs (2) are transformed to

$\displaystyle\dot{\theta}_{j}$ $\displaystyle=$ $\displaystyle
1-\cos\left(\theta_{j}\right)\;\;$ (10) $\displaystyle+$
$\displaystyle\left[1+\cos\left(\theta_{j}\right)\right]\left[\eta_{j}+Js(t)+I_{c}(t)\right].$

We integrated these equations by the Euler method using a time step of
$dt=10^{-4}$. We have generated the values of the Lorentzian distributed (3)
heterogeneous parameter deterministically using
$\eta_{j}=\bar{\eta}+\Delta\tan(\pi/2)(2j-N-1)/(N+1)]$ for $j=1,\ldots,N$. For
more details on modelling the Eqs. (10), see Ref. [42]. From the Eqs. (10), we
estimated the Kuramoto order parameter [2]

$Z=\frac{1}{N}\sum\limits_{j=1}^{N}\exp(i\theta_{j})$ (11)

and used its relation with the spiking rate $r$ and the mean membrane
potential $v$ [39]:

$r=\frac{1}{\pi}\operatorname{Re}\left(\frac{1-Z^{*}}{1+Z^{*}}\right),\quad
v=\operatorname{Im}\left(\frac{1-Z^{*}}{1+Z^{*}}\right),$ (12)

where $Z^{*}$ denots complex conjugate of $Z$.

Results derived from the microscopic model (10) for $N=10^{4}$ neurons are
presented in Fig. 3 by thin red curves. They are in good agreement with the
results obtained from the reduced mean-field Eqs. (5). Thus, the control
algorithm works well for a large population of $N=10^{4}$ neurons, and the
mean-field theory correctly predicts the dynamics of the population in the
presence of control. To demonstrate network dynamics at the microscopic level,
Fig. 3(c) shows raster plots of 200 randomly selected neurons. Without
stimulation ($t<5$), most neurons spike coherently. Turning on the control at
$t>5$ destroys the coherent spiking and stabilizes the initially unstable
incoherent state.

Although the results of the mean-field equations and the microscopic model are
very close, there is a fundamental difference in the asymptotic dynamics of
these two models. As $t\to\infty$, the dynamic variables $(r,v)$ of the mean-
field equations approach exactly the unstable fixed point $(r^{*},v^{*})$ of
the uncontrolled system, and the control perturbation vanishes $I_{c}(t)\to
0$. In the microscopic model, the variables $(r,v)$ exhibit small fluctuations
around the fixed point $(r^{*},v^{*})$, and the control perturbation
$I_{c}(t)$ fluctuates around zero. Figure 4 shows the dependence of the
variance $\mathrm{Var}(I_{c})$ of the control perturbation in the post-
transient regime on the network size $N$. The variance decreases with
increasing $N$ and vanishes at $N\to\infty$. This dependence is well described
by the power law $\mathrm{Var}(I_{c})\sim N^{-\gamma}$ with $\gamma\approx
1.3$.

Figure 4: The variance $\mathrm{Var}(I_{c})$ of the control perturbation in
the post-transient regime as a function of the network size $N$. The asterisks
show the result of the numerical simulation, and the dashed line shows the
power-law approximation $\mathrm{Var}(I_{c})=CN^{-\gamma}$ with $C=1450$ and
$\gamma\approx 1.3$.

Let us now consider the control of coherent oscillations near a homoclinic
bifurcation. We will use the following set of the parameters: $\Delta=1$,
$V_{th}=50$, $\bar{\eta}=-7$, and $J=21$. For these parameters, the free
($I_{c}=0$) system (5) has a stable limit cycle and outside it a saddle point
with coordinates $(r^{*},v^{*})\approx(0.4073,-0.3908)$ and two real
eigenvalues $\lambda_{1,2}\approx(2.5306,-3.9255)$. Stabilization of the
incoherent state associated with the saddle point cannot be attained with an
ordinary LPF and requires the use of an unstable LPF with a negative parameter
$\omega_{c}$. The eigenvalues of the saddle point in presence of the control
are determined by the Eqs. (7) and (8). The dependence of the two largest real
parts of the eigenvalues on $k$ for a fixed $\omega_{c}=-1$ is shown in Fig.
5.

Figure 5: Linear stability of a saddle incoherent state of a population of QIF
neurons controlled by an unstable controller with a negative parameter
$\omega_{c}=-1$. Dependence of two largest real parts of eigenvalues of the
closed loop system on the control gain $k$. The inset shows the root loci of
the characteristic Eq. (7) in the complex plane $\lambda$ as $k$ changes from
$0$ to $\infty$. The crosses on the real axes indicate the location of the
eigenvalues at $k=0$, and the dot at the origin shows the location of one of
the eigenvalues at $k=\infty$. Network parameters: $\Delta=1$, $V_{th}=50$,
$\bar{\eta}=-7$, and $J=21$.

The saddle point stabilization mechanism is best understood from the root loci
diagram shown in the inset. Here we show the evolution of eigenvalues in the
complex plane $\lambda$ as $k$ changes from $0$ to $\infty$. Two crosses on
the real axes determine the location of the eigenvalues at $k=0$. One of them
$\lambda=2.5306$ corresponds to a free network, and the other
$\lambda=-\omega_{c}=1$ corresponds to a disabled unstable controller. With
the increase of $k$, they approach each other on the real axes, collide and
pass to the complex plane. At $k\approx 15.3$, they cross symmetrically into
the left half-plane (Hopf bifurcation). For very large $k\approx 91.8$, we
have a collision on the real axis again, and then one of the roots goes to
infinity, while the other approaches the origin. For $k>15.3$, the closed loop
system is stable.

Figure 6 shows the results of stabilization of a saddle incoherent state with
unstable controller parameters $\omega_{c}=-1$ and $k=20$. As in Fig. 3, the
dynamics derived from the mean-field equations are shown as thick gray curves,
and the corresponding dynamics derived from the microscopic model of $10^{4}$
neurons are shown as thin red curves. Again, there is complete agreement
between the mean-field theory and the microscopic theory. For $t<10$, there is
no control, and the system is in the limit cycle mode, which is close to a
homoclinic bifurcation. For $t>10$, the control is activated and the system
approaches a stabilized incoherent saddle point. In the mean-field theory, the
control perturbation vanishes asymptotically, while in the microscopic model
it experiences small fluctuations around zero. Note that the steady-state
spiking rate in saddle equilibrium is much lower than in focus equilibrium
[cp. post-transient dynamics in Figs. 3(b) and 6(b)].

Figure 6: Suppression of coherent oscillations in a population of QIF neurons
by stabilization of a saddle incoherent state with an unstable controller at
$\omega_{c}=-1$ and $k=20$. As in Fig. 3, the dynamics derived from the mean-
field equations are shown as thick gray curves, and the corresponding dynamics
derived from the microscopic model of $10^{4}$ neurons are shown as thin red
curves. All other designations are the same as in Fig. 3. The control turns on
at $t=10$. The network parameters correspond to Fig. 5.

##  4 Controlling two interacting populations of excitatory and inhibitory QIF
neurons

Let us now consider the control of a more complex network built from two
connected populations of excitatory and inhibitory QIF neurons. We follow the
model discussed in Ref. [16] whose network architecture mimics the network
architecture used in Parkinson’s disease models. Such models are usually based
on two interacting neural populations of the subthalamic nucleus (STN)
consisting of excitatory neurons and the external segment of the globus
pallidus (GPe) consisting of inhibitory neurons (cf.,e.g., Ref. [58]). It was
shown in [16] that synchronous oscillations can be very effectively suppressed
by HF stimulation of the inhibitory population, while HF stimulation of the
excitatory population is ineffective. Here we want to test whether our control
algorithm applied to the excitatory population can suppress synchronization.

The microscopic model of the network considered here is determined by the set
of $2N$ neurons’ membrane potentials $\\{V_{j}^{(E,I)}\\}_{j=1,\ldots,N}$.
They satisfy the system of $2N$ ODEs [16]:

$\displaystyle\tau_{m}\dot{V}_{j}^{(E,I)}$ $\displaystyle=$
$\displaystyle({V}_{j}^{(E,I)})^{2}+\eta_{j}^{(E,I)}+\mathcal{I}_{j}^{(E,I)},$
(13) $\displaystyle\;\;\text{if}\;\;{V}_{j}^{(E,I)}\geq
V_{p}\;\;\text{then}\;\;{V}_{j}^{(E,I)}\leftarrow V_{r},$

where, $V_{j}^{(E,I)}$ is the membrane potential of neuron $j$ in the
excitatory (E) or the inhibitory (I) population, and $\tau_{m}$ is the
membrane time constant. The threshold potential assumption is the same as in
the previous model: $V_{p}=-V_{r}=\infty$. The heterogeneous parameters
$\eta_{j}^{(E,I)}$ for populations E and I are taken from two independent
Lorentzian distributions:

$g_{E,I}(\eta)=\frac{1}{\pi}\frac{\Delta_{E,I}}{(\eta-\bar{\eta}_{E,I})^{2}+\Delta_{E,I}^{2}},$
(14)

where $\Delta_{E,I}$ and $\bar{\eta}_{E,I}$ are respectively the width and the
center of the distribution for the populations E and I. The last term
$\mathcal{I}_{j}^{(E,I)}$ in Eqs. (2) describes synaptic coupling and external
stimulation in the respective populations:

$\displaystyle\mathcal{I}_{j}^{(E)}$ $\displaystyle=$ $\displaystyle-
J_{IE}r_{I}(t)+I_{c}(t),$ (15a) $\displaystyle\mathcal{I}_{j}^{(I)}$
$\displaystyle=$ $\displaystyle J_{EI}r_{E}(t)-J_{II}r_{I}(t).$ (15b)

Unlike the previous model, here the interaction between neurons is provided by
instantaneous pulses. Each time the potential of a given neuron reaches
$\infty$, it resets to $-\infty$, and the neuron emits a Dirac delta spike,
which contributes to the output of the network. The mean synaptic rates of E
and I populations are as follows:

$\displaystyle r_{E,I}(t)=\lim_{\tau_{s}\to
0}\frac{\tau_{m}}{\tau_{s}N}\sum_{i=1}^{N}\sum_{k}\int_{t-\tau_{s}}^{t}\delta(t^{\prime}-(t_{i}^{k})_{E,I})dt^{\prime},$
(16)

where $\delta(t)$ is the Dirac delta function and $(t_{i}^{k})_{E,I}$ is the
time of the $k$th spike of the $i$th neuron in E and I population,
respectively. Parameters $J_{EI}$, $J_{IE}$ and $J_{II}$ denote synaptic
weights. The current $J_{EI}r_{E}(t)$ excites I neurons due to the synaptic
activity of E population and the current $-J_{IE}r_{I}(t)$ inhibits E neurons
due to the synaptic activity of the I population. The current
$-J_{II}r_{I}(t)$ recurrently inhibits neurons in population I. We are
considering a stimulation protocol in which only the excitatory population is
stimulated, so the control current $I_{c}(t)$ is only included in the Eq.
(15a).

In the limit $N\to\infty$, this microscopic model reduces to an exact closed
system of four ODEs for four biophysical quantities, mean firing rates
$r_{E,I}$ and mean membrane potentials $v_{E,I}$ of populations E and I [39,
16]:

$\displaystyle\tau_{m}\dot{r}_{E}$ $\displaystyle=$
$\displaystyle\Delta_{E}/\pi+2r_{E}v_{E},$ (17a)
$\displaystyle\tau_{m}\dot{v}_{E}$ $\displaystyle=$
$\displaystyle\bar{\eta}_{E}+v_{E}^{2}-\pi^{2}r_{E}^{2}-J_{IE}r_{I}+I_{c}(t),$
(17b) $\displaystyle\tau_{m}\dot{r}_{I}$ $\displaystyle=$
$\displaystyle\Delta_{I}/\pi+2r_{I}v_{I},$ (17c)
$\displaystyle\tau_{m}\dot{v}_{I}$ $\displaystyle=$
$\displaystyle\bar{\eta}_{I}+v_{I}^{2}-\pi^{2}r_{I}^{2}+J_{EI}r_{E}-J_{II}r_{I}.$
(17d)

Bifurcation analysis of an uncontrolled ($I_{c}=0$) system (17) showed a wide
variety of different dynamic modes [16]. Here we focus on the case when the
system has a single attractor, the limit cycle. Specifically, we consider the
following set of system parameters: $\Delta_{E}=0.05$, $\bar{\eta}_{E}=0.5$,
$\Delta_{I}=0.5$, $\bar{\eta}_{I}=-4$, $J_{EI}=20$, $J_{IE}=5$, $J_{II}=0.5$,
and $\tau_{m}=14$ ms. At these parameters, the system, along with a stable
limit cycle, has an unstable fixed point, which is a high dimensional focus
with coordinates
$(r_{E}^{*},v_{E}^{*},r_{I}^{*},v_{I}^{*})\approx(0.1319,-0.0603,0.0663,-1.1990)$
and two pairs of complex conjugate eigenvalues $\lambda_{1,2}\approx(0.0448\pm
1.0304i)/\tau_{m}$ and $\lambda_{3,4}\approx(-2.5634\pm 0.8190i)/\tau_{m}$.
Our goal is to stabilize this fixed point using the control algorithm defined
by Eqs. (1), with the constraint that the available network output is the mean
membrane potential of the excitatory population, $v=v_{E}$, and the control
current $I_{c}$ is applied only to the excitatory population. Linear stability
of the fixed point in the presence of control can be analyzed in a similar way
as in the previous model. Now the characteristic equation has five
eigenvalues. The dependence of the $\max[\mathrm{Re}(\lambda)]$ on the control
gain $k$ for three different values of the cutoff frequency $\omega_{c}$ is
shown in Fig. 7. Again we see that the stability condition
$\max[\mathrm{Re}(\lambda)]<0$ is satisfied in a wide range of the control
parameters $k$ and $\omega_{c}$.

Figure 7: Linear stability of incoherent state associated with a high-
dimensional focus in a system of two interacting populations of excitatory and
inhibitory QIF neurons in the presence of control. The entire network is
controlled using the output and input of the excitatory population only. The
maximum real part of the eigenvalues as a function of the control gain $k$ is
shown for different values of the cutoff frequency $\omega_{c}$ of LPF.
Network parameters: $\Delta_{E}=0.05$, $\bar{\eta}_{E}=0.5$, $\Delta_{I}=0.5$,
$\bar{\eta}_{I}=-4$, $J_{EI}=20$, $J_{IE}=5$, $J_{II}=0.5$, and $\tau_{m}=14$
ms.

Figure 8 shows the performance of the control algorithm for fixed values of
$\omega_{c}=0.5/\tau_{m}$ and $k=0.5$. The dynamics of the free ($t<300$ ms)
and controlled ($t>300$ ms) network, obtained from the mean-field Eqs. (17),
are shown as thick gray curves. The control switches the state of the system
from coherent limit cycle oscillations to the stabilized incoherent state and
the feedback perturbation asymptotically vanishes. These results are
consistent with numerical simulations of a microscopic model with $N=10^{4}$
neurons in each excitatory and inhibitory population (thin red curves). As in
the previous case, we changed the variables

$V_{j}^{(E,I)}=\tan(\theta_{j}^{(E,I)}/2)$ (18)

to rewrite the Eqs. (13) in terms of theta neurons:

$\displaystyle\tau_{m}\dot{\theta}_{j}^{(E,I)}$ $\displaystyle=$
$\displaystyle 1-\cos\left(\theta_{j}^{(E,I)}\right)$ (19) $\displaystyle+$
$\displaystyle\left[1+\cos\left(\theta_{j}^{(E,I)}\right)\right]\left[\eta_{j}^{(E,I)}+\mathcal{I}_{j}^{(E,I)}\right].$

We integrated these equations by the Euler method with a time step of
$dt=2\times 10^{-5}$. For the numerical implementation of Eq. (16), we set
$\tau_{s}=5\times 10^{-5}\tau_{m}$. To estimate the variables of the mean-
field theory, we calculated the Kuramoto order parameters

$Z_{E,I}=\frac{1}{N}\sum\limits_{j=1}^{N}\exp(i\theta_{j}^{(E,I)})$ (20)

for each population and evaluated the mean spiking rates and mean membrane
potentials for populations E and I as [39]:

$r_{E,I}=\frac{1}{\pi}\operatorname{Re}\left(\frac{1-Z_{E,I}^{*}}{1+Z_{E,I}^{*}}\right),\
v_{E,I}=\operatorname{Im}\left(\frac{1-Z_{E,I}^{*}}{1+Z_{E,I}^{*}}\right),$
(21)

where $Z_{E,I}^{*}$ denotes complex conjugate of $Z_{E,I}$. Panels (a), (c)
and (e) in Fig. 8 show a good agreement of time traces obtained from mean
filed equations and microscopic model. Panels (b) and (d) are raster plots of
randomly selected $500$ neurons in excitatory and inhibitory populations,
respectively.

Figure 8: Suppression of coherent oscillations in a system of two interacting
populations of excitatory and inhibitory QIF neurons by stabilization of
unstable incoherent state associated with a high-dimensional focus. Dynamics
of mean spiking rate of (a) excitatory and (c) inhibitory populations, and (e)
control perturbation applied to the excitatory population. The dynamics
derived from the mean-field equations are shown as thick gray curves, and the
corresponding dynamics derived from the microscopic model with $10^{4}$
neurons in each excitatory and inhibitory population are shown as thin red
curves. (b), (d) Raster plots of 500 randomly selected neurons in E and I
populations, respectively. The control turns on at $t=300$ ms. Network
parameters as in Fig. 7. Controller parameters: $\omega_{c}=0.5/\tau_{m}$ and
$k=0.5$.

##  5 Controlling a population of chaotically spiking Hindmarsh-Rose neurons

As a final example, consider the control of synchronous oscillations in a
heterogeneous population of electrically coupled Hindmarsh-Rose neurons [55]:

$\displaystyle\dot{v}_{j}$ $\displaystyle=$ $\displaystyle
y_{j}-v_{j}^{3}+3v_{j}^{2}-z_{j}+I_{j}+K(v-v_{j})+I_{c}(t),$ (22a)
$\displaystyle\dot{y}_{j}$ $\displaystyle=$ $\displaystyle
1-5v_{j}^{2}-y_{j},$ (22b) $\displaystyle\dot{z}_{j}$ $\displaystyle=$
$\displaystyle r[\nu(v_{j}-\kappa)-z_{j}],\quad j=1,\ldots,N.$ (22c)

Here, $v_{j}$, $y_{j}$ and $z_{j}$ are the membrane potential, the spiking
variable and the adaptation current of the $j$th neuron, respectively. The
variable

$v=\frac{1}{N}\sum_{i=1}^{N}v_{i}$ (23)

is the mean membrane potential. The heterogeneity of neurons is provided by
currents $I_{j}$, which we randomly select from a Gaussian distribution with a
mean value of $3$ and a variance of $0.1$. Parameters $r=0.06$, $\nu=4$ and
$\kappa=-1.56$ are chosen such that free ($K=0$ and $I_{c}=0$) neurons
generate chaotic bursts. The term $K(v-v_{j})$ in the Eq. (22a) determines the
electrical coupling between neurons, where $K$ is the coupling strength. To
get synchronized oscillations of the uncontrolled population, we take this
parameter large enough, $K=0.1$. The last term $I_{c}(t)$ in this equation is
the control current given by Eqs. (1).

Figure 9 shows how the control with fixed parameters $\omega_{c}=0.05$ and
$k=2$ changes the dynamics of a population of $N=10^{4}$ coupled neurons.
Without control ($t<700$), synchronous oscillations of large amplitude are
observed in the dynamics of the mean membrane potential $v(t)$, and coherent
bursts are visible on the raster plot. Activation of control at $t>700$
effectively suppresses synchronous oscillations of the mean membrane
potential, and neurons demonstrate incoherent bursts. As in previous examples,
only small amplitude oscillations around zero are observed in the asymptotic
dynamics of control perturbation $I_{c}(t)$. Figure 9(d) demonstrates that
control almost does not affect the amplitude dynamics of individual neurons.
As an example, we show the time trace of the membrane potential of the first
neuron $v_{1}(t)$ before and after activation of control.

Note that, unlike the previous examples, there is no known way to reduce this
model to a low-dimensional system. Thus, here we cannot theoretically estimate
the mean value of the membrane potential of an unstable incoherent state in
the thermodynamic limit and determine whether the equilibrium is associated
with an unstable focus or saddle and how its stability changes in the presence
of control. However, our algorithm does not require such detailed knowledge
and, with an appropriate choice of control parameters, works just as well as
in previous relatively simple models that allow a low-dimensional reduction in
the thermodynamic limit. Numerical simulations of this model show that our
algorithm works only when $\omega_{c}>0$ and fails when $\omega_{c}<0$. This
allows us to conclude that the unstable equilibrium in this model is an
unstable focus.

Figure 9: Suppression of coherent oscillations in a population of electrically
coupled Hidmarsh-Rose neurons. The control is activated at the time $t=700$.
(a) Dynamics of the mean membrane potential. (b) Raster plot of $100$ randomly
selected neurons. (c) and (d) Time traces of the control perturbation and the
membrane potential of the first neuron, respectively. Network parameters:
$r=0.06$, $\nu=4$, $\kappa=-1.56$, $K=0.1$ and $N=10^{4}$. Heterogeneous
currents $I_{j}$ in Eq. (22a) are randomly selected from a Gaussian
distribution with a mean value of $3$ and a variance of $0.1$. Controller
parameters: $\omega_{c}=0.05$ and $k=2$.

##  6  Conclusions

We considered the problem of suppressing collective synchronous oscillations
in large-scale neural networks. This problem is relevant in neurology, as
excessive synchronized oscillations in certain areas of the brain are often
associated with various neurological disorders [7, 8, 9, 10]. Synchronized
oscillations usually appear when an equilibrium incoherent state of the
network becomes unstable. Information about unstable network states is
difficult to extract from experimental data. We have shown that a priory
unknown unstable incoherent states of large-scale neural networks can be
effectively stabilized using a simple first order feedback controller based on
a low-pass filter. Initially, this controller was developed for stabilization
of unknown unstable equilibrium points of low-dimensional dynamical systems
[52, 53] and has not yet been tested for high-dimensional systems such as
neural networks, consisting of a huge number of interacting neurons.

We have demonstrated the effectiveness of our control algorithm on three
examples of neural networks. The first two examples refer to QIF neurons. In
the thermodynamic limit, microscopic models of networks built from QIF neurons
can be reduced to exact low-dimensional systems of mean-field equations. This
greatly simplifies the analysis of the effect of control on network dynamics.
In the first example, we demonstrated the suppression of synchronous
oscillations in a population of excitatory QIF neurons interacting through
synaptic pulses of finite width. Here we have stabilized two types of
incoherent states associated with an unstable focus and a saddle equilibrium
point. Until now, the control of the saddle equilibrium state in neural
networks has not been considered in the literature. Here we have achieved
stabilization of the saddle state with the help of an unstable controller. In
the second example, we considered the control of a network built from two
connected populations of excitatory and inhibitory QIF neurons, whose
architecture mimics that, used in Parkinson’s disease models [58]. Previously,
it was shown that high-frequency stimulation of the inhibitory population can
effectively suppress synchronization in such a network, but stimulation of the
excitatory population is ineffective [16]. Here, our algorithm provided
effective stabilization of the incoherent state of the network by using the
output and input of the excitatory population. For both first examples, the
results derived from the mean-field equations were confirmed by numerical
simulations of the respective microscopic models. We have shown that networks
of $10^{4}$ neurons are quantitatively well described by mean-field equations.
In the third example, we demonstrated the suppression of coherent oscillations
in a population of electrically coupled Hindmarsh-Rose neurons. Low-
dimensional reduction of the equations of the microscopic model is impossible
in this case. However, the direct application of the control algorithm to the
microscopic model showed that it works just as well as in the previous two
examples. Note that successful stabilization of unstable incoherent states
makes them experimentally observable, and this can serve as a quantitative
benchmark for assessing the quality of neural network models.

Finally, we summarize the main advantages of the proposed algorithm for
suppressing coherent oscillations in large-scale neural networks: (i) the
algorithm does not require any detailed knowledge of the network model and its
unstable incoherent equilibria; (ii) the algorithm is robust to changes in
control parameters; (iii) the algorithm can stabilize not only incoherent
states associated with an unstable focus but also with a saddle equilibrium
point; (iv) for large networks the algorithm is weakly invasive: the control
perturbation decreases according to a power law with increasing network size
and vanishes as the network size tends to infinity; (v) the algorithm is
adaptive, which means that it provides tracking of the equilibrium states in
the case of slowly varying system parameters (see [52, 53] for details).

In this paper, we limited ourselves to the consideration of the simplest
first-order controller to stabilize unknown incoherent states. More complex
networks may require higher-order generalized adaptive controllers [52]. In
addition, we emphasize that mean-field equations derived from microscopic
dynamics accurately describe synchronization processes in large networks, and
these models are well suited for testing and developing various algorithms for
suppressing unwanted coherent oscillations.

## Acknowledgments

This work is supported by grant No. S-MIP-21-2 of the Research Council of
Lithuania.

## References

  * [1] A. T. Winfree, The Geometry of Biological Time, Springer, 2001.
  * [2] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence, Springer-Verlag, 2003\.
  * [3] A. Pikovsky, M. Rosenblum, J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences, Cambridge University Press, 2001.
  * [4] S. Boccaletti, A. N. Pisarchik, C. I. del Genio, A. Amann, Synchronization: From Coupled Systems to Complex Networks, Cambridge University Press, 2018.
  * [5] W. Singer, Neuronal synchrony: A versatile code for the definition of relations?, Neuron 24 (1) (1999) 49–65.
  * [6] J. Fell, N. Axmacher, The role of phase synchronization in memory processes, Nat Rev Neurosci 12 (2) (2011) 105–118.
  * [7] C. Hammond, H. Bergman, P. Brown, Pathological synchronization in Parkinson’s disease: networks, models and treatments, Trends Neurosci. 30 (7) (2007) 357–364.
  * [8] P. Jiruska, M. de Curtis, J. G. R. Jefferys, C. A. Schevon, S. J. Schiff, K. Schindler, Synchronization and desynchronization in epilepsy: controversies and hypotheses, The Journal of Physiology 591 (4) (2013) 787–797.
  * [9] M. Gerster, R. Berner, J. Sawicki, A. Zakharova, A. Škoch, J. Hlinka, K. Lehnertz, E. Schöll, FitzHugh-Nagumo oscillators on complex networks mimic epileptic-seizure-related synchronization phenomena, Chaos: An Interdisciplinary Journal of Nonlinear Science 30 (12) (2020) 123130.
  * [10] P. A. Tass, O. V. Popovych, Unlearning tinnitus-related cerebral synchrony with acoustic coordinated reset stimulation: theoretical concept and modelling, Biol. Cybern. 106 (1) (2012) 27–36.
  * [11] A. Benabid, P. Pollak, D. Hoffmann, C. Gervason, M. Hommel, J. Perret, J. de Rougemont, D. Gao, Long-term suppression of tremor by chronic stimulation of the ventral intermediate thalamic nucleus, The Lancet 337 (8738) (1991) 403–406.
  * [12] M. L. Kringelbach, N. Jenkinson, S. L. Owen, T. Z. Aziz, Translational principles of deep brain stimulation, Nat. Rev. Neurosci. 8 (2007) 623–635.
  * [13] A. M. Lozano, H. Eltahawy, How does dbs work?, in: M. Hallett, L. Phillips, D. Schomer, J. Massey (Eds.), Advances in Clinical Neurophysiology, Vol. 57 of Supplements to Clinical Neurophysiology, Elsevier, 2004, pp. 733 –736.
  * [14] A. M. Lozano, N. Lipsman, H. Bergman, P. Brown, S. Chabardes, J. W. Chang, K. Matthews, C. C. McIntyre, T. E. Schlaepfer, M. Schulder, Y. Temel, J. Volkmann, J. K. Krauss, Deep brain stimulation: current challenges and future directions, Nat Rev Neurol 15 (3) (2019) 148–160.
  * [15] K. Pyragas, V. Novičenko, P. A. Tass, Mechanism of suppression of sustained neuronal spiking under high-frequency stimulation, Biological Cybernetics 107 (6) (2013) 669–684.
  * [16] K. Pyragas, A. P. Fedaravičius, T. Pyragienė, Suppression of synchronous spiking in two interacting populations of excitatory and inhibitory quadratic integrate-and-fire neurons, Phys. Rev. E 104 (2021) 014203\.
  * [17] J. K. Krauss, N. Lipsman, T. Aziz, A. Boutet, P. Brown, J. W. Chang, B. Davidson, W. M. Grill, M. I. Hariz, A. Horn, M. Schulder, A. Mammis, P. A. Tass, J. Volkmann, A. M. Lozano, Technology of deep brain stimulation: current status and future directions, Nat Rev Neurol 17 (2) (2021) 75–87.
  * [18] K. Pyragas, A. P. Fedaravičius, T. Pyragienė, P. A. Tass, Optimal waveform for entrainment of a spiking neuron with minimum stimulating charge, Phys. Rev. E 98 (2018) 042216.
  * [19] K. Pyragas, A. P. Fedaravičius, T. Pyragienė, P. A. Tass, Entrainment of a network of interacting neurons with minimum stimulating charge, Phys. Rev. E 102 (2020) 012221.
  * [20] M. G. Rosenblum, A. S. Pikovsky, Controlling synchronization in an ensemble of globally coupled oscillators, Phys. Rev. Lett. 92 (2004) 114102.
  * [21] M. G. Rosenblum, A. S. Pikovsky, Delayed feedback control of collective synchrony: An approach to suppression of pathological brain rhythms, Phys. Rev. E 70 (2004) 041904.
  * [22] C. Hauptmann, O. Popovych, P. A. Tass, Delayed feedback control of synchronization in locally coupled neuronal networks, Neurocomputing 65-66 (2005) 759–767.
  * [23] C. Hauptmann, O. Popovych, P. A. Tass, Effectively desynchronizing deep brain stimulation based on a coordinated delayed feedback stimulation via several sites: a computational study, Biol. Cybern. 93 (2005) 463–470.
  * [24] C. A. S. Batista, S. R. Lopes, R. L. Viana, A. M. Batista, Delayed feedback control of bursting synchronization in a scale-free neuronal network, Neural Networks 23 (2010) 114–124.
  * [25] O. V. Popovych, C. Hauptmann, P. A. Tass, Effective desynchronization by nonlinear delayed feedback, Phys. Rev. Lett. 94 (2005) 164102.
  * [26] O. V. Popovych, C. Hauptmann, P. A. Tass, Control of neuronal synchrony by nonlinear delayed feedback, Biol. Cybern. 95 (2006) 69–85.
  * [27] O. V. Popovych, P. A. Tass, Synchronization control of interacting oscillatory ensembles by mixed nonlinear delayed feedback, Phys. Rev. E 82 (2010) 026204.
  * [28] N. Tukhlina, M. Rosenblum, A. Pikovsky, J. Kurths, Feedback suppression of neural synchrony by vanishing stimulation, Phys. Rev. E 75 (2007) 011918.
  * [29] L. Ming, W. Yongjun, J. Peng, Washout filter aided mean field feedback desynchronization in an ensemble of globally coupled neural oscillators, Biol. Cybern. 101 (2009) 214–246.
  * [30] G. Montaseri, M. J. Yazdanpanah, A. Pikovsky, M. Rosenblum, Synchrony suppression in ensembles of coupled oscillators via adaptive vanishing feedback, Chaos 23 (2013) 033122.
  * [31] K. Pyragas, O. V. Popovych, P. A. Tass, Controlling synchrony in oscillatory networks with a separate stimulation-registration setup, Europhys. Lett. 80 (4) (2007) 40002.
  * [32] I. Ratas, K. Pyragas, Controlling synchrony in oscillatory networks via an act-and-wait algorithm, Phys. Rev. E 90 (2014) 032914.
  * [33] I. Ratas, K. Pyragas, Eliminating synchronization in bistable networks, Nonlinear Dyn. 83 (3) (2016) 1137–1151.
  * [34] B. Rosin, M. Slovik, R. Mitelman, M. Rivlin-Etzion, S. N. Haber, Z. Israe, E. Vaadia, H. Bergman, Closed-loop deep brain stimulation is superior in ameliorating parkinsonism, Neuron 72 (2011) 370–384.
  * [35] A. Berényi, M. Belluscio, D. Mao, G. Buzsáki, Closed-loop control of epilepsy by transcranial electrical stimulation, Science 337 (2012) 735–737.
  * [36] V. H. P. Louzada, N. A. M. Araújo, J. S. Andrade Jr., H. J. Herrmann, How to suppress undesired synchronization, Scientific Reports 2 (2012) 658.
  * [37] S. Coombes, Next generation neural population models, Frontiers in Applied Mathematics and Statistics 9 (2023).
  * [38] E. Ott, T. M. Antonsen, Low dimensional behavior of large systems of globally coupled oscillators, Chaos: An Interdisciplinary Journal of Nonlinear Science 18 (3) (2008) 037113.
  * [39] E. Montbrió, D. Pazó, A. Roxin, Macroscopic description for networks of spiking neurons, Phys. Rev. X 5 (2015) 021028.
  * [40] G. Bard Ermentrout, D. H. Terman, Mathematical Foundations of Neuroscience, Springer, New York, 2010.
  * [41] D. Pazó, E. Montbrió, From quasiperiodic partial synchronization to collective chaos in populations of inhibitory neurons with delay, Phys. Rev. Lett. 116 (2016) 238101.
  * [42] I. Ratas, K. Pyragas, Macroscopic self-oscillations and aging transition in a network of synaptically coupled quadratic integrate-and-fire neurons, Phys. Rev. E 94 (2016) 032215.
  * [43] I. Ratas, K. Pyragas, Macroscopic oscillations of a quadratic integrate-and-fire neuron network with global distributed-delay coupling, Phys. Rev. E 98 (2018) 052224.
  * [44] I. Ratas, K. Pyragas, Noise-induced macroscopic oscillations in a network of synaptically coupled quadratic integrate-and-fire neurons, Phys. Rev. E 100 (2019) 052211.
  * [45] E. Montbrió, D. Pazó, Exact mean-field theory explains the dual role of electrical synapses in collective synchronization, Phys. Rev. Lett. 125 (2020) 248101.
  * [46] M. Segneri, H. Bi, S. Olmi, A. Torcini, Theta-nested gamma oscillations in next generation neural mass models, Front. Comput. Neurosci. 14 (2020) 47.
  * [47] D. S. Goldobin, M. di Volo, A. Torcini, Reduction methodology for fluctuation driven population dynamics, Phys. Rev. Lett. 127 (2021) 038301.
  * [48] D. S. Goldobin, Mean-field models of populations of quadratic integrate-and-fire neurons with noise on the basis of the circular cumulant approach, Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (8) (2021) 083112.
  * [49] V. Pyragas, K. Pyragas, Dynamics of a network of quadratic integrate-and-fire neurons with bimodal heterogeneity, Phys. Lett. A 416 (2021) 127677.
  * [50] V. Pyragas, K. Pyragas, Mean-field equations for neural populations with $q$-Gaussian heterogeneities, Phys. Rev. E 105 (2022) 044402.
  * [51] V. Pyragas, K. Pyragas, Effect of Cauchy noise on a network of quadratic integrate-and-fire neurons with non-Cauchy heterogeneities, Physics Letters A 480 (2023) 128972.
  * [52] K. Pyragas, V. Pyragas, I. Z. Kiss, J. L. Hudson, Stabilizing and tracking unknown steady states of dynamical systems, Phys. Rev. Lett. 89 (2002) 244103\.
  * [53] K. Pyragas, V. Pyragas, I. Z. Kiss, J. L. Hudson, Adaptive control of unknown unstable steady states of dynamical systems, Phys. Rev. E 70 (2004) 026215.
  * [54] I. Panagiotopoulos, J. Starke, W. Just, Control of collective human behavior: Social dynamics beyond modeling, Phys. Rev. Res. 4 (2022) 043190.
  * [55] J. L. Hindmarsh, R. M. Rose, A model of neuronal bursting using three coupled first order differential equations, Proc. R. Soc. Lond. B. 221 (1984) 87–102.
  * [56] K. Pyragas, Continuous control of chaos by self-controlling feedback, Phys. Lett. A 170 (6) (1992) 421–428.
  * [57] K. Pyragas, Control of chaos via an unstable delayed feedback controller, Phys. Rev. Lett. 86 (2001) 2265–2268.
  * [58] D. Terman, J. E. Rubin, A. C. Yew, C. J. Wilson, Activity patterns in a model for the subthalamopallidal network of the basal ganglia, The Journal of Neuroscience 22 (7) (2002) 2963.