RESEARCH ARTICLE | MAY 22 2024

Control of seizure-like dynamics in neuronal populations
with excitability adaptation related to ketogenic diet Special Collection: Regime switching in coupled nonlinear systems: sources, prediction, and control

Sebastian Eydam ; Igor Franović; Louis Kang

Chaos 34, 053128 (2024)
https://doi.org/10.1063/5.0180954


View
OnlineExportCitation

28 June 2024 05:00:50

ChaosARTICLEpubs.aip.org/aip/cha

Control of seizure-like dynamics in neuronal
populations with excitability adaptation related
to ketogenic diet
Cite as: Chaos 34, 053128 (2024); doi: 10.1063/5.0180954
Submitted: 13 October 2023 · Accepted: 6 May 2024 ·Published Online: 22 May 2024
Sebastian Eydam,1,a)Igor Franović,2,b)View Online

Export Citation

CrossMarkand Louis Kang1,c)

AFFILIATIONS
1Neural Circuits and Computations Unit, RIKEN Center for Brain Science, 2-1 Hirosawa, Wako 351-019, Japan
Scientific Computing Laboratory, Center for the Study of Complex Systems, Institute of Physics Belgrade, University of Belgrade,
Pregrevica 118, Belgrade 11080, Serbia
2

Note: This paper is part of the Focus Issue on Regime switching in coupled nonlinear systems: sources, prediction, and control.
Author to whom correspondence should be addressed: richard.eydam@riken.jpb)Electronic mail: franovic@ipb.ac.rsc)Electronic mail: louis.kang@riken.jpa)

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-fieldsystem 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 forswitching 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.
Published under an exclusive license by AIP Publishing. https://doi.org/10.1063/5.0180954

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-gatedpotassium channels in neuronal membranes. We build a modelof a heterogeneous, globally coupled population of inherently
excitable or tonic spiking excitatory quadratic integrate-and-fireChaos 34, 053128 (2024); doi: 10.1063/5.0180954
Published under an exclusive license by AIP Publishingneurons influenced by an ATP-dependent hyperpolarizing current, endowed by an additional equation describing the changes
in the total ATP concentration. The stability of solutions andtheir bifurcations are investigated by analyzing the corresponding mean-field model via the numerical path-following technique.
We disclose qualitatively different scenarios of multistability andbifurcations depending on the coupling strength. Interpreting the
low-activity asynchronous stationary state as the physiologicallydesired (homeostatic) one, and the emergence of synchronous
solutions featuring collective oscillations as the signature of
seizure-like events, we demonstrate three control strategies tosuppress the seizure dynamics by inducing a switch to the homeostatic state. Interestingly, such switches may be triggered by the34, 053128-1

28 June 2024 05:00:50

ABSTRACTChaoschanges in the ATP production rate, as a form of parametricperturbation, both for smooth or hysteretic transitions between
synchronous and asynchronous states.

I. INTRODUCTION

Chaos 34, 053128 (2024); doi: 10.1063/5.0180954
Published under an exclusive license by AIP Publishingpubs.aip.org/aip/cha

metabolism triggered by the ketogenic diet (KD).50–55 Constructing
the bifurcation diagrams to reveal critical transitions and regions ofmultistability, 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 neurologicaldisorders.49,56–58 The neuronal energy demand is typically high evenin the resting state because it requires maintaining ion fluxes againstconcentration and electric gradients.48 Such ion fluxes are mainly
perpetuated by the Na+ /K+ pump which consumes energy in theform of ATP,59,60 produced in the glycolysis process unfolding
close to the neuronal membrane. It has long been understood thatsome neurological disorders, such as epilepsy47,61–63 or Parkinson’s
disease,64,65 involve a strong metabolic component. Epilepsy is thethird most prevalent neurological disorder66,67 and is characterized
by recurring seizures. Despite the advanced anti-convulsive medications, about a third of patients still remain drug-resistant andrequire alternative therapies.67 Among the latter, KD,50,52–55 whose
origins date back to the 1920s, is nowadays often recommended tochildren 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 ofketosis, 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) targetedby medications.50
KD is believed to induce two different anti-seizure types ofeffects: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 theglycolytic ATP production results in the depletion of ATP in anear-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–71 Our approach involves a model of local neuronal dynamics supplied with the hyperpolarizing current associatedwith 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 themodel of a globally coupled heterogeneous population of intrinsically excitable or tonic spiking excitatory quadratic integrateand-fire neurons supplied by the individual terms for the ATPdependent hyperpolarizing currents and an equation for the global
dynamics of ATP concentration. Using the approach from Ref. 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

34, 053128-2

28 June 2024 05:00:50

The recent advent of the next-generation neural mass (NGNM)
models1–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 thecontrol strategies for switching between the states interpreted interms of neurophysiology as regular or pathological ones. Earlier
phenomenological neural mass models6–8 described the dynamics
of neuronal populations in terms of firing rate equations, typicallyassuming 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 featuresof local and collective dynamics. These issues have been resolvedin the NGNM models, based on applying the Ott–Antonsen9–11
or the equivalent Lorentzian Ansatz1 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 inneuroscience, is reflected in oscillations of the mean firing rate,12
the NGNM framework further establishes a relation between themean firing rate and the population synchrony described by thecomplex Kuramoto order parameter, a synchrony measure inheritedfrom the theory of coupled phase oscillators.1,4 NGNM models havethe advantage of being exact in the thermodynamic limit and are
amenable to bifurcation analysis. So far, NGNM models have beenused to address a broad range of theoretical problems, including
understanding the impact and interplay of chemical and electrical synapses,13–15 the effects of deterministic or stochastic drive,16–20
presence of quenched randomness in network topology,21 finite-size
fluctuations,22,23 interaction of neuronal populations,24–26 as well as
formation of bumps and waves in neural fields.27–29 Moreover, a wide
spectrum of applications has been considered from explaining theonset of gamma and theta-nested gamma oscillations,12,30–33 crossfrequency couplings,34 abnormal beta rebound in schizophrenia,35
the impact of deep brain stimulation,36 to studies on working
memory,25,37 propagation of epileptic seizures,38 and modeling of theentire 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.37,39–41 In general, neuronal systems feature two broadtypes of adaptive behavior.42 The first is synaptic plasticity,43,44 which
involves modification of coupling strengths due to feedback withneuronal activity. The other type comprises adaptation of local
kinetics, including spike frequency adaptation41,45,46 or changes inneuronal excitability due to metabolic constraints.47–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

ARTICLEChaosARTICLEpubs.aip.org/aip/cha

a spike whenever its membrane potential reaches the peak value Vp ,after which it is reset to −Vr . To allow for subsequent analyticaltractability, we set Vp = −Vr = ∞.
The terms Iext (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 areassumed to be instantaneous such that the mean synaptic activity isgiven by the normalized output signal of the population,
N

S(t) =

FIG. 1. Schematic of the excitability adaptation associated with the changes inthe 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 theresting potential and enabling spiking.

II. MODEL OF MICROSCOPIC DYNAMICS
We consider a large heterogeneous population of synapticallycoupled excitatory quadratic integrate-and-fire neurons72,73 whose
dynamics is given byV̇j = V2j + ηj + KS(t) + IATP(t) + Iext (t),jIATP(t) = −αVj (t)
j

C̃
,
C(t)

(1)

where Vj for j = {1, . . . , N}, N  1 are the neuronal membranepotentials. The population diversity is manifested through heterogeneity of local bifurcation parameters ηj , which follow theLorentzian density distributiong(η) =1
1
π (η − η̄)2 + 12

(2)

and may in terms of neurophysiology be interpreted as heterogeneous components of external bias current. Quadratic integrateand-fire neurons are paradigmatic for type I excitability such that an
isolated unit undergoes a transition from excitable regime (ηj . 0)√
to tonic spiking (ηj > 0) with a frequency fj = ηj /π via a SNIPER(saddle node of infinite period) bifurcation at ηj = 0. A neuron fires

Chaos 34, 053128 (2024); doi: 10.1063/5.0180954
Published under an exclusive license by AIP Publishing(3)

k\tj ≤twhere tkj is the time of the kth spike of neuron j.
The excitability adaptation triggered by the changes in theenergy metabolism due to the ketogenic diet is mediated by the ATPdependent hyperpolarizing currents IATP
, which may be regarded asj
an effective parametric perturbation to the local bifurcation parameters ηj . The hyperpolarizing currents reflect the action of the ATPgated potassium channels in neuronal membranes activated by the
diet50,69,74–76 and are assumed to have a uniform conductance α. Notethat the gating of the potassium channels is in fact a complex process governed by the exchange of sodium and ATP.70,71 Nevertheless,for simplicity, here we adopt the model for IATP
from Ref. 69 whichjomits the detailed equation for the sodium dynamics. To close thesystem (1), one requires an additional equation for the dynamics ofATP concentration C(t),Ċ =C̃ − C  S(t) C
,
−
τ
C̃

(4)

where C̃ denotes the maximal ATP concentration. The first termin the r.h.s. accounts for the ATP production with a characteristic
timescale τ , whereas the second term describes the consumption of
ATP due to spiking, with  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 frombecoming negative. In a notable contrast to Ref. 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 ofindividual neurons. Adopting such an approach facilitates the subsequent derivation of the mean-field model. Systems (1) and (4)provide a complete description of a heterogeneous population ofquadratic 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,13,24 obtained by the change of variables Vj= tan(θj /2)77,78 in system (1). The phase dynamics of theta neurons
then readsθ̇j = (1 − cos θj ) + (1 + cos θj )(ηj + KS(t) + Iext (t)) − αC̃sin θj ,C(t)(5)

which should be complemented with Eq. (4) for the dynamics ofATP concentration. The spiking events from representation (1)34, 053128-3

28 June 2024 05:00:50

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.


1 X X δ t − tkj ,N j=1 kChaosARTICLE

just correspond to the phase variables’ crossing the π value. The
Lorentzian distribution of local bifurcation parameters ηj is gen-.13 Theerated deterministically via ηj = η̄ + 1 tan (π/2) 2j−N−1
N+1
numerical simulations of the model (4) and (5) are carried out for
the fixed parameter set 1 = 1.0, α = 1.0,  = 1.0, C̃ = 1.0, N = 104
and varying η̄, K, and τ . While η̄ and K, respectively, characterize
the intrinsic population dynamics and the impact of interactions inmodifying the local dynamics, τ may be seen as the main parameterthat describes the metabolic changes associated with the ketogenic
diet, consistent with both in vivo and in vitro findings that thediet affects the ATP production capacity.50,79–83 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.84 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 agiven time step and then apply changes to the phases and the ATPconcentration according to θj (t)+ = θj (t)− + (1 + cos(θj (t)− ))KS(t)and C(t)+ = C(t)− − S(t)C(t)− /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 convolutionfunction with the time width 0.1, similar to.13

Applying the general reduction approach from Ref. 1 to microscopic models (1)–(4), in this section we derive the mean-fieldsystem describing the collective behavior of a neuronal population affected by the excitability adaptation reflecting the metabolicchanges related to ketogenic diet. In the thermodynamic limit
N → ∞, the population state can be characterized by the conditional density function ρ(V|η, t) such that ρ(V|η, t)dV gives thefraction of neurons with membrane potential within the interval(V, V + dV) and local bifurcation parameter η at time t. The density
function satisfies the continuity equation,
" (
)#
C̃
∂
∂
2
ρ V + η + KS + Iext − αV
.
(6)
ρ=−
∂t
∂V
C
According to the Lorentzian Ansatz,1,13,85 connected to the
Ott–Antonsen Ansatz via a conformal mapping,1 the solutions of (6)generically converge to a Lorentzian distribution,19
ρ(V|η, t) =x(η, t)1,
π V − y(η, t) 2 + x(η, t)2(7)

where the time-dependent parameters x(η, t) and y(η, t) capture themacroscopic dynamics in a reduced subspace. In particular, the local
firing rate r(η, t) for a fixed η is given by the probability flux throughthe infinity threshold r(η, t) = ρ(V → ∞|η, t)V̇(V → ∞|η, t). This
yields the relation to the Lorentzian half-width x(η, t) = π r(η, t)such that the population total firing rate satisfies1
Z
1 +∞r(t) =x(η, t)g(η)dη.
(8)
π −∞

Chaos 34, 053128 (2024); doi: 10.1063/5.0180954
Published under an exclusive license by AIP Publishing

Similarly, the mean membrane potential is given byZ +∞y(η, t)g(η)dη.
v(t) =(9)−∞

Introducing the complex variable w(η, t) = x(η, t) + iy(η, t)
and substituting the Lorentzian Ansatz (7) into the continuity
equation leads to C̃αw(η, t)ẇ(η, t) = i η + Kr(t) − w(η, t)2 + Iext (t) −
,
C(t)

(10)

where the total output of the finite population S(t) is replaced withthe 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, theparticular choice of Lorentzian distribution (2) allows for the maximal reduction of (10) by solving the integrals over η in Eqs. (8)and (9) via the residue theorem.1,13–15,24 By analytically continuingη into the complex plane and closing the integration contour overthe half-plane Im(η) < 0, one readily finds that the integrals in
definitions (8) and (9) only depend on the value of w at the poleη = η̄ − i1 of the g(η) diversity distribution, yielding the relation
π r(t) + iv(t) = w(η̄ − i1, t). Evaluating Eq. (10) at η = η̄ − i1 andrecalling 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),ṙ = 1/π + (2v − α C̃/C)r,v̇ = η̄ − π 2 r2 + v2 + Kr − αvC̃/C + Iext (t),(11)

C̃ − C  r C
Ċ =
.
−
τ
C̃
Note that system (11) is expected to feature a certain separationof timescales, since the ATP production, characterized by the rate
1/τ , 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/τ  1, so in the following, we consider τ 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 theonset of collective bursting in neuronal populations with short-term
plasticity86 and the mechanism behind emergent excitability.87
By carrying out the stability and bifurcation analysis of themean-field model (11), in Sec. IV, we demonstrate how the classical physical picture from Ref. 1 qualitatively changes in the presenceof excitability adaptation associated with ketogenic diet.
IV. BIFURCATION ANALYSIS OF THE MEAN-FIELD
MODEL
The stability of equilibria of the macroscopic system (11) andtheir bifurcations are analyzed by the numerical path-following
technique using the software package BifurcationKit.88 The bifurcation diagrams are presented in the (η̄, τ ) parameter plane keeping
the remaining parameters fixed. The structure of the bifurcation
diagrams qualitatively depends on the coupling strength K, with

34, 053128-4

28 June 2024 05:00:50

III. MEAN-FIELD SYSTEM

pubs.aip.org/aip/cha

ChaosARTICLEpubs.aip.org/aip/cha

FIG. 2. Stability and bifurcations of equilibria of (11) in the (η̄, τ )-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, and 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 and 6– 8.
Remaining parameters are 1 = 1.0,  = 1.0, α = 1.0, C̃ = 1.0.

π (V−α C̃/(2C )) +η+Kr −(α C̃/2C )The two asynchronous states coexist within a small wedgeshaped region of bistability organized around a co-dimension two

Chaos 34, 053128 (2024); doi: 10.1063/5.0180954
Published under an exclusive license by AIP Publishingcusp 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 thehigher-activity state r∗ > 0.3 for sufficiently small τ . 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) η̄ describe the disappearance of the higheractivity (lower-activity) state.
In contrast to the classical result for an all-to-all synaptically coupled population of quadratic integrate-and-fire neurons1
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 thegreen region in Fig. 2(a). However, unlike the scenario in mixed
populations of excitable units and oscillators with a Kuramototype coupling,89 the branch of Hopf bifurcations does not emanatefrom a Bogdanov–Takens point. Sufficiently above the Hopf bifurcation, the local dynamics of the synchronous state consists ofsynchronous 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) andv(t) are approached, see the example in Figs. 3(e)–3(h). In the neuroscience context, the states with an increased level of synchrony, suchas this one, may be interpreted as describing seizure-like dynamics.
In particular, in our model, each maximal synchronous burst can beassociated with a seizure event. Note that for decreasing adaptation

34, 053128-5

28 June 2024 05:00:50

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 theexample in Figs. 3(a)–3(d) which show the asymptotic dynamicsof r(t), v(t), C(t), and the corresponding spike raster plot for thecorresponding 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) sub 2
α C̃
population for η < 2C
− Kr∗ , characterized by the density dis∗



q
2
α C̃
α C̃
− (η + Kr∗ ) , and atribution ρ(V|η) = δ V − 2C
∗ −
2C∗
 2
α C̃
spiking subpopulation η ≥ 2C
− Kr∗ with a Lorentzian-shaped∗√2
η+Kr∗ −(α C̃/2C∗ ).
density distribution ρ(V|η) = 
∗ 2
∗
∗ 2

Chaosrate (increasing τ ), one finds a reentrant transition to asynchronous
state, in the sense that for higher τ at fixed η̄, 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 (GH) bifurcations, 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 loweractivity and the higher-activity asynchronous state may undergo atransition 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 Refs. 13 and 34, the Hopf curves exhibit a change ofcriticality in generalized Hopf (Bautin) points, turning from subcritical to supercritical. From the GH points derive curves of folds oflimit cycles, which are not explicitly shown due to additional structures associated with them. Nevertheless, for both GH points, wehave been able to numerically corroborate the coexistence between
the respective stable equilibrium and the stable limit cycle justabove the subcritical Hopf bifurcations. In the neuroscience context,
the coexistence between asynchronous and synchronous states has

Chaos 34, 053128 (2024); doi: 10.1063/5.0180954
Published under an exclusive license by AIP Publishingpubs.aip.org/aip/cha

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, thebistability domains between asynchronous and synchronous states
are confined by the branches of homoclinic bifurcations. There, thelimit 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 justmention that the homoclinic bifurcation curve associated with thedisappearance of the limit cycle born from the destabilization of thehigher-activity asynchronous state extends toward the right curve of
saddle-nodes for smaller τ . On the other hand, the homoclinic bifurcation associated with the vanishing of the limit cycle derived fromthe destabilization of the lower-activity asynchronous state remains
close to the cusp point, terminating at a fold-homoclinic point on theleft branch of the saddle-nodes near the CP. Apart from inside thewedge-shaped region, bistability between asynchronous and synchronous states is also found immediately above the right (higherη̄) 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 aco-dimension three bifurcation point. In particular, we suspect thatfor 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 coincideand subsequently disappear.
In Sec. V, we demonstrate the different strategies that may be
applied to induce controlled switches between the seizure-like andhomeostatic dynamical regimes.

V. CONTROL OF SWITCHING DYNAMICS
Developing control strategies to suppress the seizure-like statesand promote the homeostatic ones is of fundamental importance topotential applications in neuroscience. In Sec. IV, we have shownthat the presence of excitability adaptation associated with themetabolic changes triggered by ketogenic diet qualitatively enriches
the bifurcation diagram compared to the classical scenarios1,13 wheresuch 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 theexpense 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 becast into two groups: one associated with inducing critical transitions and the other related to triggering switches between coexistingstates away from criticality. In particular, a supercritical or a subcritical (hysteretic) transition may be induced by the parametricperturbation of the ATP production rate, i.e., the production time
constant τ . 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 ofstates.

34, 053128-6

28 June 2024 05:00:50

FIG. 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 example, of an asynchronous state [(a)–(d)] and synchronous state[(e)–(h)]. The two states, corresponding to a stable focus and a stable limit cycleof system (11), coexist for τ = 8.15, η̄ = −1.6 and K = 15 [star in Fig. 2(b)].
Remaining parameters: N = 104 , 1 =  = α = C̃ = 1.0.

ARTICLEChaosChaos 34, 053128 (2024); doi: 10.1063/5.0180954
Published under an exclusive license by AIP Publishingpubs.aip.org/aip/cha

28 June 2024 05:00:50

Let us first consider the scenarios concerning the variationof 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 thecontext 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 τ , 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 thetransition, the same effect can be achieved by appropriately reducing
τ . In terms of neurophysiology, the former method is still preferred,because the mean rate of the asynchronous state decreases with
τ , taking rather low values r∗ ≈ 0.2 immediately above the state’s
reappearance.
For the case of stronger couplings from Fig. 2(b), varying τ
may induce a hysteretic transition in the vicinity of subcritical Hopf
bifurcations associated with the destabilization of the lower- or thehigher-activity asynchronous state. As an example, let us focus on
the case of subcritical Hopf bifurcation controlled by reducing τ
for fixed K = 15 and η̄ = −1.6. Slightly above the bifurcation, themean-field system features bistability between the focus equilibrium and the limit cycle illustrated in Fig. 3 for τ = 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 τ
may induce a hysteretic transition between the asynchronous andsynchronous states. In the context of neuroscience, a hysteretic transition between the homeostatic and seizure-like states implies thepossibility of a strong jump in the signal amplitude for the seizure
onset or termination. Note that a strong jump associated with thehysteretic transition in general implies irreversibility, in the sensethat 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, orrather the ATP production time constant τ , is introduced as follows. Starting from the asynchronous state, a sudden decrease ofτ to the value τ = 7.65 below the subcritical Hopf bifurcation atτ = τHB , cf. the dashed line in Fig. 4(a), triggers the switch tothe synchronous state. Then, a sudden increase of τ to the valueslightly above τHB is insufficient to induce a regime shift, i.e., theseizure-like state persists. To suppress the synchronous activity andreinstate the asynchronous state, one is required to introduce anadditional parameter perturbation such that τ exceeds the fold ofcycles bifurcation that gives rise to the system’s bistability above τ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
Figs. 4(b)–4(d), and the corresponding time traces (dashed lines)
obtained for the population of theta neurons given by microscopic
model (5).
Since implementing sudden parameter changes may not be
possible in real experiments, we have also examined a hysteretic

ARTICLE

FIG. 4. Hysteretic transition between the lower-activity asynchronous state andthe synchronous state triggered by sudden changes in ATP production rate in
vicinity of subcritical Hopf bifurcation of system (11) at τHB ≈ 8.122. (a) Variationof ATP production constant τ ; (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 = 104
theta neurons; and (f) system’s orbit in phase space of macroscopic variables.
Slightly above τHB , there is bistability between synchronous and asynchronous
states, while below τHB the limit cycle associated with synchronous state is the
only attractor. Parameters: K = 15.0, η̄ = −1.6, 1 =  = α = C̃ = 1.0, andN = 104 .

34, 053128-7

Chaostransition between asynchronous and synchronous states for a gradual variation in τ , see Fig. 5. The stimulation protocol consists ingradually varying τ from τ = 8.15 to τ = 7.85 through the subcritical Hopf bifurcation and back, cf. Fig. 5(a). One observes that thesystem switches from the asynchronous to the synchronous state
and remains there as τ is increased to its initial value τ = 8.15.
Moreover, the transition to the synchronous state when decreasing τdoes not occur immediately after reaching the value τ = τHB , which
conforms to the effect of delayed (dynamic) bifurcation,90,91 classically observed for slow passages through the bifurcation threshold.
While the asymptotic dynamics is captured well by the mean-fieldmodel, the transient associated with the delayed bifurcation does
not entirely match the behavior of the microscopic system, cf.
Figs. 5(b)–5(d). Note that the similar hysteretic scenario is observed
when considering the subcritical Hopf bifurcation within the wedgeshaped region in Fig. 2(b).

Chaos 34, 053128 (2024); doi: 10.1063/5.0180954
Published under an exclusive license by AIP Publishingpubs.aip.org/aip/cha

FIG. 6. Switch from the lower-activity asynchronous state to the coexisting synchronous state and back induced by excitatory and inhibitory external pulse
currents Iext (t) (a), respectively. (b)–(d) Matching between the time traces of themacroscopic model (11) and the microscopic system (5) in terms of r(t), v(t), andC(t). (e) Spike raster plot for (5). Parameters: τ = 2.9, K = 15.0, η̄ = −2.7,
1 =  = α = C̃ = 1.0 [diamond in Fig. 2(b)] ,and N = 104 .

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 bysimply adding or washing out ATP as desired. Both methods canbe applied in parameter regions supporting bistability between thesynchronous and asynchronous states away from criticality.
In Fig. 6, it is demonstrated how the excitatory or inhibitory
external pulse currents Iext (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 thesense that the altered regimes persist after the external stimulation
is terminated. Note that the mean-field system captures both theasymptotic and transient dynamics of the microscopic system, see
Figs. 6(b)–6(d).

34, 053128-8

28 June 2024 05:00:50

FIG. 5. Continuous variation of τ from τ = 8.15 to τ = 7.85 and back (a)induces a hysteretic transition between the lower-activity asynchronous state andthe synchronous state in the vicinity of the subcritical Hopf bifurcation of system (11) at τHB ≈ 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 spikeraster plot for (5) is shown in (e). One observes a delayed (dynamic) bifurcation to synchronous state as τ is decreased. Parameters: K = 15.0, η̄ = −1.6,
1 =  = α = C̃ = 1.0, and N = 104 .

ARTICLEChaosVI. 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
homeostasis42,49,56–58,92,93 and in development of certain neurologicaldisorders, such as epilepsy52–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 afascinating ability to address the mechanisms of neurological disorders with a strong metabolic component and to even gain a footholdin understanding the impact of certain therapeutic procedures. In
the present paper, we have introduced and analyzed a model of aheterogeneous, globally coupled population of quadratic integrateand-fire neurons where the individual excitability is modified via an
ATP-sensitive potassium current. Implementing the classical reduction approach from Ref. 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 andChaos 34, 053128 (2024); doi: 10.1063/5.0180954
Published under an exclusive license by AIP Publishingpubs.aip.org/aip/cha

FIG. 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 theseizure-like state. The homeostatic state is regained following an instantaneous,small-amplitude negative ATP shock to the seizure-like state. (a)–(c) 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) and (d) spike raster plot for (5). Parameters: τ = 2.9, K = 15.0,η̄ = −2.7, 1 =  = α = C̃ = 1.0 [diamond in Fig. 2(b)] and N = 104 .

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 modelin terms of the ATP production rate and the mean neuronal
excitability. The structure of the bifurcation diagram has beenshown 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 orstronger 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.

34, 053128-9

28 June 2024 05:00:50

A similar effect can be accomplished by ATP shocks, see Fig. 7,
which comprise sharp externally induced increases (positive ATPshocks) 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 higheractivity 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 issufficient 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 thegiven 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)–7(c).
External stimulation currents and ATP shocks can also be
applied to trigger switches between the coexisting higher- and loweractivity asynchronous states. An example in Fig. 8 first shows how abrief negative ATP shock induces the regime shift from the higherto 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)–8(c).
The examples provided in Sec. V indicate a spectrum of strategies that may in general be applied to regain the homeostaticstate, and selecting an optimal strategy in any particular instance isdependent on the corresponding system parameters.

ARTICLEChaosMoreover, for intermediate production rates and a sufficiently highmean 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 τ . Compared to the classical scenario in heterogeneous populations of quadratic integrate-and-fireneurons 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-andfire neurons with finite synaptic time constants,13 populations withinstantaneous chemical and electrical synapses,14,15 as well as mixed
populations of excitable and oscillatory units with Kuramoto-type
interactions,89,94 where the transition to synchrony unfolds via thesupercritical 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

Chaos 34, 053128 (2024); doi: 10.1063/5.0180954
Published under an exclusive license by AIP Publishingpubs.aip.org/aip/cha

two points, namely, a cusp and pairs of Bogdanov–Takens,generalized Hopf and fold-homoclinic points. In contrast toRefs. 13–Montbrio2020, 94, and 89, the presence of excitability adaptation here facilitates a transition to synchrony not only
from the higher-activity asynchronous state but also from themacroscopic equilibrium corresponding to the lower-activity asynchronous state. Moreover, in both instances, the Hopf bifurcations
emanating from the Bogdanov–Takens point are subcritical andchange 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,26,31,95 orstochastic stimulation.21 In this context, our results indicate that
excitability adaptation may be an important additional source ofmultistability in the collective dynamics of neuronal populations.
In relation to explaining the impact of the ketogenic diet as atreatment 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. Notethat this interpretation is different from Ref. 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 neuronalmetabolism on the emergent dynamics. In that work, the transition
to a high-activity asynchronous state was interpreted as a seizurelike 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 applicationof 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 thelatter 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 thenormal state and back. Interestingly, earlier numerical studies on
populations of Hodgkin–Huxley neurons including a similar modelof ATP-gated potassium adaptation currents have revealed only asmooth (supercritical) transition from the synchronous (seizurelike) to the asynchronous (normal) state with the increase of ATP
production rate.69–71 Nevertheless, our study has also revealed thepossibility of hysteretic transitions facilitated by the bistability ofthese 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 ofhysteretic transitions may prove advantageous for controlling the
seizure-like states, given that the induced transition to the loweractivity 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

34, 053128-10

28 June 2024 05:00:50

FIG. 8. Switch between coexisting higher- and lower-activity asynchronous states[stable foci of system (11)] and vice versa triggered by respective negative andpositive ATP shocks [see arrows in (c)]. (a)–(c) Time traces for the mean-fieldmodel 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: τ = 2.5, K = 15.0,η̄ = −2.65, 1 =  = α = C̃ = 1.0 [triangle in Fig. 2(b)]. and N = 104 .

ARTICLEChaosChaos 34, 053128 (2024); doi: 10.1063/5.0180954
Published under an exclusive license by AIP Publishingpubs.aip.org/aip/cha

elaborately designed experiments, which could then lead to further
development and refinement of the current model.
ACKNOWLEDGMENTS
We would like to thank Professor Eckehard Schöll, Dr.
Matthias Wolfrum, and Dr. Predrag Janjić 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
No. ANSO-CR-PP-2022-05.
AUTHOR DECLARATIONSConflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Sebastian Eydam: Conceptualization (equal); Formal analysis(equal); Investigation (equal); Methodology (equal); Visualization(equal); Writing – original draft (equal); Writing – review & editing(equal). Igor Franović: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review &editing (equal). Louis Kang: Conceptualization (supporting); Formal analysis (supporting); Investigation (supporting); Methodology(supporting); Visualization (supporting); Writing – review & editing(supporting).
DATA AVAILABILITY
The data that support the findings of this study are available
from the corresponding author upon reasonable request.
REFERENCES
1
E. Montbrió, D. Pazó, and A. Roxin, “Macroscopic description for networks ofspiking neurons,” Phys. Rev. X 5, 021028 (2015).
2
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).
3
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.
4
A. Byrne, R. D. O’Dea, M. Forrester, J. Ross, and S. Coombes, “Next-generationneural mass and field modeling,” J. Neurophysiol. 123, 726 (2020).
5
S. Coombes, “Next generation neural population models,” Front. Appl. Math.
Stat. 9, 1128224 (2023).
6
S. Coombes, P. beim Graben, R. Potthast, and J. Wright, Neural Fields: Theory
and Applications (Springer Berlin, Heidelberg, 2014).
7
G. Deco, V. K. Jirsa, P. A. Robinson, M. Breakspear, and K. Friston, “Thedynamic brain: From spiking neurons to neural masses and cortical fields,” PLoS
Comput. Biol. 4, e1000092 (2008).
8
N. Deschle, J. Ignacio Gossn, P. Tewarie, B. Schelter, and A. Daffertshofer, “Onthe validity of neural mass models,” Front. Comput. Neurosci. 14, 581040 (2021).
9
E. Ott and T. M. Antonsen, “Low dimensional behavior of large systems ofglobally coupled oscillators,” Chaos 18, 037113 (2008).
10
E. Ott and T. M. Antonsen, “Long time evolution of phase oscillator systems,”
Chaos 19, 023117 (2009).

34, 053128-11

28 June 2024 05:00:50

of seizures from the electroencephalography recordings. Classically,
for the continuous transitions, one expects standard indicators ofcriticality related to critical slowing down91,96,97 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 pro98 and contra,99 one may speculate that certain seizures may in fact derive from the continuous, andsome from the discontinuous hysteretic transitions, as indicated byour model.
One should point out that the current model does not provide
an exhaustive description of the effects of epilepsy and ketogenic
diet on the local neuronal kinetics. In particular, epilepsy is known
to impact the dynamics of leaky sodium, potassium, and chloride currents by triggering changes in the corresponding reversal
potentials.100–102 Also, the current model is limited by the fact thatit does not consider the reversal potential of the ATP-dependent
potassium current and the way in which it is affected by the ketogenic diet. While this is partly related to the lack of conclusive
experimental evidence, it is also associated with the fact that thesechanges likely cannot be considered isolated from the dynamics ofthe leaky currents, including their implicit dependencies on surrounding glia, whose taking into account would make the model andits subsequent analysis substantially more complex.
We believe that our results, together with Refs. 38 and 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-generationneural mass models. This is expected to provide a valuable complement to experimental studies in at least two ways. First, it willallow for a qualitative insight into multistability which may be lacking or not be easily accessible to experiments and second, 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 unfoldingin 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 forthe 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 fromthe current mean-field model, similar to those observed in Refs. 37
and 40. Finally, the immediate impact of ketone bodies on inhibition of excitatory synaptic channels could be taken into account asa factor modulating the synaptic efficacy. More generally, an interesting challenge would be to understand how the classical physicalpicture developed for activity of neuronal populations within the
NGNM framework changes by incorporating the different formsof homeostatic plasticity,103 which involve negative feedback loops
to constrain the network activity within the desired physiologicallimit. Also, as a next step, it would be important to provide atleast a qualitative validation of our theoretical predictions through

ARTICLEChaosChaos 34, 053128 (2024); doi: 10.1063/5.0180954
Published under an exclusive license by AIP Publishingpubs.aip.org/aip/cha

36
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).
37
H. Taher, A. Torcini, and S. Olmi, “Exact neural mass model for synaptic-based
working memory,” PLoS Comput. Biol. 16, e1008533 (2020).
38
M. Gerster, H. Taher, A. Škoch, J. Hlinka, M. Guye, F. Bartolomei, V. Jirsa,
A. Zakharova, and S. Olmi, “Patient-specific network connectivity combinedwith a next generation neural mass model to test clinical hypothesis of seizure
propagation,” Front. Syst. Neurosci. 15, 675272 (2021).
39
R. Gast, H. Schmidt, and T. R. Knösche, “A mean-field description of bursting dynamics in spiking neural networks with short-term adaptation,” NeuralComput. 32, 1615 (2020).
40
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).
41
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).
42
J. Sawicki, R. Berner, S. A. M. Loos, M. Anvari, R. Bader, W. Barfuss, N. Botta,
N. Brede, I. Franović, 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).
43
L. Abbott and S. Nelson, “Synaptic plasticity: Taming the beast,” Nat. Neurosci.
3, 1178 (2000).
44
N. Caporale and Y. Dan, “Spike timing-dependent plasticity: A Hebbian learning rule,” Ann. Rev. Neurosci. 31, 25 (2008).
45
G. E. Ha and E. Cheong, “Spike frequency adaptation in neurons of the central
nervous system,” Exp. Neurobiol. 26, 179 (2017).
46
G. Fuhrmann, H. Markram, and M. Tsodyks, “Spike frequency adaptation andneocortical rhythms,” J. Neurophysiol. 88, 761 (2002).
47
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).
48
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).
49
I. Franović, S. Eydam, S. Yanchuk, and R. Berner, “Collective activity burstingin a population of excitable units adaptively coupled to a pool of resources,” Front.
Netw. Physiol. 2, 841829 (2022).
50
A. Lutas and G. Yellen, “The ketogenic diet: Metabolic influences on brain
excitability and epilepsy,” Trends Neurosci. 36, 32 (2013).
51
J. W. Wheless, “History of the ketogenic diet,” Epilepsia 49, 3 (2008).
52
Y. Fei, R. Shi, Z. Song, and J. Wu, “Metabolic control of epilepsy: A promisingtherapeutic target for epilepsy,” Front. Neurol. 11, 592514 (2020).
53
J. M. Rho, “How does the ketogenic diet induce anti-seizure effects?,” Neurosci.
Lett. 637, 4 (2017).
54
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).
55
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).
56
J. A. Roberts, K. K. Iyer, S. Vanhatalo, and M. Breakspear, “Critical rolefor resource constraints in neural models,” Front. Syst. Neurosci 8, 154(2014).
57
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).
58
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).
59
I. M. Glynn, “A hundred years of sodium pumping,” Annu. Rev. Physiol. 64, 1(2002).
60
M. D. Forrest, “The sodium-potassium pump is an information processing
element in brain computation,” Front. Physiol. 5, 472 (2014).

34, 053128-12

28 June 2024 05:00:50

11T. 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).
12
F. Devalle, A. Roxin, and E. Montbrió, “Firing rate equations require a spikesynchrony mechanism to correctly describe fast oscillations in inhibitory networks,” PLoS Comput. Biol. 13, e1005881 (2017).
13
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).
14
B. Pietras, F. Devalle, A. Roxin, A. Daffertshofer, and E. Montbrió, “Exact firing
rate model reveals the differential effects of chemical versus electrical synapses inspiking networks,” Phys. Rev. E 100, 042412 (2019).
15
E. Montbrió and D. Pazó, “Exact mean-field theory explains the dual role ofelectrical synapses in collective synchronization,” Phys. Rev. Lett. 125, 248101(2020).
16
D. S. Goldobin, “Mean-field models of populations of quadratic integrate-andfire neurons with noise on the basis of the circular cumulant approach,” Chaos 31,
083112 (2021).
17
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).
18
P. Clusella and E. Montbrió, “Exact low-dimensional description for fast neuraloscillations with low firing rates,” arXiv.2208.05515 (2022).
19
B. Pietras, R. Cestnik, and A. Pikovsky, “Exact finite-dimensional description
for networks of globally coupled spiking neurons,” Phys. Rev. E 107, 024315(2023).
20
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).
21
D. S. Goldobin, M. di Volo, and A. Torcini, “Reduction methodology forfluctuation driven population dynamics,” Phys. Rev. Lett. 127, 038301 (2021).
22
V. V. Klinshov and S. Y. Kirillov, “Shot noise in next-generation neural massmodels for finite-size networks,” Phys. Rev. E 106, L062302 (2022).
23
V. V. Klinshov, P. S. Smelov, and S. Y. Kirillov, “Constructive role of shot noisein the collective dynamics of neural networks,” Chaos 33, 061101 (2023).
24
I. Ratas and K. Pyragas, “Symmetry breaking in two interacting populations ofquadratic integrate-and-fire neurons,” Phys. Rev. E 96, 042212 (2017).
25
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).
26
K. Pyragas, A. P. Fedaravičius, and T. Pyragienė, “Suppression of synchronous
spiking in two interacting populations of excitatory and inhibitory quadratic
integrate-and-fire neurons,” Phys. Rev. E 104, 014203 (2021).
27
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).
28
H. Schmidt and D. Avitabile, “Bumps and oscillons in networks of spiking
neurons,” Chaos 30, 033133 (2020).
29
A. Byrne, J. Ross, R. Nicks, and S. Coombes, “Mean-field models for EEG/MEG:
From oscillations to waves,” Brain Topogr. 35, 36 (2022).
30
G. Dumont and B. Gutkin, “Macroscopic phase resetting-curves determine
oscillatory coherence and signal transfer in inter-coupled neural circuits,” PLoS
Comput. Biol. 15, e1007019 (2019).
31
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).
32
M. Segneri, H. Bi, S. Olmi, and A. Torcini, “Theta-nested gamma oscillations innext generation neural mass models,” Front. Comput. Neurosci. 14, 47 (2020).
33
S. Keeley, A. Byrne, A. Fenton, and J. Rinzel, “Firing rate models for gamma
oscillations,” J. Neurophysiol. 121, 2181 (2019).
34
A. Ceni, S. Olmi, A. Torcini, and D. Angulo-Garcia, “Cross frequency coupling in next generation inhibitory neural mass models,” Chaos 30, 053121(2020).
35
Á. 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.

ARTICLEChaosChaos 34, 053128 (2024); doi: 10.1063/5.0180954
Published under an exclusive license by AIP Publishingpubs.aip.org/aip/cha

82

L. B. Gano, M. Patel, and J. M. Rho, “Ketogenic diets, mitochondria, andneurological diseases,” J. Lipid Res. 55, 2211 (2014).
83
C. G. J. Saris and S. Timmers, “Ketogenic diets and ketone suplementation: A
strategy for therapeutic intervention,” Frontiers Nutr. 9, 947567 (2022).
84
C. Rackauckas and Q. Nie, “DifferentialEquations.jl—A performant andfeature-rich ecosystem for solving differential equations in Julia,” J. Open Res.
Softw. 5, 15 (2017).
85
F. Devalle, E. Montbrió, and D. Pazó, “Dynamics of a large system of spiking
neurons with synaptic delay,” Phys. Rev. E 98, 042214 (2018).
86
H. Taher, D. Avitabile, and M. Desroches, “Bursting in a next generation neuralmass model with synaptic dynamics: A slow–fast approach,” Nonlinear Dyn. 108,
4261 (2022).
87
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).
88
R. Veltz, see https://hal.archives-ouvertes.fr/hal-02902346 for “Bifurcationkit.jl”(2020).
89
V. Klinshov and I. Franović, “Two scenarios for the onset and suppression ofcollective oscillations in heterogeneous populations of active rotators,” Phys. Rev.
E 100, 062211 (2019).
90
Dynamic Bifurcations, Lecture Notes in Mathematics Vol. 1493, edited by E.
Benoît (Springer, Berlin, 1991).
91
C. Kuehn, “A mathematical framework for critical transitions: Bifurcations,fast–slow systems and stochastic dynamics,” Phys. D 240, 1020 (2011).
92
R. D. Fields, “A new mechanism of nervous system plasticity: Activitydependent myelination,” Nat. Rev. Neurosci. 16, 756 (2015).
93
D. Attwell and S. B. Laughlin, “An energy budget for signaling in the grey matterof the brain,” J. Cereb. Blood Flow Metab. 21, 1133 (2001).
94
L. F. Lafuerza, P. Colet, and R. Toral, “Nonuniversal results induced by diversity distribution in coupled excitable systems,” Phys. Rev. Lett. 105, 084101(2010).
95
S. Ostojic, N. Brunel, and V. Hakim, “Synchronization properties of networksof electrically coupled neurons in the presence of noise and heterogeneities,”
J. Comput. Neurosci. 26, 369 (2009).
96
C. Kuehn, “A mathematical framework for critical transitions: Normal forms,
variance and applications,” J. Nonlinear Sci. 23, 457 (2013).
97
C. Meisel, A. Klaus, C. Kuehn, and D. Plenz, “Critical slowing down governs thetransition to neuron spiking,” PLoS Comput. Biol. 11, e1004097 (2015).
98
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).
99
T. Wilkat, T. Rings, and K. Lehnertz, “No evidence for critical slowing downprior to human epileptic seizures,” Chaos 29, 091104 (2019).
100
J. R. Cressman, G. Ullah, J. Ziburkus, S. J. Schiff, and E. Barreto, “The influence of sodium and potassium dynamics on excitability, seizures, and the stabilityof persistent states: I. Single neuron dynamics,” J. Comput. Neurosci. 26, 159(2009).
101
M. Du, J. Li, R. Wang, and Y. Wu, “The influence of potassium concentration
on epileptic seizures in a coupled neuronal model in the hippocampus,” Cogn.
Neurodyn. 10, 405 (2016).
102
G. Ullah and S. J. Schiff, “Assimilating seizure dynamics,” PLoS Comput. Biol.
6, e1000776 (2010).
103
J. Zierenberg, J. Wilting, and V. Priesemann, “Homeostatic plasticity andexternal input shape neural network dynamics,” Phys. Rev. X 8, 031018 (2018).

34, 053128-13

28 June 2024 05:00:50

61
P. Bazzigaluppi, A. E. Amini, I. Weisspapir, B. Stefanovic, and P. L. Carlen,“Hungry neurons: Metabolic insights on seizure dynamics,” Int. J. Mol. Sci. 18,2269 (2017).
62
M. Patel, “A metabolic paradigm for epilepsy,” Epilepsy Curr. 18, 318 (2018).
63
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).
64
H. Büeler, “Impaired mitochondrial dynamics and function in the pathogenesis
of Parkinson’s disease,” Exp. Neurol. 218, 235 (2009).
65
D. Haddad and K. Nakamura, “Understanding the susceptibility of dopamine
neurons to mitochondrial stressors in Parkinson’s disease,” FEBS Lett. 589, 3702
(2015).
66
E. Beghi, “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).
67
M.-C. Picot, M. Baldy-Moulinier, J.-P. Daurs, P. Dujols, and A. Crespel, “Theprevalence of epilepsy and pharmacoresistant epilepsy in adults: A populationbased study in a western European country,” Epilepsia 49, 1230 (2008).
68
E. G. Neal and J. H. Cross, “Efficacy of dietary treatments for epilepsy,” J. Hum.
Nutr. Diet. 23, 113–119 (2010).
69
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).
70
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,” Proc. Natl. Acad. Sci. U. S. A. 103, 5597 (2006).
71
S. Ching, P. L. Purdon, S. Vijayan, N. J. Kopell, and E. N. Brown, “A neurophysiological–metabolic model for burst suppression,” Proc. Natl. Acad. Sci. U. S. A.
109, 3095 (2012).
72
G. B. Ermentrout and D. H. Terman, Mathematical Foundations of Neuroscience(Springer, New York, 2010).
73
E. M. Izhikevich, Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting (The MIT Press, Cambridge, 2007).
74
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).
75
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).
76
G. Yellen, “Ketone bodies, glycolysis, and KATP channels in the mechanism of
the ketogenic diet,” Epilepsia 49, 80 (2008).
77
G. B. Ermentrout and N. Kopell, “Parabolic bursting in an excitable system
coupled with a slow oscillation,” SIAM J. Appl. Math. 46, 233 (1986).
78
B. Ermentrout, “Type I membranes, phase resetting curves, and synchrony,”
Neural Comput. 8, 979 (1996).
79
V. J. Miller, R. A. LaFountain, E. Barnhart, T. S. Sapper, J. Short, W. D. Arnold,
P. N. Hyde, C. D. Crabtree, M. L. Kackley, W. J. Kraemer, F. A. Villamena, andJ. S. Volek, “A ketogenic diet combined with exercise alters mitochondrial function in human skeletal muscle while improving metabolic health,” Am. J. Physiol.
Endocrinol. Metab. 319, E995 (2020).
80
S. A. Masino and J. M. Rho, “Mechanisms of ketogenic diet action,” in Jasper’s
Basic Mechanisms of the Epilepsies (Oxford University Press, 2012).
81
K. J. Bough, J. Wetherington, B. Hassel, J. F. Pare, J. W. Gawryluk, J. G. Greene,
R. Shaw, Y. Smith, J. D. Geiger, and R. J. Dingledine, “Mitochondrial biogenesis in
the anticonvulsant mechanism of the ketogenic diet,” Ann. Neurol. 60, 223 (2006).

ARTICLE