Desynchronization of large-scale neural networks by stabilizing unknown
  unstable incoherent equilibrium states
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.
Tatjana Pyragiene, Kestutis Pyragas
Desynchronization of large-scale neural networks by stabilizing unknown unstable incoherent equilibrium states

###### 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.

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:

\[\dot{w} = \omega_{c}(v-w), \tag{1a}\] \[I_{c} = k(w-v), \tag{1b}\]

where \(w\) is a dynamic variable of the controller (LPF). The control algorithm has two adjustable parameters: the cut-off 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

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.
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]:

\[\dot{V}_{j} = V_{j}^{2}+\eta_{j}+Js(t)+I_{c}(t), \tag{2}\] \[\mbox{if }V_{j}\geq V_{p}\mbox{ then }V_{j}\gets 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}}, \tag{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}), \tag{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\rightarrow\infty\), the above microscopic model reduces to an exact system of two ordinary differential equations (ODEs) [42]

\[\dot{r} = \Delta/\pi+2rv, \tag{5a}\] \[\dot{v} = \bar{\eta}+v^{2}-\pi^{2}r^{2}+Js(t)+I_{c}(t) \tag{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]. \tag{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 \tag{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} \tag{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+[(V_{th}-v^{*})/(\pi r^{*})]^{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. 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\). 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), \tag{9}\]

which transforms QIF neurons into theta neurons. The advantage of theta neurons is that they avoid the discontinuy 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

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.

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\).
\[\dot{\theta}_{j} = 1-\cos\left(\theta_{j}\right) \tag{10}\] \[+ \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_{j=1}^{N}\exp(i\theta_{j}) \tag{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), \tag{12}\]

where \(Z^{*}\) denotes 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 \(\operatorname{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 \(\operatorname{Var}(I_{c})\sim N^{-\gamma}\) with \(\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. 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 4: The variance \(\operatorname{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 \(\operatorname{Var}(I_{c})=CN^{-\gamma}\) with \(C=1450\) and \(\gamma\approx 1.3\).

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\).
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)].

## 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]:

\[\tau_{m}\dot{V}_{j}^{(E,I)} = (V_{j}^{(E,I)})^{2}+\eta_{j}^{(E,I)}+\mathcal{I}_{j}^{(E,I)}, \tag{13}\] \[\text{if}\;\;V_{j}^{(E,I)}\geq V_{p}\;\;\text{then}\;\;V_{j}^{(E, I)}\gets 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}}, \tag{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:

\[\mathcal{I}_{j}^{(E)} = -J_{IE}r_{I}(t)+I_{c}(t), \tag{15a}\] \[\mathcal{I}_{j}^{(I)} = J_{EI}r_{E}(t)-J_{II}r_{I}(t). \tag{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:

\[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}, \tag{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).

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.
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]:

\[\tau_{m}\dot{r}_{E} = \Delta_{E}/\pi+2r_{E}v_{E}, \tag{17a}\] \[\tau_{m}\dot{v}_{E} = \bar{\eta}_{E}+v_{E}^{2}-\pi^{2}r_{E}^{2}-J_{IE}r_{I}+I_{c}(t),\] (17b) \[\tau_{m}\dot{r}_{I} = \Delta_{I}/\pi+2r_{I}v_{I},\] (17c) \[\tau_{m}\dot{v}_{I} = \bar{\eta}_{I}+v_{I}^{2}-\pi^{2}r_{I}^{2}+J_{EI}r_{E}-J_{II}r_{I}. \tag{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

\[\left(r_{E}^{*},v_{E}^{*},r_{I}^{*},v_{I}^{*}\right)\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 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) \tag{18}\]

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

\[\tau_{m}\dot{\theta}_{j}^{(E,I)} = 1-\cos\left(\theta_{j}^{(E,I)}\right) \tag{19}\] \[+ \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_{j=1}^{N}\exp(i\theta_{j}^{(E,I)}) \tag{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), \tag{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.

## 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]:

\[\dot{v}_{j} = y_{j}-v_{j}^{3}+3v_{j}^{2}-z_{j}+I_{j}+K(v-v_{j})+I_{c}(t), \tag{22a}\] \[\dot{y}_{j} = 1-5v_{j}^{2}-y_{j},\] (22b) \[\dot{z}_{j} = r[\nu(v_{j}-\kappa)-z_{j}],\quad j=1,\ldots,N. \tag{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

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 cut-off 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.
\[v=\frac{1}{N}\sum_{i=1}^{N}v_{i} \tag{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.

## 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

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_{\rm m}\) and \(k=0.5\).
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.

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_{\mathrm{c}}=0.05\) and \(k=2\).