# Spiral Attractors in a Reduced Mean-Field Model of Neuron-Glial Interaction

S. Olenin  Control Theory Department, Lobachevsky University, 603022, Russia,
Nizhny Novgorod, Gagarin ave, 23  S. Stasenko  Laboratory of advanced methods
for high-dimensional data analysis, Lobachevsky University, 603022, Russia,
Nizhny Novgorod, Gagarin ave, 23  T. Levanova  Control Theory Department,
Lobachevsky University, 603022, Russia, Nizhny Novgorod, Gagarin ave, 23
<EMAIL_ADDRESS>
###### Abstract

This paper investigates various bifurcation scenarios of the appearance of
bursting activity in the phenomenological mean-field model of neuron-glial
interactions. In particular, we show that the homoclinic spiral attractors in
this system can be the source of several types of bursting activity with
different properties.

††preprint: AIP/123-QED

> It is well known that bursting activity plays an important role in the
> processes of transmission of neural signals. In terms of population
> dynamics, macroscopic bursting can be described using a mean-field approach.
> Mean field theory provides a useful tool for analysis of collective behavior
> of a large populations of interacting units, allowing to reduce the
> description of corresponding dynamics to just a few equations. Recently a
> new phenomenological model was proposed that describes bursting population
> activity of a big group of excitatory neurons, taking into account short-
> term synaptic plasticity and the astrocytic modulation of the synaptic
> dynamics Olenin, Levanova, and Stasenko (2023). The purpose of the present
> study is to investigate various bifurcation scenarios of the appearance of
> bursting activity in the phenomenological model. We show that the birth of
> bursting population pattern can be connected both with the cascade of period
> doubling bifurcations and further development of chaos according to the
> Shilnikov scenario, which leads to the appearance of a homoclinic attractor
> containing a homoclinic loop of a saddle-focus equilibrium with the two-
> dimensional unstable invariant manifold. We also show that the homoclinic
> spiral attractors observed in the system under study generate several types
> of bursting activity with different properties.

##  I Introduction

One of the most interesting and important patterns of population activity is
the bursting population activity Chen  _et al._ (2022), which consists of two
stages. The first one includes oscillations with high frequency, repeated at
relatively short intervals of time, and the second one is the state of
"quiescence". Bursting activity can be observed in many perceptual and
behavioral experiments and underlie both normal physiological Meister  _et
al._ (1991); Krahe and Gabbiani (2004); Salinas and Sejnowski (2001); Axmacher
_et al._ (2006) and pathophysiological processes (e.g., epilepsy Avoli and
Olivier (1987); Hofer  _et al._ (2022); Stasenko, Hramov, and Kazantsev
(2023)). Among the various informational and physiological processes related
to bursting activity the following ones can be mentioned: restoration of
synaptic transmission after its disruption Lisman (1997), neurotransmitter
release Lisman (1997), formation of long-term potentiation Lisman (1997),
selectivity of synaptic communication Izhikevich  _et al._ (2003), sensory
inputs and increasing the information capacity Gabbiani  _et al._ (1996);
Oswald  _et al._ (2004); Lesica and Stanley (2004); Reinagel  _et al._ (1999);
Segev  _et al._ (2004) and many others Hulata  _et al._ (2004); Stasenko and
Kazantsev (2023a, b).

The mechanisms underlying the appearance of bursting population activity are
not yet fully understood Gast, Schmidt, and Knösche (2020). To date, most
computational models deal with the study of single cell dynamics, see, e.g.
Izhikevich  _et al._ (2003); Guckenheimer  _et al._ (1997). To reproduce
bursting activity at the population level one can use implicit properties of
interacting excitatory and inhibitory populations Kudela, Franaszczuk, and
Bergey (2003); Zeldenrust, Wadman, and Englitz (2018) or explicitly add burst-
provoking factors (e.g. action of a neuromodulator Marder and Thirumalai
(2002), feedforward inhibition Zeldenrust and Wadman (2013), spike-frequency
adaptation Vreeswijk and Hansel (2001) or their combinations Markram and
Tsodyks (1996); Masquelier and Deco (2013); Rozhnova  _et al._ (2021);
Stasenko and Kazantsev (2023c)). Another valuable features that should be
taken into account when creating such models are astrocytic modulation of
neuronal activity and short-term synaptic plasticity Barabash, Levanova, and
Stasenko (2021, 2023). Adding these factors to the model can lead to rich and
nontrivial dynamics Stasenko, Lazarevich, and Kazantsev (2020); Stasenko and
Kazantsev (2022). In particular, in Cortes  _et al._ (2013) it was shown that
the addition of short-term synaptic plasticity to the Tsodyks-Markram model
leads to the emergence of chaotic bursting related to spiral chaos, which
occurs according to the Shilnikov scenario Gonchenko  _et al._ (2022).

Shilnikov given a phenomenological description of his scenario in the paper
Shilnikov (1984) for the case of one-parameter families of multidimensional
systems. This scenario is extremely simple (by modulo of certain intermediate
details that themselves can be very complicated), and therefore, it should
come as no surprise that it is often seen in many models. In the three-
dimensional case, the essence of the Shilnikov scenario is as follows.

At the beginning, the attractor of a system is a stable equilibrium. Then it
loses stability under a supercritical (soft) Andronov-Hopf bifurcation: the
equilibrium itself becomes a saddle-focus of type (1,2), i.e. with one-
dimensional stable and two-dimensional unstable invariant manifolds, and a
stable limit cycle is born, which becomes the attractor. With further changes
in parameter, this cycle becomes focal, i.e. it takes complex conjugate
multipliers, and the so-called Shilnikov funnel is created, into which all
trajectories from the absorbing region are drawn, with the exception of one
stable separatrix of the saddle-focus. This funnel persists even when the
limit cycle loses its stability and the attractor becomes strange. Moreover,
the Shilnikov attractor can be formed here, when a homoclinic loop of the
saddle-focus arises (one of its stable separatrices falls on the two-
dimensional unstable manifold). In this case, the famous Shilnikov criterion
of chaos Shil’nikov (1965) is certainly satisfied when the divergence at the
saddle-focus is negative.

It is worth noting that Shilnikov chaos has been observed in various systems,
both in experiments and in models. It has been experimentally detected in
various systems, including laser systems Pisarchik, Meucci, and Arecchi
(2000); Pisarchik  _et al._ (2001); Pisarchik, Meucci, and Arecchi (2001);
Arecchi, Meucci, and Gadomski (1987); Viktorov, Klemer, and Karim (1995),
optically bistable devices Farjas  _et al._ (1996), electronic circuits Kingni
_et al._ (2013), memristive circuits Njitacke  _et al._ (2016), and numerous
other systems Dangoisse  _et al._ (1988); Hennequin  _et al._ (1989); de
Tomasi  _et al._ (1989). Additionally, numerical investigations have revealed
the appearance of Shilnikov chaos in mathematical models of diverse scientific
fields, such as biology and medicine Griffith, Parthimos, and Edwards (2009);
Cortes  _et al._ (2013); Korotkov, Kazakov, and Levanova (2019); Shilnikov and
Kolomiets (2008), ecology Bakhanova  _et al._ (2018); Barrio  _et al._ (2011);
Rosenzweig (1973); Hastings and Powell (1991), chemistry Gaspard and Nicolis
(1983), economics Barnett  _et al._ (2022); Bella, Mattana, and Venturi
(2017), engineering Wei, Qi-Chang, and Rui-Lan (2010); Yao and Zhang (2005),
and many others, see, e.g., papers Borisov, Kazakov, and Sataev (2016);
Bakhanova  _et al._ (2020, 2023); Malykh  _et al._ (2020); Karatetskaia,
Shykhmamedov, and Kazakov (2021) and references therein.

One of the main features of Shilnikov chaos is that one can distinguish in it
two rather different stages: the stage during which the trajectory passes near
the equilibrium point (the stage of small amplitude), and the stage, when it
passes along a global piece of the homoclinic loop (the stage of large
amplitude). These stages are the source of specific bursting patterns (so-
called homoclinic burstings) that differ from typical patterns of bursting
activity when oscillations with high frequencies alternate with stages of
quiescence (non-homoclinic burstings). In Bakhanova  _et al._ (2018), using
the Rosenzweig-MacArthur system as an example, it was shown that homoclinic
attractors, which appear due to Shilnikov bifurcation scenario, can be a
source of such types of bursting activity.

In this paper, we focus on similar scenarios of the appearance of different
types of regular and chaotic bursting population activity in a recently
proposed mean-field model of neuron-glial interaction Olenin, Levanova, and
Stasenko (2023). This model is a simplification of the previously proposed
4-dimensional model of population activity Barabash, Levanova, and Stasenko
(2021, 2023), which, in its turn, is based on the Tsodyks-Markram mean-field
model Mongillo, Barak, and Tsodyks (2008) and takes into account the main
features of neuron-glial interaction Lazarevich, Stasenko, and Kazantsev
(2017). We show that mixed-mode oscillations arising from a Shilnikov’s
homoclinic loop, together with a slow-fast behavior, lead to the appearance of
homoclinic bursting population activity in the system.

The paper is organized as follows. In Section II we introduce and briefly
describe the reduced model of population activity. We present the results of
two-parameter bifurcation analysis in Section III. We describe the onset of
spiral chaos according to Shilnikov scenario in Section IV. In Section V we
compare different types of bursting activity based on its features and
mathematical images. In Section VI we discuss our results before draw our
conclusions in Section VII.

##  II The model

To investigate mechanisms underlying population bursting activity in the
presence of astrocytic influence, in our previous paper Olenin, Levanova, and
Stasenko (2023) we proposed a new reduced mean-field model of neuron-glial
interactions. This model is a combination of two models, the Wilson-Cowan
model Wilson and Cowan (1972) describing the dynamics of populations of
excitatory and inhibitory neurons, and the Tsodyks-Markram model Mongillo,
Barak, and Tsodyks (2008); Tsodyks, Pawelzik, and Markram (1998) of short-term
synaptic plasticity, to which we added the equation of neuron-glial
interaction, phenomenologically described using the mean- field approach
Stasenko, Lazarevich, and Kazantsev (2020); Gordleeva (2012). For simplicity,
we focus only on the population of excitatory neurons, as done in Wilson and
Cowan (1972).

The model describes the synchronous dynamics of a population of excitatory
neurons (e.g., pyramidal neurons) with a constant inhibitory input from
another population of neurons (e.g., interneurons) in the presence of
astrocytic influence. Mathematically, the reduced model can be written in the
form of the following three-dimensional system of ODEs:

$\left\\{\begin{array}[]{l}\tau\dot{E}=-E+\alpha\ln{(1+\exp{\frac{JU(y)xE+I_{0}}{\alpha}})}\\\
\dot{x}=\frac{1-x}{\tau_{D}}-U(y)xE\\\
\dot{y}=\frac{-y}{\tau_{y}}+\beta\sigma_{y}(x)\\\ \end{array}\right.$ (1)

Here the phase variables $E,x$ and $y$ are as follows: $E$ is the value of
synchronous activity of population of excitatory neurons, $x$ and $y$ are,
respectively, probabilities of neurotransmitter and gliotransmitter release.
Note that variable $E$ is not a membrane potential, but only reflects the rate
of population activity of a large population of neurons. System (1) depends
also on a number of parameters that are certain characteristics of
phenomenological functions of phase variables used for the modelling of
synaptic processes. So, the parameter $I_{0}$ is a constant that defines the
inhibitory input from the population of interneurons applied to the population
of excitatory neurons. The parameter $\alpha$ determines the threshold for
increasing the synchronous activity of the population of excitatory neurons.
Here the parameters $\tau$, $\tau_{D}$ and $\tau_{y}$ are characteristic time
constants for population activity, synaptic depression and gliotransmitter
relaxation, respectively.

Function $U(y)$ determines the probability of neurotransmitter release, while
the combination $JU(y)xE$ introduces a positive feedback:

$U(y)=U_{0}+\frac{\Delta U_{0}}{1+e^{-50(y-y_{thr})}},$ (2)

where the parameter $U_{0}$ determines an averaged probability of
neurotransmitter release, the parameter $\Delta U_{0}$ corresponds to the
change in $U_{0}$ induced by the astrocytic influence, the parameter $y_{thr}$
is the activation threshold for the gliotransmitter that sets a level of
gliotransmission required to manifest an effect of neurotransmitter release.

Function $\sigma_{y}(x)$ is an activation function for astrocyte that depends
on the state of neuron:

$\sigma_{y}(x)=\frac{1}{1+e^{-20(x-x_{thr})}},$ (3)

with activation threshold $x_{thr}$ determines the level of neural activity at
which the gliotransmitter is released by astrocyte. In system (1) parameter
$\beta$ preceding the function $\sigma_{y}(x)$ sets the released part of the
gliotransmitter.

The interaction between neurons and astrocytes can be simplistically described
as follows. Being activated, excitatory neurons release the neurotransmitter
$x$, which, in turn, activates astrocytes via a cascade of biochemical
reactions. As a result, a gliotransmitter $y$ is released. The
gliotransmitter, in its turn, can change the probability of neurotransmitter
release. Thus, a feedback loop is formed Araque  _et al._ (1999). This
phenomenon is depicted in the model in such a way that the function $U$
(activation for neuron) depends on $y$, and the function $\sigma$ (activation
for astrocyte) depends on $x$.

It should be noted that the proposed model does not incorporate the mechanism
of synaptic facilitation. The model includes short-term synaptic plasticity
only in the form of synaptic depression. The regulation of different patterns
of population activity in the model is solely controlled by the dynamics of
the astrocytes through astrocytic facilitation of synaptic transmission. Our
model has all the limitations inherent in a Wilson-Cowan type model.

We choose $I_{0}$ and $U_{0}$ as control parameters of the model as it was
done in Olenin, Levanova, and Stasenko (2023). All other parameters were fixed
as follows: $\Delta U_{0}=0.305$, $\tau_{D}=0.07993$, $\tau_{y}=3.3$,
$\tau=0.013$, $\beta=0.3$, $x_{thr}=0.75$, $y_{thr}=0.4$, $J=3.07$,
$\alpha=1.58$.

##  III Two-parameter bifurcation analysis

In our recent paper Olenin, Levanova, and Stasenko (2023) it was shown that
the dynamics of system (1) can be chaotic. In this paper we go beyond the
above-mentioned studies and explore peculiarities of chaotic bursting
population activity in system (1).

Figure 1: Chart of MLE $\lambda$ on the parameter plane ($I_{0}$, $U_{0}$) and
its enlarged fragments. Red marker corresponds to the regions of chaotic
dynamics ($\lambda>0$), blue marker – to the regions of quiescence
($\lambda<0$), white marker – to the regions with periodic dynamics
($\lambda=0$). Here $ah$ label marks Andronov-Hopf bifurcation, $f_{1}$ and
$f_{2}$ – fold bifurcations, $h_{1}$ and $h_{2}$ – homoclinic bifurcations.
Point $BT$ marks the Bogdanov-Takens bifurcation, point $GH$ – generalized
Hopf bifurcation, $cusp$ – cusp point. Letters $A$ and $B$ in right lower
panel and $C$ and $D$ in left lower panel mark combinations of parameters
$I_{0}$ and $U_{0}$ that correspond to different patterns of bursting
population activity. See detailed description in the text.

In order to do this, we conducted a two-parameter bifurcation analysis of
system (1). First of all, we constructed the chart of maximal Lyapunov
exponent (MLE) $\lambda$ on the plane of control parameters $I_{0}$ and
$U_{0}$, see fig. 1. This chart consists of different colored regions whose
colors mean the following. The red colored regions correspond to those
parameter values at which system (1) has chaotic dynamics (strange attractor,
$\lambda>0$). Regions with shades of blue correspond to those parameter values
at which system (1) has a stable equilibrium (a "quiescence" state,
$\lambda<0$). Finally, white regions correspond to those parameter values at
which system (1) has a stable limit cycle (a periodic oscillations,
$\lambda=0$ and other Lyapunov exponents are negative).

As can be seen from fig. 1, chaotic regions include stability windows, i.e.
regions in the parameter space which correspond to the emergence of stable
periodic orbits (here, multi-turn limit cycles). Such windows of stability are
typical for chaotic non-hyperbolic systems and have universal form called
’saddle area’ (or ’squid’), see Gonchenko, Li, and Turaev (2022) for details.

Also, we performed a numerical bifurcation analysis of system (1). Namely,
using Bifurcationkit package Veltz (2020) and Matcont software Dhooge,
Govaerts, and Kuznetsov (2003), we plotted in the MLE chart bifurcation curves
both for equilibria of system (1) and for homolinic loops to these equilibria
of saddle and saddle-focus type.

Three of these curves $ah$, $f_{1}$ and $f_{2}$ relate to bifurcations of
equilibria.

Curves $f_{1}$ and $f_{2}$ relate to the fold bifurcations. Note that in the
region between the curves $f_{1}$ and $f_{2}$ (let us call this region
$Q_{3}$), which converge at the $cusp$ point, there are three equilibria
$O_{1},O_{2}$ and $O_{3}$, and outside this region there is only one. The
equilibria $O_{2}$ and $O_{3}$ merge and disappear as a result of fold
bifurcation when passing through the curve $f_{1}$. Here this bifurcation is
always of a saddle-node type: upon entering the region $Q_{3}$, a stable
equilibrium $O_{3}$ and a saddle $O_{2}$ of type (2,1) are born. On the curve
$f_{2}$, where the equilibria $O_{1}$ and $O_{2}$ merge, there is a
bifurcation point $BT$ of codimension 2, the so-called Bogdanov-Takens point,
that divides the curve $f_{2}$ into two segments, from the point $BT$ to the
$cusp$ point, where the equilibrium $O_{12}$ is a saddle-node, and from the
point $BT$ to the left, where $O_{12}$ is a saddle-saddle.

The red curve $ah$ corresponds to an Andronov-Hopf bifurcation of the
equilibrium $O_{1}$. On this curve there is the bifurcation point $GH$ of
codimension 2, the so-called generalized Hopf (or Bautin Kuznetsov, Kuznetsov,
and Kuznetsov (1998)) bifurcation point, which divides the curve $ah$ into two
segments, above the point $GH$, where the bifurcation is supercritical, and
below the point $GH$, where the bifurcation is subcritical. During the
supercritical Andronov-Hopf bifurcation, the stable point $O_{1}$ becomes a
saddle-focus (1,2) and a stable limit cycle is born. During the subcritical
Andronov-Hopf bifurcation, the stable point $O_{1}$ is surrounded by a saddle
limit cycle, which disappears as a result of sticking into the point $O_{1}$,
which becomes a saddle-focus (1,2).

It is worth to note some features of dynamics of system (1) for parameter
values from the region $Q_{3}$. Firstly, there is a phenomenon of
multistability for some parts of this region. Thus, near the $cusp$ point the
equilibria $O_{1}$ and $O_{3}$ are attractors; in some subregion above the
curve $ah$, attractors are $O_{3}$ and the limit cycle $C_{1}$ that was born
from $O_{1}$. Also here one can, in principle, observe such types of
multistability when the attractor $O_{3}$ coexists with attractors generated
by bifurcations of the cycle $C_{1}$ (including chaotic ones after an infinite
cascade of bifurcations). In addition, it should also be noted that when
leaving the region $Q_{3}$ through a part of the curve $f_{1}$ above the $FP$
point (Feigenbaum point), one can observe the scenario of transition to chaos
via intermittency Pomeau and Manneville (1980).

There are also several bifurcation points of codimension two presented in fig.
1. The point $GH$ on the curve $ah$ marks a generalized Hopf bifurcation. The
higher branch of $ah$ curve corresponds to supercritical Andronov-Hopf
bifurcation, while the lower branch – to subcritical Andronov-Hopf
bifurcation. At the bifurcation point $GH$ the equilibrium $O_{1}$ has a pair
of purely imaginary eigenvalues and a zero first Lyapunov coefficient. Point
$BT$ on the curve $f_{2}$ marks the Bogdanov-Takens bifurcation. At $BT$
bifurcation point, three bifurcation curves intersect: the Andronov-Hopf
bifurcation curve $ah$, the fold bifurcation curve $f_{2}$ and the homoclinic
bifurcation curve $h_{2}$. At the bifurcation point $BT$ the equilibrium
$O_{1}$ has the double zero eigenvalue. In the $cusp$ bifurcation point, where
equilibrium state has zero eigenvalue and the quadratic coefficient for the
saddle-node bifurcation vanishes, two fold bifurcation curves $f_{1}$ and
$f_{2}$ converge.

(a)  (b)

Figure 2: Homoclinic loop to (a) saddle-focus $O_{1}$ for
$I_{0}\approx-1.7187$, $U_{0}\approx 0.2659$, (b) saddle $O_{2}$ for
$I_{0}\approx-1.3032$, $U_{0}\approx 0.1799$.

Of particular interest for further research are two green curves $h_{1}$ and
$h_{2}$, which mark homoclinic bifurcations. On the curve $h_{1}$ a homoclinic
loop to saddle-focus $O_{1}$ appears (see fig. 2(a)), while on the curve
$h_{2}$ homoclinic loop to a saddle $O_{2}$ arises (see fig. 2(b)).

Now let us study in detail the parameter region near $h_{1}$ homoclinic curve,
since it can be a source of spiral homoclinic attractors born due to Shilnikov
scenario and, therefore, a source of different patterns of bursting activity
associated with it.

##  IV Shilnikov scenario and homoclinic attractors

(a) (b) (c)
(d) (e) (f)

Figure 3: Shilnikov’s scenario in system (1): (a) stable equilibrium $O_{1}$,
$I_{0}=-0.9$; (b) stable limit cycle, $I_{0}=-1.6$; (c) two-period stable
limit cycle, $I_{0}=-1.66$; (d) 4-period stable limit cycle, $I_{0}=-1.695$;
(e) attractor after a cascade of period doubling bifurcations,
$I_{0}=-1.6972$; and (f) a homoclinic spiral attractor, $I_{0}=-1.706$, the
pink curve corresponds to the homoclinic trajectory, the dark grey curves – to
the trajectories on the attractor. In all cases $U_{0}=0.265$. In cases
(b)-(f) saddle-focus $SF(1,2)$ is depicted near attractor.

Let us describe the main stages of Shilnikov’s scenario that is observed in
system (1) for fixed $U_{0}=0.265$ and decreasing $I_{0}$, see fig. 3. For
$I_{0}>I_{ah}\approx-1.1065337$, there is an asymptotically stable equilibrium
state $O_{1}$ (fig. 3(a)). For $I_{0}=I_{ah}$ it undergoes a supercritical
Andronov-Hopf bifurcation, as a result of which the equilibrium $O_{1}$
becomes a saddle-focus $(1,2)$ and a stable limit cycle $C_{1}$ is born (fig.
3(b)). Then, for further decrease in the value of $I_{0}$ this limit cycle
becomes focal and the two-dimensional unstable manifold $W^{u}(O_{1})$ begins
to wrap around it, forming a funnel-like configuration. The boundary of this
funnel consists of a saddle-focus equilibrium and its two-dimensional unstable
invariant manifold, which has the shape of a bowl, the edges of which are bent
inside the funnel. Besides, all trajectories from the absorbing region are
drawn into the funnel, except for one stable separatrix of $O_{1}$ that tends
to $O_{1}$ from the outer side of the funnel. With further varying $I_{0}$,
the limit cycle $C_{1}$ undergoes a cascade of period-doubling bifurcations
(fig. 3(c)-(d)), while the size of the funnel grows. As a result, a strange
attractor of Feigenbaum type appears (fig. 3(e)). Finally, for
$I_{0}\approx-1.7062$ a homoclinic loop to saddle-focus $O_{1}$ appears. In
fig. 3(f) one can see a strange attractor containing this loop. Strange
attractors of such type, containing either an equilibrium state in the case of
flows or a fixed point in the case of mappings, are generally called
homoclinic attractors Gonchenko, Gonchenko, and Shilnikov (2012); Gonchenko
_et al._ (2014). In this particular case, the attractor has name the Shilnikov
attractor Gonchenko  _et al._ (2019).

The structure of the set of trajectories of spiral attractors can be quite
diverse. The basic elements of the geometry of such attractors form mainly
those trajectories that pass near a global piece of the unstable manifold
$W^{u}(O1)$, which, in turn, is the boundary of the Shilnikov funnel. These
trajectories can make one, two, three or more turns before returning to the
vicinity of the saddle-focus. In papers devoted to the study of the topology
of attractors, e.g., in Rössler (1979); Argoul, Arneodo, and Richetti (1987);
Letellier, Dutertre, and Maheu (1995), such attractors are often called screw,
funnel and multifunnel attractors. The type of spiral attractor (connected
with the a number of turns) affects the number of oscillations in the
additional high-amplitude stage of bursting Bakhanova  _et al._ (2018).

Note that spiral attractors with different topology appear in models of single
neurons Hindmarsh and Rose (1984); Shilnikov and Kolomiets (2008), as well as
in the classical mean-filed models Cortes  _et al._ (2013). Moreover, in all
these cases the onset of chaos occurs due to the Shilnikov scenario.

##  V Different types of bursting activity in a reduced model of neuron-glial
interaction

System (1) belongs to the class of slow-fast systems ($E$ is the fast
variable, while $x$ and $y$ are slow ones), which also have spiral attractors.
It is therefore not surprising that system (1) can exhibit various types of
bursting population activity. Since spiral attractors belong to the class of
quasiattractors, they either contain stable periodic orbits (as a rule, of big
periods and with thin absorbing domains) or such orbits arise under
arbitrarily small perturbations. These periodic orbits repeat a global
behavior of homoclinic loops in the attractors, so one can often observe near
a chaotic spiral attractor its "a regular version" (a multi-round stable limit
cycle existing for parameter values from the corresponding windows of
stability). In contrast to spiral attractors, that exist on certain
bifurcation line, stable limit cycles persist under small perturbations and
exist for some open regions of the parameter space.

(a)
(b)

Figure 4: Chaotic non-homoclinic bursting activity (point A, $I_{0}=-1.723$,
$U_{0}=0.282$). (a) The phase portrait of the attractor. (b) Time series
$E(t)$.

According to this remark, we can now divide all observed patterns of bursting
population activity associated with spiral homoclinic attractors in system (1)
into four types, depending on which attractor corresponds to them. Below we
will give examples of such patterns, which are observed, in particular, at the
values of the control parameters $I_{0}$ and $U_{0}$ corresponding to the
selected points $A,B,C$ and $D$ on the MLE chart in fig. 1. Here we follow the
classification introduced in Bakhanova  _et al._ (2018).

###  V.1 Chaotic non-homoclinic bursting activity (point A)

Such type of bursting is generated by non-homoclinic multi-funnel attractors,
which appear in chaotic regions far from the homoclinic bifurcation curve
$h_{1}$. Examples of phase portrait of this attractor can be seen in fig.
4(a), while time series $E(t)$ are presented in fig. 4(b). The main feature of
this type of bursting population activity is a random number of fast
oscillations (stage 1, red marker) that alternates with motions along the
stable slow manifolds (stage 2, blue marker) of system (1).

###  V.2 Regular bursting activity (point B)

(a)
(b)

Figure 5: Regular non-homoclinic bursting activity (point B, $I_{0}=-1.728$,
$U_{0}=0.282$). (a) The phase portrait of the attractor. (b) Time series
$E(t)$.

This type of bursting population activity is generated by the stable limit
cycles, which appear in the stability windows of non-homoclinic multi-funnel
attractors, see fig. 5. Such limit cycles repeat the behavior of trajectories
on the corresponding attractors, and have the same number of turns. Comparing
to the previous case (point A), the length of both stages is fixed.

###  V.3 Chaotic homoclinic bursting activity (point C)

(a)
(b)

Figure 6: Chaotic homoclinic bursting activity (point C,
$I_{0}\approx-1.7071$, $U_{0}\approx 0.26506$). (a) The phase portrait of the
attractor. (b) Time series $E(t)$.

This type of population bursting activity is associated with the co-existence
of slow-fast behavior in system (1) with a drift near a homoclinic orbit of a
multi-funnel attractor. It appears in the chaotic regions along the homoclinic
bifurcation curve $h_{1}$. As one can see in fig. 6, a random number of fast
oscillations (stage 1, red marker) alternates with two types of slow motions:
previously observed motions along the stable slow manifold of the fast
subsystem (stage 2, blue marker) and small amplitude oscillations near a
saddle-focus equilibrium (stage 3, green marker). Slow motions associated with
trajectories passing near the saddle-focus $O_{1}$ appear when phase point
returns to the neighbourhood of the saddle-focus after producing fast
oscillations. The length of stage 3 depends on how close a trajectory passes
near the saddle-focus $O_{1}$. After the slow motions at stage 3 phase
trajectory oscillates near the saddle-focus (stage 4, yellow marker) and then
goes to the stable slow manifold (stage 2, blue marker).

###  V.4 Regular near homoclinic bursting activity (point D)

(a)
(b)

Figure 7: Regular near homoclinic bursting activity (point C, $I_{0}=-1.714$,
$U_{0}=0.2657$). (a) The phase portrait of the attractor. (b) Time series
$E(t)$.

This type of population bursting activity is associated with the multi-round
limit cycles which appear close enough to the saddle-focus $O_{1}$, see 7(a).
The motions along the corresponding part of these cycles are slow, see 7(b)
(stage 3). Fast oscillations (stage 1) alternate with regular slow motions
near the saddle-focus $O_{1}$ (stage 3) and with regular motions along the
stable slow manifold of the fast subsystem (stage 2). The transition between
stage 3 and stage 2 corresponds to the passing from a neighbourhood of the
saddle-focus $O_{1}$ to the slow stable manifold. The lengths of all stages
are fixed.

(a)
(b)  (c)
(d)  (e)

Figure 8: Examples of regular near homoclinic bursting activity. (a) Part of
MLE chart in the vicinity of the intersection of $f_{1}$ fold bifurcation
curve and $h_{1}$ homoclinic bifurcation curve. Points G and K mark
combinations of governing parameters $I_{0}$ and $U_{0}$ that correspond to
regular population bursting patterns with different length of quiescence
stages. (b) The phase portrait and (c) time series $E(t)$ for
$I_{0}\approx-1.7088$, $U_{0}\approx 0.2651$ (point K). (d) The phase portrait
and (e) time series $E(t)$ for $I_{0}\approx-1.7236$, $U_{0}\approx 0.2663$
(point G).

As one can see in fig. 8, the interburst interval increases while travelling
along the homoclinic curve $h_{1}$ from right to left. At the left borderline
of white region (fig. 8(a)) saddle-node bifurcation on the homoclinic curve
takes place. Before it the called bottleneck occur – a region in the phase
space of system (1) where the phase point slows down along the trajectory,
which is a sign of upcoming bifurcation. A part of attractor (stable slow
manifold, stage 2) stretches towards a future equilibrium. A similar mechanism
for increasing the interburst interval is observed outside the homoclinic
bifurcation curve $h_{1}$ near the saddle-node bifurcation on the circle on
the left border of the white region.

##  VI Discussion

An important challenge in neuroscience is to understand how complex patterns
of neural activity and brain functions are shaped by the collective dynamics
of big populations of neurons. Most of the research on this problem has
focused on large-scale numerical simulations Brette  _et al._ (2007);
Izhikevich and Edelman (2008); Izhikevich (2007). An alternative is mean field
theory, which provides insight into the macroscopic states of large neural
networks in terms of collective neuronal activity or firing rates Montbrió,
Pazó, and Roxin (2015). This kind of analysis allows one to draw conclusions
about a set of key macroscopic parameters for a neuronal population. Knowledge
of a set of such most important averaged characteristics of ensembles for
describing their dynamics allows one to formulate phenomenological self-
consistent models based on observations.

Our mean-field model of neuron-glial interactions is based on the well-known
Wilson-Cowan model, which considers integral neuronal dynamics without
reference to individual spikes. A distinctive feature of our model is the
modeling of three components: the averaged population activity of neurons
($E$), the concentration of neuro- and gliotransmitters ($x$ and $y$,
correspondingly), taking into account synaptic plasticity in the form of
synaptic depression and astrocytic potentiation. The model allows one to study
patterns of synchronized activity in large neuronal populations regardless of
the specific number of neurons, local dynamics of neurons and glial cells, or
the precise synaptic weights. In this approach all properties are averaged
across the distribution of possible weights, and in the limit of infinite
network size. All observed effects of neuron-like dynamics and neuron-glial
interaction are general and are determined only by the presence of a feedback
loop between the presynaptic neuron and the glial cell.

The proposed model allows us to reproduce and analyze various temporal
patterns observed at the population level in big neuronal ensembles in
biological experiments Lebedeva  _et al._ (2023); Wagenaar, Nadasdy, and
Potter (2006). In particular, in Wagenaar, Nadasdy, and Potter (2006) it was
shown that neural activity observed in neuronal cultures using multielectrode
recordings may include bursting and superbursting activity (microscopic
level), which in our mean-field model correspond to oscillatory and bursting
activity (macroscopic level). In this case, chaotic dynamics can be considered
as one of the mechanisms of synchronization and switching between states of
neural activity Rabinovich and Abarbanel (1998).

##  VII Conclusion

In this paper we have studied scenarios for the birth of spiral homoclinic
attractors in system (1) that is a mean-field model of neuron-glial
interaction. This model was proposed in Barabash, Levanova, and Stasenko
(2023) and it describes the synchronous neuronal activity of a population of
excitatory neurons under astrocytic modulation. We have shown that spiral
attractors in system (1) arise according to the Shilnikov scenario, that leads
to the appearance of a homoclinic attractor containing a homoclinic loop to
the saddle-focus equilibrium with a two-dimensional unstable manifold. As a
result, the system under study can generate several types of bursting
population activity with different properties.

The obtained results can help to gain new insights into the nature of some
specific patterns of activity that may arise in models of neuron-glial
interactions. We believe that the obtained results are important not only from
the point of view of nonlinear dynamics, but also from the point of view of
neurodynamics and biology.

###### Acknowledgements.

We are grateful to E.A. Grines, H.G.E. Meijer, R. Veltz, C. Rackauckas and G.
Datseris for their helpful advice with numerical calculations. We also thank
S.V. Gonchenko for careful reading. This study was financially supported by
the Russian Science Foundation grant no. 19-72-10128.

## Data Availability Statement

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

## References

  * Olenin, Levanova, and Stasenko (2023) S. M. Olenin, T. A. Levanova, and S. V. Stasenko, “Dynamics in the reduced mean-field model of neuron–glial interaction,” Mathematics 11 (2023).
  * Chen  _et al._ (2022) Y. Chen, H. Douglas, B. J. Medina, M. Olarinre, J. H. Siegle, and R. E. Kass, “Population burst propagation across interacting areas of the brain,” Journal of Neurophysiology 128, 1578–1592 (2022).
  * Meister  _et al._ (1991) M. Meister, R. O. Wong, D. A. Baylor, and C. J. Shatz, “Synchronous bursts of action potentials in ganglion cells of the developing mammalian retina,” Science 252, 939–943 (1991).
  * Krahe and Gabbiani (2004) R. Krahe and F. Gabbiani, “Burst firing in sensory systems,” Nature Reviews Neuroscience 5, 13–23 (2004).
  * Salinas and Sejnowski (2001) E. Salinas and T. J. Sejnowski, “Correlated neuronal activity and the flow of neural information,” Nature reviews neuroscience 2, 539–550 (2001).
  * Axmacher  _et al._ (2006) N. Axmacher, F. Mormann, G. Fernández, C. E. Elger, and J. Fell, “Memory formation by neuronal synchronization,” Brain research reviews 52, 170–182 (2006).
  * Avoli and Olivier (1987) M. Avoli and A. Olivier, “Bursting in human epileptogenic neocortex is depressed by an n-methyl-d-aspartate antagonist,” Neuroscience letters 76, 249–254 (1987).
  * Hofer  _et al._ (2022) K. T. Hofer, Á. Kandrács, K. Tóth, B. Hajnal, V. Bokodi, E. Z. Tóth, L. Erőss, L. Entz, A. G. Bagó, D. Fabó, _et al._ , “Bursting of excitatory cells is linked to interictal epileptic discharge generation in humans,” Scientific Reports 12, 6280 (2022).
  * Stasenko, Hramov, and Kazantsev (2023) S. V. Stasenko, A. E. Hramov, and V. B. Kazantsev, “Loss of neuron network coherence induced by virus-infected astrocytes: A model study,” Scientific Reports 13, 6401 (2023).
  * Lisman (1997) J. E. Lisman, “Bursts as a unit of neural information: making unreliable synapses reliable,” Trends in neurosciences 20, 38–43 (1997).
  * Izhikevich  _et al._ (2003) E. M. Izhikevich, N. S. Desai, E. C. Walcott, and F. C. Hoppensteadt, “Bursts as a unit of neural information: selective communication via resonance,” Trends in neurosciences 26, 161–167 (2003).
  * Gabbiani  _et al._ (1996) F. Gabbiani, W. Metzner, R. Wessel, and C. Koch, “From stimulus encoding to feature extraction in weakly electric fish,” Nature 384, 564–567 (1996).
  * Oswald  _et al._ (2004) A.-M. M. Oswald, M. J. Chacron, B. Doiron, J. Bastian, and L. Maler, “Parallel processing of sensory input by bursts and isolated spikes,” Journal of Neuroscience 24, 4351–4362 (2004).
  * Lesica and Stanley (2004) N. A. Lesica and G. B. Stanley, “Encoding of natural scene movies by tonic and burst spikes in the lateral geniculate nucleus,” Journal of Neuroscience 24, 10731–10740 (2004).
  * Reinagel  _et al._ (1999) P. Reinagel, D. Godwin, S. M. Sherman, and C. Koch, “Encoding of visual information by lgn bursts,” Journal of neurophysiology 81, 2558–2569 (1999).
  * Segev  _et al._ (2004) R. Segev, I. Baruchi, E. Hulata, and E. Ben-Jacob, “Hidden neuronal correlations in cultured networks,” Physical Review Letters 92, 118102 (2004).
  * Hulata  _et al._ (2004) E. Hulata, I. Baruchi, R. Segev, Y. Shapira, and E. Ben-Jacob, “Self-regulated complexity in cultured neuronal networks,” Physical review letters 92, 198105 (2004).
  * Stasenko and Kazantsev (2023a) S. V. Stasenko and V. B. Kazantsev, “Information encoding in bursting spiking neural network modulated by astrocytes,” Entropy 25, 745 (2023a).
  * Stasenko and Kazantsev (2023b) S. V. Stasenko and V. B. Kazantsev, “Dynamic image representation in a spiking neural network supplied by astrocytes,” Mathematics 11, 561 (2023b).
  * 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 computation 32, 1615–1634 (2020).
  * Guckenheimer  _et al._ (1997) J. Guckenheimer, R. Harris-Warrick, J. Peck, and A. Willms, “Bifurcation, bursting, and spike frequency adaptation,” Journal of computational neuroscience 4, 257–277 (1997).
  * Kudela, Franaszczuk, and Bergey (2003) P. Kudela, P. J. Franaszczuk, and G. K. Bergey, “Changing excitation and inhibition in simulated neural networks: effects on induced bursting behavior,” Biological cybernetics 88, 276–285 (2003).
  * Zeldenrust, Wadman, and Englitz (2018) F. Zeldenrust, W. J. Wadman, and B. Englitz, “Neural coding with bursts—current state and future perspectives,” Frontiers in computational neuroscience 12, 48 (2018).
  * Marder and Thirumalai (2002) E. Marder and V. Thirumalai, “Cellular, synaptic and network effects of neuromodulation,” Neural Networks 15, 479–493 (2002).
  * Zeldenrust and Wadman (2013) F. Zeldenrust and W. J. Wadman, “Modulation of spike and burst rate in a minimal neuronal circuit with feed-forward inhibition,” Neural Networks 40, 1–17 (2013).
  * Vreeswijk and Hansel (2001) C. v. Vreeswijk and D. Hansel, “Patterns of synchrony in neural networks with spike adaptation,” Neural computation 13, 959–992 (2001).
  * Markram and Tsodyks (1996) H. Markram and M. Tsodyks, “Redistribution of synaptic efficacy between neocortical pyramidal neurons,” Nature 382 (1996).
  * Masquelier and Deco (2013) T. Masquelier and G. Deco, “Network bursting dynamics in excitatory cortical neuron cultures results from the combination of different adaptive mechanism,” PloS one 8, e75824 (2013).
  * Rozhnova  _et al._ (2021) M. A. Rozhnova, E. V. Pankratova, S. V. Stasenko, and V. B. Kazantsev, “Bifurcation analysis of multistability and oscillation emergence in a model of brain extracellular matrix,” Chaos, Solitons & Fractals 151, 111253 (2021).
  * Stasenko and Kazantsev (2023c) S. V. Stasenko and V. B. Kazantsev, “Bursting dynamics of spiking neural network induced by active extracellular medium,” Mathematics 11, 2109 (2023c).
  * Barabash, Levanova, and Stasenko (2021) N. Barabash, T. Levanova, and S. Stasenko, “Stsp model with neuron-glial interaction produced bursting activity,” in  _2021 Third International Conference Neurotechnologies And Neurointerfaces (CNN)_ (IEEE, 2021) pp. 12–15.
  * Barabash, Levanova, and Stasenko (2023) N. Barabash, T. Levanova, and S. Stasenko, “Rhythmogenesis in the mean field model of the neuron–glial network,” The European Physical Journal Special Topics , 1–6 (2023).
  * Stasenko, Lazarevich, and Kazantsev (2020) S. V. Stasenko, I. A. Lazarevich, and V. B. Kazantsev, “Quasi-synchronous neuronal activity of the network induced by astrocytes,” Procedia Computer Science 169, 704–709 (2020).
  * Stasenko and Kazantsev (2022) S. Stasenko and V. Kazantsev, “3d model of bursting activity generation,” in  _2022 Fourth International Conference Neurotechnologies and Neurointerfaces (CNN)_ (IEEE, 2022) pp. 176–179.
  * Cortes  _et al._ (2013) J. M. Cortes, M. Desroches, S. Rodrigues, R. Veltz, M. A. Muñoz, and T. J. Sejnowski, “Short-term synaptic plasticity in the deterministic Tsodyks–Markram model leads to unpredictable network dynamics,” Proceedings of the National Academy of Sciences 110, 16610–16615 (2013).
  * Gonchenko  _et al._ (2022) S. Gonchenko, A. Kazakov, D. Turaev, and A. L. Shilnikov, “Leonid shilnikov and mathematical theory of dynamical chaos,” Chaos: An Interdisciplinary Journal of Nonlinear Science 32 (2022).
  * Shilnikov (1984) L. P. Shilnikov, “Bifurcation theory and turbulence,” Nonlinear and Turbulent Processes in Physics , 1627 (1984).
  * Shil’nikov (1965) L. Shil’nikov, “A case of the existence of a countable number of periodic motions,” in  _Sov. Math. Doklady_ , Vol. 6 (1965) pp. 163–166.
  * Pisarchik, Meucci, and Arecchi (2000) A. Pisarchik, R. Meucci, and F. Arecchi, “Discrete homoclinic orbits in a laser with feedback,” Physical Review E 62, 8823 (2000).
  * Pisarchik  _et al._ (2001) A. Pisarchik, F. Arecchi, R. Meucci, and A. DiGarbo, “Synchronization of shilnikov chaos in a co~ 2 laser with feedback,” LASER PHYSICS-LAWRENCE- 11, 1235–1239 (2001).
  * Pisarchik, Meucci, and Arecchi (2001) A. Pisarchik, R. Meucci, and F. Arecchi, “Theoretical and experimental study of discrete behavior of shilnikov chaos in a co 2 laser,” The European Physical Journal D-Atomic, Molecular, Optical and Plasma Physics 13, 385–391 (2001).
  * Arecchi, Meucci, and Gadomski (1987) F. Arecchi, R. Meucci, and W. Gadomski, “Laser dynamics with competing instabilities,” Physical review letters 58, 2205 (1987).
  * Viktorov, Klemer, and Karim (1995) E. A. Viktorov, D. R. Klemer, and M. A. Karim, “Shil’nikov case of antiphase dynamics in a multimode laser,” Optics communications 113, 441–448 (1995).
  * Farjas  _et al._ (1996) J. Farjas, J. Rosell, R. Herrero, R. Pons, F. Pi, and G. Orriols, “Equivalent low-order model for a nonlinear diffusion equation,” Physica D: Nonlinear Phenomena 95, 107–127 (1996).
  * Kingni  _et al._ (2013) S. T. Kingni, L. Keuninckx, P. Woafo, G. Van der Sande, and J. Danckaert, “Dissipative chaos, shilnikov chaos and bursting oscillations in a three-dimensional autonomous system: theory and electronic implementation,” Nonlinear Dynamics 73, 1111–1123 (2013).
  * Njitacke  _et al._ (2016) Z. Njitacke, H. Fotsin, A. N. Negou, D. Tchiotsop, _et al._ , “Coexistence of multiple attractors and crisis route to chaos in a novel memristive diode bidge-based jerk circuit,” Chaos, Solitons & Fractals 91, 180–197 (2016).
  * Dangoisse  _et al._ (1988) D. Dangoisse, A. Bekkali, F. Papoff, and P. Glorieux, “Shilnikov dynamics in a passive q-switching laser,” Europhysics Letters 6, 335 (1988).
  * Hennequin  _et al._ (1989) D. Hennequin, M. Lefranc, A. Bekkali, D. Dangoisse, and P. Glorieux, “Characterization of shilnikov chaos in a co 2 laser containing a saturable absorber,” Measures of Complexity and Chaos , 299–302 (1989).
  * de Tomasi  _et al._ (1989) F. de Tomasi, D. Hennequin, B. Zambon, and E. Arimondo, “Instabilities and chaos in an infrared laser with saturable absorber: experiments and vibrorotational model,” JOSA B 6, 45–57 (1989).
  * Griffith, Parthimos, and Edwards (2009) T. M. Griffith, D. Parthimos, and D. H. Edwards, “Vasomotion: the case for chaos,” Journal of biorheology 23, 11–23 (2009).
  * Korotkov, Kazakov, and Levanova (2019) A. G. Korotkov, A. O. Kazakov, and T. A. Levanova, “Effects of memristor-based coupling in the ensemble of fitzhugh–nagumo elements,” The European Physical Journal Special Topics 228, 2325–2337 (2019).
  * Shilnikov and Kolomiets (2008) A. Shilnikov and M. Kolomiets, “Methods of the qualitative theory for the Hindmarsh–Rose model: A case study–a tutorial,” International Journal of Bifurcation and chaos 18, 2141–2168 (2008).
  * Bakhanova  _et al._ (2018) Y. V. Bakhanova, A. O. Kazakov, A. G. Korotkov, T. A. Levanova, and G. V. Osipov, “Spiral attractors as the root of a new type of “bursting activity” in the Rosenzweig–MacArthur model,” The European Physical Journal Special Topics 227, 959–970 (2018).
  * Barrio  _et al._ (2011) R. Barrio, F. Blesa, S. Serrano, and A. Shilnikov, “Global organization of spiral structures in biparameter space of dissipative systems with shilnikov saddle-foci,” Physical Review E 84, 035201 (2011).
  * Rosenzweig (1973) M. L. Rosenzweig, “Exploitation in three trophic levels,” The American Naturalist 107, 275–294 (1973).
  * Hastings and Powell (1991) A. Hastings and T. Powell, “Chaos in a three-species food chain,” Ecology 72, 896–903 (1991).
  * Gaspard and Nicolis (1983) P. Gaspard and G. Nicolis, “What can we learn from homoclinic orbits in chaotic dynamics?” Journal of statistical physics 31, 499–518 (1983).
  * Barnett  _et al._ (2022) W. A. Barnett, G. Bella, T. Ghosh, P. Mattana, and B. Venturi, “Shilnikov chaos, low interest rates, and new keynesian macroeconomics,” Journal of Economic Dynamics and Control 134, 104291 (2022).
  * Bella, Mattana, and Venturi (2017) G. Bella, P. Mattana, and B. Venturi, “Shilnikov chaos in the lucas model of endogenous growth,” Journal of Economic Theory 172, 451–477 (2017).
  * Wei, Qi-Chang, and Rui-Lan (2010) W. Wei, Z. Qi-Chang, and T. Rui-Lan, “Shilnikov sense chaos in a simple three-dimensional system,” Chinese Physics B 19, 030517 (2010).
  * Yao and Zhang (2005) M. Yao and W. Zhang, “Multipulse shilnikov orbits and chaotic dynamics for nonlinear nonplanar motion of a cantilever beam,” International Journal of Bifurcation and Chaos 15, 3923–3952 (2005).
  * Borisov, Kazakov, and Sataev (2016) A. V. Borisov, A. O. Kazakov, and I. R. Sataev, “Spiral chaos in the nonholonomic model of a chaplygin top,” Regular and Chaotic Dynamics 21, 939–954 (2016).
  * Bakhanova  _et al._ (2020) Y. V. Bakhanova, A. O. Kazakov, E. Y. Karatetskaya, A. D. Kozlov, and K. Saphonov, “On homoclinic attractors of three-dimensional flows,” Izvestiya VUZ. Applied Nonlinear Dynamics 28, 231–258 (2020).
  * Bakhanova  _et al._ (2023) Y. V. Bakhanova, S. Gonchenko, A. Gonchenko, A. Kazakov, and E. Samylina, “On shilnikov attractors of three-dimensional flows and maps,” Journal of Difference Equations and Applications 29, 1184–1201 (2023).
  * Malykh  _et al._ (2020) S. Malykh, Y. Bakhanova, A. Kazakov, K. Pusuluri, and A. Shilnikov, “Homoclinic chaos in the rössler model,” Chaos: An Interdisciplinary Journal of Nonlinear Science 30 (2020).
  * Karatetskaia, Shykhmamedov, and Kazakov (2021) E. Karatetskaia, A. Shykhmamedov, and A. Kazakov, “Shilnikov attractors in three-dimensional orientation-reversing maps,” Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (2021).
  * Mongillo, Barak, and Tsodyks (2008) G. Mongillo, O. Barak, and M. Tsodyks, “Synaptic theory of working memory,” Science 319, 1543–1546 (2008).
  * Lazarevich, Stasenko, and Kazantsev (2017) I. A. Lazarevich, S. Stasenko, and V. B. Kazantsev, “Synaptic multistability and network synchronization induced by the neuron–glial interaction in the brain,” JETP letters 105, 210–213 (2017).
  * Wilson and Cowan (1972) H. R. Wilson and J. D. Cowan, “Excitatory and inhibitory interactions in localized populations of model neurons,” Biophysical journal 12, 1–24 (1972).
  * Tsodyks, Pawelzik, and Markram (1998) M. Tsodyks, K. Pawelzik, and H. Markram, “Neural networks with dynamic synapses,” Neural computation 10, 821–835 (1998).
  * Gordleeva (2012) S. Y. Gordleeva, “Bi-directional astrocytic regulation of neuronal activity within a network,” Frontiers in computational neuroscience 6, 92 (2012).
  * Araque  _et al._ (1999) A. Araque, V. Parpura, R. P. Sanzgiri, and P. G. Haydon, “Tripartite synapses: glia, the unacknowledged partner,” Trends in neurosciences 22, 208–215 (1999).
  * Gonchenko, Li, and Turaev (2022) S. Gonchenko, D. Li, and D. Turaev, “Universal bifurcation patterns in the unfolding of a pair of homoclinic tangencies,” arXiv preprint arXiv:2212.03964 (2022).
  * Veltz (2020) R. Veltz, “BifurcationKit.jl,”  (2020).
  * Dhooge, Govaerts, and Kuznetsov (2003) A. Dhooge, W. Govaerts, and Y. A. Kuznetsov, “Matcont: a matlab package for numerical bifurcation analysis of odes,” ACM Transactions on Mathematical Software (TOMS) 29, 141–164 (2003).
  * Kuznetsov, Kuznetsov, and Kuznetsov (1998) Y. A. Kuznetsov, I. A. Kuznetsov, and Y. Kuznetsov,  _Elements of applied bifurcation theory_ , Vol. 112 (Springer, 1998).
  * Pomeau and Manneville (1980) Y. Pomeau and P. Manneville, “Intermittent transition to turbulence in dissipative dynamical systems,” Communications in Mathematical Physics 74, 189–197 (1980).
  * Gonchenko, Gonchenko, and Shilnikov (2012) A. Gonchenko, S. Gonchenko, and L. Shilnikov, “Towards scenarios of chaos appearance in three-dimensional maps, nelin,” (2012).
  * Gonchenko  _et al._ (2014) A. Gonchenko, S. Gonchenko, A. Kazakov, and D. Turaev, “Simple scenarios of onset of chaos in three-dimensional maps,” International Journal of Bifurcation and Chaos 24, 1440005 (2014).
  * Gonchenko  _et al._ (2019) S. V. Gonchenko, A. S. Gonchenko, A. O. Kazakov, A. D. Kozlov, and Y. V. Bakhanova, “Mathematical theory of dynamical chaos and its applications: Review. part 2. spiral chaos of three-dimensional flows,” Izvestiya VUZ. Applied Nonlinear Dynamics 27, 7–52 (2019).
  * Rössler (1979) O. E. Rössler, “Continuous chaos—four prototype equations,” Annals of the New York Academy of Sciences 316, 376–392 (1979).
  * Argoul, Arneodo, and Richetti (1987) F. Argoul, A. Arneodo, and P. Richetti, “Experimental evidence for homoclinic chaos in the Belousov-Zhabotinskii reaction,” Physics Letters A 120, 269–275 (1987).
  * Letellier, Dutertre, and Maheu (1995) C. Letellier, P. Dutertre, and B. Maheu, “Unstable periodic orbits and templates of the Rössler system: toward a systematic topological characterization,” Chaos: An Interdisciplinary Journal of Nonlinear Science 5, 271–282 (1995).
  * Hindmarsh and Rose (1984) J. L. Hindmarsh and R. Rose, “A model of neuronal bursting using three coupled first order differential equations,” Proceedings of the Royal society of London. Series B. Biological sciences 221, 87–102 (1984).
  * Brette  _et al._ (2007) R. Brette, M. Rudolph, T. Carnevale, M. Hines, D. Beeman, J. M. Bower, M. Diesmann, A. Morrison, P. H. Goodman, F. C. Harris, _et al._ , “Simulation of networks of spiking neurons: a review of tools and strategies,” Journal of computational neuroscience 23, 349–398 (2007).
  * Izhikevich and Edelman (2008) E. M. Izhikevich and G. M. Edelman, “Large-scale model of mammalian thalamocortical systems,” Proceedings of the national academy of sciences 105, 3593–3598 (2008).
  * Izhikevich (2007) E. M. Izhikevich,  _Dynamical systems in neuroscience_ (MIT press, 2007).
  * Montbrió, Pazó, and Roxin (2015) E. Montbrió, D. Pazó, and A. Roxin, “Macroscopic description for networks of spiking neurons,” Physical Review X 5, 021028 (2015).
  * Lebedeva  _et al._ (2023) A. Lebedeva, K. Maltseva, N. Barabash, R. Sokolov, A. Rozov, and T. Levanova, “The effect of astrocytes in synaptic transmission and rhythmogenesis in hippocampus,” in  _2023 Fifth International Conference Neurotechnologies and Neurointerfaces (CNN)_ (IEEE, 2023) pp. 42–46.
  * Wagenaar, Nadasdy, and Potter (2006) D. A. Wagenaar, Z. Nadasdy, and S. M. Potter, “Persistent dynamic attractors in activity patterns of cultured neuronal networks,” Physical Review E 73, 051907 (2006).
  * Rabinovich and Abarbanel (1998) M. Rabinovich and H. Abarbanel, “The role of chaos in neural systems,” Neuroscience 87, 5–14 (1998).
  * Bakhanova  _et al._ (2022) Y. V. Bakhanova, S. Gonchenko, A. Gonchenko, A. Kazakov, and E. Samylina, “On Shilnikov attractors of three-dimensional flows and maps,” Journal of Difference Equations and Applications , 1–18 (2022).
  * Perea, Navarrete, and Araque (2009) G. Perea, M. Navarrete, and A. Araque, “Tripartite synapses: astrocytes process and control synaptic information,” Trends in neurosciences 32, 421–431 (2009).
  * Agulhon  _et al._ (2008) C. Agulhon, J. Petravicz, A. B. McMullen, E. J. Sweger, S. K. Minton, S. R. Taves, K. B. Casper, T. A. Fiacco, and K. D. McCarthy, “What is the role of astrocyte calcium in neurophysiology?” Neuron 59, 932–946 (2008).
  * Barres (2008) B. A. Barres, “The mystery and magic of glia: a perspective on their roles in health and disease,” Neuron 60, 430–440 (2008).
  * Braun, Lisboa, and Gallas (1992) T. Braun, J. A. Lisboa, and J. A. C. Gallas, “Evidence of homoclinic chaos in the plasma of a glow discharge,” Phys. Rev. Lett. 68, 2770–2773 (1992).