# Control of seizure-like dynamics in neuronal populations with excitability
adaptation related to ketogenic diet

Sebastian Eydam<EMAIL_ADDRESS>Neural Circuits and Computations Unit,
RIKEN Center for Brain Science, 2-1 Hirosawa, 351-019 Wako, Japan  Igor
Franović<EMAIL_ADDRESS>Scientific Computing Laboratory, Center for the
Study of Complex Systems, Institute of Physics Belgrade, University of
Belgrade, Pregrevica 118, 11080 Belgrade, Serbia  Louis Kang
<EMAIL_ADDRESS>Neural Circuits and Computations Unit, RIKEN Center for
Brain Science, 2-1 Hirosawa, 351-019 Wako, Japan

###### Abstract

We consider a heterogeneous, globally coupled population of excitatory
quadratic integrate-and-fire neurons with excitability adaptation due to a
metabolic feedback associated with ketogenic diet, a form of therapy for
epilepsy. Bifurcation analysis of a three-dimensional mean-field system
derived in the framework of next-generation neural mass models allows us to
explain the scenarios and suggest control strategies for the transitions
between the neurophysiologically desired asynchronous states and the
synchronous, seizure-like states featuring collective oscillations. We reveal
two qualitatively different scenarios for the onset of synchrony. For weaker
couplings, a bistability region between the lower- and the higher-activity
asynchronous states unfolds from the cusp point, and the collective
oscillations emerge via a supercritical Hopf bifurcation. For stronger
couplings, one finds seven co-dimension two bifurcation points, including
pairs of Bogdanov-Takens and generalized Hopf points, such that both lower-
and higher-activity asynchronous states undergo transitions to collective
oscillations, with hysteresis and jump-like behavior observed in vicinity of
subcritical Hopf bifurcations. We demonstrate three control mechanisms for
switching between asynchronous and synchronous states, involving parametric
perturbation of the adenosine triphosphate (ATP) production rate, external
stimulation currents, or pulse-like ATP shocks, and indicate a potential
therapeutic advantage of hysteretic scenarios.

> Adaptation has a profound impact on the dynamics of coupled systems across a
> wide variety of fields. In neuronal systems, it has two facets: one
> associated with synaptic plasticity, and the other related to modification
> of local kinetics, such as the spiking frequency or the excitability
> feature. Here we introduce a model of excitability adaptation reflecting the
> changes in neuronal metabolism linked to the ketogenic diet, a dietary
> treatment for epilepsy based on replacing carbohydrates with fat. The diet
> induces a shift in the main mechanism of ATP production, which in turn leads
> to activation of ATP-gated potassium channels in neuronal membranes. We
> build a model of a heterogeneous, globally coupled population of inherently
> excitable or tonic spiking excitatory quadratic integrate-and-fire neurons
> influenced by an ATP-dependent hyperpolarizing current, endowed by an
> additional equation describing the changes in the total ATP concentration.
> The stability of solutions and their bifurcations are investigated by
> analyzing the corresponding mean-field model via the numerical path-
> following technique. We disclose qualitatively different scenarios of
> multistability and bifurcations depending on the coupling strength.
> Interpreting the low-activity asynchronous stationary state as the
> physiologically desired (homeostatic) one, and the emergence of synchronous
> solutions featuring collective oscillations as the signature of seizure-like
> events, we demonstrate three control strategies to suppress the seizure
> dynamics by inducing a switch to the homeostatic state. Interestingly, such
> switches may be triggered by the changes in the ATP production rate, as a
> form of parametric perturbation, both for smooth or hysteretic transitions
> between synchronous and asynchronous states.

##  I Introduction

The recent advent of the next-generation neural mass (NGNM) models [1, 2, 3,
4, 5] has allowed for a deeper understanding of multistability and critical
transitions in the dynamics of neuronal populations, and has provided us with
the means to systematically explore the control strategies for switching
between the states interpreted in terms of neurophysiology as regular or
pathological ones. Earlier phenomenological neural mass models [6, 7, 8]
described the dynamics of neuronal populations in terms of firing rate
equations, typically assuming a sigmoid-like transfer function between the
mean membrane potential and the mean firing rate. Such models relied on the
assumption of uncorrelated local dynamics and could neither capture the
scenarios where the state of synchrony undergoes qualitative change nor
establish an exact relationship between the features of local and collective
dynamics. These issues have been resolved in the NGNM models, based on
applying the Ott-Antonsen [9, 10, 11] or the equivalent Lorentzian Ansatz [1]
to populations of theta neurons and quadratic integrate-and-fire neurons,
respectively. NGNM models characterize population dynamics in terms of mean-
field equations for the mean firing rate and the mean membrane potential.
While the notion of spike synchrony, classically invoked in neuroscience, is
reflected in the increase of the mean firing rate, the NGNM framework further
establishes a conformal mapping [1] between the mean firing rate and the
population synchrony described by the Kuramoto order parameter, a synchrony
measure inherited from the theory of coupled phase oscillators. Given that the
mapping is a real valued function of the complex argument, different values of
the Kuramoto order parameter may correspond to the same mean firing rate. NGNM
models have the advantage of being exact in the thermodynamic limit and are
amenable to bifurcation analysis. So far, NGNM models have been used to
address a broad range of theoretical problems, including understanding the
impact and interplay of chemical and electrical synapses [12, 13, 14], the
effects of deterministic or stochastic drive [15, 16, 17, 18, 19], presence of
quenched randomness in network topology [20], finite-size fluctuations [21,
22], interaction of neuronal populations [23, 24, 25], as well as formation of
bumps and waves in neural fields [26, 27, 28]. Moreover, a wide spectrum of
applications has been considered, from explaining the onset of gamma and
theta-nested gamma oscillations [29, 30, 31, 32, 33], cross-frequency
couplings [34], abnormal beta rebound in schizophrenia [35], the impact of
deep brain stimulation [36], to studies on working memory [24, 37],
propagation of epileptic seizures [38], and modeling of the entire brain
functional networks [4].

A corpus of problems that have so far been less explored within the formalism
of NGNM models concerns the impact of neuronal adaptation [39, 37, 40, 41]. In
general, neuronal systems feature two broad types of adaptive behavior [42].
The first is synaptic plasticity [43, 44], which involves modification of
coupling strengths due to feedback with neuronal activity. The other type
comprises adaptation of local kinetics, including spike frequency adaptation
[45, 46, 41] or changes in neuronal excitability due to metabolic constraints
[47, 48, 49]. Our study addresses the latter problem. In particular, we use
the NGNM framework to investigate how the population dynamics is affected by
the excitability adaptation due to a severe change of neuronal energy
metabolism triggered by the ketogenic diet (KD) [50, 51, 52, 53, 54, 55].
Constructing the bifurcation diagrams to reveal critical transitions and
regions of multistability, we demonstrate different control strategies that
may be used to induce switches between asynchronous (”homeostatic”, i.e.
neurophysiologically desired equilibrium states with low firing rates) and
synchronous, seizure-like regimes.

Within recent years, there has been a growing interest in the impact of energy
metabolism constraints, both on maintaining neuronal homeostasis and the
emergence of neurological disorders [56, 57, 58, 49]. The neuronal energy
demand is typically high even in the resting state because it requires
maintaining ion fluxes against concentration and electric gradients [48]. Such
ion fluxes are mainly perpetuated by the $Na^{+}/K^{+}$ pump which consumes
energy in the form of ATP [59, 60], produced in the glycolysis process
unfolding close to the neuronal membrane. It has long been understood that
some neurological disorders, such as epilepsy [47, 61, 62, 63] or Parkinson’s
disease [64, 65], involve a strong metabolic component. Epilepsy is the third
most prevalent neurological disorder [66, 67], and is characterized by
recurring seizures. Despite the advanced anti-convulsive medications, about a
third of patients still remain drug-resistant and require alternative
therapies [67]. Among the latter, KD [50, 52, 53, 54, 55], whose origins date
back to the 1920s, is nowadays often recommended to children with
pharmacology-resistant seizures [68]. KD is based on a 4:1 intake ratio of fat
vs. carbohydrates plus proteins, which gives rise to a drastic metabolic
change called ketosis. In the state of ketosis, the glycolytic ATP production
reduces, becoming effectively replaced by an alternative oxidation mechanism
in mitochondria, fueled by the ketone bodies generated in the liver using
diet-related fatty acids. The KD’s effectiveness apparently derives from
leveraging anti-convulsive mechanisms that are not (sufficiently) targeted by
medications [50].

Figure 1: Schematic of the excitability adaptation associated with the changes
in the energy metabolism due to ketogenic diet. The diet causes the main
mechanism of ATP production to switch from sub-membrane glycolysis to an
oxidation process in mitochondria. In parallel, activation of the ATP-gated
potassium channels in the neuronal membrane gives rise to an ATP-dependent
hyperpolarizing current. ATP balance is further influenced by consumption due
to maintaining the resting potential and enabling spiking.

KD is believed to induce two different anti-seizure types of effects [50]: one
concerns a direct impact of ketone bodies to inhibition of excitatory synaptic
channels, and the other is related to the metabolic switch from glycolytic to
mitochondria-based ATP production. Our study focuses on the latter group of
effects, whose schematic is provided in Fig. 1. In brief, abolishing the
glycolytic ATP production results in the depletion of ATP in a near-membrane
region, which causes activation of ATP-sensitive potassium channels in the
membrane. Their activation has long been known to induce hyperpolarizing
currents that reduce neuronal excitability [69, 70, 71]. Our approach involves
a model of local neuronal dynamics supplied with the hyperpolarizing current
associated with ATP-gated potassium channels, as well as the ATP dynamics
describing its production and utilization.

The paper is organized as follows. In Sec. II, we introduce the model of a
globally coupled heterogeneous population of intrinsically excitable or tonic
spiking excitatory quadratic integrate-and-fire neurons supplied by the local
terms for the ATP-dependent hyperpolarizing currents and an equation for the
global dynamics of ATP concentration. Using the approach from [1], in Sec. III
we derive the corresponding mean-field model, which compared to the classical
results, contains additional dependencies of the mean firing rate and mean
membrane potential on the ATP concentration level. Section IV provides the
bifurcation analysis of the reduced system, revealing two qualitatively
different bifurcation scenarios in terms of transitions to synchrony and the
system’s multistability depending on the coupling strength. In Sec. V we
demonstrate mechanisms for controlled switching between the homeostatic and
seizure-like states, including external stimulation, parametric and dynamical
perturbation. Section VI provides a summary and discussion of our results.

##  II Model of microscopic dynamics

We consider a large heterogeneous population of synaptically coupled
excitatory quadratic integrate-and-fire neurons [72, 73] whose dynamics is
given by

$\displaystyle\dot{V}_{j}$
$\displaystyle=V_{j}^{2}+\eta_{j}+KS(t)+I_{j}^{\mathrm{ATP}}(t)+I_{\mathrm{ext}}(t)$
$\displaystyle I_{j}^{\mathrm{ATP}}(t)$ $\displaystyle=-\alpha
V_{j}(t)\frac{\tilde{C}}{C(t)}$ (1)

where ${V_{j}}$ for $j=\\{1,…,N\\},\,N\gg 1$ are the neuronal membrane
potentials. The population diversity is manifested through heterogeneity of
local bifurcation parameters $\eta_{j}$, which follow the Lorentzian density
distribution

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

and may in terms of neurophysiology be interpreted as heterogeneous components
of external bias current. Quadratic integrate-and-fire neurons are
paradigmatic for type I excitability, such that an isolated unit undergoes a
transition from excitable regime ($\eta_{j}\lesssim 0$) to tonic spiking
($\eta_{j}>0$) with a frequency $f_{j}=\sqrt{\eta_{j}}/\pi$ via a SNIPER
(saddle-node of infinite period) bifurcation at $\eta_{j}=0$. A neuron fires a
spike whenever its membrane potential reaches the peak value $V_{p}$, after
which it is reset to $-V_{r}$. To allow for subsequent analytical
tractability, we set $V_{p}=-V_{r}=\infty$.

The terms $I_{\mathrm{ext}}(t)$ and $KS(t)$ describe a common external
stimulation and the total synaptic current, respectively. The chemical
synapses are characterized by a uniform coupling strength $K$ and are assumed
to be instantaneous, such that the mean synaptic activity is given by the
normalized output signal of the population

$\displaystyle S(t)=\frac{1}{N}\sum_{j=1}^{N}\sum_{k\backslash
t_{j}^{k}<t}\delta(t-t_{j}^{k}),$ (3)

where $t_{j}^{k}$ is the time of the k-th spike of neuron $j$. The
excitability adaptation triggered by the changes in the energy metabolism due
to the ketogenic diet is mediated by the ATP-dependent hyperpolarizing
currents $I_{j}^{\mathrm{ATP}}$. The latter reflects the action of the ATP-
gated potassium channels in neuronal membranes activated by the diet [69, 74,
75], and are assumed to have a uniform conductance $\alpha$. Note that the
gating of the potassium channels is in fact a complex process governed by the
exchange of sodium and ATP [71, 70]. Nevertheless, for simplicity, here we
adopt the model for $I_{j}^{\mathrm{ATP}}$ from [69] which omits the detailed
equation for the sodium dynamics. To close the system (1), one requires an
additional equation for the dynamics of ATP concentration $C(t)$

$\displaystyle\dot{C}$
$\displaystyle=\frac{\tilde{C}-C}{\tau}-\frac{\epsilon\,S(t)\,C}{\tilde{C}},$
(4)

where $\tilde{C}$ denotes the maximal ATP concentration. The first term in the
r.h.s. accounts for the ATP production with a characteristic timescale $\tau$,
whereas the second term describes the consumption of ATP due to spiking, with
$\epsilon$ being the amount of ATP consumed per each spike [69]. The
consumption part is assumed to be proportional to the ATP concentration to
keep the concentration from becoming negative. In a notable contrast to [69],
the dynamics of ATP concentration is introduced via a global variable rather
than the set of microscopic variables that depend on the spike trains of
individual neurons. Adopting such an approach facilitates the subsequent
derivation of the mean-field model. System (1), (4) provides a complete
description of a heterogeneous population of quadratic integrate-and-fire
neurons subjected to excitability adaptation reflecting the metabolic feedback
associated with ketogenic diet.

The difficulties with numerically simulating the reset of membrane potentials
in system (1) are resolved in the standard way by performing simulations of an
equivalent model of a population of theta neurons [12, 23], obtained by the
change of variables $V_{j}=\tan(\theta_{j}/2)$ [76, 77] in system (1). The
phase dynamics of theta neurons then reads

$\displaystyle\dot{\theta}_{j}$
$\displaystyle=(1-\cos\theta_{j})+(1+\cos\theta_{j})(\eta_{j}+KS(t)+I_{\mathrm{ext}}(t))$
$\displaystyle-\alpha\frac{\tilde{C}}{C(t)}\sin{\theta_{j}},$ (5)

which should be complemented with Eq. (4) for the dynamics of ATP
concentration. The spiking events from the representation (1) just correspond
to the phase variables’ crossing the $\pi$ value. The Lorentzian distribution
of local bifurcation parameters $\eta_{j}$ is generated deterministically via
$\eta_{j}=\bar{\eta}+\Delta\tan((\pi/2)\frac{2j-N-1}{N+1})$ [12]. The
numerical simulations of the model (5), (4) are carried out for the fixed
parameter set $\Delta=1.0,\alpha=1.0,\epsilon=1.0,\tilde{C}=1.0,N=10^{4}$ and
varying $\bar{\eta}$, $K$ and $\tau$. While $\bar{\eta}$ and $K$ respectively
characterize the intrinsic population dynamics and the impact of interactions
in modifying the local dynamics, $\tau$ may be seen as the main parameter that
describes the metabolic changes associated with the ketogenic diet. For
numerical integration of the model (5) and (4), we use an error-controlled 5/4
Runge-Kutta scheme called Tsit5 implemented in the Julia package
DifferentialEquations.jl [78]. During the simulation, the impact of the spikes
and the ATP consumption are handled separately. In particular, we determine
$S(t)$ by counting the number of spikes within a given time step and then
apply changes to the phases and the ATP concentration according to
$\theta_{j}(t)^{+}=\theta_{j}(t)^{-}+(1+\cos(\theta_{j}(t)^{-}))KS(t)$ and
$C(t)^{+}=C(t)^{-}-S(t)\epsilon C(t)^{-}/\tilde{C}$, where the superscripts
$-$ and $+$ denote the values before and after the impact of the spikes,
respectively. In a subsequent comparison with the mean-field model, we present
the mean firing rate of the population of theta neurons as a time-smoothed
signal of $S(t)$ using a rectangular convolution function with the time width
$0.1$, similar to [12].

##  III Mean-field system

Applying the general reduction approach from [1] to the microscopic model (1)-
(4), in this Section we derive the mean-field system describing the collective
behavior of a neuronal population affected by the excitability adaptation
reflecting the metabolic changes related to ketogenic diet. In the
thermodynamic limit $N\rightarrow\infty$, the population state can be
characterized by the conditional density function $\rho(V|\eta,t)$, such that
$\rho(V|\eta,t)dV$ gives the fraction of neurons with membrane potential
within the interval $(V,V+dV)$ and local bifurcation parameter $\eta$ at time
$t$. The density function satisfies the continuity equation

$\displaystyle\frac{\partial}{\partial t}\rho$
$\displaystyle=-\frac{\partial}{\partial
V}\left[\rho\left\\{V^{2}+\eta+KS+I_{\mathrm{ext}}-\alpha
V\frac{\tilde{C}}{C}\right\\}\right].$ (6)

According to the Lorentzian Ansatz [1, 12, 79], connected to the Ott-Antonsen
Ansatz via a conformal mapping [1], the solutions of (6) generically converge
to a Lorentzian distribution [18]

$\displaystyle\rho(V|\eta,t)$
$\displaystyle=\frac{1}{\pi}\frac{x(\eta,t)}{\left[V-y(\eta,t)\right]^{2}+x(\eta,t)^{2}},$
(7)

where the time-dependent parameters $x(\eta,t)$ and $y(\eta,t)$ capture the
macroscopic dynamics in a reduced subspace. In particular, the local firing
rate $r(\eta,t)$ for a fixed $\eta$ is given by the probability flux through
the infinity threshold
$r(\eta,t)=\rho(V\to\infty|\eta,t)\dot{V}(V\to\infty|\eta,t)$. This yields the
relation to the Lorentzian half-width $x(\eta,t)=\pi r(\eta,t)$, such that the
population total firing rate satisfies [1]

$\displaystyle r(t)$
$\displaystyle=\frac{1}{\pi}\int_{-\infty}^{+\infty}x(\eta,t)g(\eta)d\eta.$
(8)

Similarly, the mean membrane potential is given by

$\displaystyle v(t)$
$\displaystyle=\int_{-\infty}^{+\infty}y(\eta,t)g(\eta)d\eta.$ (9)

Introducing the complex variable $w(\eta,t)=x(\eta,t)+iy(\eta,t)$ and
substituting the Lorentzian Ansatz (7) into the continuity equation leads to

$\displaystyle\dot{w}(\eta,t)$
$\displaystyle=i\left(\eta+Kr(t)-w(\eta,t)^{2}+I_{\mathrm{ext}}(t)\right)-\frac{\tilde{C}\alpha
w(\eta,t)}{C(t)},$ (10)

where the total output of the finite population $S(t)$ is replaced with the
mean rate $r(t)$, the population’s output in the thermodynamic limit. The
infinite set of integro-differential equations (10) holds for arbitrary choice
of the diversity distribution. Nevertheless, the particular choice of
Lorentzian distribution (2) allows for the maximal reduction of (10) by
solving the integrals over $\eta$ in Eq. (8) and Eq. (9) via the residue
theorem [1, 12, 23, 13, 14]. By analytically continuing $\eta$ into the
complex plane and closing the integration contour over the half-plane
$\operatorname{Im}(\eta)<0$, one readily finds that the integrals in
definitions (8) and (9) only depend on the value of $w$ at the pole
$\eta=\bar{\eta}-i\Delta$ of the $g(\eta)$ diversity distribution, yielding
the relation $\pi r(t)+iv(t)=w(\bar{\eta}-i\Delta,t)$. Evaluating Eq. (10) at
$\eta=\bar{\eta}-i\Delta$ and recalling Eq. (4), one arrives at the mean-field
model describing the macroscopic dynamics of the population in terms of the
mean firing rate $r(t)$, mean membrane potential $v(t)$ and the total ATP
concentration $C(t)$

$\displaystyle\dot{r}$ $\displaystyle=\Delta/\pi+(2v-\alpha\tilde{C}/C)r,$
$\displaystyle\dot{v}$ $\displaystyle=\bar{\eta}-\pi^{2}r^{2}+v^{2}+Kr-\alpha
v\tilde{C}/C+I_{\mathrm{ext}}(t),$ $\displaystyle\dot{C}$
$\displaystyle=\frac{\tilde{C}-C}{\tau}-\frac{\epsilon\,r\,C}{\tilde{C}}.$
(11)

Note that the system (11) is expected to feature a certain separation of
timescales, since the ATP production, characterized by the rate $1/\tau$,
should take place on a timescale slower than that of the local neuronal
dynamics. Nevertheless, there is no reason to a priori assume a strong
separation of timescales $1/\tau\ll 1$, so in the following we consider $\tau$
as a free control parameter. One should however bear in mind that the slow-
fast decomposition within the NGNM framework has so far provided important
insights concerning the onset of collective bursting in neuronal populations
with short-term plasticity [80] and the mechanism behind emergent excitability
[81].

By carrying out the stability and bifurcation analysis of the mean-field model
(11), in Sec. IV we demonstrate how the classical physical picture from [1]
qualitatively changes in the presence of excitability adaptation associated
with ketogenic diet.

##  IV Bifurcation analysis of the mean-field model

Figure 2: Stability and bifurcations of equilibria of (11) in the
$(\bar{\eta},\tau)$-plane for coupling strengths $K=10$ and $K=15$ in a) and
b), respectively. Shading in different colors, cf. the colorbar, indicates the
type of equilibrium (stable node, stable focus, saddle focus, saddle node)
obtained by simulating (11) from random initial conditions. Bullets denote co-
dimension two bifurcations: Bogdanov-Takens (BT), cusp (CP), and generalized
Hopf (GH). Lines indicate Hopf (HB) and saddle-node (SN) bifurcation curves
(supercritical HB – solid blue lines, subcritical HB – dashed blue lines).
Other symbols (star, triangle, diamond) respectively indicate parameter sets
in Fig. 3, Fig. 6, Fig. 7 and Fig. 8. Remaining parameters are
$\Delta=1.0,\epsilon=1.0,\alpha=1.0,\tilde{C}=1.0$.

The stability of equilibria of the macroscopic system (11) and their
bifurcations are analyzed by the numerical path-following technique using the
software package BifurcationKit [82]. The bifurcation diagrams are presented
in the $(\bar{\eta},\tau)$ parameter plane keeping the remaining parameters
fixed. The structure of the bifurcation diagrams qualitatively depends on the
coupling strength $K$, with the physical picture gaining complexity for larger
$K$. To illustrate this, we provide two examples of bifurcation diagrams,
namely for $K=10$ in Fig. 2(a) and $K=15$ in Fig. 2(b). On top of the
bifurcation diagrams are overlaid the fixed-point solutions of (11) obtained
numerically starting from random initial conditions. The classification of
macroscopic equilibria is indicated by the color-coding scheme. One finds four
different types of equilibria, including two stable ones (stable node and
stable focus) and two unstable ones (saddle and saddle focus).

The bifurcation diagram for $K=10$ in Fig. 2(a) features two types of stable
equilibria, namely a stable node and a stable focus. The microscopic structure
of these states involves asynchronous dynamics and is qualitatively the same
for arbitrary $K$, see the example in Figs. 3(a)-(d) which show the asymptotic
dynamics of $r(t)$, $v(t)$, $C(t)$ and the corresponding spike raster plot for
the corresponding microscopic model (5), respectively. From Eq. (6), it
follows that for an equilibrium $(r,v,C)=(r^{*},v^{*},C^{*})$, the population
is divided into two parts: an excitable (non-spiking) subpopulation for
$\eta<(\frac{\alpha\tilde{C}}{2C^{*}})^{2}-Kr^{*}$, characterized by the
density distribution
$\rho(V|\eta)=\delta(V-(\frac{\alpha\tilde{C}}{2C^{*}}-\sqrt{\frac{\alpha\tilde{C}}{2C^{*}}^{2}-(\eta+Kr^{*})}))$,
and a spiking subpopulation
$\eta\geq(\frac{\alpha\tilde{C}}{2C^{*}})^{2}-Kr^{*}$ with a Lorentzian-shaped
density distribution
$\rho(V|\eta)=\frac{\sqrt{\eta+Kr^{*}-(\alpha\tilde{C}/2C^{*})^{2}}}{\pi((V-\alpha\tilde{C}/(2C^{*}))^{2}+\eta+Kr^{*}-(\alpha\tilde{C}/2C^{*})^{2})}$.

The two asynchronous states coexist within a small wedge-shaped region of
bistability organized around a co-dimension two cusp point (CP), where two
branches of saddle-node bifurcations (SN) meet. Within that region, the stable
node (focus) corresponds to the lower (higher)-activity state. Outside of the
bistability domain, the mean rate of the nodal equilibrium remains generally
low, whereas the focus equilibrium corresponds to the higher-activity state
$r^{*}>0.3$ for sufficiently small $\tau$. From the neuroscience perspective,
the lower-activity state indicates the normal, i.e. homeostatic activity of
the population. The saddle-node curves for smaller (larger) $\bar{\eta}$
describe the disappearance of the higher-activity (lower-activity) state.

Figure 3: Comparison between the mean-field model (11) (dashed lines) and
microscopic (solid lines) system (5) in terms of $r(t)$, $v(t)$ and $C(t)$
show a good agreement, as well as spike raster plots obtained by simulating
the microscopic system for examples of an asynchronous state (a)-(d) and
synchronous state (e)-(h). The two states, corresponding to a stable focus and
a stable limit cycle of system (11), coexist for $\tau=8.15$,
$\bar{\eta}=-1.6$ and $K=15$ (star in Fig. 2(b)). Remaining parameters:
$N=10^{4}$, $\Delta=\epsilon=\alpha=\tilde{C}=1.0$.

In contrast to the classical result for an all-to-all synaptically coupled
population of quadratic integrate-and-fire neurons [1], in the presence of
excitability adaptation one finds a transition to synchrony mediated by the
supercritical Hopf bifurcation (HB), which destabilizes the higher-activity
state, cf. the boundaries of the green region in Fig. 2(a). However, unlike
the scenario in mixed populations of excitable units and oscillators with a
Kuramoto-type coupling [83], the branch of Hopf bifurcations does not emanate
from a Bogdanov-Takens point. Sufficiently above the Hopf bifurcation, the
local dynamics of the synchronous state consists of synchronous population
bursts. Within the macroscopic oscillation cycle, one observes an accumulation
of synchrony where more units are recruited to emit synchronous spikes as the
maxima of $r(t)$ and $v(t)$ are approached, see the example in Figs. 3(e)-(h).
In the neuroscience context, the states with an increased level of synchrony,
such as this one, may be interpreted as describing seizure-like dynamics. In
particular, in our model, each maximal synchronous burst can be associated
with a seizure event. Note that for decreasing adaptation rate (increasing
$\tau$), one finds a reentrant transition to asynchronous state, in the sense
that for higher $\tau$ at fixed $\bar{\eta}$, the focus equilibrium regains
stability via an inverse supercritical Hopf bifurcation.

The physical picture substantially changes for a larger coupling strength
$K=15$, see Fig. 2(b). In particular, the bifurcation diagram is organized
around seven co-dimension two bifurcation points, including a cusp and pairs
of Bogdanov-Takens points (BT), generalized Hopf bifurcations (GH) and fold-
homoclinic points (not shown). The most important novelty is that in the
presence of excitability adaptation and for sufficiently large $K$, both the
lower-activity and the higher-activity asynchronous state may undergo a
transition to synchrony. In contrast to Fig. 2(a), these transitions are
associated with branches of Hopf bifurcation curves emanating from the
Bogdanov-Takens points. Unlike the classical scenario described in [12, 34],
the Hopf curves exhibit a change of criticality in generalized Hopf (Bautin)
points, turning from subcritical to supercritical. From the GH points derive
curves of folds of limit cycles, which are not explicitly shown due to
additional structures associated with them. Nevertheless, for both GH points,
we have been able to numerically corroborate the coexistence between the
respective stable equilibrium and the stable limit cycle just above the
subcritical Hopf bifurcations. In the neuroscience context, the coexistence
between asynchronous and synchronous states has interesting implications in
the context of potential control of seizure dynamics, as elaborated in Sec. V.

Within the wedge-shaped region organized around the CP, the bistability
domains between asynchronous and synchronous states are confined by the
branches of homoclinic bifurcations. There, the limit cycles derived from HBs
associated with the BTs vanish. Nevertheless, the homoclinic bifurcation
curves are not shown to avoid overburdening an already complicated bifurcation
diagram. We just mention that the homoclinic bifurcation curve associated with
the disappearance of the limit cycle born from the destabilization of the
higher-activity asynchronous state extends toward the right curve of saddle-
nodes for smaller $\tau$. On the other hand, the homoclinic bifurcation
associated with the vanishing of the limit cycle derived from the
destabilization of the lower-activity asynchronous state remains close to the
cusp point, terminating at a fold-homoclinic point on the left branch of the
saddle-nodes near the CP. Apart from inside the wedge-shaped region,
bistability between asynchronous and synchronous states is also found
immediately above the right (higher $\bar{\eta}$) branch of subcritical Hopf
bifurcations, cf. the right dashed line in Fig. 2(b).

We believe that the qualitative change of the bifurcation diagram observed for
varying $K$ is due to the system being close to a co-dimension three
bifurcation point. In particular, we suspect that for increasing $K$, there is
a scenario of unfolding of a degenerate BT bifurcation, after which the
supercritical Hopf curve and the bistability tongue detach from each other.
Consequently, there should exist a critical coupling strength, where the two
BT points coincide and subsequently disappear.

In Sec. V, we demonstrate the different strategies that may be applied to
induce controlled switches between the seizure-like and homeostatic dynamical
regimes.

##  V Control of switching dynamics

Developing control strategies to suppress the seizure-like states and promote
the homeostatic ones is of fundamental importance to potential applications in
neuroscience. In Sec. IV, we have shown that the presence of excitability
adaptation associated with the metabolic changes triggered by ketogenic diet
qualitatively enriches the bifurcation diagram compared to the classical
scenarios [1, 12] where such adaptation mechanism is absent. This naturally
raises the question of whether the diet-related emergent dynamics may
facilitate new means to promote the lower-activity asynchronous states at the
expense of regimes featuring excessive synchrony.

It turns out that there are indeed several control strategies available to
induce the desired regime switches. These strategies may be cast into two
groups: one associated with inducing critical transitions, and the other
related to triggering switches between coexisting states away from
criticality. In particular, a supercritical or a subcritical (hysteretic)
transition may be induced by the parametric perturbation of the ATP production
rate, i.e. the production time constant $\tau$. On the other hand, within the
domains supporting bistability between the synchronous and asynchronous
states, one may introduce an external stimulation current or apply a dynamic
perturbation to ATP concentration level to induce the desired switch of
states.

Let us first consider the scenarios concerning the variation of ATP production
rate. For lower $K$, the bifurcation diagram in Fig. 2(a) features a confined
synchrony domain (green region) bounded by a reentrant transition to the
asynchronous state. In the context of neuroscience, the fact that the
transition is mediated by the supercritical Hopf bifurcation implies less
sharp onset or termination of seizure-like states. In order to induce a
critical transition from synchrony to the asynchronous state, one is required
to systematically increase $\tau$, i.e. decrease the ATP production rate,
which could in principle be achieved by modulating any number of steps along
the metabolic pathway, including the carbohydrate source via the ketogenic
diet. Nevertheless, given the reentrant character of the transition, the same
effect can be achieved by appropriately reducing $\tau$. In terms of
neurophysiology, the former method is still preferred, because the mean rate
of the asynchronous state decreases with $\tau$, taking rather low values
$r^{*}\approx 0.2$ immediately above the state’s reappearance.

Figure 4: Hysteretic transition between the lower-activity asynchronous state
and the synchronous state triggered by sudden changes in ATP production rate
in vicinity of subcritical Hopf bifurcation of system (11) at
$\tau_{\mathrm{HB}}\approx 8.122$. (a) variation of ATP production constant
$\tau$; (b)-(d) evolution of macroscopic variables $r(t)$, $v(t)$ and $C(t)$
from (11) (dashed lines) and their counterparts from the microscopic system
(5) (solid lines); (e) raster plot of spike times for population of $N=10^{4}$
theta neurons, and (f) system’s orbit in phase space of macroscopic variables.
Slightly above $\tau_{\mathrm{HB}}$, there is bistability between synchronous
and asynchronous states, while below $\tau_{\mathrm{HB}}$ the limit cycle
associated with synchronous state is the only attractor. Parameters: $K=15.0$,
$\bar{\eta}={-}1.6$, $\Delta=\epsilon=\alpha=\tilde{C}=1.0$, and $N=10^{4}$.

For the case of stronger couplings from Fig. 2(b), varying $\tau$ may induce a
hysteretic transition in the vicinity of subcritical Hopf bifurcations
associated with the destabilization of the lower- or the higher-activity
asynchronous state. As an example, let us focus on the case of subcritical
Hopf bifurcation controlled by reducing $\tau$ for fixed $K=15$ and
$\bar{\eta}=-1.6$. Slightly above the bifurcation, the mean-field system
features bistability between the focus equilibrium and the limit cycle
illustrated in Fig. 3 for $\tau=8.15$, whereas below the bifurcation, the
limit cycle remains the system’s only attractor. In Fig. 4, it is demonstrated
how sudden shifts in $\tau$ may induce a hysteretic transition between the
asynchronous and synchronous states. In the context of neuroscience, a
hysteretic transition between the homeostatic and seizure-like states implies
the possibility of a strong jump in the signal amplitude for the seizure onset
or termination. Note that a strong jump associated with the hysteretic
transition in general implies irreversibility, in the sense that once the jump
has occurred, a small reverse of the parameter value cannot induce the
transition back to the initial state. In terms of controlling the seizure-like
activity, such irreversibility may be favorable because the metabolic
interventions that recover the homeostatic state could be persistent and have
a lasting therapeutic effect.

The parametric perturbation of the ATP production rate, or rather the ATP
production time constant $\tau$, is introduced as follows. Starting from the
asynchronous state, a sudden decrease of $\tau$ to the value $\tau=7.65$ below
the subcritical Hopf bifurcation at $\tau=\tau_{\mathrm{HB}}$, cf. the dashed
line in Fig. 4(a), triggers the switch to the synchronous state. Then, a
sudden increase of $\tau$ to the value slightly above $\tau_{\mathrm{HB}}$ is
insufficient to induce a regime shift, i.e. the seizure-like state persists.
To suppress the synchronous activity and reinstate the asynchronous state, one
is required to introduce an additional parameter perturbation such that $\tau$
exceeds the fold of cycles bifurcation that gives rise to the system’s
bistability above $\tau_{\mathrm{HB}}$. One observes a good matching between
the time series of macroscopic variables $r(t)$, $v(t)$ and $C(t)$, indicated
by the solid lines in Fig. 4(b)-(d), and the corresponding time traces (dashed
lines) obtained for the population of theta neurons given by the microscopic
model (5).

Figure 5: Continuous variation of $\tau$ from $\tau=8.15$ to $\tau=7.85$ and
back (a) induces a hysteretic transition between the lower-activity
asynchronous state and the synchronous state in the vicinity of the
subcritical Hopf bifurcation of system (11) at $\tau_{\mathrm{HB}}\approx
8.122$ described in Fig. 4. Time traces of the macroscopic variables $r(t)$,
$v(t)$ and $C(t)$ from (11) (dashed lines) are compared to their counterparts
from the microscopic system (5) (solid lines) in (b)-(d), whereas the spike
raster plot for (5) is shown in (e). One observes a delayed (dynamic)
bifurcation to synchronous state as $\tau$ is decreased. Parameters: $K=15.0$,
$\bar{\eta}={-}1.6$, $\Delta=\epsilon=\alpha=\tilde{C}=1.0$, and $N=10^{4}$.

Since implementing sudden parameter changes may not be possible in real
experiments, we have also examined a hysteretic transition between
asynchronous and synchronous states for a gradual variation in $\tau$, see
Fig. 5. The stimulation protocol consists in gradually varying $\tau$ from
$\tau=8.15$ to $\tau=7.85$ through the subcritical Hopf bifurcation and back,
cf. Fig. 5(a). One observes that the system switches from the asynchronous to
the synchronous state and remains there as $\tau$ is increased to its initial
value $\tau=8.15$. Moreover, the transition to the synchronous state when
decreasing $\tau$ does not occur immediately after reaching the value
$\tau=\tau_{\mathrm{HB}}$, which conforms to the effect of delayed (dynamic)
bifurcation [84, 85], classically observed for slow passages through the
bifurcation threshold. While the asymptotic dynamics is captured well by the
mean-field model, the transient associated with the delayed bifurcation does
not entirely match the behavior of the microscopic system, cf. Figs. 5(b)-(d).
Note that the similar hysteretic scenario is observed when considering the
subcritical Hopf bifurcation within the wedge-shaped region in Fig. 2(b).

In the cases when the controlled manipulation of ATP production rate is
inaccessible or difficult to achieve, the switches between the seizure-like
and homeostatic states may be triggered by introducing appropriate external
stimulation currents, or by inducing ATP shocks. The latter conform to the
dynamic perturbations of the ATP concentration level, which can be realized in
experiments by simply adding or washing out ATP as desired. Both methods can
be applied in parameter regions supporting bistability between the synchronous
and asynchronous states away from criticality.

In Fig. 6, it is demonstrated how the excitatory or inhibitory external pulse
currents $I_{\mathrm{e}xt}(t)$ may cause the switches in population dynamics.
In particular, the inhibitory pulse of sufficient duration can suppress
synchrony by inducing the transition to the coexisting asynchronous state. The
regime shifts are permanent in the sense that the altered regimes persist
after the external stimulation is terminated. Note that the mean-field system
captures both the asymptotic and transient dynamics of the microscopic system,
see Figs. 6(b)-(d).

Figure 6: Switch from the lower-activity asynchronous state to the coexisting
synchronous state and back induced by excitatory and inhibitory external pulse
currents $I_{\mathrm{ext}}(t)$ (a), respectively. (b)-(d) Matching between the
time traces of the macroscopic model (11) and the microscopic system (5) in
terms of $r(t)$, $v(t)$ and $C(t)$; (e) Spike raster plot for (5). Parameters:
$\tau=2.9$, $K=15.0$, $\bar{\eta}={-}2.7$,
$\Delta=\epsilon=\alpha=\tilde{C}=1.0$ (diamond in Fig. 2(b)) and $N=10^{4}$.

A similar effect can be accomplished by ATP shocks, see Fig. 7, which comprise
sharp externally-induced increases (positive ATP shocks) or decreases
(negative ATP shocks) of the ATP concentration level. ATP shocks may
permanently alter the system dynamics by switching the system between the
attraction basins of the coexisting states. For the parameters in Fig. 7, the
asynchronous state is the lower-activity equilibrium (stable focus) of system
(11), whereas the synchronous state is born from the destabilization of the
higher-activity asynchronous state. One observes that applying a positive ATP
shock of finite width to the asynchronous state is required to trigger the
switch to the synchronous state. On the other hand, even an instantaneous
negative ATP shock of small amplitude is sufficient to suppress the seizure-
like state featuring excessive synchrony and promote the transition back to
the homeostatic regime. In principle, the sensitivity of an attractor to
switch under dynamical perturbation to a coexisting attractor depends on how
close the given attractor lies to the manifold separating between the
respective basins of attraction. Note that the mean-field model captures well
the impact of the dynamic perturbation associated with ATP shocks, both in
terms of transients and the long-term dynamics, see Figs. 7(a)-(c).

Figure 7: Switches between the coexisting seizure-like and homeostatic states
triggered by sudden changes in ATP level, indicated by arrows in (c). A
positive ATP shock of finite duration induces the switch from the homeostatic
to the seizure-like state. The homeostatic state is regained following an
instantaneous, small-amplitude negative ATP shock to the seizure-like state.
(b)-(d) Comparison between the time traces $r(t)$, $v(t)$ and $C(t)$ for the
mean-field model (11) (dashed lines) and the corresponding variables of the
microscopic system (5) (solid lines); (e) spike raster plot for (5).
Parameters: $\tau=2.9$, $K=15.0$, $\bar{\eta}={-}2.7$,
$\Delta=\epsilon=\alpha=\tilde{C}=1.0$ (diamond in Fig. 2(b)) and $N=10^{4}$.

External stimulation currents and ATP shocks can also be applied to trigger
switches between the coexisting higher- and lower-activity asynchronous
states. An example in Fig. 8 first shows how a brief negative ATP shock
induces the regime shift from the higher- to the lower-activity equilibrium.
Then, the higher-activity asynchronous state is reinstated by applying a brief
positive ATP shock to the lower-activity equilibrium. Note that both the long-
term dynamics and the damped oscillations to equilibrium are well described by
the mean-field system, see Figs. 8(a)-(c).

Figure 8: Switch between coexisting higher- and lower-activity asynchronous
states (stable foci of system (11)) and vice versa triggered by respective
negative and positive ATP shocks (see arrows in (c)). (a)-(c) Time traces for
the mean-field model $r(t)$, $v(t)$ and $C(t)$ from (11) (dashed lines) and
their counterparts for (5) (solid lines). (d) Spike raster plot for (5).
Parameters: $\tau=2.5$, $K=15.0$, $\bar{\eta}={-}2.65$,
$\Delta=\epsilon=\alpha=\tilde{C}=1.0$ (triangle in Fig. 2(b)) and $N=10^{4}$.

##  VI Summary and discussion

Within the last two decades, there has been a growing awareness that the
feedback of neuronal activity with energy metabolism plays a pivotal role both
in maintaining the neuronal homeostasis [56, 86, 57, 58, 87, 49, 42] and in
development of certain neurological disorders, such as epilepsy [52, 53, 54,
55] or Parkinson’s disease [64, 65]. Nevertheless, many earlier computational
studies have suffered from improperly including the energy constraints, using
oversimplified models of local dynamics or implementing too small networks
with rather strong finite-size effects. The recent theoretical breakthrough
associated with the advent of NGNM models has provided us with a fascinating
ability to address the mechanisms of neurological disorders with a strong
metabolic component, and to even gain a foothold in understanding the impact
of certain therapeutic procedures. In the present paper, we have introduced
and analyzed a model of a heterogeneous, globally coupled population of
quadratic integrate-and-fire neurons where the local excitability is modified
via an ATP-sensitive potassium current. Implementing the classical reduction
approach from [1], we have derived a three-dimensional mean-field model of
population dynamics. The model is intended as a first step toward
theoretically explaining the mechanism and apparent effectiveness of ketogenic
diet, a long-standing treatment for epilepsy, which nowadays receives a
revived interest for helping patients with a drug-resistant form of the
disease [50, 52].

Using the numerical path-following technique, we have analyzed the stability
and bifurcations of the derived mean-field model in terms of the ATP
production rate and the mean neuronal excitability. The structure of the
bifurcation diagram has been shown to qualitatively depend on the coupling
strength, gaining complexity for stronger couplings. We have found two
qualitatively different types of bifurcation diagrams associated with weaker
or stronger couplings.

For weaker couplings, the diagram features two stable macroscopic equilibria
corresponding to the lower- and the higher-activity asynchronous state, which
coexist within a small wedge-shaped region bounded by fold bifurcation curves
that meet at the cusp. Moreover, for intermediate production rates and a
sufficiently high mean excitability, the higher-activity equilibrium undergoes
a transition to synchrony via a supercritical Hopf bifurcation. Nevertheless,
the asynchronous state is regained for sufficiently small ATP production
rates, i.e. for larger $\tau$. Compared to the classical scenario in
heterogeneous populations of quadratic integrate-and-fire neurons with
instantaneous chemical synapses [1], our weak-coupling scenario is richer in
the sense that the excitability adaptation gives rise to an additional
Kuramoto-like transition to synchrony via supercritical Hopf bifurcation.
Nevertheless, it is also simpler than the scenario reported for populations of
quadratic integrate-and-fire neurons with finite synaptic time constants [12],
populations with instantaneous chemical and electrical synapses [13, 14], as
well as mixed populations of excitable and oscillatory units with Kuramoto-
type interactions [88, 83], where the transition to synchrony unfolds via the
supercritical Hopf bifurcation emanating from a co-dimension two Bogdanov-
Takens point.

For stronger couplings, we have revealed a highly complex bifurcation scenario
organized around seven co-dimension two points, namely a cusp and pairs of
Bogdanov-Takens, generalized Hopf and fold-homoclinic points. In contrast to
[12, 13, 14, 88, 83], the presence of excitability adaptation here facilitates
a transition to synchrony not only from the higher-activity asynchronous state
but also from the macroscopic equilibrium corresponding to the lower-activity
asynchronous state. Moreover, in both instances, the Hopf bifurcations
emanating from the Bogdanov-Takens point are subcritical and change character
in generalized Hopf points for higher mean population excitability.
Remarkably, subcritical transitions to collective oscillations have so far not
been observed in massively coupled populations of quadratic integrate-and-fire
neurons with only excitatory couplings. Such behavior has rather been found in
models involving some form of inhibition and/or connection sparseness [25, 31,
89], or stochastic stimulation [20]. In this context, our results indicate
that excitability adaptation may be an important additional source of
multistability in the collective dynamics of neuronal populations.

In relation to explaining the impact of the ketogenic diet as a treatment for
epilepsy, we have interpreted the macroscopic equilibrium corresponding to the
lower-activity asynchronous state as the “normal”, i.e. the desired state,
while the synchronous states have been interpreted as pathological, i.e.
seizure-like states. Note that this interpretation is different from [38], the
only study so far that has considered the onset and propagation of epileptic
seizures within the framework of next-generation neural mass models but
focuses on the impact of connectivity rather than neuronal metabolism on the
emergent dynamics. In that work, the transition to a high-activity
asynchronous state was interpreted as a seizure-like event, but this was
likely a consequence of the lack of the transition to synchrony in their
model. We have unveiled three different stimulation protocols that may induce
controlled switches between the seizure-like and homeostatic states. These
methods include parametric perturbation of the ATP production rate, the
application of an external (excitatory or inhibitory) pulse current, and a
brief pulse-like perturbation of ATP concentration (ATP-shock). The first
method consists of inducing a critical transition between the seizure-like and
the lower-activity asynchronous state, whereas the latter two methods apply
for parameter domains where the population is bistable. The effectiveness of
all three methods is corroborated both for the switching scenarios from the
seizure-like to the normal state and back. Interestingly, earlier numerical
studies on populations of Hodgkin-Huxley neurons including a similar model of
ATP-gated potassium adaptation currents have revealed only a smooth
(supercritical) transition from the synchronous (seizure-like) to the
asynchronous (normal) state with the increase of ATP production rate [70, 71,
69]. Nevertheless, our study has also revealed the possibility of hysteretic
transitions facilitated by the bistability of these states in the vicinity of
subcritical Hopf bifurcations, which evinces yet another advantage of
introducing the theoretical framework of NGNM models. Note that the
irreversible character of hysteretic transitions may prove advantageous for
controlling the seizure-like states, given that the induced transition to the
lower-activity asynchronous state may then have a persistent therapeutic
effect.

Viewed from a different perspective, the possibility of featuring either a
continuous or hysteretic transition to seizure-like states may have
implications with regard to the ability to predict the onset of seizures from
the electroencephalography recordings. Classically, for the continuous
transitions, one expects standard indicators of criticality related to
critical slowing down [85, 90, 91] to apply, whereas such methods may not hold
for hysteretic transitions. In light of the recently sparked controversy as to
whether signatures of critical slowing down can indeed be observed prior to
epileptic seizures, with arguments provided both pro [92] and contra [93], one
may speculate that certain seizures may in fact derive from the continuous,
and some from the discontinuous hysteretic transitions, as indicated by our
model.

We believe that our results, together with [38, 35], have set the stage for
opening a new direction of research aimed at theoretical understanding of
neurological disorders and associated therapeutic procedures by implementing
the next-generation neural mass models. This is expected to provide a valuable
complement to experimental studies in at least two ways. Firstly, it will
allow for a qualitative insight into multistability which may be lacking or
not be easily accessible to experiments; and secondly, it may help in
developing effective strategies to control the switching dynamics, especially
in terms of suppressing the pathological states. Concerning the research on
the ketogenic diet, the next step could involve incorporating some other
ingredients that may be related to its therapeutic effects. In particular, one
may attempt to derive a “compartmental” model of local dynamics that could
distinguish between the two different ATP production processes, namely the
classical, glycolytic one near the membrane and the diet-related one unfolding
in mitochondria. Also, one may attempt to circumvent the limitations of the
current model associated with the assumed uniformity of ATP concentration over
the population, aiming to accommodate for the neuron-specific ATP
concentration variables that depend on the spike trains of individual neurons
[69]. The microscopic system introduced in the latter way would likely show
significant deviations from the current mean-field model, similar to those
observed in [37, 40]. Finally, the immediate impact of ketone bodies on
inhibition of excitatory synaptic channels could be taken into account as a
factor modulating the synaptic efficacy. More generally, an interesting
challenge would be to understand how the classical physical picture developed
for activity of neuronal populations within the NGNM framework changes by
incorporating the different forms of homeostatic plasticity [94], which
involve negative feedback loops to constrain the network activity within the
desired physiological limit.

###### Acknowledgements.

We would like to thank Prof. Eckehard Schöll, Dr. Matthias Wolfrum, and Dr.
Predrag Janjić for useful discussions. I.F. acknowledges the funding from the
Institute of Physics Belgrade through grant by the Ministry of Science,
Technological Development and Innovation of the Republic of Serbia and the
partial support by the ANSO - Alliance of International Science Organizations
Collaborative Research Projects and Training Projects, grant number ANSO-CR-
PP-2022-05.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Montbrió, Pazó, and Roxin [2015] E. Montbrió, D. Pazó, and A. Roxin, “Macroscopic description for networks of spiking neurons,” Phys. Rev. X 5, 021028 (2015).
  * Bick  _et al._ [2020] C. Bick, M. Goodfellow, C. R. Laing, and E. A. Martens, “Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: a review,” J. Math. Neurosci. 10, 9 (2020).
  * Coombes, , and Byrne [2019] S. Coombes, , and Á. Byrne, “Next generation neural mass models,” in  _Nonlinear Dynamics in Computational Neuroscience_, edited by F. Corinto and A. Torcini (Springer International Publishing, Cham, 2019) pp. 1–16.
  * Byrne  _et al._ [2020] A. Byrne, R. D. O’Dea, M. Forrester, J. Ross, and S. Coombes, “Next-generation neural mass and field modeling,” J. Neurophysiol. 123, 726 (2020).
  * Coombes [2023] S. Coombes, “Next generation neural population models,” Front. Appl. Math. Stat. 9, 1128224 (2023).
  * Coombes  _et al._ [2014] S. Coombes, P. beim Graben, R. Potthast, and J. Wright,  _Neural Fields:Theory and Applications_ (Springer Berlin, Heidelberg, 2014).
  * Deco  _et al._ [2008] G. Deco, V. K. Jirsa, P. A. Robinson, M. Breakspear, and K. Friston, “The dynamic brain: From spiking neurons to neural masses and cortical fields,” PLOS Comput. Biol. 4, e1000092 (2008).
  * Deschle  _et al._ [2021] N. Deschle, J. Ignacio Gossn, P. Tewarie, B. Schelter, and A. Daffertshofer, “On the validity of neural mass models,” Front. Comput. Neurosci. 14, 581040 (2021).
  * Ott and Antonsen [2008] E. Ott and T. M. Antonsen, “Low dimensional behavior of large systems of globally coupled oscillators,” Chaos 18, 037113 (2008).
  * Ott and Antonsen [2009] E. Ott and T. M. Antonsen, “Long time evolution of phase oscillator systems,” Chaos 19, 023117 (2009).
  * Luke, Barreto, and So [2013] T. B. Luke, E. Barreto, and P. So, “Complete Classification of the Macroscopic Behavior of a Heterogeneous Network of Theta Neurons,” Neural Comput. 25, 3207 (2013).
  * Ratas and Pyragas [2016] I. Ratas and K. Pyragas, “Macroscopic self-oscillations and aging transition in a network of synaptically coupled quadratic integrate-and-fire neurons,” Phys. Rev. E 94, 032215 (2016).
  * Pietras  _et al._ [2019] B. Pietras, F. Devalle, A. Roxin, A. Daffertshofer, and E. Montbrió, “Exact firing rate model reveals the differential effects of chemical versus electrical synapses in spiking networks,” Phys. Rev. E 100, 042412 (2019).
  * Montbrió and Pazó [2020] E. Montbrió and D. Pazó, “Exact mean-field theory explains the dual role of electrical synapses in collective synchronization,” Phys. Rev. Lett. 125, 248101 (2020).
  * Goldobin [2021] 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 31, 083112 (2021).
  * di Volo  _et al._ [2022] M. di Volo, M. Segneri, D. S. Goldobin, A. Politi, and A. Torcini, “Coherent oscillations in balanced neural networks driven by endogenous fluctuations,” Chaos 32, 023120 (2022).
  * Clusella and Montbrió [2022] P. Clusella and E. Montbrió, “Effect of cauchy noise on a network of quadratic integrate-and-fire neurons with non-cauchy heterogeneities,”  (2022), https://doi.org/10.48550/arXiv.2208.05515.
  * Pietras, Cestnik, and Pikovsky [2023] B. Pietras, R. Cestnik, and A. Pikovsky, “Exact finite-dimensional description for networks of globally coupled spiking neurons,” Phys. Rev. E 107, 024315 (2023).
  * Pyragas and Pyragas [2023] V. Pyragas and K. Pyragas, “Effect of cauchy noise on a network of quadratic integrate-and-fire neurons with non-cauchy heterogeneities,” Phys. Lett. A 480, 128972 (2023).
  * Goldobin, di Volo, and Torcini [2021] D. S. Goldobin, M. di Volo, and A. Torcini, “Reduction methodology for fluctuation driven population dynamics,” Phys. Rev. Lett. 127, 038301 (2021).
  * Klinshov and Kirillov [2022] V. V. Klinshov and S. Y. Kirillov, “Shot noise in next-generation neural mass models for finite-size networks,” Phys. Rev. E 106, L062302 (2022).
  * Klinshov, Smelov, and Kirillov [2023] V. V. Klinshov, P. S. Smelov, and S. Y. Kirillov, “Constructive role of shot noise in the collective dynamics of neural networks,” Chaos 33, 061101 (2023).
  * Ratas and Pyragas [2017] I. Ratas and K. Pyragas, “Symmetry breaking in two interacting populations of quadratic integrate-and-fire neurons,” Phys. Rev. E 96, 042212 (2017).
  * Schmidt  _et al._ [2018] H. Schmidt, D. Avitabile, E. Montbrió, and A. Roxin, “Network mechanisms underlying the role of oscillations in cognitive tasks,” PLOS Comput. Biol. 14, e1006430 (2018).
  * Pyragas, Fedaravičius, and Pyragienė [2021] K. Pyragas, A. P. Fedaravičius, and T. Pyragienė, “Suppression of synchronous spiking in two interacting populations of excitatory and inhibitory quadratic integrate-and-fire neurons,” Phys. Rev. E 104, 014203 (2021).
  * Byrne, Avitabile, and Coombes [2019] A. Byrne, D. Avitabile, and S. Coombes, “Next-generation neural field model: The evolution of synchrony within patterns and waves,” Phys. Rev. E 99, 012313 (2019).
  * Schmidt and Avitabile [2020] H. Schmidt and D. Avitabile, “Bumps and oscillons in networks of spiking neurons,” Chaos 30, 033133 (2020).
  * Byrne  _et al._ [2022] A. Byrne, J. Ross, R. Nicks, and S. Coombes, “Mean-Field Models for EEG/MEG: From Oscillations to Waves,” Brain Topogr. 35, 36 (2022).
  * Devalle, Roxin, and Montbrió [2017] F. Devalle, A. Roxin, and E. Montbrió, “Firing rate equations require a spike synchrony mechanism to correctly describe fast oscillations in inhibitory networks,” PLOS Comput. Biol. 13, e1005881 (2017).
  * Dumont and Gutkin [2019] G. Dumont and B. Gutkin, “Macroscopic phase resetting-curves determine oscillatory coherence and signal transfer in inter-coupled neural circuits,” PLOS Computational Biology 15, e1007019 (2019).
  * Bi  _et al._ [2020] H. Bi, M. Segneri, M. di Volo, and A. Torcini, “Coexistence of fast and slow gamma oscillations in one population of inhibitory spiking neurons,” Phys. Rev. Res. 2, 013042 (2020).
  * Segneri  _et al._ [2020] M. Segneri, H. Bi, S. Olmi, and A. Torcini, “Theta-nested gamma oscillations in next generation neural mass models,” Front. Comput. Neurosci. 14, 47 (2020).
  * Keeley  _et al._ [2019] S. Keeley, A. Byrne, A. Fenton, and J. Rinzel, “Firing rate models for gamma oscillations,” J. Neurophysiol. 121, 2181 (2019).
  * Ceni  _et al._ [2020] A. Ceni, S. Olmi, A. Torcini, and D. Angulo-Garcia, “Cross frequency coupling in next generation inhibitory neural mass models,” Chaos 30, 053121 (2020).
  * Byrne, Coombes, and Liddle [2019] Á. Byrne, S. Coombes, and P. F. Liddle, “A neural mass model for abnormal beta-rebound in schizophrenia,” in  _Handbook of Multi-Scale Models of Brain Disorders_, edited by V. Cutsuridis (Springer, Cham., 2019) pp. 21–27.
  * Weerasinghe  _et al._ [2019] G. Weerasinghe, B. Duchet, H. Cagnan, P. Brown, C. Bick, and R. Bogacz, “Predicting the effects of deep brain stimulation using a reduced coupled oscillator model,” PLOS Comput. Biol. 15, e1006575 (2019).
  * Taher, Torcini, and Olmi [2020] H. Taher, A. Torcini, and S. Olmi, “Exact neural mass model for synaptic-based working memory,” PLOS Comput. Biol. 16, e1008533 (2020).
  * Gerster  _et al._ [2021] M. Gerster, H. Taher, A. Škoch, J. Hlinka, M. Guye, F. Bartolomei, V. Jirsa, A. Zakharova, and S. Olmi, “Patient-specific network connectivity combined with a next generation neural mass model to test clinical hypothesis of seizure propagation,” Front. Syst. Neurosci. 15, 675272 (2021).
  * Gast, Schmidt, and Knösche [2020] R. Gast, H. Schmidt, and T. R. Knösche, “A Mean-Field Description of Bursting Dynamics in Spiking Neural Networks with Short-Term Adaptation,” Neural Comput. 32, 1615 (2020).
  * Gast, Knösche, and Schmidt [2021] R. Gast, T. R. Knösche, and H. Schmidt, “Mean-field approximations of networks of spiking neurons with short-term synaptic plasticity,” Phys. Rev. E 104, 044310 (2021).
  * Ferrara  _et al._ [2023] A. Ferrara, D. Angulo-Garcia, A. Torcini, and S. Olmi, “Population spiking and bursting in next-generation neural masses with spike-frequency adaptation,” Phys. Rev. E 107, 024311 (2023).
  * Sawicki  _et al._ [2023] J. Sawicki, R. Berner, S. A. M. Loos, M. Anvari, R. Bader, W. Barfuss, N. Botta, N. Brede, I. Franović, D. J. Gauthier, S. Goldt, A. Hajizadeh, P. Hövel, O. Karin, P. Lorenz-Spreen, C. Miehl, J. Mölter, S. Olmi, E. Schöll, A. Seif, P. A. Tass, G. Volpe, S. Yanchuk, and J. Kurths, “Perspectives on adaptive dynamical systems,” Chaos 33, 071501 (2023).
  * Abbott and Nelson [2000] L. Abbott and S. Nelson, “Synaptic plasticity: taming the beast,” Nat. Neurosci. 3, 1178 (2000).
  * Caporale and Dan [2008] N. Caporale and Y. Dan, “Spike timing–dependent plasticity: A hebbian learning rule,” Annu. Rev. Neurosci. 31, 25 (2008).
  * Ha and Cheong [2017] G. E. Ha and E. Cheong, “Spike frequency adaptation in neurons of the central nervous system,” Exp. Neurobiol. 26, 179 (2017).
  * Fuhrmann, Markram, and Tsodyks [2002] G. Fuhrmann, H. Markram, and M. Tsodyks, “Spike frequency adaptation and neocortical rhythms,” J. Neurophysiol. 88, 761 (2002).
  * Katsu-Jiménez, Alves, and Giménez-Cassina [2017] Y. Katsu-Jiménez, R. M. P. Alves, and A. Giménez-Cassina, “Food for thought: Impact of metabolism on neuronal excitability,” Exp. Cell Res. 360, 41–46 (2017).
  * Fardet and Levina [2020] T. Fardet and A. Levina, “Simple models including energy and spike constraints reproduce complex activity patterns and metabolic disruptions,” PLOS Comput. Biol. 16, e1008503 (2020).
  * Franović  _et al._ [2022] I. Franović, S. Eydam, S. Yanchuk, and R. Berner, “Collective activity bursting in a population of excitable units adaptively coupled to a pool of resources,” Front. Netw. Physiol. 2, 841829 (2022).
  * Lutas and Yellen [2013] A. Lutas and G. Yellen, “The ketogenic diet: metabolic influences on brain excitability and epilepsy,” Trends Neurosci. 36, 32 (2013).
  * Wheless [2008] J. W. Wheless, “History of the ketogenic diet,” Epilepsia 49, 3 (2008).
  * Fei  _et al._ [2020] Y. Fei, R. Shi, Z. Song, and J. Wu, “Metabolic control of epilepsy: A promising therapeutic target for epilepsy,” Front. Neurol. 11, 592514 (2020).
  * Rho [2017] J. M. Rho, “How does the ketogenic diet induce anti-seizure effects?” Neurosci. Lett. 637, 4 (2017).
  * D’Andrea Meira  _et al._ [2019] I. D’Andrea Meira, T. T. Romão, H. J. P. do Prado, L. T. Krüger, M. E. P. Pires, and P. O. da Conceição, “Ketogenic diet and epilepsy: What we know so far,” Front. Neurosci. 13, 5 (2019).
  * El-Rashidy  _et al._ [2023] O. F. El-Rashidy, M. F. Nassar, W. A. Shokair, and Y. G. A. El Gendy, “Ketogenic diet for epilepsy control and enhancement in adaptive behavior,” Sci. Rep. 13, 2102 (2023).
  * Roberts  _et al._ [2014] J. A. Roberts, K. K. Iyer, S. Vanhatalo, and M. Breakspear, “Critical role for resource constraints in neural models,” Front. Syst. Neurosci 8, 154 (2014).
  * Virkar  _et al._ [2016] Y. S. Virkar, W. L. Shew, J. G. Restrepo, and E. Ott, “Feedback control stabilization of critical dynamics via resource transport on multilayer networks: How glia enable learning dynamics in the brain,” Phys. Rev. E 94, 042310 (2016).
  * Kroma-Wiley, Mucha, and Bassett [2021] K. A. Kroma-Wiley, P. J. Mucha, and D. S. Bassett, “Synchronization of coupled kuramoto oscillators under resource constraints,” Phys. Rev. E 104, 014211 (2021).
  * Glynn [2002] I. M. Glynn, “A hundred years of sodium pumping,” Annu. Rev. Physiol. 64, 1 (2002).
  * Forrest [2014] M. D. Forrest, “The sodium-potassium pump is an information processing element in brain computation,” Front. Physiol. 5, 472 (2014).
  * Bazzigaluppi  _et al._ [2017] P. Bazzigaluppi, E. A. A., I. Weisspapir, B. Stefanovic, and P. L. Carlen, “Hungry neurons: Metabolic insights on seizure dynamics,” Int. J. Mol. Sci. 18, 2269 (2017).
  * Patel [2018] M. Patel, “A metabolic paradigm for epilepsy,” Epilepsy Curr. 18, 318 (2018).
  * Kovács  _et al._ [2018] R. Kovács, Z. Gerevich, A. Friedman, J. Otáhal, O. Prager, S. Gabriel, and N. Berndt, “Bioenergetic mechanisms of seizure control,” Front. Cell. Neurosci. 12, 335 (2018).
  * Büeler [2009] H. Büeler, “Impaired mitochondrial dynamics and function in the pathogenesis of parkinson’s disease,” Exp. Neurol. 218, 235 (2009).
  * Haddad and Nakamura [2015] D. Haddad and K. Nakamura, “Understanding the susceptibility of dopamine neurons to mitochondrial stressors in parkinson’s disease,” FEBS Lett. 589, 3702 (2015).
  * Beghi  _et al._ [2019] E. Beghi _et al._ , “Global, regional, and national burden of epilepsy, 1990-2016: a systematic analysis for the global burden of disease study 2016,” Lancet Neurol. 18, 357 (2019).
  * Picot  _et al._ [2008] M.-C. Picot, M. Baldy-Moulinier, J.-P. Daurs, P. Dujols, and A. Crespel, “The prevalence of epilepsy and pharmacoresistant epilepsy in adults: A population-based study in a western european country,” Epilepsia 49, 1230 (2008).
  * Neal and Cross [2010] E. G. Neal and J. H. Cross, “Efficacy of dietary treatments for epilepsy,” J. Hum. Nutr. Diet. 23, 113–119 (2010).
  * Joo  _et al._ [2021] P. Joo, H. Lee, S. Wang, S. Kim, and A. G. Hudetz, “Network model with reduced metabolic rate predicts spatial synchrony of neuronal activity,” Front. Comput. Neurosci. 15, 738362 (2021).
  * Cunningham  _et al._ [2006] M. O. Cunningham, D. D. Pervouchine, C. Racca, N. J. Kopell, C. H. Davies, R. S. G. Jones, R. D. Traub, and M. A. Whittington, “Neuronal metabolism governs cortical network response state,” PNAS 103, 5597 (2006).
  * Ching  _et al._ [2012] S. Ching, P. L. Purdon, S. Vijayan, N. J. Kopell, and E. N. Brown, “A neurophysiological–metabolic model for burst suppression,” PNAS 109, 3095 (2012).
  * Ermentrout and Terman [2010] G. B. Ermentrout and D. H. Terman,  _Mathematical Foundations of Neuroscience_ (Springer, New York, 2010).
  * Izhikevich [2007] E. M. Izhikevich,  _Dynamical systems in neuroscience: the geometry of excitability and bursting_ (The MIT Press, Cambridge, 2007).
  * Huang  _et al._ [2007] C.-W. Huang, C.-C. Huang, J.-T. Cheng, J.-J. Tsai, and S.-N. Wu, “Glucose and hippocampal neuronal excitability: Role of atp-sensitive potassium channels,” J. Neurosci. Res. 85, 1468 (2007).
  * Martinet  _et al._ [2017] L. E. Martinet, G. Fiddyment, J. R. Madsen, E. N. Eskandar, W. Truccolo, U. T. Eden, S. S. Cash, and M. A. Kramer, “Human seizures couple across spatial scales through travelling wave dynamics,” Nat. Commun. 8, 14896 (2017).
  * Ermentrout and Kopell [1986] G. B. Ermentrout and N. Kopell, “Parabolic bursting in an excitable system coupled with a slow oscillation,” SIAM J. Appl. Math. 46, 233 (1986).
  * Ermentrout [1996] B. Ermentrout, “Type I Membranes, Phase Resetting Curves, and Synchrony,” Neural Comput. 8, 979 (1996).
  * Rackauckas and Nie [2017] C. Rackauckas and Q. Nie, “DifferentialEquations.jl–a performant and feature-rich ecosystem for solving differential equations in Julia,” Journal of Open Research Software 5 (2017).
  * Devalle, Montbrió, and Pazó [2018] F. Devalle, E. Montbrió, and D. Pazó, “Dynamics of a large system of spiking neurons with synaptic delay,” Phys. Rev. E 98, 042214 (2018).
  * Taher, Avitabile, and Desroches [2022] H. Taher, D. Avitabile, and M. Desroches, “Bursting in a next generation neural mass model with synaptic dynamics: a slow–fast approach,” Nonlinear Dyn. 108, 4261 (2022).
  * Avitabile, Desroches, and Ermentrout [2022] D. Avitabile, M. Desroches, and G. B. Ermentrout, “Cross-scale excitability in networks of quadratic integrate-and-fire neurons,” PLOS Comput. Biol. 18, e1010569 (2022).
  * Veltz [2020] R. Veltz, “Bifurcationkit.jl,”  (2020), https://hal.archives-ouvertes.fr/hal-02902346.
  * Klinshov and Franović [2019] V. Klinshov and I. Franović, “Two scenarios for the onset and suppression of collective oscillations in heterogeneous populations of active rotators,” Phys. Rev. E 100, 062211 (2019).
  * Benoît [1991] E. Benoît, ed.,  _Dynamic Bifurcations_, Lecture Notes in Mathematics, vol. 1493 (Springer, Berlin Heidelberg, 1991).
  * Kuehn [2011] C. Kuehn, “A mathematical framework for critical transitions: Bifurcations, fast–slow systems and stochastic dynamics,” Physica D: Nonlinear Phenomena 240, 1020 (2011).
  * Fields [2015] R. D. Fields, “A new mechanism of nervous system plasticity: activity-dependent myelination,” Nat. Rev. Neurosci. 16, 756 (2015).
  * Attwell and Laughlin [2001] D. Attwell and S. B. Laughlin, “An energy budget for signaling in the grey matter of the brain,” J. Cereb. Blood Flow Metab. 21, 1133 (2001).
  * Lafuerza, Colet, and Toral [2010] L. F. Lafuerza, P. Colet, and R. Toral, “Nonuniversal results induced by diversity distribution in coupled excitable systems,” Phys. Rev. Lett. 105, 084101 (2010).
  * Ostojic, Brunel, and Hakim [2009] S. Ostojic, N. Brunel, and V. Hakim, “Synchronization properties of networks of electrically coupled neurons in the presence of noise and heterogeneities,” J. Comput. Neurosci. 26, 369 (2009).
  * Kuehn [2013] C. Kuehn, “A mathematical framework for critical transitions: Normal forms, variance and applications,” J. Nonlinear Sci. 23, 457 (2013).
  * Meisel  _et al._ [2015] C. Meisel, A. Klaus, C. Kuehn, and D. Plenz, “Critical slowing down governs the transition to neuron spiking,” PLOS Comput. Biol. 11, e1004097 (2015).
  * Maturana  _et al._ [2020] M. I. Maturana, C. Meisel, K. Dell, P. J. Karoly, W. D’Souza, D. B. Grayden, A. N. Burkitt, P. Jiruska, J. Kudlacek, J. Hlinka, M. J. Cook, L. Kuhlmann, and D. R. Freestone, “Critical slowing down as a biomarker for seizure susceptibility,” Nat. Commun. 11, 2172 (2020).
  * Wilkat, Rings, and Lehnertz [2019] T. Wilkat, T. Rings, and K. Lehnertz, “No evidence for critical slowing down prior to human epileptic seizures,” Chaos 29, 091104 (2019).
  * Zierenberg, Wilting, and Priesemann [2018] J. Zierenberg, J. Wilting, and V. Priesemann, “Homeostatic plasticity and external input shape neural network dynamics,” Phys. Rev. X 8, 031018 (2018).