# Bursting in a next generation neural mass model with synaptic dynamics: a
slow-fast approach

Halgurd Taher1 1Inria Sophia Antipolis Méditerranée Research Centre, 2004
Route des Lucioles, 06902 Valbonne, France<EMAIL_ADDRESS>Daniele
Avitabile1,2 2Department of Mathematics, Faculteit der Exacte Wetenschappen,
Vrije Universiteit (VU University Amsterdam), De Boelelaan 1081a, 1081 HV
Amsterdam, Netherlands<EMAIL_ADDRESS>Mathieu Desroches1
<EMAIL_ADDRESS>
###### Abstract

We report a detailed analysis on the emergence of bursting in a recently
developed neural mass model that takes short-term synaptic plasticity into
account. Neural mass models are capable of mimicking the collective dynamics
of large scale neuronal populations in terms of a few macroscopic variables
like mean membrane potential and firing rate. The one being used here
particularly important, as it represents an exact meanfield limit of
synaptically coupled quadratic integrate & fire neurons, a canonical model for
type I excitability. In absence of synaptic dynamics, a periodic external
current with a slow frequency $\varepsilon$ can lead to burst-like dynamics.
The firing patterns can be understood using techniques of singular
perturbation theory, specifically slow-fast dissection. In the model with
synaptic dynamics the separation of timescales leads to a variety of slow-fast
phenomena and their role for bursting is rendered inordinately more intricate.
Canards are one of the main slow-fast elements on the route to bursting. They
describe trajectories evolving nearby otherwise repelling locally invariant
sets of the system and are found in the transition region from subthreshold
dynamics to bursting. For values of the timescale separation nearby the
singular limit $\varepsilon=0$, we report peculiar jump-on canards, which
block a continuous transition to bursting. In the biologically more plausible
regime of $\varepsilon$ this transition becomes continuous and bursts emerge
via consecutive spike-adding transitions. The onset of bursting is of complex
nature and involves mixed-type like torus canards, which form the very first
spikes of the burst and revolve nearby fast-subsystem repelling limit cycles.
We provide numerical evidence for the same mechanisms to be responsible for
the emergence of bursting in the quadratic integrate & fire network with
plastic synapses. The main conclusions apply for the network, owing to the
exactness of the meanfield limit.

  * Date:

##  1 Introduction

In the past decade a novel approach in meanfield theory, the so-called Ott-
Antonson (OA) ansatz [1, 2], has received major attention. The ansatz serves
as a recipe to perform an exact reduction from a high dimensional network of
interacting phase oscillators towards a low dimensional dynamical system, that
describes the macroscopic behavior. It was firstly applied to the prototypical
model for synchronization, namely the Kuramoto model [3]. The term exact
refers to the property of the meanfield limit to rely on first principles.
Thus, in the thermodynamic limit one obtains an exact agreement of the
meanfield system with respect to the macroscopic dynamics of the underlying
network.

Remarkably, the OA ansatz is applicable to a broader class of phase oscillator
networks and thereby has found its way into the field of computational
neuroscience. One system of interest is the Quadratic Integrate & Fire (QIF)
neuron, which represents a canonical model for Saddle-Node on Invariant Circle
(SNIC) bifurcations and so-called type I excitability. Given this canonical
role of the QIF neuron, it is standing to reason to investigate the dynamics
emerging on macroscopic scale. Here the OA ansatz comes into play. Under
certain assumptions the QIF neuron is rendered equivalent to the Ermentrout-
Kopell model, also referred to as $\theta$-model [4]. It represents a phase
oscillator model for neuronal dynamics, which indeed fits into the class of
problems that can be treated within the OA framework. The OA ansatz applied to
an ensemble of QIF neurons, in the work of Montbrió, Pazó and Roxin (MPR), has
lead to an upsurge of next generation neural mass models [5].

These models try to capture the macroscopic dynamics of networks of spiking
neurons, using just a few variables, like here, the population firing rate and
mean membrane potential. They hence contribute crucially to studies of
collective phenomena in these otherwise high dimensional dynamical systems and
allow for simple analytical and numerical treatment. Various applications of
the MPR model has been studied in recent years. They range from the inclusion
of delayed synaptic interactions [6], giving rise to chaos, to studies of
cortical oscillations in multipopulation models [7] and cross-frequency
coupling [8]. While the original MPR model accounts for chemical synapses, the
methodology can straightforwardly been applied to include as well electrical
synapses, formed by gap junctions between neurons [9, 10]. Extensions of the
MPR model towards sparse networks [11] and fluctuation driven dynamics have
been proposed [12].

These examples of QIF networks and their neural mass counterparts can give
rise to interesting dynamical regimes, typically evoked by bistability in the
system. Indeed, the original MPR model exhibits a parameter regime where a
stable node and a focus coexist [5]. This is particularly relevant for the
macroscopic response of neuronal ensembles when they are subject to an
external current: bistability implies that a time dependent external drive can
lead to interesting firing rhythms.

A recent study takes into account synaptic dynamics in form of exponentially
decaying action potentials [13]. In the present work however, we want to
explore a QIF network that accounts for synaptic dynamics in form of short-
term synaptic plasticity (STP), thus adding to the biological plausibility.
According to the phenomenological STP model of Tsodyks and Markram two
opposing effects must be distinguished: depression, i.e. weakening, and
facilitation, i.e. strengthening, of synaptic connections [14]. To our
knowledge, the role of STP in next generation neural mass models has received
just little attention, despite being highly relevant in neuroscience. Previous
macroscopic models of STP typically make use of the Wilson-Cowan (WC) model,
hence are of heuristic nature [15, 16]. Nevertheless they helped to develop a
novel synaptic theory of working memory [17], a cognitive system for short-
term information storage and manipulation in the brain. The synaptic theory of
WM has triggered various theoretical studies, which use the WC model with STP
in multipopulation topologies in order to understand basic WM operations like
information loading and recall, as well as to estimate the maximum WM capacity
[18, 19].

In a recent study an extension of the MPR firing rate equations towards STP
was proposed, in order to model WM [20]. The meanfield limit, in presence of
STP, remains exact. Therefore one can exploit this limit, in order to get
insight into the emergence of firing patterns in the network. An aspect that
can lead to complex behavior is the timescale separation, which comes along
with STP. Depression and facilitation might indeed act on different
timescales. As an example, measurements in the prefrontal cortex suggest that
the facilitation of synapses can be maintained for seconds, while depression
decays within a few hundred milliseconds [21].

Synaptic dynamics and additional timescales enrich the dynamical landscape, by
giving rise to bistability involving limit cycles [20]. This is the foundation
for bursting rhythms to emerge. Bursting refers to dynamics that alternates
between a quiescent phase and rapid oscillations. When slowly forcing the
population of QIF neurons, by virtues of a slowly drifting external current,
the system can transit from a quasi-static motion to rapid oscillations
associated with the presence of stable cycles in the system with constant
external current.

Bursting has been found in various experimental studies in neuroscience [22,
23, 24, 25, 26, 27, 28, 29] and theoretical approaches [30, 31, 32, 33, 34,
35] not only aim at classifying the observed dynamics, but also mimicking and
revealing the mechanisms responsible for the emergence of bursting. While
bursting in spiking neural networks is subject of recent studies [36, 37], the
mechanisms responsible for their emergence often remain unclear: exploring the
state space of large scale networks is tedious and the addition of slow-fast
aspects complicates the problem. The exactness of the MPR model helps to
overcome this limitation: analytic tools and bifurcation analysis applied to
the neural mass model allow to draw conclusions for the microscopic network.

The main results of this work are related to the emergence of bursting in a
QIF network with STP. In particular, we investigate the transition from
subthreshold oscillations to bursting in presence of an external slow and
periodic current. The forcing introduces a clear timescale separation into the
problem, giving rise to intricate slow-fast phenomena and allowing for the
application of slow-fast dissection methods, to be described later. As an
outlook, the findings comprise a differentiation of the route to bursting,
depending on the timescale separation. For strongly separated timescales, far
away from biologically plausible scenarios, the route is complicated, possibly
discontinuous in parameter space and it is related to a certain type of so-
called canards. However, moderate timescale separation reveals a number of
intermingled slow-fast mechanisms that lead to a continuous transition from
subthreshold oscillations to bursting and are related to different types of
canards. Our results are supported by slow-fast arguments and numerical
evidence.

A first illustration of the dynamical regime of interest in this work is
displayed in Fig. 1. Panel (a) depicts the response of a large scale network
consisting of $N=10^{5}$ QIF neurons to a slow external sinusoidal current.
The second panel (b) shows the firing rate of the QIF network, as well as the
firing rate of the meanfield limit. Both systems undergo a quiescence phase of
low firing activity. When the external current exceeds a certain level, the
systems start to burst, characterized by a rapid series of synchronized firing
at high rates.

Figure 1: Spiking neuron network and meanfield limit: (a) Scatter plot of a
subset of 2000 representative neurons out of $10^{5}$. Each dot represent a
spike. (b) Firing rate of the network (black) and of the corresponding
meanfield limit (red). The blue curve shows a time dependent external current
applied to the two systems and is of sinusoidal form.

In order to understand how these bursts emerge, we have to encapsulate two
main aspects. First, in the upcoming Section 2, we will introduce the QIF
network model with STP as well as the corresponding meanfield limit, and we
will analyze the state space structure and dynamics. Second, the presence of a
slow external drive calls for the application of slow-fast dissection. For
this it is essential to introduce the general slow-fast framework and revise
well known slow-fast mechanisms, which play a role for our model, in Section
3. As a next step, this generic methodology for timescale separated problems
is applied to the present model in Section 4. Dissection is crucial for
understanding the results in Section 5, where canards and in particular jump-
on canards are studied using slow-fast arguments. This paves the way to
investigate the mechanisms responsible for the emergence of bursting, as done
in Section 6. Finally, in Section 7 the initially posed problem of bursting on
macroscopic scale is approached by a comparison of meanfield dynamics versus
QIF network dynamics in the bursting regime.

##  2 A next-generation neural mass model with synaptic dynamics

Spiking neuron models can be characterized by their response to the injection
of a current, which is often measured in terms of the $f$-$I$ curve,
determining the relation of firing frequency $f$ versus input current $I$. The
dynamics of Hodgkin-Huxley type neurons can either be in the excitable or
tonic regime [38]. Excitable neurons in absence of input approach an
equilibrium. However, sufficient input can excite the membrane potential
beyond a threshold leading to the firing of a single action potential, before
going back to the rest state. Tonic neurons on the other hand fire
periodically with a frequency $f$. Based on the behaviour at the transition
from excitable to tonic dynamics one can distinguish (at least) two classes of
membranes. For class I neurons the $f$-$I$ curve is continuous and transitions
from quiescence ($f=0$) to repetitive firing at arbitrarily slow frequencies
($f>0$). Typically it occurs at a SNIC bifurcation. Class II neurons on the
other hand exhibit a discontinuous $f$-$I$ curve, leading to finite firing
rates at the onset of tonic firing, and they are usually associated with a
Hopf bifurcation. This Hopf bifurcation is often subcritical, for example in
the Hodgkin-Huxley, FitzHugh-Nagumo [39, 40] and Morris-Lecar model [41].

###  2.1 Network of spiking neurons and meanfield limit

The model of interest for our work the canonical model for type I
excitability: the QIF neuron. In a network of $N$ synaptically coupled neurons
the membrane potentials $V_{i}(t)$ obey Eq. 1.

$\displaystyle\dot{V}_{i}$ $\displaystyle=V_{i}^{2}+\eta_{i}+Jr(t)+I_{1}(t)$
(1a) if $V_{i}>V_{\rm thresh}$: $V\leftarrow V_{\rm r}$ $\displaystyle r(t)$
$\displaystyle=\frac{1}{N}\sum_{j=1}^{N}\sum_{k:t_{j}^{(k)}<t}\delta(t-t_{j}^{(k)})$
(1b)

The total current applied to the neuron is a sum of the constant component
$\eta_{i}$, the synaptic current $Jr(t)$, with synaptic weight $J$ and an
external, possibly time-dependent, current $I_{1}(t)$. Variable $r$ denotes
the instantaneous firing rate and is composed of the single neuron spike
trains $\sum_{k:t_{j}^{(k)}<t}\delta(t-t_{j}^{(k)})$, where $t_{j}^{(k)}$
denotes the $k$-th spike time of neuron $j$ entering into the Dirac $\delta$
function. Whether a QIF neuron is excitable or tonic depends on $\eta_{i}$.
Given $r=0$, the neuron is excitable for $\eta_{i}<0$ and in tonic firing
state for $\eta_{i}>0$. Firing occurs whenever $V_{i}$ exceeds the threshold
$V_{\rm thresh}$ at which the reset rule applies leading to a reset of the
potential to $V_{\rm r}$.

When performing the thermodynamic limit $N\to\infty$ and imposing certain
conditions, the mean dynamics of the above microscopic model leads to a
reduced macroscopic description in terms of the mean membrane potential $v(t)$
and firing rate $r(t)$, namely the MPR model [5]. The derivation is based on
the OA ansatz and yields an exact reduction [1]. Thus the collective behavior
of the QIF network, aside from finite size fluctuations, is in perfect
agreement with the MPR model. Following the OA Ansatz and the derivation of
the MPR model the following assumptions have to be made in order to obtain an
exact firing rate formulation.

(i) The threshold and reset voltage have to be considered in the limit $V_{\rm
thresh}=-V_{\rm r}\rightarrow\infty$, rendering the QIF neuron identical to
the $\theta$-model [4].

(ii) The excitabilities $\eta_{i}$ are drawn from a Lorentzian distribution
$g(\eta)=\frac{1}{\pi}\frac{\Delta}{(\eta-\bar{\eta})^{2}+\Delta^{2}}$,
centred at $\bar{\eta}$ and with the width parameter $\Delta$.

(iii) Neurons are all-to-all coupled. This way each QIF neuron receives
identical synaptic current $Jr(t)$.

(iv) The QIF network has to be considered in the thermodynamic limit
$N\rightarrow\infty$.

The resulting MPR model consists of two ordinary differential equations for
$r(t)$ and $v(t)$ given in Eq. 2.

$\displaystyle\dot{r}$ $\displaystyle=\frac{\Delta}{\pi}+2rv$ (2a)
$\displaystyle\dot{v}$ $\displaystyle=v^{2}-(\pi
r)^{2}+Jr+\bar{\eta}+I_{1}(t)$ (2b)

Despite having a rather simple state-space structure, the MPR model can give
rise to interesting periodic patterns when externally forced. In the case of
constant $I_{1}$ periodic solutions are absent and one can find node, focus
and saddle equilibria. However, there are regions of parameter space in which
bistability between the node and focus appears. Hence slow periodic forcing,
for example given by $I_{1}=A\sin(\varepsilon t),\quad 0<\varepsilon\ll 1$,
can lead to a hysteretic loop in these regions, as shown in [5]. In this case
the trajectories consist of a low firing rate segment and a high firing rate
segment with damped oscillations, related to the presence of foci in the
system with constant $I_{1}$. Orbits of this type can already be seen as
bursting patterns characterized by an alternation of slow drifts and fast
oscillations. However, in the limit $\varepsilon=0$ of infinitely slow
forcing, the fast oscillations vanish. In that case, the resulting cycles can
be classified as relaxation oscillations, which are introduced in Section 3.

###  2.2 Neural mass model with short-term plasticity

Our present study investigates bursting patterns in an extended version of Eq.
2 that accounts for STP as described in the phenomenological model of Tsodyks
and Markram [14]. Short-term synaptic depression is related to
neurotransmitter depletion. Each neuron $i$ has a limited amount
$X_{i}(t)\in[0,1]$ of resources (i.e vesicles ready to be released). Spiking
is followed by the emittance of presynaptic action potentials. Upon their
arrival at the synaptic terminal a fraction $U_{i}(t)\in[U_{0},1]$ of the
neurotransmitters is released into the synaptic cleft, resulting in the
generation of postsynaptic potentials (PSPs). Therefore, each presynaptic
spike is linked to the utilization and reduction of resources available for
the generation of upcoming PSPs, consequently leading to a decrease of future
postsynaptic excitations. The resource $X_{i}$ exponentially recovers to its
base value of $X_{i}=1$ on a timescale
$\tau_{\mathrm{d}}=$200\text{\,}\mathrm{m}\mathrm{s}$$ (depression timescale).

Facilitation, as opposed to depression, leads to enhanced PSPs and is related
the neurotransmitter release probability at the synaptic terminals, which is
modelled by the utilization factor $U_{i}$. The release probability (and
therefore $U_{i}$) depends on the intracellular calcium concentration. The
neurotransmitter release is associated with the accumulation of calcium ions
in the presynaptic terminal, hence each spike leads to an increase of $U_{i}$.
Calcium concentration and the utilization factor decay to the base level
$U_{i}=U_{0}$ on the facilitation timescale
$\tau_{\mathrm{f}}=$1500\text{\,}\mathrm{m}\mathrm{s}$$.

We will focus on an implementation of STP into the model on macroscopic level
(m-STP) as suggested in [20], in order maintain the exactness of the firing
rate model. In other words, depression and facilitation, accounted for by
$X_{i}$ and $U_{i}$ respectively, will not be treated on single neurons level,
but rather on population level, with the depression and facilitation variables
$x(t)$ and $u(t)$, respectively. This results in $N$ membrane potential
equations and two synaptic equations for the QIF network, as given in Eq. 3.

$\displaystyle\dot{V}_{i}$ $\displaystyle=V_{i}^{2}+\eta_{i}+Juxr+I_{1}(t)$
(3a) $\displaystyle\dot{x}$ $\displaystyle=\frac{1-x}{\tau_{\mathrm{d}}}-uxr$
(3b) $\displaystyle\dot{u}$
$\displaystyle=\frac{U_{0}-u}{\tau_{\mathrm{f}}}+U_{0}(1-u)r$ (3c)

The amount of resources $x$ of the QIF network reduces when the population
firing rate $r$ increases, while at the same time the utilization $u$
increases. Both quantities enter into the effective synaptic weight $Jux$. The
extended system given in Eq. 4, in the following referred to as neural mass
with STP (NMSTP), represents an exact meanfield limit of the QIF network given
in Eq. 3. The state variables are the firing rate $r$, mean membrane potential
$v$, amount of resources $x$ and utilization factor $u$.

$\displaystyle\dot{r}$ $\displaystyle=\frac{\Delta}{\pi}+2rv$ (4a)
$\displaystyle\dot{v}$ $\displaystyle=v^{2}-(\pi
r)^{2}+Juxr+\bar{\eta}+I_{1}(t)$ (4b) $\displaystyle\dot{x}$
$\displaystyle=\frac{1-x}{\tau_{\mathrm{d}}}-uxr$ (4c) $\displaystyle\dot{u}$
$\displaystyle=\frac{U_{0}-u}{\tau_{\mathrm{f}}}+U_{0}(1-u)r$ (4d)

We note that Eq. 4 will evolve on the fastest timescale of our problem. This
holds despite the fact that it already possesses multiple timescale via
$\tau_{\mathrm{d}}$ and $\tau_{\mathrm{f}}$. However, as we will discuss
later, this inherent timescale separation of the NMSTP is subtle and not
observable everywhere in state space. Nevertheless, it has significant impact
on how the transition from subthreshold (non-bursting) behaviour to bursting
occurs (see also Sections 5.2 and 6.3).

###  2.3 Dynamics under constant forcing

Most of the parameters values used for Eq. 4 will remain fixed and, if not
stated differently, given in Table 1. Note that time is measured in units of
the membrane time constant $\tau_{\mathrm{m}}$. For more details on the
numerical methods, we refer to the supplementary material of this work.

Table 1: Parameters and their values, which are fixed throughout this work, if not stated differently. Symbol | Description | Value
---|---|---
$\Delta$ | Width of Lorentzian | 0.5
$\bar{\eta}$ | Centre of Lorentzian | -1.7
$J$ | Synaptic weight | 30
$U_{0}$ | Baseline utilization | 0.1
$\tau_{\mathrm{m}}$ | Membrane time constant | $20\text{\,}\mathrm{m}\mathrm{s}$
$\tau_{\mathrm{d}}$ | Depression timescale | $$200\text{\,}\mathrm{m}\mathrm{s}$/\tau_{\mathrm{m}}$
$\tau_{\mathrm{f}}$ | Facilitation timescale | $$1500\text{\,}\mathrm{m}\mathrm{s}$/\tau_{\mathrm{m}}$

We will outline the different dynamical regimes in presence of a constant
current $I_{1}(t)=\text{const.}$, using the above parameter values. In Fig.
2(a) the resulting bifurcation diagram is displayed.

For currents $I_{1}\lesssim 0.25$ we find a branch of stable node equilibria
at low firing rates. The branch develops into a family of foci and
destabilizes around $I_{1}\approx 0.25$ via a subcritical Hopf bifurcation
($H_{\rm L}$) followed by two saddle-node (fold) bifurcation at $F_{\rm L}$
and $F_{\rm U}$ (black dots), where
$F_{k}=(r_{k},v_{k},x_{k},u_{k},I_{k}),\quad k\in\\{\mathrm{L},\mathrm{U}\\}$,
denotes the equilibrium and parameter values of the bifurcations. These folds
can also be found in absence of STP, in which case the upper branch is stable.
However, in Fig. 2(a) the instability persists throughout the S-shaped curve
up until the upper supercritical Hopf bifurcation $H_{\rm U}$. The lower Hopf
bifurcation $H_{\rm L}$ generates a family of unstable limit cycles that
undergoes a fold bifurcation of cycles, giving rise to stable periodic
solutions.

Figure 2: Solution families of the system with constant forcing: (a)
Bifurcation diagram $r$ versus $I_{1}$ of Eq. 4. For $I_{1}\lesssim 0.25$
there exists only one fixed point (FP, solid black line). At $I_{1}\approx
0.25$ the FP destabilizes via a subcritical Hopf-Bifurcation ($H_{\rm L}$,
lower orange dot), creating a branch of unstable limit cycles (LC, purple
dashed line). Two saddle-node bifurcations (black dots) of the unstable FP
branch (dashed black line) occur in a narrow regime of $I_{1}$, folding the
branch twice. Stability is regained for $I_{1}\gtrsim 0.7$ at a supercritical
Hopf-Bifurcation ($H_{\rm U}$, upper orange dot). The unstable LC stabilizes
(solid purple line) via a saddle-node bifurcation of cycles (purple dot) and
vanishes at the second Hopf Bifurcation. The purple line marks the maximum
firing rates of the LC branch. (c - e) Periodic solution
$(r(t),v(t),x(t),u(t))$ vs. time $t$ for $I_{1}=0.4$ marked in panel (a) by a
dashed red line. The red curves show simulations results of the NMSTP, the
black ones of the network Eq. 3. (b) Spike scatter plot for 20000
representative neurons in the network out of $N=100000$.

One of these trajectories $(r(t),v(t),x(t),u(t))$ is presented in Fig. 2(b-f)
as a function of time. It is superimposed onto the corresponding variables
calculated by simulating a QIF network governed by Eq. 3 and consisting of
$N=100000$ neurons. For this network, the firing rate is estimated via binning
of time, i.e, by counting the number of spikes per time bin of width $\Delta
t=10^{-2}$. The average membrane potential for the network reads
$v(t)=\frac{1}{N}\sum_{j=1}^{N}V_{j}(t)$.

The primary mechanism driving the oscillations is an interplay of so-called
populations bursts and the ensuing synaptic depression and facilitation. At
the microscopic scale, population bursts are emitted via a cascade of spikes
throughout the network, which as a consequence leads to the facilitation of
synapses, leveraging the firing activity further; see Fig. 2(b,c,f). The
consequent depression suppresses the activity, but recovers on the timescale
$\tau_{\mathrm{d}}$ allowing for the emittance of population bursts in a
periodic manner.

Notably, in the $I_{1}$-interval depicted in the inset of Fig. 2(a) we find
bistability between equilibria and limit cycles. We can therefore predict that
a time dependent slow current $I_{1}(t)$, evolving across this region, will
lead to a dynamic transition from the equilibrium branch to the stable limit
cycles, giving rise to bursting. This exact example can be found in Fig. 1.

Overall, in contrast to the QIF network without STP and original MPR model,
where no limit cycles exist, STP gives rise to bistability among equilibria
and cycles. In Ref. [5] a slow periodic currents leads to the emergence of
macroscopic relaxation-type oscillations in the network. We want to
investigate how the presence of STP impacts the response of the system towards
such an input. Simulations of QIF networks are difficult computationally.
However, the expected agreement of QIF network results and the NMSTP depicted
in Fig. 2(b-f) justifies to perform the upcoming analysis using solely the
NMSTP. We will return to the implication of NMSTP dynamics for the network in
Section 7.

###  2.4 Dynamics under slow periodic forcing

Owing to the previous observations in the system with constant forcing, we
will introduce a slow periodic drive into the model via the external current
$I_{1}$. We impose that it evolves periodically and on a timescale
considerably larger than the slowest timescale of the neural mass, namely the
facilitation decay time $\tau_{\mathrm{f}}$. In order to remain in a general
framework, $I_{1}(t)$ will be sinusoidal, given by
$I_{1}(t)=A\sin{(\varepsilon t)}$, with period
$T=\frac{2\pi}{\varepsilon}\gg\tau_{\mathrm{f}}$ and amplitude $A$. Throughout
this work we set
$\tau_{\mathrm{f}}=$1500\text{\,}\mathrm{m}\mathrm{s}$/\tau_{\mathrm{m}}=75$,
therefore the separation between forcing and slowest intrinsic timescale of
the fast subsystem is calculated as
$\tau_{\mathrm{f}}/T=\varepsilon\frac{\tau_{\mathrm{f}}}{2\pi}\approx
10\varepsilon$.

Through the choice of $I_{1}$ to be explicitly time dependent, the system
given in Eq. 4 becomes non-autonomous. This in turn comes along with hurdles
in the application of slow-fast dissection. Thus, in order to retrieve an
autonomous system, a second forcing variable $I_{2}$ is introduced. The
dynamics of ($I_{1},I_{2})$ follows a Hopf normal form as given below.

$\displaystyle\dot{I_{1}}$ $\displaystyle=\varepsilon
g_{1}(I_{1},I_{2})=\varepsilon\left[I_{1}(a-I_{1}^{2}-I_{2}^{2})+I_{2}\right]$
(5a) $\displaystyle\dot{I_{2}}$ $\displaystyle=\varepsilon
g_{2}(I_{1},I_{2})=\varepsilon\left[I_{2}(a-I_{1}^{2}-I_{2}^{2})-I_{1}\right]$
(5b)

The Hopf bifurcation at $a=0$ gives rise to stable limit cycles of the form
$(I_{1},I_{2})=A\cdot(\sin{\varepsilon t},\cos{\varepsilon t})$, with
amplitude $A=\sqrt{a}$ and angular frequency $\varepsilon$, in the following
referred to as forcing cycle. To assure equivalence of the explicitly defined
$I_{1}(t)=A\sin(\varepsilon t)$ and the one generated by the Hopf form Eq. 5,
the initial conditions $\left((I_{1}(t_{0}),I_{2}(t_{0})\right)$ will lie on
$(I_{1},I_{2})=A\cdot(\sin{\varepsilon t},\cos{\varepsilon t})$. The full
system is given by the NMSTP in presence of slow external forcing, i.e, Eqs. 4
to 5.

To understand the impact of this slow forcing, it is advantageous to
superimpose solutions of the full problem on the $r$ vs. $I_{1}$ bifurcation
diagram of the unforced system, this is at the core of the slow-fast
dissection introduced by J. Rinzel [32, 33, 34]. An example of a purely slow
trajectory is shown in Fig. 3(a1-c1) and labeled $\boldsymbol{\gamma}_{0}(t)$.
The firing rate $r(t)$, shown in panel (b1), increases and decreases following
the same pattern as the forcing $I_{1}(t)$ in panel (a1). Moreover, in
$r-I_{1}$ projection, shown in panel (c1), it becomes clear that the dynamics
takes place nearby the equilibrium branch of the unforced system. The forcing
introduces a drift of the equilibrium, slow enough to be followed by the
dynamics in an $\mathcal{O(\varepsilon)}$ neighbourhood of the branch.

While this example can be understand as a quasi-static motion, the more
complex solutions $\boldsymbol{\gamma}_{1}(t)$ to $\boldsymbol{\gamma}_{4}(t)$
in columns 2 and 3 of Fig. 3, exhibit canard dynamics and bursting,
respectively. A more rigorous analysis is required, including a slow-fast
dissection of the model. In the next sections a generic framework in which
timescale separated problems can be treated is introduced and applied to the
NMSTP.

Figure 3: Typical solutions $\boldsymbol{\gamma}_{0}(t)$ to
$\boldsymbol{\gamma}_{4}(t)$ of the full system: (a) Periodic forcing current
$I_{1}(t)$ and (b) firing rate $r(t)$ vs. time $t$. (c) Same trajectories
superimposed on the bifurcation diagram of the unforced system in $r$-$I_{1}$.
The parameter values are as follows: $\varepsilon=10^{-5}$;
$\boldsymbol{\gamma}_{0}$: $A\approx 0.2487$; $\boldsymbol{\gamma}_{1}$ to
$\boldsymbol{\gamma}_{3}$: in increasing order exponentially close to
$A\approx 0.2507$; $\boldsymbol{\gamma}_{4}$: $\varepsilon=10^{-3}$, $A\approx
0.2553$.

##  3 Slow-fast framework and state of the art

The dynamics of slow-fast systems can be regarded in terms of fast variables
$\mathbf{X}_{\mathrm{f}}(t)\in\mathbb{R}^{k}$ and slow variables
$\mathbf{X}_{\mathrm{s}}(t)\in\mathbb{R}^{l}$. Their dynamics is governed by
the differential equations given in Eq. 6 and here referred to as full system,

$\displaystyle\dot{\mathbf{X}}_{\mathrm{f}}$
$\displaystyle=~{}\mathbf{F}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}})$
(6a) $\displaystyle\dot{\mathbf{X}}_{\mathrm{s}}$
$\displaystyle=\varepsilon\mathbf{G}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}})$
(6b)

with fast-time parametrisation $t$ (the overdot denoting differentiation with
respect to $t$),
$\mathbf{F}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}}):\mathbb{R}^{k}\times\mathbb{R}^{l}\;\to\;\mathbb{R}^{k}$
and
$\mathbf{G}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}}):\mathbb{R}^{k}\times\mathbb{R}^{l}\;\to\;\mathbb{R}^{l}$.
Here the separation of timescales is reflected by a small parameter
$0<\varepsilon\ll 1$. We will refer to this type of system as $k$-fast
$l$-slow system.

A different formulation of the full system is obtained in Eq. 7 by
parametrizing it in slow time $\tau:=\varepsilon t$.

$\displaystyle\varepsilon\mathbf{X}^{\prime}_{\mathrm{f}}$
$\displaystyle=\mathbf{F}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}})$
(7a) $\displaystyle\mathbf{X}^{\prime}_{\mathrm{s}}$
$\displaystyle=\mathbf{G}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}})\quad.$
(7b)

The derivative with respect to slow time $\tau$ is denoted
$(\\_)^{\prime}:=\differential/\differential\tau(\\_)=\frac{1}{\varepsilon}\dot{(\\_)}$.
The two representations Eq. 6 and Eq. 7 are equivalent, but they allow to
exploit the premise of slow-fast systems, namely the timescale separation
given by a small value of $\varepsilon$, in different ways. It is natural to
consider the singular limit $\varepsilon=0$ and take the fast and slow time
parametrizations into account. One obtains two different subsystems, which
represent a dissection of slow and fast dynamics of the full system.

In the first case we obtain the fast subsystem Eq. 8. This limit can be used
to understand dynamics of the full system for which $\mathbf{X}_{\mathrm{f}}$
evolves fast and results in the following $k+l$ ODEs, $l$ of which being
trivial:

$\displaystyle\dot{\mathbf{X}}_{\mathrm{f}}$
$\displaystyle=\mathbf{F}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}})$
(8a) $\displaystyle\dot{\mathbf{X}}_{\mathrm{s}}$ $\displaystyle=\mathbf{0}$
(8b)

Indeed, the dynamics of the slow variables $\mathbf{X}_{\mathrm{s}}$ is
trivial and their value does not change in time. As a matter of fact they can
be treated as parameters entering into the dynamics of
$\mathbf{X}_{\mathrm{f}}$.

The second limit $\varepsilon\rightarrow 0$, now done in the slow-time
parametrization Eq. 7, yields the slow subsystem, namely:

$\displaystyle\mathbf{0}$
$\displaystyle=\mathbf{F}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}})$
(9a) $\displaystyle\mathbf{X}^{\prime}_{\mathrm{s}}$
$\displaystyle=\mathbf{G}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}})\quad.$
(9b)

Eq. 9 is also referred to as the reduced system and it is represented by a
differential-algebraic system, in which the dynamics of the slow variables
remains unchanged with respect to the full system and is governed by
$\mathbf{X}^{\prime}_{\mathrm{s}}=\mathbf{G}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}})$.
The dynamics of the fast variables on the other hand are hidden within the $k$
algebraic constraints Eq. 9a. They define the critical manifold:

$\displaystyle
S_{0}=\\{(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}})~{}|~{}\mathbf{F}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}})=0\\}\quad,$
(10)

usually a $l$-dimensional manifold embedded in $\mathbb{R}^{(k+l)}$.

In the slow subsystem the dynamics of the fast variables is slaved to the slow
variables, their relation is given by the critical manifold’s equations, which
defines the state space of this limiting problem: motion of the slow subsystem
takes place on $S_{0}$. At the same time points of $S_{0}$ correspond to
equilibria of the fast subsystem, as is clear from equation Eq. 8a. By joining
solutions of the different subsystems at specific points, singular orbits can
be constructed. They are trajectories resulting from the concatenation of slow
and fast segments, for which the dynamics is determined by the respective
subsystem.

Solutions of systems of the type of Eq. 6 with both slow and fast segments are
$\varepsilon$-perturbations of singular orbits and the prototypical slow-fast
cycles are relaxation oscillations [42]. The way these cycles emerge in
parameter space is rather peculiar and involve the famous canard solutions,
which we review in the context of the classical van der Pol (VdP) system
below.

### Classical canards in the van der Pol oscillator.

Here we briefly present the essentials of classical canards in the
prototypical VdP system given in Eq. 11, consisting of one fast variable $x$
and one slow variable $y$. Once written in first-order form, the system’s
equations read:

$\displaystyle\dot{x}$ $\displaystyle=y-\frac{x^{3}}{3}+x$ (11a)
$\displaystyle\dot{y}$ $\displaystyle=\varepsilon(a-x)\quad.$ (11b)

The critical manifold $S_{0}$ is one dimensional and S-shaped, parametrized by
$x$ as $y(x)=\frac{x^{3}}{3}-x$. In VdP the critical manifold
$\\{(x,y);\;y=y(x)\\}$ has a local maximum and a local minimum, at $x=\pm 1$,
respectively. Therefore $S_{0}$ folds twice and has three branches, of which
the middle one is repelling and the other two attracting. They correspond to
unstable and stable equilibria of the fast subsystem, respectively. The slow
nullcline, determined by $\dot{y}=0$, is the straight line $\\{x=a\\}$ and, at
$a=\pm 1$, it intersects $S_{0}$ at the fold points perpendicularly.

Overall, the VdP system, for $0<\varepsilon\ll 1$ and with $a$ as a
bifurcation parameter, displays four regimes, from subthreshold oscillations
to relaxation oscillations. The former is characterized by a quasi-static
motion that takes place entirely nearby one of the stable branches of $S_{0}$
with a slow frequency and low amplitude. Relaxation oscillations on the other
hand comprise a quasi-static motion nearby both attracting branches of
$S_{0}$, joined via fast transitions: when the branch becomes repelling at one
of the fold points, also called jump-off point in this context, a rapid jump
to the opposed attracting branch of $S_{0}$ occurs.

The subthreshold regime terminates at a supercritical Hopf bifurcation, at
which stable small Hopf cycles emanate, which are not yet of relaxation type
and only exist in an $O(\varepsilon)$ distance from the Hopf point. This is
followed by an exponentially narrow parameter interval, for which the orbits
grow in an explosive manner when the parameter is varied. This phenomenon is
known as a canard explosion [43] and the associated canards separate the small
Hopf cycles from relaxation oscillations [44]. They evolve for some time near
the repelling sheet of $S_{0}$, before jumping to one of the attracting
sheets. They represent a mechanism which allows to connect the attracting
sheet of the critical manifold to the repelling one. In our work, as in many
problems involving bursting solutions, canard dynamics play an essential role.
They separate the parameter regime for subthreshold oscillations from the one
where bursting can occur.

The classical canard can be understood by taking into account the slow flow,
i.e, the flow of Eq. 9, which describes the slow dynamics on $S_{0}$, and will
be derived later. For VdP it reads
$(x^{\prime},y^{\prime})=(\frac{a-x}{x^{2}-1},a-x)$ and is in general
undefined at the folds ($x=\pm 1$). However, its $x$-component reduces to
$x^{\prime}=-\frac{1}{x\pm 1}$ for $a=\pm 1$. In this case the flow remains
defined at the corresponding fold $x=\pm a$, but undefined at the other fold
$x=\mp a$. Despite the intersection of the slow nullcline $x=a$ with the cubic
nullcline, no equilibrium exists. Instead a turning point forms, which allows
a continuous passage through the corresponding fold without any obstruction.
This passage from the attracting to the repelling sheet of $S_{0}$ occurs in
finite time and is the basis of singular canard orbits. Therefore the turning
point is also referred to as canard point.

#### Folded-saddle canards.

The understanding of canards based on the VdP system can be extended towards
1-fast 2-slow systems. In the most simple case the parameter $a$ of the VdP
system is subject to a slow drift given by $\dot{a}=\varepsilon\mu$, with a
constant speed $\mu$. Naturally, the slow nullcline is absent in this case and
the intersection with the fast subsystem’s fold, which is forming the turning
point, can not occur. Nevertheless canards can be found and are a result of a
slow passage effect: the parameter $a$ dynamically transitions through the
canard explosion of the original problem.

In VdP the existence of turning points is conditioned by the fact that the
slow variable receives feedback from the fast variable , leading to a “turn”
of the slow flow direction, depending on the value of $x$. In the neural mass
model with STP in presence of external periodic forcing, however, this
feedback is absent. The turn of the slow flow appears naturally via the form
of the external forcing. In both cases the passage through the fold underlies
the same mechanisms and it is well understood by making use of an auxiliary
system, called the desingularized reduced system (DRS), which is introduced in
introduced in Section Desingularization.

Making use of the DRS, the formation of canards as described above, is
reflected by so-called folded-saddle singularities, which allow a passage from
the attracting to the repelling sheet of $S_{0}$. Ultimately, this leads to
the existence of canard orbits in the extended VdP system, as well as in the
neural mass with STP in presence of periodic forcing.

#### Torus and mixed-type canards.

The term canard is not restricted to dynamics taking place in the vicinity of
(or on) attracting and repelling manifolds, which represent equilibria. In
general, it refers to any type of solution evolving near attracting and
repelling invariant sets associated with the fast subsystem. These invariant
sets can correspond to equilibria but also to limit cycles. Following this
definition, a particular type of canard can be found in elliptic bursters
[45], which require at least a 2-fast 1-slow system. Here elliptic bursting
can arise due to a subcritical Hopf-Bifurcation (in the fast subsystem) giving
rise to unstable limit cycles, which stabilize via a fold bifurcation of
cycles. The Hopf bifurcation initiates the burst, while the fold of cycles
marks their termination. Usually in elliptic bursters the full dynamics follow
the family of stable limit cycles of the fast subsystem. However so-called
torus canards can be found for small enough $\varepsilon$ [46, 47]. They
describe orbits following a stable family of fast subsystem cycles and
switching to the unstable one past the fold.

A hybrid of classical canards and torus canards, so-called mixed-type canards,
were reported in [48]. They describe trajectories that spend time near
repelling branches of equilibria as well as limit cycles, and can therefore be
seen as a mix of classical canards and torus canards. Segments of these
solutions evolve nearby unstable equilibria and connect to unstable limit
cycles of the fast subsystem. These last types of canards are found in
bursting systems, which are ubiquitous in the modeling of neural activity at
both single-cell and population level.

##  4 Slow-fast dissection of the model

We will start a systematic investigation of the full system by dissecting it
into a slow and fast subsystem. The full problem represents a 4-fast 2-slow
system with $\mathbf{X}_{\mathrm{f}}=(r,v,x,u)$ and
$\mathbf{X}_{\mathrm{s}}=(I_{1},I_{2})$. Their dynamics is governed by the
right hand sides $\mathbf{F}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}})$
and $\mathbf{G}(\mathbf{X}_{\mathrm{s}})$, recalled below.

$\displaystyle\mathbf{F}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}})$
$\displaystyle=\begin{pmatrix}\frac{\Delta}{\tau_{\mathrm{m}}\pi}+2rv\\\
v^{2}+J\tau_{\mathrm{m}}uxr-(\pi\tau_{\mathrm{m}}r)^{2}+I_{1}\\\
(1-x)/\tau_{\mathrm{d}}-uxr\\\
(U_{0}-u)/\tau_{\mathrm{f}}+U_{0}(1-u)r\end{pmatrix}$ (12)
$\displaystyle\mathbf{G}(\mathbf{X}_{\mathrm{s}})$
$\displaystyle=\varepsilon\begin{pmatrix}g_{1}(I_{1},I_{2})\\\
g_{2}(I_{1},I_{2})\end{pmatrix}=\varepsilon\begin{pmatrix}\left(I_{1}(a-I_{1}^{2}-I_{2}^{2})+I_{2}\right)\\\
\left(I_{2}(a-I_{1}^{2}-I_{2}^{2})-I_{1}\right)\end{pmatrix}$ (13)

The equilibrium branches of the fast subsystem, shown in Fig. 2(a), are
defined via $\\{\mathbf{X}_{\mathrm{f}}|\dot{\mathbf{X}_{\mathrm{f}}}=0\\}$,
with the slow variable coordinates acting as bifurcation parameters. This
naturally coincides with the definition Eq. 14 of the critical manifold
$S_{0}$

$\displaystyle\mathbf{F}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}})$
$\displaystyle=\mathbf{0}$ (14)

We can therefore already infer the shape of $S_{0}$. It corresponds to the
cartersian product $S^{\ast}\times\\{I_{2}|I_{2}\in\mathbb{R}\\}$, where
$S^{\ast}$ denotes the S-shaped branch of equilibria of the fast subsystem and
is given in Eq. 15.

$\displaystyle
S^{\ast}:=\\{(r,v,x,u,I_{1})|\mathbf{F}(\mathbf{X}_{\mathrm{f}},\mathbf{X}_{\mathrm{s}})=0\\}$
(15)

Hence the associated fold set $\mathcal{F}$ has two 1D connected components,
namely the two lines $\mathcal{F}_{\rm U}:=\\{F_{\rm
U}\\}\times\\{I_{2}|I_{2}\in\mathbb{R}\\}$ and $\mathcal{F}_{\rm L}:=\\{F_{\rm
L}\\}\times\\{I_{2}|I_{2}\in\mathbb{R}\\}$. This means that the two fold lines
are parametrised by $I_{2}$. These two fold lines of $S_{0}$, along
$\mathcal{F}_{\rm L}$ and $\mathcal{F}_{\rm U}$, can be seen in Fig. 4.

Figure 4: Slow-fast dissection and critical manifold: Solutions
$\boldsymbol{\gamma}_{0}(t)$ to $\boldsymbol{\gamma}_{4}(t)$ of the full
system superimposed on the critical manifold $S_{0}$ in
$(I_{1},I_{2},r)$-space. The parameter values are identical to those in Fig.
3.

Points of $S_{0}$ are equilibria of the fast subsystem. This mean that, in
addition, local stability properties of the fast subsystem can be associated
with points of the critical manifold. Stable (unstable) parts of the
equilibrium branch in Fig. 2 will become attracting (repelling) sheets of the
critical manifold. This property can be seen as an indicator for fast flow in
the full system, distant from $S_{0}$. Here fast dynamics act and $(r,v,x,u)$
will be repelled and attracted accordingly. The stability changes along the
set of Hopf bifurcation given by $\mathcal{H}_{\rm L}:=\\{H_{\rm
L}\\}\times\\{I_{2}|I_{2}\in\mathbb{R}\\}$ and $\mathcal{H}_{\rm U}:=\\{H_{\rm
U}\\}\times\\{I_{2}|I_{2}\in\mathbb{R}\\}$ . Attracting (repelling) parts of
$S_{0}$ are marked as green (light green) surfaces throughout the present
work.

###  4.1 Singular dynamics: fast subsystem

The extended model Eq. 4 is able to generate periodic oscillations due to
plastic synapses even in absence of time dependent forcing
($I_{1}=\text{const.}$), as outlined in Section 2.2. Their existence depends
on the exact choice of parameters values, one of the important ones being the
total nonsynaptic current given by $\bar{\eta}+I_{1}$. Limit cycles can arise
via a plethora of bifurcation scenarios. Here we considered the case of a
subcritical Hopf bifurcation followed by a fold of limit cycles, giving rise
to stable oscillatory behaviour, when considering $I_{1}$ as a bifurcation
parameter (see Fig. 2). Since $S_{0}$ consists of fast subsystem equilibria,
the bifurcation diagram is essentially a projection of $S_{0}$ onto the
$I_{1}$-$r$ plane.

###  4.2 Singular dynamics: slow subsystem

#### Slow flow.

We make further use of the dissection by studying the slow flow on the
critical manifold $S_{0}$. In the slow subsystem Eq. 9 the state space is
reduced to $S_{0}$, described by four algebraic conditions in Eq. 14, the
solutions of which depend on the slow variable $I_{1}$ entering into the
membrane potential equation. The state variables in this limit are subject to
the slow flow
$(\mathbf{X}_{\mathrm{f}}^{\prime},\mathbf{X}_{\mathrm{s}}^{\prime})$
describing their dynamics on $S_{0}$. For
$\mathbf{X}_{\mathrm{s}}=(I_{1},I_{2})$ this is explicitly given via the Hopf
normal form Eq. 13. For the fast variables $\mathbf{X}_{\mathrm{f}}=(r,v,x,u)$
however, the algebraic constraints define $\mathbf{X}_{\mathrm{f}}$ as well as
$\mathbf{X}_{\mathrm{f}}^{\prime}$ on $S_{0}$ implicitly. In this case the
flow can be obtained by taking the total (slow) time derivative of Eq. 9a as
done in Eq. 16,

$\displaystyle
0=\derivative{\tau}\mathbf{F}(\mathbf{X}_{\mathrm{f}}(\tau),\mathbf{X}_{\mathrm{s}}(\tau))=\partialderivative{\mathbf{F}}{\mathbf{X}_{\mathrm{f}}}\derivative{\mathbf{X}_{\mathrm{f}}}{\tau}+\partialderivative{\mathbf{F}}{\mathbf{X}_{\mathrm{s}}}\derivative{\mathbf{X}_{\mathrm{s}}}{\tau}\quad,$
(16)

where $\partialderivative{(\cdot)}{\boldsymbol{a}}$ is the Jacobian of
$(\cdot)$ with respect to $\mathbf{a}$. If the Jacobian
$\partialderivative{\mathbf{F}}{\mathbf{X}_{\mathrm{f}}}$ is invertible, i.e.,
$\det(\partialderivative{\mathbf{F}}{\mathbf{X}_{\mathrm{f}}})\neq 0$, then
the slow flow of $\mathbf{X}_{\mathrm{f}}$ can be calculated and results in
Eq. 17.

$\displaystyle\derivative{\mathbf{X}_{\mathrm{f}}}{\tau}\equiv\mathbf{X}^{\prime}_{\mathrm{f}}=-\left(\partialderivative{\mathbf{F}}{\mathbf{X}_{\mathrm{f}}}\right)^{-1}\left(\partialderivative{\mathbf{F}}{\mathbf{X}_{\mathrm{s}}}\mathbf{X}^{\prime}_{\mathrm{s}}\right)$
(17)

This slow flow is only defined on $S_{0}$ and represents a system of ODEs
capturing the dynamics of $\mathbf{X}_{\mathrm{f}}$ and
$\mathbf{X}_{\mathrm{s}}$ on the manifold.

In the particular case of the neural mass we have
$\mathbf{X}_{\mathrm{f}}=(r,v,x,u)$, $\mathbf{X}_{\mathrm{s}}=(I_{1},I_{2})$
and $\mathbf{F}$, $\mathbf{G}$ given in Eqs. 12 and 13, respectively. Values
of $\mathbf{X}_{\mathrm{f}}$ on $S_{0}$ are determined via the algebraic
conditions Eq. 14, which entangle the fast variables to each other and to
$I_{1}$. Given either $r$, $v$, $x$ or $u$, the other components can be
calculated straightforwardly. In other words, it is sufficient to consider the
slow flow of one the fast variables, here $r$, to understand the slow
dynamics. Taken this into account we finally obtain the slow flow
$(r^{\prime},I_{1}^{\prime},I_{2}^{\prime})$ given in Eq. 18

$\displaystyle r^{\prime}$
$\displaystyle=g_{1}(I_{1},I_{2})Ar(1+ru\tau_{\mathrm{d}})(1+rU_{0}\tau_{\mathrm{f}})/D$
(18a) $\displaystyle I_{1}^{\prime}$ $\displaystyle=g_{1}(I_{1},I_{2})$ (18b)
$\displaystyle I_{2}^{\prime}$ $\displaystyle=g_{2}(I_{1},I_{2})$ (18c)

We can see the relevance of the denominator $D$, given in Eq. 19, by noting
that $D=0$ defines the singular points of $S_{0}$ at which the fast subsystem
Jacobian $\partialderivative{\mathbf{F}}{\mathbf{X}_{\mathrm{f}}}$ is
singular.

$\displaystyle D=2\left[(\pi
r)^{2}+v^{2}\right](r\tau_{\mathrm{d}}u+1)(r\tau_{\mathrm{f}}U_{0}+1)-Jxr(r\tau_{\mathrm{f}}U_{0}+u)$
(19)

Saddle-node bifurcations of the fast subsystem are characterized by the
identical condition
$\det(\partialderivative{\mathbf{F}}{\mathbf{X}_{\mathrm{f}}})=0$ and singular
points are equivalent to the fold curves $\mathcal{F}_{\rm L}$ and
$\mathcal{F}_{\rm U}$ shown in Fig. 4. Therefore, the slow flow is undefined
along these lines and the slow subsystem fails to describe the slow dynamics
for trajectories intersecting with $\mathcal{F}_{\rm L}$ or $\mathcal{F}_{\rm
U}$.

#### Desingularization.

This limitation can be mitigated by introducing an auxiliary system and
desingularizing Eq. 18 through the application of a nonlinear time rescaling
$\tau\mapsto D\cdot\tau$. One obtains the desingularized reduced system (DRS)
given in Eq. 20 with $\hat{\tau}:=D\tau$.

$\displaystyle\frac{\differential r}{d\hat{\tau}}$
$\displaystyle=g_{1}(I_{1},I_{2})Ar(1+ru\tau_{\mathrm{d}})(1+rU_{0}\tau_{\mathrm{f}})$
(20a) $\displaystyle\frac{\differential I_{1}}{d\hat{\tau}}$
$\displaystyle=g_{1}(I_{1},I_{2})\cdot D$ (20b)
$\displaystyle\frac{\differential I_{2}}{d\hat{\tau}}$
$\displaystyle=g_{2}(I_{1},I_{2})\cdot D$ (20c)

The DRS benefits from the fact that the singularities are resolved, allowing
to investigate the slow dynamics near and on the fold lines $\mathcal{F}_{\rm
L}$ and $\mathcal{F}_{\rm U}$. At the same time new equilibria are introduced
satisfying $D=0$. Additionally, as a consequence of the employed non-linear
time rescaling, the flow direction is not preserved. At the fold curves, with
$D=0$, a change of sign of $D$ takes place. Hence, between $\mathcal{F}_{\rm
L}$ and $\mathcal{F}_{\rm U}$, i.e, on the middle sheet of $S_{0}$, the flow
of the DRS is opposite to the slow flow.

Slow trajectories which entirely remain on the same sheet of $S_{0}$ can be
easily understood using the slow flow Eq. 18. Canard orbits evolve on
attracting as well as repelling sheets of $S_{0}$ and therefore require a view
on the dynamics near the folds. For this, we will determine the equilibria of
the DRS Eq. 20 in the following and analyze their invariant manifolds.
Equilibria of the DRS which satisfy $D=0$ necessarily coincide with a fold of
$S_{0}$. As we will show in the following section, this gives rise to so-
called folded singularities.

#### Folded saddle and folded homoclinics.

Up to three focus equilibria are located at
$p_{1}=(I_{1},I_{2},r)=(0,0,r_{k})$, where the $r_{k}$ are points of $S_{0}$
given $(I_{1},I_{2})=(0,0)$. In this work we will remain at $\bar{\eta}$
values for which only one equilibrium $p_{0}=(0,0,r_{0})$ exists. This point
is the only equilibrium of the slow flow Eq. 18. In the DRS (Eq. 20) however,
the condition $\\{g_{1}(I_{1},I_{2})=0,D=0\\}$ yields additional fixed points
$p_{1}$ and $p_{2}$ located on the fold lines at $D=0$.

The three equilibria of the DRS are displayed in Fig. 5(a) on $S_{0}$ in
$r-I_{2}$ projection. In the full system, for sufficiently small $A$ and
$\varepsilon$, solutions lie close to the bottom, attracting, sheet of
$S_{0}$. When the amplitude is increased, these cycles can pass very close to
$p_{1}$ and start to follow the middle, repelling sheet. One way to understand
this canard dynamics is to make use of the properties of the $p_{k}$ in the
DRS and their role for the slow subsystem.

Figure 5:  Folded-saddle and jump-on canards: (a) Critical manifold $S_{0}$ in
$r-I_{2}$ projection superimposed with the slow flow [green arrows, see Eq.
18]. The black dots $p_{0}$ (unstable focus), $p_{1}$ (saddle), $p_{2}$
(center) denote equilibria of the DRS [Eq. 20]. The point $p_{1}$ denotes a
folded-saddle equilibrium with the associated stable (unstable) eigendirection
indicated by a solid (dashed) arrow along the slow flow. The orange curves
$\mathcal{M}_{\rm FS}$ mark the stable and unstable manifolds of $p_{1}$,
forming heteroclinic connections through $p_{1}$. (b) $S_{0}$ in
$(I_{1},I_{2},r)$-space. The curves $\boldsymbol{\gamma}_{2}(t)$ and
$\boldsymbol{\gamma}_{3}(t)$ are solutions of the full system. The objects
$p_{1},p_{2}$ and $\mathcal{M}_{\rm FS}$ depend on the choice of $A$, here
they correspond to the value used to obtain $\boldsymbol{\gamma}_{3}(t)$.
Other parameters values are as in Figs. 3 and 4.

Located on the bottom sheet of $S_{0}$, $p_{0}$ results from the Hopf form
given in Eq. 5 and is an unstable focus at ($I_{1},I_{2})=(0,0)$. On the other
hand, $p_{2}$ lies on the upper fold line $\mathcal{F}_{\rm U}$ and denotes a
center, i.e., it has purely imaginary complex conjugate eigenvalues. The
equilibrium $p_{1}$ can be found on the lower fold line $\mathcal{F}_{\rm L}$
and is of saddle type. At a specific value of $A$, namely, when the forcing
cycle intersects with $p_{1}$, an 8-shaped double homoclinic connection
$\mathcal{M}_{\rm FS}=\mathcal{M}_{\rm FS}^{\rm L}\cup\mathcal{M}_{\rm
FS}^{\rm U}$ forms, consisting of two parts, which are connected via $p_{1}$;
see the orange curve in Fig. 5(a,b). The connection $\mathcal{M}_{\rm FS}^{\rm
L}$ is located on the lower sheet, while $\mathcal{M}_{\rm FS}^{\rm U}$ spans
the middle and upper sheet of $S_{0}$. They revolve around the unstable focus
$p_{0}$ and center $p_{2}$, respectively, and are the stable and unstable
manifolds of $p_{1}$.

The points $p_{2}$, $p_{1}$, and in particular the invariant manifolds
associated with $p_{1}$, play an important role for the slow subsystem. Due to
the negative sign $D<0$ on the middle sheet of $S_{0}$ the slow flow is
reversed with respect to the DRS. As a consequence the DRS saddle $p_{1}$ and
center $p_{2}$ become folded singularities of the slow subsystem. These folded
saddle ($p_{1}$) and folded center ($p_{2}$) are not equilibria of the slow
flow. However, for the slow dynamics they have similar impact on the dynamics
as there unfolded counterparts, but with the crucial difference of reversed
flow direction between $\mathcal{F}_{\rm L}$ and $\mathcal{F}_{\rm U}$.
Accordingly, the folded saddle $p_{1}$ has significant influence on the
dynamics of the slow subsystem along $\mathcal{M}_{\rm FS}$, as described in
the following.

(i) First of all, trajectories in the DRS evolving on $\mathcal{M}_{\rm FS}$
necessarily approach $p_{1}$ asymptotically from the direction of the stable
eigenvector, but can never pass through the saddle.

(ii) In the slow subsystem however, the folded saddle $p_{1}$ allows a pass-
through along this direction. Below $\mathcal{F}_{\rm L}$ trajectories are
attracted to and above repelled from $p_{1}$.

(iii) The double homoclinic connection $\mathcal{M}_{\rm FS}$ of the DRS is
referred to as a folded homoclinic in the slow subsystem. For this solution
the passage of trajectories through $p_{1}$ occurs in finite time [49].

(iv) Using the same type of argument, the invariant manifold $\mathcal{M}_{\rm
FS}^{\rm U}$ around the folded center $p_{2}$ becomes disconnected at the two
intersections with $\mathcal{F}_{\rm U}$, due to a reversal of the slow flow
direction. Solutions of the slow subsystem on $\mathcal{M}_{\rm FS}^{\rm U}$
can not cross this line.

As a consequence of the previous properties (i)-(iv), a singular canard exists
in the slow subsystem, given a specific value of $A$. It evolves along the
folded homoclinic $\mathcal{M}_{\rm FS}$ below the bottom fold line
$\mathcal{F}_{\rm L}$ and extends, while remaining on $S_{0}$, beyond the
folded-saddle $p_{1}$ until the upper fold $\mathcal{F}_{\rm U}$.

#### Singular canard orbits.

For the construction of singular orbits, we note once more that the middle
sheet of $S_{0}$ is repelling while the bottom sheet is attracting.
Accordingly, a continuum of fast segments emerging from the middle sheet and
connecting to the bottom sheet exist in the fast subsystem. This family of
fast orbits collides with $\mathcal{M}_{\rm FS}$ and likewise with the
singular canard described above. As a result, infinitely many singular orbits
can be constructed, by merging the singular canard at arbitrary positions on
the middle sheet of $S_{0}$, with fast segments.

These singular orbits evolve on the bottom sheet of $S_{0}$, continue through
the folded-saddle $p_{1}$, while following the folded homoclinic, and jump at
different heights, in terms of the coordinate $r$, from the middle to the
bottom sheet of $S_{0}$. The full system solution $\boldsymbol{\gamma}_{1}(t)$
in Figs. 3 and 4 displays this type of dynamics. The singular canard, hence
also the family of singular canard orbits, can at most reach the upper fold
line $F_{\rm U}$. Here the slow flow is undefined and the reduction of state
space to $S_{0}$ fails to describe the dynamics. This is additionally
reflected by the fact that $\mathcal{M}_{\rm FS}^{\rm U}$ is disconnected at
the intersections with $\mathcal{F}_{\rm U}$, due to the folded property of
$p_{2}$; see (iv) above. The singular canard orbit which reaches up until this
point is the maximal canard.

##  5 Full system dynamics: beyond singular orbits and classical canards

The solutions $\boldsymbol{\gamma}_{0}(t)$ to $\boldsymbol{\gamma}_{4}(t)$ are
results of numerical computations for $\varepsilon>0$. As such, their slow
segments evolve not on, but in an $\mathcal{O}(\varepsilon)$ neighborhood of
$S_{0}$. This is to some extent an implication of Fenichel’s theory [50]. For
$0<\varepsilon\ll 1$, it guarantees the existence of a slow manifold
$S_{\varepsilon}$, that is in an $\mathcal{O}(\varepsilon)$ neighborhood of
$S_{0}$, if $S_{0}$ is normally hyperbolic (see below). Additionally,
$S_{\varepsilon}$ is locally invariant under the flow of the full system. The
theorem also states that stable and unstable manifolds associated to $S_{0}$
persist as $\mathcal{O}(\varepsilon)$ perturbations. In other words, the flow
on $S_{\varepsilon}$ can be seen as a perturbation of the flow on $S_{0}$; and
the flow perpendicular to $S_{0}$ as a perturbation of the fast subsystem’s
flow. For normally hyperbolic critical manifolds, one can deduce that singular
orbits persist for $\varepsilon$ and perturb into an
$\mathcal{O}(\varepsilon)$ neighborhood.

Normal hyperbolicity requires all eigenvalues of the Jacobian
$\partialderivative{\mathbf{F}}{\mathbf{X}_{\mathrm{f}}}|_{S_{0}}$ to have
non-zero real part [51]. The theorem can therefore not be applied on
$\mathcal{F}_{\rm L}$ and $\mathcal{F}_{\rm U}$, given that they describe
lines of saddle-node bifurcations. However, one can consider the three sheets
of $S_{0}$ separately, each one up to an $\mathcal{O}(\varepsilon)$
neighborhood of the folds. From this we can conclude the persistence of slow
segments nearby $S_{0}$ for $\varepsilon>0$, including the repelling segments
within a canard solution, until close to the folds. Fenichel’s theory does not
encompass whether or not a connection of these segments exists. Yet,
$\boldsymbol{\gamma}_{1}(t)$ to $\boldsymbol{\gamma}_{3}(t)$ exemplify such
connection numerically: segments that evolve close to the bottom sheet of
$S_{0}$ connect to segments close the middle sheet. A detailed treatment of
these connections nearby the non-hyperbolic points, exceeds the scope of this
work. For more rigorous approaches, we refer to non-standard analysis [44],
matched asymptotics [52] and so-called blow-up techniques [53]. With these
advanced methods, it is possible to show that orbits with canard segments of
different length perturb at different parameter values within an exponentially
narrow regime, thus leading to the canard explosion in the full system.

###  5.1 Jump-on canards

The construction of singular canard orbits linked with Fenichel’s theorem
explains the dynamics of $\boldsymbol{\gamma}_{1}(t)$ in Figs. 3 and 4, where
the forcing amplitude $A$ is large enough to surpass the lower Hopf
bifurcation $H_{\rm L}$ and the lower fold $F_{\rm L}$. Headless canards, like
this one, have a jump to the bottom branch in common. Slow-fast systems may
also have canard solutions with a head. They usually appear in systems, which
have two folds: the first one destabilizing, the second one stabilizing the
branch. Headed canards jump onto this stable upper part. In 3D systems with an
1D or 2D S-shaped critical manifold, the upper sheet is typically unstable
near the upper fold, thus preventing the existence of headed canards. Instead,
past the maximal canard, fast oscillations related to the existence of limit
cycles develop and lead to bursting solutions; see Section 6.

The trajectory $\boldsymbol{\gamma}_{2}(t)$ in Figs. 3 and 4 has the
characteristic dynamics of a headed canard. However, it is peculiar for
various reasons. First of all, we can classify this type of dynamics as jump-
on canard, since the trajectory lands (after the fast jump) on a seemingly
repelling slow manifold. Moreoever, the jump-on dynamics can occur after a
regular canard segment, as for $\boldsymbol{\gamma}_{2}$ or independent of
that, like for $\boldsymbol{\gamma}_{3}$. As a matter of fact, trajectories of
the latter type resemble relaxation oscillations, like in VdP, despite the
repelling upper sheet of $S_{0}$. This has an additional consequence: for
small enough $\varepsilon$ a continuous transition from subthreshold
oscillations to bursting is blocked by jump-on canards. The fact that fast
oscillations other than relaxation oscillations are absent beyond the maximal
canard are novel and unexpected phenomena. In the following we will address
how these solutions emerge. Their impact on the route towards bursting is
discussed in details in Section 6.3.

The two solutions $\boldsymbol{\gamma}_{2}(t)$ and
$\boldsymbol{\gamma}_{3}(t)$ are shown in Fig. 5(b) and Fig. 6 using two
different projections. They exemplify two types of jump-on canards, that have
an approach towards globally repelling equilibria of the fast subsystem in
common. The cyan solution $\boldsymbol{\gamma}_{3}$ represents an orbit which
does not interact with the folded-saddle $p_{1}$. It slowly evolves on the
lower sheet of $S_{0}$, crosses the curve of Hopf bifurcations
$\mathcal{H}_{\rm L}$ and reaches the lower fold curve $\mathcal{F}_{\rm L}$
at which the slow subsystem is singular. Fast dynamics come into play and
expectedly the dynamics will approach attractors of the fast subsystem. Such
attractors for the considered $I_{1}$ value at which the curve escapes $S_{0}$
are solely the stable limit cycles displayed in Fig. 2(a). Instead of entering
a period of bursting, the cyan trajectories approaches the upper branch of
unstable equilibria. As soon as it jumps onto on $S_{0}$, the slow subsystem
becomes a valid limit anew and the curve remains on $S_{0}$ until it reaches
$\mathcal{F}_{\rm U}$, where it jumps down to the stable sheet.

In the second case, $\boldsymbol{\gamma}_{2}$ in Fig. 5(b,c), the orbit
possesses a canard segment and jumps from the repelling middle sheet to the
repelling upper sheet of $S_{0}$. Similar to the previous case, it evolves on
$S_{0}$ until $\mathcal{F}_{\rm U}$ and finally jumps down. The global motion
is identical to that of a canard with a head in VdP.

Multiple elements of these singular cycles have to be understood. First of
all, both cases have a similar slow segment in common, namely the part of the
trajectory on the upper sheet of $S_{0}$. They can be approximated by
solutions of the reduced problem Eq. 9 and are enforced by the presence of
$p_{2}$, around which the trajectory evolves. Since the center $p_{2}$ is
folded, full rotations around it are not possible and the slow parts terminate
at $\mathcal{F}_{\rm U}$, where the slow flow is undefined. Here the
trajectory can be joined to a fast bit which connects from $\mathcal{F}_{\rm
U}$ to the attracting sheet of $S_{0}$. After this part, the dynamics on
$S_{0}$ is again governed by the slow subsystem and depending on the forcing
amplitude $A$, the two orbits take different paths. The solution
$\boldsymbol{\gamma}_{3}$ crosses $\mathcal{F}_{\rm L}$ far from $p_{1}$ and
the slow segment on the attracting sheet stops; $\boldsymbol{\gamma}_{2}$
passes through a neighbourhood of $p_{1}$ and exhibits canard dynamics before
a jump occurs. These parts of the orbits are entirely described within the
scope of the slow subsystem.

###  5.2 Nested timescale separation

Understanding the remaining segment that leads $\boldsymbol{\gamma}_{2}$ and
$\boldsymbol{\gamma}_{3}$ towards the repelling sheet of $S_{0}$ requires a
more detailed analysis. The mechanism is the same for both cases and will be
discussed in the following. Since these pieces of the orbits evolve on the
fast timescale we present a visualization of the curves in $I_{1}-v-r$ space
in Fig. 6. In this projection the critical manifold $S_{0}$ is shown as a
green curve $r(I_{1}),v(I_{1})$ with attracting (repelling) parts as a solid
(dashed) line. In the singular limit, the jump-on points
$\boldsymbol{\gamma}^{\ast}_{2,3}=(r_{2,3}^{\ast},v_{2,3}^{\ast},x_{2,3}^{\ast},u_{2,3}^{\ast})$
of the jump-on canard solutions $\boldsymbol{\gamma}_{2}$ and
$\boldsymbol{\gamma}_{3}$, marked by the red and cyan dots in Fig. 6, are of
saddle-focus type. Linearization of the dynamics reveals a weakly and strongly
attracting direction in a neighbourhood of the jump-on points together with a
repelling direction with complex conjugate eigenvalues. This suggests a 2D
stable manifold leading to $\boldsymbol{\gamma}^{\ast}_{2,3}$. We will
simplify the problem further by noting that the entire dynamics of jump-on
canards takes place close to the surface $\mathcal{M}_{rv}$ defined in Eq. 21.
On the one hand, for the slow pieces of the curve this observation is as
expected, since by definition this conditions holds on $S_{0}$.

Figure 6:  Folded-saddle and jump-on canards: $S_{0}$ in $I_{1}$-$v$-$r$ space
(green curve). The blue surface $\mathcal{M}_{rv}$ is defined by
$r=-\frac{\Delta}{\pi\tau_{\mathrm{m}}v}$. The curves
$\boldsymbol{\gamma}_{2,3}(t)$ are solutions of the full system and their
jump-on points $\boldsymbol{\gamma}_{2,3}^{\ast}$ on the upper repelling sheet
of $S_{0}$ are marked by dots. Parameters values are as in Fig. 5.

$\displaystyle\dot{r}=0\Rightarrow
r=-\frac{\Delta}{2\pi\tau_{\mathrm{m}}v}\text{\quad for }v\neq 0$ (21)

On the other hand, also the fast parts of the orbits remain on
$\mathcal{M}_{rv}$. This implies a reduction of the fast dynamics to
$\mathcal{M}_{rv}$ for the part of state space in which jump-on canards can be
found. We will exploit this reduction and investigate a secondary
differential-algebraic system resulting from the fast subsystem Eq. 8 in which
the flow $\dot{r}$ is equilibrated.

This latter step, without having a complete picture of the time scaling, is
analogous to an additional dissection of the fast subsystem. In other words,
the full system exhibits three timescales for what concerns jump-on canards:
the dynamics of $r$ takes place on a fast, that of $(v,x,u)$ on an
intermediate and that of $(I_{1},I_{2})$ on a slow timescale. The
equilibration $\dot{r}=0$ eliminates the fastest of these scales and
approximates the intermediate scale dynamics of jump-on canards, which takes
place in the vicinity of $\mathcal{M}_{rv}$.

In this new framework $\mathcal{M}_{rv}$ describes a manifold on which the
dynamics of $(v,x,u)$ take place, while $I_{1},I_{2}$ remain frozen. The jump-
on points $(v^{\ast},x^{\ast},u^{\ast})$ [red and cyan dots in Fig. 6] are
mutual points of $\mathcal{M}_{rv}$ and the upper sheet of $S_{0}$. In the
reduced problem on $\mathcal{M}_{rv}$, $(v^{\ast},x^{\ast},u^{\ast})$ are of
saddle type with eigenvalues $\lambda_{1}\gg\lambda_{2}\gg-\lambda_{3}>0$ and
therefore have associated 1D stable manifolds $\mathcal{M}_{\rm
JO}\subset\mathcal{M}_{rv}$ [orange curves in Fig. 6]. They exists for any
value of $I_{1}$ after the upper fold. For values of $I_{1}$ beyond the lower
fold, they extend down to $r=0$. The family of 1D stable manifolds associated
with the jump-on points can guide trajectories towards the upper repelling
sheet of $S_{0}$ and this way leads to the existence of jump-on canards.

So far we have discussed singular orbits, for which the fast segments connect
different sheets of the critical manifold $S_{0}$. In the case of regular
canards, for $\varepsilon>0$ small enough, these correspond to stable
equilibria of the fast subsystem. For jump-on canards they might be unstable,
but possess a stable direction, allowing to reach and stay on the repelling
sheet of $S_{0}$. The main dynamics of these singular orbits takes place on
$S_{0}$ and does not display phases of fast oscillations, as for bursting
solutions, but only single fast jumps.

##  6 Slow-fast transition to bursting: a tale of two routes

Opposed to that, the solution $\boldsymbol{\gamma}_{4}(t)$ [see Figs. 3 and 4]
and the case initially illustrated in Fig. 1 exhibit bursts: a slow segment is
followed by fast oscillations. Periodic solutions of the fast subsystem are
the underpinning elements of bursting, such that a classification in terms of
the fast subsystem’s bifurcations appears appropriate. As shown in the
bifurcation diagram Fig. 2(a), limit cycles originate and terminate at Hopf
bifurcations, and change stability at a fold of cycles. Strictly following the
classification of Izhikevich, bursting solutions in this system are of
subcritical Hopf/fold cycle type. However, the subcritical Hopf bifurcation is
closely followed by a fold of the underlying equilibrium branch. Due to a
delay effect when surpassing the subcritical Hopf bifurcation [54, 55, 56],
bursting can effectively be initiated at that fold; see e.g. Fig. 7(b4). We
will restrict our analysis to these cases. Here a more aptly description of
the bursting type is fold/fold cycle, which corresponds to elliptic bursting
in the classification of Rinzel [33].

In order to understand bursting solutions of the full system, we want to
remain in the slow-fast dissection. However, the neural mass with STP in
presence of periodic forcing turns out to be a peculiar system and numerically
difficult to handle. We are constrained by two main factors. First of all, as
soon as we leave the singular limit, i.e, for $\varepsilon>0$, slow segments
of trajectories diverge sensitively from the critical manifold. In other
words, $\varepsilon$ is required to be remarkably small to maintain a good
agreement between between full system trajectories and singular orbits.
Secondly, numerical simulation as well as numerical continuation of the full
system for small enough $\varepsilon$ are challenging, since the dynamics
appears to be stiff and require high accuracy.

Therefore a clear view on the emergence of bursting can not be gained easily
in this framework. For this reason we will provide, additionally to
geometrical arguments, numerical evidence on how bursting forms in the present
system, either via direct simulation or continuation using the full system. In
general, bursts might emerge via a spike-adding mechanism, that is, the
consecutive addition of spikes into the orbit, when varying a parameter (e.g.,
the forcing amplitude). For parabolic bursting this spike-adding is mediated
by folded-saddle canards [49]; in square-wave bursters on the other hand,
passages through a fast-subsystem saddle-homoclinic bifurcation and a folded
node determine the number of large-amplitude oscillations in the burst and
small-amplitude oscillations before the burst, respectively [57]. In the
following, we report the spike-adding mechanism for the NMSTP. At its basis is
an interaction of the canards dynamics invoked by the presence of the folded
saddle $p_{1}$, as well as unexpected torus-canard dynamics. Moreover, we will
point out the role of jump-on canards for this spike-adding transition.

###  6.1 Canard explosion and spike-adding

Figure 7: Emergence of bursting: (a) Bifurcation diagram of the full system.
The black curve shows the $L_{2}$-norm of a family of periodic solutions vs.
the forcing amplitude $A$. The blue dots show the number of spikes
$n_{\mathrm{s}}$, defined as the local maxima of for which
$r(t)>$0.21\text{\,}\tau_{\mathrm{m}}^{-1}$$. The dashed vertical line is
located at $A=A^{\ast}=0.25531851205$. (b) Solutions
$(r(t),v(t),x(t),v(t),I_{1}(t),I_{2}(t))$ in $r-I$ projection superimposed on
the bifurcation diagram of the fast subsystem. The insets show the solution in
time. (c - d) Same solutions as in (b) and critical manifold $S_{0}$ in
$(I_{1},I_{2},r)$-space (c) and in $(I_{1},v,x)$-space (d). In (d) attracting
(repelling) sheets of $S_{0}$ are visualized as a green solid (dotted) line;
the purple surface (wireframe) represents the family of stable (unstable)
limit cycles of the fast subsystem [purple branch in column (b)]. Note that
here the $x$-axis is inverted. In (b) the black dots denote $\mathcal{F}_{\rm
L}$ and $\mathcal{F}_{\rm U}$, the orange dot $\mathcal{H}_{\rm L}$, while the
black dots in (c) show the folded singularity $p_{1}$, assuming $A=A^{\ast}$.
Spikes contributing to $n_{\mathrm{s}}$ are marked with blue dots. The dashed
red curve shows the solution of the panels above. The $A$ values in (b - d)
are in increasing order from top to bottom and near $A^{\ast}$. Full system
solutions obtained at $\varepsilon=10^{-3}$.

To start with, we consider the case $\varepsilon=10^{-3}$ and investigate the
full system dynamics by performing continuation with the forcing amplitude $A$
as a parameter. We want to stress that this $\varepsilon$ value, although very
small, proves to be rather distant from the singular limit and slow-fast
dissection arguments have to be taken with caution. The initial solution is
for an $A$ value corresponding to subthreshold oscillations, like
$\boldsymbol{\gamma}_{0}(t)$ in Figs. 3 and 4, and is continued towards larger
amplitudes.

As a solution measure the $L_{2}$-norm of this family is plotted vs. $A$ in
Fig. 7(a). The first part until $A=A^{\ast}\approx 0.25531851205$ is in the
subthreshold regime. Around $A^{\ast}$ a very sharp transition occurs,
resembling a canard explosion. In this transition region the orbits already
exhibit first spikes, here defined as the number $n_{\mathrm{s}}$ of local
maxima of $r(t)$ for which $r(t)>0.21$. This is followed up by a series of
arches (on the solution branch) at $A>A^{\ast}$.

The arches are clearly related to the addition of new spikes to the orbit:
with every termination of an arch, by means of a vertical dip of the curve,
the number of spikes $n_{\mathrm{s}}$ increases. This behaviour can be better
observed for larger $A$, for which each arch is related to the adding of
exactly one spike. Prior, the arches lay more dense along $A$ and the spike
adding appears to be of more complex nature. Despite the fact that the points
of $n_{\mathrm{s}}$ vs. $A$ depend on the choice of the $r$-threshold for
which a spike is counted, it is evident that spikes are added consecutively.
It is also clear that these bursting solutions emerge, in a continuous manner,
from subthreshold oscillations.

A more detailed view on the full system dynamics near the explosive transition
around $A=A^{\ast}$ is given in the columns (b - d) of Fig. 7. Column (b)
shows the solution in $r-I_{1}$ projection superimposed on the fast
subsystem’s bifurcation diagram; column (c) in $(I_{1},I_{2},r)$-space
together with $S_{0}$; column (d) in $(I_{1},v,x)$-space.

In projection onto the ($I_{1},v,x$)-space, the critical manifold appears as a
curve ($I_{1},v(I_{1}),x(I_{1}))$). Additionally we show the family of fast
subsystem limit cycles, which emerge at the lower subcritical Hopf bifurcation
$H_{\rm L}$. Unstable periodic solutions of this branch will be denoted by
$\Gamma^{\rm r}$, stable ones after the fold of cycles by $\Gamma^{\rm a}$.
Embedding these solutions $\Gamma^{\rm a,r}(I_{1},t)=(r,v,x,u)(I_{1},t)$ into
the state space of the full system one obtains the surface
$\mathcal{P}=\mathcal{P}^{\rm a}\cup\mathcal{P}^{\rm r}$, which consists of
attracting and repelling parts $\mathcal{P}^{\rm a,r}=\\{(\Gamma^{\rm
a,r}(I_{1},t),I_{1})|t\in[0,T(I_{1})]\\}\times\\{I_{2}|I_{2}\in\mathbb{R}\\}$,
corresponding to stable and unstable branches of the solution family,
respectively. The period $T(I_{1})$ of these cycles depends on $I_{1}$.

#### Onset of fast oscillations $\boldsymbol{({\rm b}_{1},{\rm c}_{1},{\rm
d}_{1})}$.

At the smallest of the chosen $A$ values near $A^{\ast}$ one can already
observe fast oscillations, consisting of five not fully developed spikes. They
occur after the trajectory has turned around the folded saddle $p_{1}$, marked
by a black dot in $({\rm c}_{1})$. It is this motion around $p_{1}$, taking
place in the vicinity of the critical manifold $S_{0}$, which has signs of a
turning point, guiding the trajectory along the repelling sheet of $S_{0}$.
Taking the rather large $\varepsilon$ value into account, this turn hints at a
canard segment arising due to the presence of the folded saddle $p_{1}$. After
this segment, in $(I_{1},I_{2},r)$-space (panel ${\rm c}_{1}$), the trajectory
pierces through the repelling sheet of $S_{0}$ and fast oscillations set in.
These results suggest that bursting is initiated at the termination of a
canard segment, closely following the repelling middle sheet of $S_{0}$.

Taking the fast subsystem LC family $\Gamma^{\rm a,r}$ into account, a
remarkable feature of the dynamics can be seen in $(I_{1},v,d)$ space (panel
(d1)). The spikes of the burst appear to follow the family of unstable limit
cycles, thus evolving near the repelling surface $\mathcal{P}^{\rm r}$.
Unexpectedly, after piercing through the critical manifold, the bursting
solution stays in the proximity of $\mathcal{P}^{\rm r}$, instead of being
repelled from it. As $I_{1}$ slowly drifts towards smaller values, the
trajectory remains close to $\mathcal{P}^{\rm r}$. These windings around
$\mathcal{P}^{\rm r}$ correspond to the first not fully grown spikes of the
full system solution. The spikes increase in amplitude, as $I_{1}$ decreases,
but remain small; see panel (b1) and inset. An enlargement of two full system
trajectories in $(I_{1},v,x)$-space is shown in Fig. 8, with two exponentially
close $A$ values near $A^{\ast}$.

Finally, the fast oscillations terminate via an escape from $\mathcal{P}^{\rm
r}$. By approaching the bottom sheet of $S_{0}$, the dynamics change to that
of a drifting equilibrium, passing from burst to quiescence.

Figure 8: Emergence of bursting: Enlargement of Fig. 7(${\rm d}_{2}$).
Critical manifold in $(I_{1},v,x)$-space with attracting (repelling) sheets of
$S_{0}$ visualized as a green solid (dotted) line; the purple surface
(wireframe) represents the family of stable (unstable) limit cycles of the
fast subsystem [purple branch in column (b)], emerging from the lower Hopf
bifurcation $H_{\rm L}$ (orange dot). The blue and red curves show solutions
of the full system near the canard explosion. Forcing amplitude $A$ for the
blue trajectory is larger than for red and both exponentially close the canard
explosion at $A^{\ast}$. Solutions obtained for $\varepsilon=10^{-3}$.

#### Explosivity and spike-adding $\boldsymbol{({\rm b}_{2},{\rm c}_{2},{\rm
d}_{2})}$.

At a slightly larger $A$ value, close to $A^{\ast}$, a majority of the
trajectory remains essentially unchanged, with respect to the previous $A$.
The slowly drifting part along the bottom equilibrium branch of the fast
subsystem and the canard segment, as well as the first oscillations, appear
frozen. This is clear by comparing the blue trajectory with the red dashed
curves in panels (${\rm b}_{2},{\rm c}_{2},{\rm d}_{2}$), as well as in Fig.
8. The fact that part of the trajectory near the fold freezes, while the
following part changes significantly, is a strong indication for explosivity
of the solution, when varying the parameter. This strong sensitivity towards
parameter changes is typical for canard dynamics. It is caused by the presence
of repelling objects in the fast subsystem, typically, but not exclusively,
equilibria, like the middle sheet of $S_{0}$.

Indeed, in panels (${\rm b}_{2},{\rm c}_{2},{\rm d}_{2}$), the full system
trajectory possesses a canard segment staying near the middle sheet of
$S_{0}$, when it turns around the folded saddle $p_{1}$. Hence the sensitivity
is to be expected. Moreover, during the fast oscillations, the bursting
solution evolves close to $\mathcal{P}^{\rm r}$, thus adding an additional
layer of sensitivity.

Compared to the previous case in Fig. 7(${\rm b}_{1},{\rm c}_{1},{\rm
d}_{1}$), the full system solution winds around $\mathcal{P}^{\rm r}$ more
often until smaller values of $I_{1}$, before jumping back to $S_{0}$. The
trajectory essentially remains close to $\mathcal{P}^{\rm r}$, but reaches up
higher. This way, by passing through $A^{\ast}$, more and more spikes with
increasing amplitude are added to the burst. These spikes are yet not fully
grown to the amplitude of the stable limit cycles present in this $I_{1}$
region. Furthermore, we note that distance of the blue trajectory to
$\mathcal{P}^{\rm r}$ increases, as it winds around it, indicating some extent
of repulsion near the surface; see also panel (b2). This way, the full
dynamics starts to escape from $\mathcal{P}^{\rm r}$ and gets attracted to
$\mathcal{P}^{\rm a}$.

#### Emergence of bursting $\boldsymbol{({\rm b}_{3},{\rm c}_{3},{\rm
d}_{3})}$.

As $A$ increases further, the point at which the trajectory starts to escape
from $\mathcal{P}^{\rm r}$ shifts towards the lower Hopf bifurcation $H_{\rm
L}$. In fact, the last two spikes are already repelled sufficiently to evolve
close to the attracting surface $\mathcal{P}^{\rm a}$; see panels (b3,d3). In
other words, the number of revolutions near $\mathcal{P}^{\rm r}$ reduces,
while the ones around $\mathcal{P}^{\rm a}$ increases. These oscillations near
$\mathcal{P}^{\rm a}$ are of large amplitude and mark the start of a dense
burst following for larger $A$.

#### Bursting $\boldsymbol{({\rm b}_{4},{\rm c}_{4},{\rm d}_{4})}$.

At the next step in panels ${\rm b}_{4}$ \- ${\rm d}_{4}$, most of the
windings around $\mathcal{P}^{\rm r}$ have vanished and the fast oscillations
take place in the proximity of $\mathcal{P}^{\rm a}$. Additionally, the burst
consists of more spikes in total, with respect to the first considered $A$
value (panels (${\rm b}_{1},{\rm c}_{1},{\rm d}_{1}$)).

Here, an additional spike-adding mechanism beyond the critical value
$A^{\ast}$ acts and is related to the period $T(I_{1})$ of $\mathcal{P}^{\rm
a}$ as a function of $I_{1}$. When $A$ is beyond the canard explosion at
$A^{\ast}$ the full system dynamics approach $\mathcal{P}^{\rm a}$ right after
surpassing the lower fold $F_{\rm L}$, without the excursion on
$\mathcal{P}^{\rm r}$. As the amplitude $A$ is increased, this attraction to
$\mathcal{P}^{\rm a}$ occurs at larger $I_{1}$ values. The period $T(I_{1}$)
of the fast subsystem limit cycles decreases with increasing $I_{1}$ and this
finally leads to more windings around $\mathcal{P}^{\rm a}$. On top of this,
the fraction of time for which the trajectory remains near $\mathcal{P}^{\rm
a}$ increases. Both effects add more spikes to the burst and result in the
spike-adding arches observed in Fig. 7(a).

The presented results already show the complexity of how bursts are generated,
that is, via a transition through the canard explosion at $A=A^{\ast}$, which
rather surprisingly leads the canard segment to evolve around the repelling
object $\mathcal{P}^{\rm r}$. This is followed by spike growth via repulsion
from $\mathcal{P}^{\rm r}$ and attraction to $\mathcal{P}^{\rm a}$, until all
oscillations evolve near $\mathcal{P}^{\rm a}$.

###  6.2 Continuous route to bursting: spike-adding via mixed-type-like torus
canards

Before the bursting transition, the full system dynamics can be described by a
single slow frequency, determined by $\varepsilon$. After the transition, a
full cycle consists of a slow phase followed by fast spiking. It is therefore
characterized by the slow frequency and a fast one, given by properties of the
fast subsystem. In the context of bursting and in slow-fast systems, whose
fast subsystem has both stable and unstable cycles, this change of the
dynamics, from one to two frequencies, hints at a torus (Neimark-Sacker)
bifurcation in the full system.

This can indicate the existence of mixed type-torus canards (MTTCs). Indeed,
the full system dynamics nearby the canard explosion not only exhibits a
canard segment along the repelling sheet of $S_{0}$, but as well a canard
segment on the repelling higher dimensional invariant set $\mathcal{P}^{\rm
r}$. This clearly resembles mixed-type canards as described in [48]. Very
similar MTTCs have been reported in [45], where they evidently arise in the
singular limit $\varepsilon\rightarrow 0$. In particular Fig. 2 of [45]
reports dynamics where a quasi-static motion of the full system along the
attracting sheet of $S_{0}$ connects to a repelling set of limit cycles
created by a subcritical Hopf bifurcation. In the NMSTP however, the
understanding of mixed-type torus canards is more complex for various reasons
and as we will show, only observed for small, but large enough $\varepsilon$.

First of all, in the $\varepsilon$ regime for which we observe mixed-type
torus canards, the timescale separation is small enough for the canard segment
on the middle sheet of $S_{0}$ to persist as a strongly perturbed version of
its singular counterpart. One can observe a turn around the lower fold
$\mathcal{F}_{\rm U}$, mediated by the folded-saddle singularity $p_{1}$. It
forces the trajectory to pierce through $S_{0}$, bringing it very close to
$\mathcal{P}^{\rm r}$; see Fig. 8.

Secondly, the solutions $\Gamma^{r}\subset\mathcal{P}^{\rm r}$, despite being
globally repelling, possess two stable Floquet multipliers. We argue that the
associated stable directions, similar to the jump-on canard case (see Section
5.2), form due to the intrinsic timescales of the fast subsystem.
Consequently, this means that stable manifolds can be associated with
$\mathcal{P}^{\rm r}$: it can attract trajectories along certain directions.

Thirdly, for large enough $\varepsilon$, one can expect dynamics very distinct
from singular orbits. In particular, solutions may not only evolve nearby, but
also switch between different attracting branches of the fast subsystem. We
observe this transition from the middle sheet of $S_{0}$ to $\mathcal{P}^{\rm
r}$. Despite both being repulsive, the reasoning holds: the folded-saddle
canard dynamics enforces the full system to stay close to the middle sheet of
$S_{0}$; stable directions allow an attraction towards $\mathcal{P}^{\rm r}$;
the large $\varepsilon$ allows to bridge $S_{0}$ to $\mathcal{P}^{\rm r}$ and
finally torus canard dynamics allow to follow $\mathcal{P}^{\rm a}$ closely,
adding more and more spikes in the transition region around $A^{\ast}$.

###  6.3 Discontinuous route to bursting: block evoked by jump-on canards

The transition from subthreshold oscillations to bursting when increasing the
forcing amplitude $A$ explains the emergence of the very first spikes in the
burst, which occur for forcing amplitudes exponentially close to $A^{\ast}$.
They also show how the subsequent spike-adding process, reflected by the
arches in Fig. 7(a), occurs. In the following we will extend the analysis of
this transition, taking into account different values of the parameter
$\varepsilon$. As an outlook, we will describe how, for $\varepsilon$ values
closer to the singular limit, jump-on canards interfere and block the
transition.

Figure 9: Families of bursting solutions and canards: The bifurcation diagram
shows solution families of the full system in terms of the $L_{2}$-norm vs.
the shifted forcing amplitude $A-A^{\ast}$. For all branches (br.) apart from
br. 4, the $A^{\ast}$ denotes the location of the canard explosion. For br. 4
on the other hand it marks the termination of continuation due to insufficient
accuracy. There are two types of solution families: the ones which undergo a
continuous transition from subthreshold oscillations to bursting (red curves)
and those which transition from subthreshold to jump-on canard dynamics (cyan
curves). Br. 4 and 5 have identical $\varepsilon\approx 2.667521298\cdot
10^{-4}$.

In Fig. 9 solution families of periodic orbits are displayed, obtained via
numerical continuation of the full system Eqs. 4 and 5. The figure comprises
seven branches for $\varepsilon$ values ranging from $\varepsilon=5\cdot
10^{-3}$ to $\varepsilon=1\cdot 10^{-5}$ and they are aligned to the canard
explosion at $A=A^{\ast}$ (apart from branch 4). Branch 3 is identical to the
one shown in Fig. 7 and undergoes a continuous transition from subthreshold
oscillations to bursting.

As a general result, we find two types of solution families, which, beyond the
canard explosion, i.e, for $A>A^{\ast}$, take different paths. For
$\varepsilon\gtrapprox 2.7\cdot 10^{-4}$ the branches display a continuous
transition from subthreshold oscillations to bursting; see red curves in Fig.
9. For $\varepsilon\lessapprox 2.7\cdot 10^{-4}$ on the other hand the
branches evolve into families of jump-on canards (curves in Fig. 9). A value
of $\varepsilon\approx 2.667521298\cdot 10^{-4}$ separates the two
$\varepsilon$ regimes.

However, it is clear by considering branches 4 and 5 that bursting solutions
do not cease to exist for smaller $\varepsilon$ values. Instead they coexist
with the jump-on type solutions.

It remains unclear how bursting forms for $\varepsilon\lessapprox 2.7\cdot
10^{-4}$. Nevertheless the role of jump-on canards for the emergence of
bursting becomes evident: for singular or small enough $\varepsilon\lessapprox
2.7\cdot 10^{-4}$ regular canards are observed, evolving on the middle sheet
of $S_{0}$. Beyond the maximal canard, the intrinsic timescales of the fast
subsystem come into play and lead to the emergence of jump-on canards (see
Section 5.2). They block a transition towards $\mathcal{P}$ and as a
consequence, bursting remains absent for these solutions families.

For $\varepsilon\gtrapprox 2.7\cdot 10^{-4}$ however, jump-on canards cease to
exist. In the amplitude regime where they would be expected, the system
approaches $\mathcal{P}^{\rm a}$ and bursts instead. The region around the
canard explosion is populated by MTTCs and separates subthreshold oscillations
from the bursting regime. It remains an open question for future work how the
differentiation depending on $\varepsilon$ occurs and in particular which
possible bifurcations of the full dynamics result in the distinct regimes of
continuous bursting transition and blocking jump-on canards.

##  7 Network behaviour

The neural mass model with STP Eq. 4 is an exact limit of the underlying QIF
network Eq. 3 as $N\rightarrow\infty$. We want to emphasize the benefits of
neural mass models and their capability of describing neuronal dynamics at a
macroscopic scale. In Fig. 2(b-f) results using the fast subsystem and the
corresponding QIF network have been shown. For the original neural mass model
[5] the exactness of the meanfield limit has been exploited in various studies
in order to understand the collective dynamics of large neuronal populations
[6, 58, 7, 59, 8, 20]. We want to note that STP, as opposed to exponential
synapses used in Ref. [13], results in a substantially higher sensitivity
towards finite size fluctuations and numerical errors, rendering the agreement
of network and neural mass less clear, in particular in the canard regime.

Here we want to assess if the mechanisms leading to bursting in the neural
mass model persists in a finite-sized network. To our knowledge such analysis,
in particular in presence of short-term synaptic plasticity, has not been
performed. For this we introduce the external periodic forcing
$I_{1}=A\sin(\varepsilon t)$ also into the network Eq. 3 and investigate the
QIF population dynamics nearby the canard explosion of the neural mass.

Figure 10: Bursting in the QIF network: (a) The bifurcation diagram shows a
family of periodic orbits of the full neural mass system in terms of the
$L_{2}$-norm vs. the forcing amplitude $A-A^{\ast}$. It corresponds to br. 2
in Fig. 9. The black dashed lines mark $A=0.265$ and $A=0.270$ respectively.
(b - c) Trajectories of the full system (red) superimposed on results obtained
by simulating the QIF network Eq. 3 with $N=100000$ neurons and in presence of
sinusoidal forcing $I_{1}(t)=A\sin(\varepsilon t)$. In (b1,c1) the forcing
amplitude is given by $A=0.265$; in (b2,c2) by $A=0.270$. Row (b) shows the
time series of the solution in terms of $x(t)$ vs. time $t$. In row (c) they
are shown in $(I_{1},v,x)$-space together with the critical manifold $S_{0}$,
with attracting (repelling) sheets as a solid (dashed) green line, and the
invariant manifold $\mathcal{P}$ of the fast subsystem (purple wireframe and
surface). The orange dot denotes $\mathcal{H}_{\rm L}$, the black dots
$\mathcal{F}_{\rm L}$ and $\mathcal{F}_{\rm U}$. The $x$-axis is inverted in
(c).

In Fig. 10(a) the family of periodic solutions transitioning from subthreshold
oscillation to bursting is shown for $\varepsilon=2\cdot 10^{-3}$. Two dashed
black lines mark the values $A=0.265$ and $A=0.27$, respectively, for which
neural mass (red) as well as QIF network trajectories (black) are depicted in
Fig. 10(b,c). Row (b) shows the time series $x(t)$ vs. $t$, row (c) the
trajectories in $(I_{1},v,x)$-space.

The QIF network consists of $N=100000$ neurons and $\varepsilon=2\cdot
10^{-3}$ is chosen to maintain reasonable computation times. In the network,
the initial conditions are chosen according to fixed point values
$(r^{\ast},v^{\ast},x^{\ast},u^{\ast})$ obtained from the neural mass at
$I_{1}=0$. However, initializing the network firing rate and average membrane
potential at a given value is a non-trivial problem. To overcome this, the
initial conditions of the QIF network are set to $V_{i}=v^{\ast}\quad\text{for
all }i=1,\dots,N$, $x=x^{\ast}$ and $u=u^{\ast}$. After a short transient, at
time $t=0$, the network has reached equilibrium and the forcing
$I=A\sin(\varepsilon t)$ sets in.

We start the analysis by considering $A=0.265$, marked by the left most dashed
line in Fig. 10(a). It lays within the spike-adding arches of the neural mass,
which displays bursting with a multitude of fully grown spikes; see red curves
in Fig. 10(b1, c1). They mainly evolve in the proximity of the invariant set
$\mathcal{P}^{\rm a}$. This periodic solution of the full problem does not
highlight a canard segment around the folded singularity $p_{1}$, which is to
be expected, since $A$ is not close enough to the canard explosion at
$A^{\ast}$.

On the other hand, the macroscopic state of the QIF network is found to be
just at the start of the spike-adding process. The orbit possesses a canard
segment: the turn around the lower fold $F_{\rm L}$ is evident. This is clear
by looking at the black curve in Fig. 10(c1). More strikingly, a few low-
amplitude oscillations are picked up in the vicinity of the repelling
invariant set $\mathcal{P}^{\rm r}$. This means that, despite the
discrepancies, the mixed-type torus canards observed in the neural mass are
also found in the QIF network. It is plausible to assume that spike-adding in
the network follows the same mechanism as described in Section 6.

At the larger amplitude $A=0.270$, shown in Fig. 10(b2, c2), the agreement
between neural mass and network trajectory is more distinct. In both systems
the canard segment is absent and they both exhibit bursting with a comparable
number of large amplitude spikes. Taking also into account the previous case
$A=0.265$, it appears that the bifurcation structure of the network is shifted
towards larger amplitudes with respect to the neural mass.

##  8 Discussion

The results of this work lay within the intersection of various fields
including meanfield theory and slow-fast dynamics. First of all, we make use
of recent developments of meanfield theory, namely the powerful OA ansatz, in
order to understand the emergent dynamics of spiking neuronal networks on
macroscopic scale.

Secondly, we aim at more biologically plausible models in this work, through
the inclusion of synaptic dynamics, in form of STP, while preserving the
exactness of the meanfield limit with respect to the QIF network. This adds
more realism to the QIF network and neural mass, but at the same time adds to
the complexity of the collective dynamics, even in absence of external
forcing. In the part of parameter space chosen here, STP allows for the
existence limit cycles, that are absent in the original model. As we have
shown via extensive numerical evidence, STP leads to various peculiarities,
when forcing the system slowly and periodically.

That leads to the third point, which we put our focus on: the treatment of the
forced neural mass with STP using methods of singular perturbation theory, in
particular slow-fast dissection. Without STP, slow external forcing already
gives rise to bursts, as shown in [5]. The fast oscillations of these orbits
vanish in the singular limit and relaxation oscillations remain. Accounting
for STP, however, leads to more intricate dynamics.

One of the fundamental elements are canards. Expectedly, due to the slow
harmonic passage through a fold of the S-shaped critical manifold, they appear
as folded-saddle canards, which play a role for spike-adding in parabolic
bursters. Here, despite the fact that the equilibrium destabilizing
bifurcation is a subcritical Hopf bifurcation, the observed bursts are
reminiscent of elliptic bursting (i.e., subcritical fold/fold of cycles
bursting), due the slow passage through the Hopf bifurcation.

Concerning the full system dynamics, a folded saddle can be found, with
associated explosive canard solutions. It separates quiescent orbits with
purely slow dynamics from bursting ones. In this transition region we find an
intriguing interplay of slow-fast effects: jump-on canards exist close enough
to the singular limit and are associated with a subtle timescale separation of
the neural mass, allegedly invoked by STP. They connect two repelling sheets
of the critical manifold and more strikingly, block a continuous transition
from quiescence to bursting.

Jump-on canards are one of the surprising elements of this work. However, they
vanish when considering biologically more plausible frequencies of the
periodic external current. Despite the rather insufficient timescale
separation in this parameter regime, we find orbits which display a strongly
perturbed canard segment. It is remarkable that these canards, as opposed to
the jump-on case, do not block a continuous transition towards bursting, but
much more promote it. Specifically they guide the trajectories into the
proximity of unstable limit cycles of the fast subsystem. Here once again, an
unexpected mechanisms sets in and allows attraction towards these globally
repelling cycles. The final element to the spike-adding mechanisms in this
region are rapidly oscillating segments that revolve around the branch of
unstable limit cycles and consecutively add more spikes to the orbit.

Overall, the full dynamics nearby the canard explosion can be seen as mixed-
type-like canards of peculiar nature: orbits evolve nearby the repelling
middle sheet of the critical manifolds, as well nearby unstable limit cycles.
Typically mixed-type canards are observed towards the singular limit, as they
represent a slow-fast effect. The neural mass with STP exhibits mixed-type-
like canards only for large enough forcing frequency; the phenomenon
disappears for too slow forcing and is blocked by jump-on canards.

By virtue of the NMSTP model being the main subject of this work, we want to
emphasize that it represents an exact limit of the QIF network with STP.
Discrepancies between network and NMSTP arise due to numerical errors and
finite-size fluctuations, especially under slow forcing. Nevertheless, our
results clearly show a good agreement between the two models. In particular,
the network simulations display the mixed-type-like torus canard dynamics,
hence this can be regarded as a strong evidence for the same mechanisms to be
responsible for burst spike-adding in the network.

In summary, the NMSTP turns out to be a useful approach in order to
investigate the ensemble dynamics of neuronal populations in presence of STP.
Slow-fast dissection reveals the mechanisms underlying burst generation on
population level. As is it turns out, synaptic dynamics indeed enriches the
complexity of the problem, by giving rise to peculiar jump-on and mixed-type
like torus canards, which both appear due to the timescales associated with
STP.

Finally, while this work is in the scope of neural mass models and STP, the
methodology of OA reduction and slow-fast dissection, coupled with numerical
bifurcation analysis, can be applied to a much broader class of phase
oscillator systems, like the Kuramoto model. This can lead to a better
understanding of emerging collective slow-fast dynamics of large scale
networks.

## References

  * [1] Ott E and Antonsen T M 2008 Chaos: An Interdisciplinary Journal of Nonlinear Science 18 037113 ISSN 1054-1500, 1089-7682
  * [2] Ott E and Antonsen T M 2009 Chaos: An Interdisciplinary Journal of Nonlinear Science 19 023117 ISSN 1054-1500
  * [3] Kuramoto Y 1975 Lecture notes in Physics 30 420
  * [4] Ermentrout G B and Kopell N 1986 SIAM Journal on Applied Mathematics 46 233–253
  * [5] Montbrió E, Pazó D and Roxin A 2015 Physical Review X 5 021028
  * [6] Pazó D and Montbrió E 2016 Physical review letters 116 238101
  * [7] Schmidt H, Avitabile D, Montbrió E and Roxin A 2018 PLOS Computational Biology 14 e1006430 ISSN 1553-7358
  * [8] Segneri M, Bi H, Olmi S and Torcini A 2020 Frontiers in computational neuroscience 14 47
  * [9] Pietras B, Devalle F, Roxin A, Daffertshofer A and Montbrió E 2019 Physical Review E 100 042412
  * [10] Montbrió E and Pazó D 2020 Physical Review Letters 125 248101
  * [11] di Volo M and Torcini A 2018 Physical Review Letters 121 128301
  * [12] Goldobin D S, di Volo M and Torcini A 2021 Physical Review Letters 127 038301
  * [13] Avitabile D, Desroches M and Ermentrout G B 2021 Cross-scale excitability in networks of quadratic integrate-and-fire neurons https://hal.inria.fr/hal-03326530
  * [14] Tsodyks M V and Markram H 1997 Proceedings of the national academy of sciences 94 719–723
  * [15] Wilson H R and Cowan J D 1972 Biophysical journal 12 1–24
  * [16] Tsodyks M, Pawelzik K and Markram H 1998 Neural computation 10 821–835
  * [17] Mongillo G, Barak O and Tsodyks M 2008 Science 319 1543–1546
  * [18] Mi Y, Katkov M and Tsodyks M 2017 Neuron 93 323–330
  * [19] Trübutschek D, Marti S, Ojeda A, King J R, Mi Y, Tsodyks M and Dehaene S 2017 eLife 6 e23871 ISSN 2050-084X
  * [20] Taher H, Torcini A and Olmi S 2020 PLoS Computational Biology 16 e1008533
  * [21] Wang Y, Markram H, Goodman P H, Berger T K, Ma J and Goldman-Rakic P S 2006 Nature Neuroscience 9 534–542 ISSN 1546-1726
  * [22] Adams W B and Benson J A 1985 Progress in Biophysics and Molecular Biology 46 1–49 ISSN 0079-6107
  * [23] Connors B W and Gutnick M J 1990 Trends in Neurosciences 13 99–104 ISSN 0166-2236
  * [24] Gray C M and McCormick D A 1996 Science 274 109–113 ISSN 0036-8075, 1095-9203
  * [25] Schwindt P, O’Brien J A and Crill W 1997 Journal of Neurophysiology 77 2484–2498 ISSN 0022-3077
  * [26] Su H, Alroy G, Kirson E D and Yaari Y 2001 Journal of Neuroscience 21 4173–4182 ISSN 0270-6474, 1529-2401
  * [27] Amir R, Liu C N, Kocsis J D and Devor M 2002 Brain 125 421–435 ISSN 0006-8950
  * [28] Wellmer J, Su H, Beck H and Yaari Y 2002 European Journal of Neuroscience 16 259–266 ISSN 1460-9568
  * [29] Womack M and Khodakhah K 2002 Journal of Neuroscience 22 10603–10612 ISSN 0270-6474, 1529-2401
  * [30] Plant R E and Kim M 1975 Mathematical Biosciences 26 357–375 ISSN 0025-5564
  * [31] Hindmarsh J L and Rose R M 1984 Proceedings of the Royal society of London. Series B. Biological sciences 221 87–102
  * [32] Rinzel J 1985 Bursting oscillations in an excitable membrane model Ordinary and Partial Differential Equations (Proceedings of the Eighth Conference Held at Dundee, Scotland, June 25-29, 1984) (Lecture Notes in Mathematics vol 1511) ed Sleeman B D and Jarvis R J (Springer) pp 304–316
  * [33] Rinzel J 1987 A formal classification of bursting mechanisms in excitable systems Mathematical Topics in Population Biology, Morphogenesis and Neurosciences (Proceedings of an International Symposium Held in Kyoto, November 10-15, 1985) (Lecture Notes in Biomathematics vol 71) ed Termato E and Yumaguti M (Springer) pp 267–281
  * [34] Rinzel J 1987 A Formal Classification of Bursting Mechanisms in Excitable Systems International Congress of Mathematicians, Berkeley, California, USA, August 3-11, 1986 vol II (American Mathematical Society) pp 1578–1593
  * [35] Izhikevich E M 2000 International Journal of Bifurcation and Chaos 10 1171–1266 ISSN 0218-1274
  * [36] Bacak B J, Kim T, Smith J C, Rubin J E and Rybak I A 2016 elife 5 e13403
  * [37] Ersöz E K, Desroches M, Guillamon A, Rinzel J and Tabak J 2020 Journal of mathematical biology 80 2075–2107
  * [38] Hodgkin A L 1948 The Journal of physiology 107 165–181
  * [39] FitzHugh R 1961 Biophysical journal 1 445–466
  * [40] Nagumo J, Arimoto S and Yoshizawa S 1962 Proceedings of the IRE 50 2061–2070
  * [41] Morris C and Lecar H 1981 Biophysical Journal 35 193–213 ISSN 00063495
  * [42] Van der Pol B 1926 The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 2 978–992
  * [43] Brøns M 1988 Math. Eng. Ind 2 51–63
  * [44] Benoît E, Callot J L, Diener F and Diener M 1981 Collectanea Mathematica 32 37–119
  * [45] Baspinar E, Avitabile D and Desroches M 2021 Chaos: An Interdisciplinary Journal of Nonlinear Science 31 063129 ISSN 1054-1500
  * [46] Benes G N, Barry A M, Kaper T J, Kramer M A and Burke J 2011 Chaos: An Interdisciplinary Journal of Nonlinear Science 21 023131
  * [47] Burke J, Desroches M, Barry A M, Kaper T J and Kramer M A 2012 The Journal of Mathematical Neuroscience 2 1–30
  * [48] Desroches M, Burke J, Kaper T J and Kramer M A 2012 Physical Review E 85 021920
  * [49] Desroches M, Krupa M and Rodrigues S 2016 Physica D: Nonlinear Phenomena 331 58–70
  * [50] Fenichel N 1979 Journal of differential equations 31 53–98
  * [51] Hek G 2010 Journal of mathematical biology 60 347–386
  * [52] Eckhaus W 1983 Relaxation oscillations including a standard chase on French ducks Asymptotic Analysis II— (Springer) pp 449–497
  * [53] Krupa M and Szmolyan P 2001 SIAM Journal on Mathematical Analysis 33 286–314 ISSN 0036-1410
  * [54] Neishtadt A I 1987 Differentsial’nye Uravneniya 23 2060–2067, 2204 ISSN 0374-0641
  * [55] Neishtadt A I 1988 Differentsial’nye Uravneniya 24 226–233, 364 ISSN 0374-0641
  * [56] Baer S M, Erneux T and Rinzel J 1989 SIAM Journal on Applied mathematics 49 55–71
  * [57] Desroches M, Kaper T J and Krupa M 2013 Chaos: An Interdisciplinary Journal of Nonlinear Science 23 046106
  * [58] Devalle F, Roxin A and Montbrió E 2017 PLOS Computational Biology 13 e1005881 ISSN 1553-7358
  * [59] Ceni A, Olmi S, Torcini A and Angulo-Garcia D 2020 Chaos: An Interdisciplinary Journal of Nonlinear Science 30 053121