# Canard solutions in neural mass models: consequences on critical regimes

Elif Köksal Ersöz, Fabrice Wendling

<figure><img src="image_1.png" /></figure>

## To cite this version:

Elif Köksal Ersöz, Fabrice Wendling. Canard solutions in neural mass models: consequences on critical regimes. Journal of Mathematical Neuroscience, 2021, 11 (1), pp.11. 10.1186/s13408-021-00109-z. hal-03366319

HAL Id: hal-03366319

https://hal.science/hal-03366319v1

Submitted on 5 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Canard solutions in neural mass models: consequences on critical regimes

Elif Köksal Ersöz¹ and Fabrice Wendling¹*

fabrice.wendling@inserm.fr

¹Univ Rennes, INSERM, LTSI-U1099, Campus de Beaulieu, F - 35000 Rennes, France

**Abstract**

Mathematical models at multiple temporal and spatial scales can unveil the fundamental mechanisms of critical transitions in brain activities. Neural mass models (NMMs) consider the average temporal dynamics of interconnected neuronal subpopulations without explicitly representing the underlying cellular activity. The mesoscopic level offered by the neural mass formulation has been used to model electroencephalographic (EEG) recordings and to investigate various cerebral mechanisms, such as the generation of physiological and pathological brain activities. In this work, we consider a NMM widely accepted in the context of epilepsy, which includes four interacting neuronal subpopulations with different synaptic kinetics. Due to the resulting three-time-scale structure, the model yields complex oscillations of relaxation and bursting types. By applying the principles of geometric singular perturbation theory, we unveil the existence of the canard solutions and detail how they organize the complex oscillations and excitability properties of the model. In particular, we show that boundaries between pathological epileptic discharges and physiological background activity are determined by the canard solutions. Finally we report the existence of canard-mediated small-amplitude frequency-specific oscillations in simulated local field potentials for decreased inhibition conditions. Interestingly, such oscillations are actually observed in intracerebral EEG signals recorded in epileptic patients during pre-ictal periods, close to seizure onsets.

**MSC:** 34E15; 34E17; 37N25; 92B25

**Keywords:** Multiple time-scale systems; Canards; Bursting; Excitability; Epilepsy; Neural mass model

## 1 Introduction

Brain dynamics emerges from neural entities interacting at different levels, from single neurons to large-scale neural networks. At each level, transitions between different regimes, such as firing/resting states in single neurons and up/down states in neural networks, are associated with both physiological functions and pathological activity [1–3]. One of the features of the system that determines how these transitions would occur is *excitability*. The concept of neural excitability for single neurons was introduced initially by Louis Lapique in 1907 [4, 5]. Alan Hodgkin, who then *re-introduced* the concept [6], classified the neural excitability with respect to the firing rate of neurons in response to injected steps of currents. Excitability properties of neural systems can vary with internal

<figure><img src="image_5.png" /></figure>

<footnote>*Correspondence:</footnote>
dynamics, leading to different physiological and pathological behavior [7–10]. At the cor-
tical scale, for instance, variations in excitability [11] and loss of network resilience [12]
are associated with epileptic seizures. Yet, what may be as important as a transition itself
is the dynamics preceding the transition. In the context of epilepsy, for example, identi-
fication of the dynamic features along the path to a transition is crucial for intervention
and prevention of seizures.

Mathematical models of brain activity range from microscopic level of single cell dy-
namics to macroscopic level of interactions between large scale neural systems. Neural
mass models (NMMs) consider the average temporal dynamics of interconnected neural
subpopulations without explicitly representing the underlying mechanisms at the level
of single cells. The mesoscopic level offered by the neural mass formulation has been
used to model brain signals, from local field potentials (LFPs) to global electroencephalo-
graphic (EEG) recordings, and to investigate various cerebral rhythms [13–15]. NMMs
have also been used extensively to study pathological dynamics such as in epilepsy [16–
19], Alzheimer’s disease [20] and Parkinson’s disease [21, 22].

Interactions between slow and fast components of neural systems, hence, of their math-
ematical models, result in multiple time-scale complex oscillations, such as relaxation,
bursting and mixed-mode oscillations. Geometric singular perturbation theory (GSPT) is
a key tool for understanding the interaction between the geometry of the system and the
emerging multiple time-scale dynamics. In particular, canard solutions, which can exist in
multiple time-scale systems with a folded geometry, appear as building blocks of complex
oscillations in both phenomenological and neurophysiologically plausible models ranging
from single cell [23–26] to neural networks [27, 28]. The canard phenomenon in such sys-
tems has been related to neural excitability [29], excitability thresholds [23, 30–34], and
boundaries between different type of solutions, such as subthreshold oscillations and large
amplitude spiking/bursting oscillations [24, 28, 35–43]. While such canard-organized fine
structures have been shown in a wide range of two-time-scale models, recent studies
started to explore canard-mediated processes in systems with three or more time-scales
[44–46].

In this study we investigate critical regimes in the NMM initially presented in [16].
This physiologically-grounded model has been extensively used for modeling structural
and functional changes leading to epileptic activity observed in intracranial (stereoelec-
troencephalography, SEEG) signals. The model includes four interacting neuronal sub-
populations: two interconnected subpopulations of glutamatergic pyramidal neurons and
GABAergic inhibitory interneurons (somatostatin positive (SOM+), and parparvalbunim
positive (PV+), also called dendrite-projecting slow and soma-projecting fast interneu-
rons, respectively). Although the model was introduced for the CA1 region of the hip-
pocampus, implementation of these four subpopulations mediating glutamatergic and
GABAergic signaling makes it generic enough to be considered for many other cortical
regions [47]. Activity of each subpopulation is given by the corresponding average post-
synaptic potential (PSP) that is determined by two functions: 1) a “pulse to wave” function,
$S(v) = 5/(1 + \exp(0.56(6 - v)))$, transforming the incoming synaptic potentials into a firing
rate; and 2) a “wave to pulse” converting the input average firing rate into a mean PSP at
the input of each subpopulation, that is, $h(t) = Wt/\tau_w \exp(-t/\tau_w)$, where W represents the
average synaptic gain and $\tau_w$ is the average synaptic time constant mimicking the rise and
decay of actual PSPs. The system schematized in Fig. 1a reads

$$ \dot{y}_0 = y_5, \qquad (1a) $$

$$ \dot{y}_5 = \frac{A}{\tau_a} S[y_1 - y_2 - y_3] - \frac{2}{\tau_a} y_5 - \frac{1}{\tau_a^2} y_0, \qquad (1b) $$

$$ \dot{y}_1 = y_6, \qquad (1c) $$

$$ \dot{y}_6 = \frac{A}{\tau_a} \{p(t) + C_2 S[C_1 y_1]\} - \frac{2}{\tau_a} y_6 - \frac{1}{\tau_a^2} y_1, \qquad (1d) $$

$$ \dot{y}_2 = y_7, \qquad (1e) $$

$$ \dot{y}_7 = \frac{B}{\tau_b} C_4 S[C_3 y_0] - \frac{2}{\tau_b} y_7 - \frac{1}{\tau_b^2} y_2, \qquad (1f) $$

$$ \dot{y}_3 = y_8, \qquad (1g) $$

$$ \dot{y}_8 = \frac{G}{\tau_g} C_7 S \left[ C_5 y_0 - \frac{C_6}{C_4} y_4 \right] - \frac{2}{\tau_g} y_8 - \frac{1}{\tau_a^2} y_3. \qquad (1h) $$

The variables $y_{0,1}$ stand for the excitatory PSPs mediated by two pyramidal neuron subpopulations, $y_2$ and $y_3$ are inhibitory PSPs mediated by the SOM+ and PV+ interneuron subpopulations, respectively. Variables $y_j$ ($j \in \{5, 6, 7, 8\}$) are the auxiliary variables that are introduced to convert the second-order differential equations describing the wave to pulse functions to first-order differential equations [13]. The parameters A, B, G are the synaptic gains, the $C_i$ are the connectivity constants representing the average number of synaptic contacts, $p(t)$ is the external (noisy) cortical input ($p(t) = p + \xi$, where p is the mean of external input $\xi$ is a random variable following a normal distribution with $N(0, \sigma^2)$). The synaptic time constants are given by $\tau_a$, $\tau_b$, $\tau_g$. The major contribution to LFPs (as recorded by intracranial electrodes in patients candidate to surgery) corresponds to the PSPs summarized at the level of pyramidal neurons, which are geometrically aligned “in palisades”, i.e. one relative to the other and perpendicular to the plane of the cortical layers. In the model, the LFP is given by the sum of excitatory PSP (EPSP) and inhibitory PSPs (IPSPs) received by the glutamatergic pyramidal cells, hence $LFP = y_1 - y_2 - y_3$.

As introduced in [48], under the following variable conversion:

$$ \left( \frac{t}{\tau_g}, \frac{y_0}{\tau_a}, \frac{y_1}{\tau_a}, \frac{y_2}{\tau_b}, \frac{y_3}{\tau_g}, y_5, y_6, y_7, y_8 \right) \Rightarrow (\tilde{t}, \nu_0, \nu_1, \nu_2, \nu_3, \nu_5, \nu_6, \nu_7, \nu_8), $$

with $\delta = \tau_g/\tau_a$ and $\varepsilon = \tau_a/\tau_b$, system (1a)–(1h) can be written in the following deterministic ($\sigma = 0$) slow–fast form:

$$ \frac{dv_3}{dt} = y_8 := F_3(y_8), \qquad (2a) $$

$$ \frac{dy_8}{dt} = GS[C_5 \tau_a \nu_0 - C_6 \tau_b \nu_2] - \nu_3 - 2y_8 := F_8(\nu_0, \nu_2, \nu_3, y_8), \qquad (2b) $$

$$ \frac{dv_0}{dt} = \delta y_5 := \delta F_0(y_5), \qquad (2c) $$
$$ \begin{aligned} \frac{dy_5}{dt} &= \delta(AS[A\tau_a p + C_2\tau_a v_1 - C_4\tau_b v_2 - C_7\tau_g v_3] - v_0 - 2y_5) \quad (2d) \\ &:= \delta F_5(v_0, v_1, v_2, v_3, y_5), \end{aligned} $$

$$ \frac{dv_1}{dt} = \delta y_6 := \delta F_1(y_6), \quad (2e) $$

$$ \frac{dy_6}{dt} = \delta(AS[C_1\tau_a v_0] - v_1 - 2y_6) := \delta F_6(v_0, v_1, y_6), \quad (2f) $$

$$ \frac{dv_2}{dt} = \delta\epsilon y_7 := \delta\epsilon F_2(y_7), \quad (2g) $$

$$ \frac{dy_7}{dt} = \delta\epsilon(BS[C_3\tau_a v_0] - v_2 - 2y_7) := \delta\epsilon F_7(v_0, v_2, y_7). \quad (2h) $$

In this manuscript, we will consider system (2a)–(2h) for the slow–fast analysis and (1a)–(1h) for simulations under the stochastic input. We will be using the parameter set given in Table 1, unless otherwise stated, for which $\delta = 0.3$ and $\epsilon = 0.2$. Numerical bifurcation analysis is performed in AUTO-07p software [49]. The stochastic differential equations were integrated using the Euler–Murayama method with a step size $dt = 10e^{-4}$ second in XPPAUT software [50].

As noticed in [48], system (2a)–(2h) is a three-time-scale system written in *fast form* with ($v_3, y_8$) being fast, ($v_0, y_5, v_1, y_6$) slow and ($v_2, y_7$) super-slow variables. Köksal Ersöz et al. [48] have focused on electrophysiological pre-ictal bursting patterns recorded in human patients just before the onset of seizure. Pre-ictal bursting patterns are characterized by fast oscillatory discharges (which will be referred as *spikes*) followed by a slower oscillation (a simulated pattern with the parameter set in Table 1 is exemplified in Fig. 1b). The authors have reproduced pre-ictal bursting and unveiled the mechanism yielding these solutions by decorticating the three-time-scale structure of the model. They have discussed appropriate stimulation strategies for aborting of the pre-ictal bursting, hence, for preventing a subsequent epileptic seizure. However, they did not focus on possible slow–fast transitions. Here we extend the slow–fast analysis initiated in [48] by investigating the role of slow-manifolds in transitions to relaxation and bursting type of solutions. We will focus on how canard trajectories shape the different routes from physiological to pathological brain activity. In what follows, we will go briefly through the multiple-time-scale analysis presented in [48], and then show different canard structures present in the model and how they take a part in critical transitions. Finally, we will see the system’s response to stochastic inputs near critical regimes, and make a remark on the slow oscillations observed along the path to seizure in SEEG signals recorded during pre-surgical evaluation of two patients with drug-resistant epilepsy.

<table><caption>**Table 1** Parameter values for the bursting-type discharges</caption><thead><tr><th>A (mV)</th><th>B (mV)</th><th>G (mV)</th><th>ρ (Hz)</th><th>C<sub>1</sub></th><th>C<sub>2</sub></th><th>C<sub>3</sub></th><th>C<sub>4</sub></th><th>C<sub>5</sub></th><th>C<sub>6</sub></th><th>C<sub>7</sub></th><th>τ<sub>a</sub> (s)</th><th>τ<sub>b</sub> (s)</th><th>τ<sub>g</sub> (s)</th></tr></thead><tbody><tr><td>5</td><td>5</td><td>35</td><td>90</td><td>135</td><td>108</td><td>80</td><td>25</td><td>450</td><td>121</td><td>121</td><td>0.01</td><td>0.05</td><td>0.003</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Figure 1** Model diagram and time series of a bursting solution. (a) Model diagram showing excitatory (red connections) and inhibitory (blue connections) interactions between subpopulations of pyramidal neurons (PYR and PYR') and inhibitory interneurons (PV+ and SOM+). The post-synaptic potential of a subpopulation, which is the output of $h(t)$, is multiplied by a synaptic coefficient $C_i$ before being received by another subpopulation. (b) Time series of a bursting solution for the parameter set in Table 1. The panels from the top to the bottom show the time courses of post-synaptic potentials of PV+ ($v_3$), PYR ($v_0$), SOM+ ($v_2$) and the local field potential (LFP), i.e. $A\tau_a p + C_2\tau_a v_1 - C_4\tau_b v_2 - C_7\tau_g v_3$</figcaption></figure>

## 2 Preliminaries

System (2a)–(2h) expressed in the fast time $\tilde{t}$ is called a *fast system*. The *slow system* is obtained by defining $\tilde{t}_s = \delta\tilde{t}$,

$$ \delta \frac{dv_3}{d\tilde{t}_s} = F_3(y_8), \qquad (3a) $$

$$ \delta \frac{dy_8}{d\tilde{t}_s} = F_8(v_0, v_2, v_3, y_8), \qquad (3b) $$

$$ \frac{dv_0}{d\tilde{t}_s} = F_0(y_5), \qquad (3c) $$

$$ \frac{dy_5}{d\tilde{t}_s} = F_5(v_0, v_1, v_2, v_3, y_5), \qquad (3d) $$

$$ \frac{dv_1}{d\tilde{t}_s} = F_1(y_6), \qquad (3e) $$

$$ \frac{dy_6}{d\tilde{t}_s} = F_6(v_0, v_1, y_6), \qquad (3f) $$

$$ \frac{dv_2}{d\tilde{t}_s} = \varepsilon F_2(y_7), \qquad (3g) $$

$$ \frac{dy_7}{d\tilde{t}_s} = \varepsilon F_7(v_0, v_2, y_7), \qquad (3h) $$

where the functions $F_i(\cdot)$ are as defined in (2a)–(2h). The *super-slow system* is obtained by defining $\tilde{t}_{ss} = \varepsilon\tilde{t}_s = \varepsilon\delta\tilde{t}$:

$$ \delta\varepsilon\frac{dv_3}{d\tilde{t}_{ss}} = F_3(y_8), \qquad (4a) $$

$$ \delta\varepsilon\frac{dy_8}{d\tilde{t}_{ss}} = F_8(v_0, v_2, v_3, y_8), \qquad (4b) $$
$$ \varepsilon \frac{dv_0}{dt_{ss}} = F_0(y_5), \qquad (4c) $$

$$ \varepsilon \frac{dy_5}{dt_{ss}} = F_5(v_0, v_1, v_2, v_3, y_5), \qquad (4d) $$

$$ \varepsilon \frac{dv_1}{dt_{ss}} = F_1(y_6), \qquad (4e) $$

$$ \varepsilon \frac{dy_6}{dt_{ss}} = F_6(v_0, v_1, y_6), \qquad (4f) $$

$$ \frac{dv_2}{dt_{ss}} = F_2(y_7), \qquad (4g) $$

$$ \frac{dy_7}{dt_{ss}} = F_7(v_0, v_2, y_7). \qquad (4h) $$

Systems (2a)–(2h), (3a)–(3h) and (4a)–(4h) describe different dynamics in the singular limits $\varepsilon \to 0$ and/or $\delta \to 0$, although they are equivalent for $\varepsilon \neq 0$ and $\delta \neq 0$. Letting $\delta \to 0$ in (2a)–(2h) yields the *fast layer problem* (2a)–(2b) which describes the dynamics of the fast variables $(v_3, y_8)$ for fixed values of the slow $(v_0)$ and super-slow $(v_2)$ variables. The critical manifold is defined by the equilibrium points of the fast layer problem, that is,

$$ S^0 = \{(v_3, y_8, v_0, y_5, v_1, y_6, v_2, y_7) \in \mathbb{R}^8 \mid GS[C_5\tau_a v_0 - C_6\tau_b v_2] - v_3 = 0\}, \qquad (5) $$

which is eventually in the $(y_8 = 0)$-space. Since the eigenvalues of the Jacobian matrix of the fast layer problem defined by (2a)–(2b) with respect to $(v_3, y_8)$ are $\lambda_{1,2} = -1$, the 6-dimensional critical manifold $S^0$ is normally hyperbolic and stable, thus, it is perturbed to local slow manifolds for sufficiently small $\delta > 0$. Therefore, the fast dynamics can be approximated by slow dynamics as suggested by the Fenichel theorem [51].

Setting $\delta \to 0$ in (3a)–(3h) gives an algebraic-differential slow reduced problem,

$$ 0 = F_3(y_8), \qquad (6a) $$

$$ 0 = F_8(v_0, v_2, v_3, y_8), \qquad (6b) $$

$$ \frac{dv_0}{dt_s} = F_0(y_5), \qquad (6c) $$

$$ \frac{dy_5}{dt_s} = F_5(v_0, v_1, v_2, v_3, y_5), \qquad (6d) $$

$$ \frac{dv_1}{dt_s} = F_1(y_6), \qquad (6e) $$

$$ \frac{dy_6}{dt_s} = F_6(v_0, v_1, y_6), \qquad (6f) $$

$$ \frac{dv_2}{dt_s} = \varepsilon F_2(y_7), \qquad (6g) $$

$$ \frac{dy_7}{dt_s} = \varepsilon F_7(v_0, v_2, y_7), \qquad (6h) $$

which describes the slow dynamics restricted to $S^0$. System (6a)–(6h) is a two-time-scale system of 4 slow/2 super-slow variables with $\varepsilon$ being the time-scaling parameter. The equilibria of the slow layer problem in the $\varepsilon \to 0$ limit defines the super-slow manifold $L^0$, which
is reduced to

$$L^0 = \{(v_3, y_8, v_0, y_5, v_1, y_6 v_2, y_7) \in S^0 | \\ \qquad \qquad \qquad \qquad \qquad \qquad AS[A \tau_a p + C_2 \tau_a AS[C_1 \tau_a v_0] - C_4 \tau_b v_2 - C_7 \tau_g v_3] - v_0 = 0\}, \quad (7)$$

and restricted to $S^0$ by the algebraic condition $v_3 = GS[C_5 \tau_a v_0 - C_6 \tau_b v_2] = K(v_0, v_2)$ in (6a)-(6b). The super-slow dynamics restricted to the 2-dimensional manifold $L^0$, hence to $S^0$, are given by the super-slow reduced system in the $\varepsilon \to 0$ limit of (4a)-(4h).

In order to investigate the super-slow flow on $L^0$, we consider the two-time-scale system (6a)-(6h) with the fast variable $v_3$ on $S^0$, i.e. $v_3 = K(v_0, v_2)$, and rewrite the slow reduced system (6a)-(6h) as

$$\frac{dv_0}{d\tilde{t}_s} = F_0(y_5), \quad (8a)$$

$$\frac{dy_5}{d\tilde{t}_s} = F_5(v_0, v_1, v_2, K(v_0, v_2), y_5), \quad (8b)$$

$$\frac{dv_1}{d\tilde{t}_s} = F_1(y_6), \quad (8c)$$

$$\frac{dy_6}{d\tilde{t}_s} = F_6(v_0, v_1, y_6), \quad (8d)$$

$$\frac{dv_2}{d\tilde{t}_s} = \varepsilon F_2(y_7), \quad (8e)$$

$$\frac{dy_7}{d\tilde{t}_s} = \varepsilon F_7(v_0, v_2, y_7). \quad (8f)$$

Applying the time-scaling $\tilde{t}_s = \varepsilon \tilde{t}_s$ and taking the singular limit $\varepsilon \to 0$ give the algebraic-differential system

$$0 = F_0(y_5), \quad (9a)$$

$$0 = F_5(v_0, v_1, v_2, K(v_0, v_2), y_5), \quad (9b)$$

$$0 = F_1(y_6), \quad (9c)$$

$$0 = F_6(v_0, v_1, y_6), \quad (9d)$$

$$\frac{dv_2}{d\tilde{t}_s} = F_2(y_7), \quad (9e)$$

$$\frac{dy_7}{d\tilde{t}_s} = F_7(v_0, v_2, y_7). \quad (9f)$$

The algebraic conditions (9a)-(9d) define the 'critical manifold' of (8a)-(8f) which is equivalent to $L^0$ given by (7). Notice that $L^0$ is restricted in the zero plane of the $(y_5, y_6)$-space. Assuming that $v_2$ is some function of $v_0$ on $L^0$, i.e. $v_2 = M(v_0)$, the fold points on $L^0$ are defined by

$$\mathcal{F} = \{(v_0, v_1, v_2, v_3, y_5, y_6, y_7, y_8) \in L^0 \mid v_2 = M(v_0), \frac{\partial M(v_0)}{\partial v_0} = 0\}. \quad (10)$$
<figure><img src="image_7.png" /><figcaption>**Figure 2** Critical manifold, slow manifold and folded singular points. (a) Critical manifold $S^0$ (green surface), super-slow manifold $L^0$ (red curve) and a bursting orbit in the $(v_0, v_2, v_3)$-space. The curve $L^0$ is divided into three branches at $\mathcal{F}_1$ and $\mathcal{F}_2$ (red dots) where it changes stability. The middle branch of the $L^0$ $(L_m^0)$ curve between $\mathcal{F}_1$ and $\mathcal{F}_2$ is unstable (dashed). The stable left-hand side branch $(L_l^0)$, $\mathcal{F}_1$, $L_m^0$ and $\mathcal{F}_2$ are entirely on the almost horizontal part of $S^0$ (approximately on the $(v_3 \approx 0)$-plane). The stable right-hand side branch of $L^0$ $(L_r^0)$ expands both on the horizontal and vertical parts of $S^0$. Arrows indicate the corresponding time-scale (single-headed for super-slow, double-headed for slow dynamics). (b) Super-slow manifold $L^0$ (red surface), fold curves $\mathcal{F}_{1,2}$ (black lines) and folded singular points $p_{1,2}$ (red dots) in the $(y_7, v_2, v_0)$-space. Arrows indicate the corresponding time-scale</figcaption></figure>

Figure 2a shows $S^0$ and $L^0$ in the $(v_0, v_2, v_3)$-space, and Fig. 2b $L^0$ in the $(y_7, v_2, v_0)$-space. The super-slow manifold $L^0$ expands between the lower horizontal and vertical planes of $S^0$. The part of curve $L^0$ on the lower horizontal plane of $S^0$ is folded with respect to $v_2$ at along the fold curves $\mathcal{F}_1$ and $\mathcal{F}_2$ defined by (10), i.e. $\mathcal{F} = \mathcal{F}_1 \cup \mathcal{F}_2$. In this projection, the 1-D fold curves divide $L^0$ into two stable (left-hand side $L_l^0$ and right-hand side $L_r^0$) and one unstable (middle $L_m^0$) branches on the $(v_0, v_2, v_3)$-space. We also verify that four eigenvalues of $L^0$ (two real and two complex conjugate) have negative real parts along the stable parts of $L^0$. One of the real eigenvalues changes sign along $\mathcal{F}_{1,2}$, hence the unstable middle branch is of saddle type. Along the stable and unstable branches $L^0$ is normally hyperbolic, so $L^0$ is perturbed to local super-slow manifolds for small values of $\epsilon > 0$ within (6a)–(6h); see the extension of Fenichel theory for systems with more than two time-scales [52]. On the other hand, the dynamics near the non-hyperbolic fold curves $\mathcal{F}_{1,2}$ should be investigated by using the elements of GSPT.

As being the usual strategy, we consider the desingularized version of the super-slow dynamics on $L^0$ is given by the desingularized slow reduced system (DSRS), reading

$$ \frac{dv_0}{d\hat{\tau}_s} = -y_7, \qquad (11a) $$

$$ \frac{dy_7}{d\hat{\tau}_s} = -\frac{\partial \mathcal{M}(v_0)}{\partial v_0} F_7(v_0, \mathcal{M}(v_0), y_7), \qquad (11b) $$

where $\tilde{\tau}_s = \frac{\partial \mathcal{M}(v_0)}{\partial v_0} \hat{\tau}_s$. The equilibria of (11a)–(11b) on the fold set $\mathcal{F}$ are located at $(v_0^{p_1}, y_7^{p_1}) = (1.2343, 0)$ and $(v_0^{p_2}, y_7^{p_2}) = (9.9976, 0)$ for the parameter set given in Table 1. These equilibrium points, which are not generally the true equilibria of (8a)–(8f), are related to the folded singularities of (8a)–(8f), hence of (2a)–(2h). On the other hand, equilibrium points $(v_0^{F_7}, y_7^{F_7})$, i.e. $F_7(v_0^{F_7}, y_7^{F_7}) = 0$, are ordinary singularities since they are also equilibria of (8a)–(8f), hence of (2a)–(2h). Figure 2b shows $L^0$, fold curves $\mathcal{F}_{1,2}$ and folded singular points $p_{1,2}$ in the $(y_7, v_2, v_0)$-space.
Stability of the equilibrium points of the desingularized (slow) reduced system on the fold set determines the type of the folded singularities of the original system. Classification of these equilibrium points is based on the linear stability analysis. When the desingularized (slow) reduced system is planar, this analysis can be done using the trace and the determinant of the Jacobian matrix at the fold equilibrium. If both are different from zero, the fold equilibrium can be a folded saddle, a folded node or a folded focus. If the determinant equals zero but not the trace, then the desingularized flow has a degenerate equilibrium point, which is a folded saddle-node. A folded saddle-node is either related to a saddle-node bifurcation of the folded equilibria or a transcritical bifurcation of a folded equilibrium with an ordinary equilibrium. The latter case refers to the *folded saddle-node type II* (FSN II) singularity [53, 54], where a folded node becomes a folded saddle and a regular saddle becomes a regular node. The original system exhibits a singular Hopf bifurcation close to a FSN II singularity [55, 56].

The Jacobian matrix of (11a)–(11b) has the following general form:

$$J = \begin{bmatrix} 0 & -1 \\ -\frac{\partial^2 \mathcal{M}(v_0^*)}{\partial v_0^2} F_7(v_0^*, y_7^*) - \frac{\partial \mathcal{M}(v_0^*)}{\partial v_0} \frac{\partial F_7(v_0^*, y_7^*)}{\partial v_0} & 2 \frac{\partial \mathcal{M}(v_0^*)}{\partial v_0} \end{bmatrix}, \quad (12)$$

where $(v_0^*, y_7^*)$ stands for the equilibrium point of interest. Since on the folded equilibria $2\frac{\partial \mathcal{M}(v_0^{p1,p2})}{\partial v_0} = 0$, the trace and determinant of (12) on the folded equilibria read

$$\mathbf{tr}(J^{p1,p2}) = 0, \quad \det(J^{p1,p2}) = -\frac{\partial^2 \mathcal{M}(v_0^{p1,p2})}{\partial v_0^2} F_7(v_0^{p1,p2}, y_7^{p1,p2}). \quad (13)$$

The trace and determinant of (12) on the regular equilibria read

$$\mathbf{tr}(J^{F_7}) = 2 \frac{\partial \mathcal{M}(v_0^{F_7})}{\partial v_0}, \quad \det(J^{F_7}) = - \frac{\partial \mathcal{M}(v_0^{F_7})}{\partial v_0} \frac{\partial F_7(v_0^{F_7}, y_7^{F_7})}{\partial v_0}. \quad (14)$$

Notice that the generic folded singularity condition is violated due to the fact that $\frac{\partial(F_7 v_2 + F_7 y_7)}{\partial v_0} = 0$ in (12), and $\mathbf{tr}(J^{p1,p2}) = 0$ in (13). Therefore, the folded singularities determined by (11a)–(11b) are not generic and a folded equilibrium is one of the following types: a saddle for $\det(J^{p1,p2}) < 0$, a center for $\det(J^{p1,p2}) > 0$, a nilpotent for $\det(J^{p1,p2}) = 0$. The latter degenerate type corresponds to a point in the parameter space at which a folded singularity and a regular singularity meet, i.e. $\mathbf{tr}(J^{F_7}) = 0$ and $\det(J^{F_7}) = 0$ in (14). Consequently, the equilibrium points of (11a)–(11b) related to the folded and regular singularities undergo degenerate transcritical bifurcations where (12) has two zero-eigenvalues.

Figure 3 shows the bifurcation diagram of (11a)–(11b) with respect to B in the region of interest. Two straight lines of the equilibria $(v_0^{p1}, y_7^{p1})$ and $(v_0^{p2}, y_7^{p2})$ intersect with the regular equilibria curve $F_7(v_0, y_7)$ at two bifurcation points, $BP_1$ at $B_{BP1} \approx 16.7817$ and $BP_2$ at $B_{BP2} \approx 5.4817$, which are degenerate transcritical bifurcations. For $B < B_{BP1}$, the equilibrium $(v_0^{p1}, y_7^{p1})$ is a center with two complex conjugate eigenvalues. After the bifurcation at $BP_1$, $(v_0^{p1}, y_7^{p1})$ becomes a saddle. Consequently, the system (8a)–(8f) (and (2a)–(2h)) has a folded-saddle singularity near $p_1$ for $B > B_{BP1}$. The equilibrium $(v_0^{p2}, y_7^{p2})$ is of a saddle type for $B < B_{BP2}$ and becomes a center with two complex conjugate eigenvalues at $B = B_{BP2}$.
<figure><img src="image_8.png" /><figcaption>**Figure 3** Bifurcation diagram of (11a)–(11b) with respect to *B*. Equilibrium points $(v_0^{p1}, v_0^{p2})$ lie on the lower red horizontal line $(v_0^{p1}, y_7)$, and $(v_0^{p2}, v_7^{p2})$ on the upper red horizontal line $(v_0^{p2}, y_7^{p2})$. Dashed parts of the red lines represent saddle type, solid parts represent center type solutions. The true equilibrium points line on the black curve, $F_7$. Dashed part of the black curve represents saddle type, solid parts stable focus type solutions. The saddle type equilibrium points along $(v_0^{p1}, y_7^{p1})$ and $(v_0^{p2}, y_7^{p2})$ change to center at the intersections with $F_7$ at $BP_1$ and $BP_2$, respectively</figcaption></figure>

Hence, the system (8a)–(8f) (and (2a)–(2h)) has a folded-saddle singularity near $p_2$ for $B < B_{BP2}$. Finally, in a neighborhood of $BP_1$, the equilibrium points along the $F_7(v_0, y_7)$ curve are of saddle type for $B < B_{BP1}$ and stable focus for $B > B_{BP1}$. Similarly, in a neighborhood of $BP_2$, the equilibrium points along the $F_7(v_0, y_7)$ curve are of stable focus type for $B < B_{BP2}$ and of saddle type for $B > B_{BP2}$.

As mentioned above, a generic transcritical bifurcation of regular and folded singularities is related to a FSN II singularity. In our case, a folded saddle becomes a folded center and a stable focus becomes a saddle at the *degenerate* transcritical bifurcation points $BP_1$ and $BP_2$. Furthermore, system (2a)–(2h) can undergo the (singular) Hopf bifurcations close to $BP_2$ and $BP_1$ in the parameter space (see Fig. 4d), as it will be detailed in the following sections. Hence, the interaction of regular and non-generic folded singularities can be referred as a *degenerate FSN II singularity*. A degenerate FSN II singularity in (2a)–(2h) stems from the structure of the NMM, which is defined as a second-order system that violates the generic folded singularity condition $\operatorname{tr}(J) \neq 0$.

Folded singularities can lead to canard solutions in the original system. In a planar slow-fast system, a singular Hopf bifurcation can occur near a folded singularity, which is then called a *canard point*. In such case, the amplitude of the periodic orbits bifurcated at the singular Hopf point increases stiffly in a narrow interval of the parameter (scaled by the time-scale separation parameter) that controls the transition from small amplitude to relaxation oscillations [57]. This phenomenon is known as *canard explosion* [26, 58]. A canard-explosive branch hosts small canards following the unstable branch of the critical manifold and one stable branch (so-called *canard-without-head solutions*), large canards following the unstable branch of the critical manifold and two stable branches (so-called *canard-with-head solutions*), and a *maximal canard solution* that follows the longest repelling branch. In planar multiple time-scale systems, canard solutions are tightly connected to excitability and firing thresholds [30, 31]. In higher dimensional multiple time-scale systems with at least two slow variables, the folded-singularities are generic, hence they are robust to small parameter perturbations, and canard solutions associated with folded singularities connect stable and unstable branches of a folded critical manifold [36, 53, 59–61]. Canards of folded node and FSN II singularities support mixed mode oscillations [27, 36, 44]. FSN II singularities have been identified in neuronal models where the exchange from an excitable to a relaxation oscillatory state is accompanied by subthreshold oscillations [24, 28, 42, 62]. Folded-saddle canards have been shown to be sculpting firing threshold manifolds, as well [33, 34, 63–65].

In our problem, the critical manifold $S^0$ (5) is hyperbolic, whereas the super-slow manifold $L^0$ (8a)–(8f) has a folded structure. Thus, the critical transitions occur mainly in the 6-dimensional reduced system given by (8a)–(8f). As the analysis above have shown, (8a)–
<figure><img src="image_9.png" /><figcaption>**Figure 4** Bifurcation diagrams of (2a)–(2h) with a function of ($B, C_3$). (a) Bifurcation diagram of (2a)–(2h) on the ($B, C_3$) plane. Curves are named, respectively, after the limit point (LP, black curves), Hopf (H, red curves) and homoclinic (HOM, blue curves) bifurcations in panels (b–f). Only the LP bifurcations interacting with canard solutions are plotted. Black squares indicate cusp (CP), red circles indicate Bogdanov–Takens (BT) and red squares indicate generalized Hopf (GF) bifurcations. The regions marked by black, green and purple boxes are zoomed in black, green and purple framed insets. The region where the homoclinic curve tips to the $LP_1$ is zoomed inside the green inset. (b–f) Bifurcation diagrams of (2a)–(2h) as a function of $B$ for different values of $C_3$. The limit point bifurcations of interest are marked by black squares, Hopf bifurcations by red circles, and homoclinic connections by blue stars. Stable and unstable solutions are represented by continuous and dashed curves, respectively. Along the curves of equilibrium points, (2a)–(2h) undergoes four Hopf bifurcations ($H_{1,2,3,4}$) for $C_3 = \{50, 80, 145\}$ (c, d, e) and three Hopf bifurcations ($H_{1,2,3}$) for $C_3 = 15$ (b) and ($H_{1,2,4}$) for $C_3 = 200$ (f).</figcaption></figure>

(8f) has degenerate folded singularities along the fold curve at $p_1$ and $p_2$. Notice that, since the system has neither a folded node nor a FSN II, small amplitude oscillations do not exist near $p_1$ or $p_2$. But the folded saddle, degenerate FSN II and singular Hopf bifurcations can lead to canard solutions governing the critical transitions in (8a)–(8f) (hence in (2a)–(2h)). On the other hand, the bursting behavior cannot be captured by (8a)–(8f) because (8a)–(8f) is restricted in the critical manifold $S^0$, whereas the fast oscillations of the bursting orbits leave off $S^0$. So the bursting solutions exist in the full system (2a)–(2h) (see [48] for a detailed analysis of the bursting solutions). As a result, our problem yields both three and two time-scaled behaviors. In the next section, we investigate canard dynamics near $p_1$ and $p_2$.

### 3 Multiple time-scale oscillations and canard transitions

Transitions near the folded singularities of (2a)–(2h) which lead to canard solutions depend on the system parameters. The reader may refer to Table 1 for the parameter values,
unless otherwise stated. The connectivity strength from the pyramidal cell population on the subpopulation of the SOM+ interneurons, C<sub>3</sub>, and their synaptic gain, B, appear as two crucial parameters controlling the transitions by affecting the curve F<sub>7</sub> in (11a)–(11b) (see Fig. 3). Figure 4a shows a 2-parameter bifurcation diagram in the plane (B, C<sub>3</sub>). Depending on C<sub>3</sub>, system (2a)–(2h) undergoes several Hopf bifurcations as a function of B. The first two Hopf bifurcations, H<sub>1</sub> and H<sub>2</sub>, yield harmonic oscillations, whereas the periodic branches appearing at H<sub>3</sub> and/or H<sub>4</sub> connect to multiple time-scale oscillations. Under the variations in (B, C<sub>3</sub>), H<sub>1</sub> and H<sub>2</sub> persist; and the emerging periodic orbits do not change qualitatively. On the other hand, H<sub>3</sub> and H<sub>4</sub> undergo Bogdanov–Takens (BT) bifurcations BT<sub>1,2</sub> and the corresponding periodic branches vary qualitatively. The periodic orbits emerging at H<sub>3</sub> and H<sub>4</sub> can end on homoclinic connections, namely HOM<sub>1,2,3,4</sub>.

Figures 4b–f exemplify qualitative variations in (2a)–(2h) as a function of B for different values of C<sub>3</sub>. For C<sub>3</sub> < C<sub>3,BT1</sub> ≈ 18.9, the system undergoes three Hopf bifurcations, for instance in Fig. 4b for C<sub>3</sub> = 15. The branch of periodic solutions starting at H<sub>3</sub> terminates at a homoclinic connection, HOM<sub>1</sub>. As C<sub>3</sub> increases, HOM<sub>1</sub> and LP<sub>1</sub> get closer while the amplitude and the number of spikes of the periodic orbits increase. The spike adding occurs as the HOM<sub>1</sub> curve folds back and forth in the (B, C<sub>3</sub>)-space (see the black framed inset in Fig. 4a for an example folding). At C<sub>3</sub> = C<sub>3,BT1</sub> another Hopf bifurcation, H<sub>4</sub>, appears yielding a new branch of periodic orbits making a second homoclinic connection, HOM<sub>2</sub> (green framed inset in Fig. 4a). Consequently, HOM<sub>1</sub> and HOM<sub>2</sub> points coexist in a narrow range of (B, C<sub>3</sub>). The HOM<sub>1</sub> curve touches the LP<sub>1</sub> curve at (B, C<sub>3</sub>) ≈ (23.98, 22.43), folds back and continues in the parameter space, which then we name as HOM<sub>3</sub> curve (dashed zone in the green framed inset in Fig. 4a). The curves HOM<sub>1</sub> and HOM<sub>3</sub> stay very close to each other in (23.98 < B < 24.46, 21.94 < C<sub>3</sub> < 22.43), before HOM<sub>3</sub> bends in the C<sub>3</sub> direction at (B, C<sub>3</sub>) ≈ (24.46, 22.43). As it happens, the branch of periodic orbits curls below LP<sub>1</sub> in the B-space and eventually connects to HOM<sub>3</sub>. With increasing C<sub>3</sub>, this branch of periodic orbits advances further towards the stable equilibrium points while introducing a region of multi-attractors of nodes, saddles, unstable small oscillations and stable large amplitude bursting oscillations (see Fig. 4c for C<sub>3</sub> = 50, dynamics will be detailed in Sect. 3.2). Concurrently, H<sub>4</sub> moves away from LP<sub>1</sub> and HOM<sub>2,3</sub> approach LP<sub>2</sub>. In (20.09 < B < 20.26, 54.08 < C<sub>3</sub> < 54.43), HOM<sub>2</sub> and HOM<sub>3</sub> curves are connected by a section that is parallel to the LP<sub>2</sub> curve (purple framed inset in Fig. 4a). For C<sub>3</sub> > 54.43, the branch of periodic orbits initiated at H<sub>4</sub> connects to the branch of large amplitude multiple time-scale oscillations (Fig. 4d).

System (2a)–(2h) does not have any LPs between H<sub>3</sub> and H<sub>4</sub> for C<sub>3,CP1</sub> ≤ C<sub>3</sub> ≤ C<sub>3,CP2</sub>. At C<sub>3</sub> = C<sub>3,CP2</sub> ≈ 141.4, as the lower branch of equilibrium points curls below H<sub>3</sub>, the connection between the large amplitude orbits and H<sub>3</sub> is broken up on a saddle-saddle homoclinic bifurcation (the equilibrium points in a neighborhood of LP<sub>3</sub> for B ≥ B<sub>LP3</sub> are saddles [66, 67]). As a consequence, the branch of periodic orbits starting from H<sub>3</sub> terminates on a homoclinic connection, HOM<sub>4</sub> (for C<sub>3</sub> = 145 in Fig. 4e). This homoclinic connection remains until H<sub>3</sub> disappears at C<sub>3</sub> = C<sub>3,BT2</sub> ≈ 157. Beyond C<sub>3</sub> > C<sub>3,BT2</sub>, the large amplitude bursting orbits introduced by H<sub>4</sub> terminate on a saddle-node homoclinic connection (for instance at C<sub>3</sub> = 200 in Fig. 4f).

The Hopf bifurcations H<sub>3</sub> and H<sub>4</sub> occur close to the folded singularities p<sub>2</sub> and p<sub>1</sub>, respectively. System (2a)–(2h) can yield canard solutions close to these points in the parameter space of B, such as B ≈ B<sub>H3</sub> and B ≈ B<sub>H4</sub>. In the following section, we will show the
canard-mediated transition from sinusoidal oscillations initiated by $H_3$ to large amplitude bursting/relaxation type solutions. Subsequently, Sect. 3.2 will detail the canard dynamics and related excitability near $H_4$, in particular, the type-I excitability for $C_3 = 50$ and type-II excitability for $C_3 = 80$.

## 3.1 Canard-mediated transitions between sinusoidal and multiple time-scale oscillations

Köksal Ersöz et al. [48] have realized that the number of spikes of a bursting solution of (2a)–(2h) depends on the amount of the PSP received by the PV+ interneuron subpopulation, hence on the EPSP coming from the pyramidal cell subpopulation and the IPSP from the SOM+ interneurons. For instance, increasing the IPSP on the PV+ interneurons by increasing $B$ decreases the number of spikes while driving the oscillations one peak to the next one in the parameter space (see Fig. 4c–f and 5). The connectivity constant from the pyramidal cell subpopulation to the PV+ interneuron subpopulation, $C_5$, directly scales the EPSP on this subpopulation, therefore determines the maximum number of fast spikes of the bursting oscillations, or more generally, the type of the multiple time-scale oscillations.

Figure 5 exemplifies how $C_5$ modulates the large amplitude oscillations between $H_3$ and $H_4$ on the bifurcation diagram of (2a)–(2h) for $C_3 = 80$. Increasing $C_5$ decreases the amplitude of the oscillations, moves $H_1$ and $H_2$ slightly to the right, but does not affect considerably the locations $H_3$ and $H_4$ with respect to $B$ (Fig. 5a). The supercritical Hopf bifurcation at $H_3$ yields a branch of sinusoidal periodic oscillations (Fig. 5b) that folds back and forth as $B$ varies and enters in a regime of multiple time-scale periodic oscillations. These oscillations are of relaxation type for small values of $C_5$, and of bursting type for large values of $C_5$. Furthermore, the stable sinusoidal and multiple time-scale periodic oscillations can coexist depending on the values of $C_5$ (see Fig. 5b at $B \approx 4.8$).

The form of the branch of periodic solutions between $H_3$ and $H_4$ in Fig. 5a indicates the type of the multiple time-scale oscillations for a certain parameter combination. For $C_5 = 80$ (black diagram in Fig. 5) the smoothly decreasing amplitude of $v_0$ with $B$ indicates that the corresponding orbits are of relaxation type (exemplified in Fig. 6). The

<figure><img src="image_10.png" /><figcaption>**Figure 5** Variation of large amplitude solutions with respect to ($B, C_5$). (a) Bifurcation diagrams of (2a)–(2h) with respect to $B$ for $C_3 = 80$ and different values of $C_5$. Curves and Hopf bifurcations ($H_{1,2,3,4}$, dots) are colored with respect to the color codes of $C_5$ values. Stable and unstable solutions are represented by continuous and dashed curves, respectively. For the sake of simplicity, the periodic solutions between $H_1$ and $H_2$ are not shown. (b) Zoom into the region of transitions between sinusoidal and large amplitude multiple-time-scale solutions in $B \in [4.7, 4.9]$</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 6** Example canard orbits along the transition from sinusoidal oscillations to relaxation oscillations. Zoom near the bifurcation diagram for $C_3 = 80$, $C_5 = 80$ and $B \in [4.7, 4.9]$ (see Fig. 5a for the whole diagram). Continuous and dash curves represent stable and unstable solutions, respectively. Hopf ($H_3$, red dot) is marked on the diagram. Numbered orbits from 1-7 are given in panels (**b-d**). The orange curve traces the frequency of the oscillations emerging at $H_3$. The orange curve traces the frequency of the oscillations emerging at $H_3$. (**b**) Periodic orbits marked in panel (**a**), $L^0$ (red curve), fold curves $F_{1,2}$ (red points) and the critical surface $S^0$ (green surface) are projected on the ($v_0, v_2, v_3$)-space. Arrows indicate the corresponding time-scale (single-headed for super-slow, double-headed for slow dynamics). (**c**) Periodic orbits marked in panel (**a**), $L^0$ (red surface), fold curves $F_{1,2}$ (black curves) and folded singular points $p_{1,2}$ (red dots) are projected on the ($y_7, v_2, v_0$)-space. Arrows indicate the corresponding time-scale. (**d**) Time series of the periodic orbits on panels (**b**, **c**) with respective color codes. Period is normalized to 1 ($\tilde{t}/\tilde{T} = 1$, where $\tilde{T}$ represents period of a cycle)</figcaption></figure>

horizontal zigzags along the upper part of the periodic branches obtained for greater values of $C_5$ indicate the presence of bursting solutions along these periodic branches and the number of their fast spikes. For instance, the 5 peaks that we count between $H_3$ and $H_4$ for $C_5 = 350$ (blue diagram in Fig. 5) signify that the maximum number of fast spikes for $C_5 = 350$ is 4. Such a bursting orbit is obtained for sufficiently small values of $B$ ($B = 5$, for instance). Then as $B$ increases, the bursting orbits lose their fast spikes one by one through the peaks of the horizontal branch. They become relaxation cycles ($B = 16$, for instance), before shrinking and disappearing via a subcritical Hopf bifurcation at $H_4$.
The periodic solutions connected to $H_1$ and $H_2$ do not interact with the singular fold points of $L^0$, $p_1$ and $p_2$, but the ones near $H_3$ and $H_4$ do because $H_3$ and $H_4$ take place close to the degenerate FSN II singularities on the fold curve $\mathcal{F}$. As a consequence, the multiple time-scale orbits emanating at singular Hopf bifurcations $H_3$ and $H_4$ can undergo canard explosion along which canard trajectories sculpt the periodic oscillations. Figure 6 shows example orbits along the periodic branch that follow $H_3$ for $C_3 = 80$ and $C_5 = 80$. As the periodic branch folds with respect to $B$ and becomes unstable (Fig. 6a), the sinusoidal orbits of 4.5–6 Hz start to interact with $p_2$. In particular, they move along the unstable branch of $L^0$, $L_m^0$, before jumping back to the stable branch $L_r^0$. Hence, the periodic orbits become canard orbits (the 1st orbit). As $B$ varies along the periodic branch in the parameter space, the canard orbits grow in amplitude along $L_m^0$ (the second, third and fourth orbits) until they stretch out between $\mathcal{F}_1$ and $\mathcal{F}_2$ (the fifth orbit). The canard orbits that oscillate between $L_m^0$ and $L_r^0$ can be interpreted as canard-without-head orbits, and the fifth orbit as the maximal canard since it has the largest period of the canard family of the periodic branch under consideration. Soon after the fifth orbit, the trajectories jump to the attracting branch $L_l^0$, get a shape of canard-with-head solutions, and become stable (the sixth orbit). As $B$ increases, the relaxation cycles appear with parts exclusively following the attracting branches of $L^0$ and jumping close to the fold points.

As mentioned in the introduction, canard solutions play a fundamental role in separating different dynamical regimes. The unstable canard orbits in Fig. 6 (from the first to the fifth) appear as an other example of this phenomenon by accompanying the transition from sinusoidal oscillations to relaxation oscillations. For instance, sinusoidal oscillations and large amplitude canard-with-head cycles coexist for $B \in (4.75, 4.81)$ and the canard-without-head cycles form the boundary between them, as seen clearly in Fig. 6a. While increasing $C_5$ introduces bursting type of solutions, it can preserve the bistability between the bursting and sinusoidal oscillations, for example for $C_5 = \{250, 350, 425\}$ in Fig. 5. Notice that with increasing $C_5$, the initially smooth branch of periodic orbits becomes steeper, gains vertical zigzags that move to the right along the $B$-axis, and the region of bistability decreases.

Figure 7 zooms into the region of canard orbits following the sinusoidal solutions of 4.8–6 Hz started at $H_3$ for $C_5 = 450$. As the stable sinusoidal oscillations grow in amplitude with increasing $B$, they start to interact with $p_2$ and to follow the bits of $L_m^0$ (the first orbit). Soon after, the orbits undergo a LP bifurcation (where the branch of the periodic orbits folds back at $B \approx 4.821$) and become unstable. As $B$ varies in the parameter space along the periodic branch, the orbits moving along $L_m^0$ in the super-slow time-scale grow in amplitude and they start to interact with the vertical panel of $S^0$ as they jump to $L_r^0$ in the slow time-scale. So, the orbits become canard orbits.

As the part of the trajectory along $L_m^0$ grows in amplitude, the trajectory gets attracted by $L_r^0$ along the vertical panel of $S^0$ and it spirals around $L_r^0$ before landing on the horizontal plane of $S^0$. This interaction with the vertical panel of $S^0$ occurs in the fast time-scale, and eventually yields fast spikes, i.e., bursting-type canard oscillations. For instance, the second orbit in Fig. 7c and 7d has one fast spike. The number of spikes increases as the trajectory stays longer and longer along $L_m^0$ while $B$ varies. More precisely, the number of spikes changes by one as we pass from one fold to another on the same side of along the snaking periodic branch with respect to $B$ (Fig. 7b). For instance, solution-2 has 1 spike, solution-3, which is two fold below, has 3 spikes and so on. The spike adding continues
<figure><img src="image_12.png" /><figcaption>**Figure 7** Example canard orbits along the transition from sinusoidal oscillations to bursting oscillations (a) Bifurcation diagram for $C_3 = 80$, $C_5 = 450$ and $B \in [4.7, 4.9]$ (see Fig. 5a for the whole diagram). Stable and unstable solutions are represented by continuous and dashed curves, respectively. Hopf ($H_3$, red dot) is marked on the diagram. The rectangular region is zoomed in panel (b). Numbered orbits from 1-7 are given in panels (c-d). The orange curve traces the frequency of the oscillations emerging at $H_3$. (c) Periodic orbits marked in panels (a, b), $L^0$ (red curve), fold curves $F_{1,2}$ (red points) and the critical surface $S^0$ (green surface) are projected on the ($v_0, v_2, v_3$)-space. Arrows indicate the flow direction and its time-scale (single-head for super-slow, double-head for slow, triple-head for fast). (d) Periodic orbits marked in panel (a), $L^0$ (red surface), fold curves $F_{1,2}$ (black curves) and folded singular points $p_{1,2}$ (red dots) are projected on the ($y_7, v_2, v_0$)-space. Arrows indicate the corresponding time-scale. (e) Time series of the periodic orbits on panels are (b, c) shown with respective color codes. Period is normalized to 1 ($\tilde{t}/\tilde{T} = 1$, where $\tilde{T}$ represents period of a cycle)</figcaption></figure>

until the canard orbits (analogous to canard-without-head orbits) expands between two folded singularities, hence till the occurrence of the maximal canard of the family (approximated by the sixth orbit). After the maximal canard, the canard cycles start to follow $L_l^0$ in the super-slow time-scale (analogous to canard-with-head orbits) and they become stable (the seventh orbit). The part of the trajectory along $L_m^0$ decreases as $B$ increases further.
For $C_3 = 80$, the system undergoes a complete canard explosion along the periodic branch following $H_3$ since it visits the whole canard family from small to large (Figs. 6 and 7). However, what we may observe for small values of $C_3$ is an incomplete canard explosion terminating at a homoclinic connection. For instance, for $C_3 = 15$ (see Fig. 4b), the sinusoidal oscillations along the periodic branch initiated at $H_3$ change qualitatively by interacting with $L_p^0$ as $B$ varies and we observe homoclinic canard-without-head orbits at $HOM_1$. The orbits terminating on $HOM_3$ are homoclinic canard-with-head orbits surrounded by stable large amplitude oscillations. Increasing $C_3$ completes the canard explosion and the system enters into an excitable regime which will be detailed in the following section.

## 3.2 Canard-mediated transitions and excitability

According to Hodgkin's [6] classification of neural excitability, type-I excitable neurons have continuous frequency-injected current curves, whereas type-II excitable neurons have discontinuous frequency-injected current curves. Rinzel and Ermentrout [68, 69] linked the type-I excitability to a SNIC bifurcation and the type-II excitability to a Hopf bifurcation. De Maesschalck and Wechselberger [29] explained the transition between the two excitability types via an intermediate regime of type-I excitability associated with a codimension-2 Bogdanov-Takens (BT) bifurcation in a planar system. They showed the existence of incomplete canard transitions in this transitory regime. Later on, transitions between the neuronal excitability types was shown to be induced by the inhibitory and excitatory autapse in the Morris-Lecar model [70]. Folded singularities and corresponding canard solutions in higher dimensional systems also have been shown to be shaping systems' excitability properties [24, 28, 33, 34, 63–65].

System (2a)–(2h) can yield large amplitude oscillations in response to certain forms of stimulation (due to stochastic inputs, for instance) after being initiated from an equilibrium point for a $B$ value close to $H_4$, $LP_1$ and $LP_2$ in Figs. 4c–4f. Hence, system (2a)–(2h) is excitable in these regions and the excitability properties of (2a)–(2h) determined by the parameter $C_3$ (see Fig. 4). Indeed, the local pictures in these regions are similar to the ones investigated in [29, 70]. In particular, system (2a)–(2h) is type-I excitable for $C_3 \in (22.43, 54.43)$, basically between the homoclinic/saddle-saddle interactions near $LP_1$ and $LP_2$. In this parameter region, the large amplitude oscillations terminate on a homoclinic orbit for which the firing frequency is zero. System (2a)–(2h) is type-II excitable for $C_3 > 54.43$ for which the termination is issued via a Hopf bifurcation. In both cases, canard solutions shape the resulting dynamics.

Figure 8 zooms in near the excitable region for $C_3 = 50$ (see Fig. 4c for the whole diagram). For a particular value of $B$ for $B < B_{LP1}$, the only attractor is the large amplitude bursting oscillation (the 1st orbit). In $B_{LP1} < B < B_{H4}$ the unstable attractors of the equilibrium points appear. The subcritical Hopf bifurcation at $B = B_{H4}$ initiates a branch of periodic orbits that terminates on the homoclinic point $HOM_2$, which bounds the canard explosion near $H_4$. For $B_{HOM2}$, a homoclinic canard-without-head orbit (the second orbit) coexists with a large stable bursting orbit of canard type (the third orbit). At $B_{HOM3}$ a homoclinic canard-with-head orbit (the fifth orbit) appears together with an outer large amplitude canard cycle (the fourth orbit). The large amplitude canard cycle grows in amplitude and disappears on a saddle-node of periodic orbits (SNPO) at $B = B_{SNPO}$ (the sixth orbit). We also notice that, as $HOM_3$ gets closer to $LP_2$ for $C_3 \approx 54.4$ and $B \approx B_{LP2}$, the canard orbits on the $HOM_3$ become of without-head type.
<figure><img src="image_13.png" /><figcaption>**Figure 8** Example canard orbits near the type-I excitable regime. (a) Bifurcation diagram for $C_3 = 50$, $C_5 = 450$ and $B \in [18, 22.5]$ (see Fig. 4c for the whole diagram). Stable and unstable solutions are represented by continuous and dashed curves, respectively. Limit point ($LP_{1,2}$, black squares), Hopf ($H_4$, red dot), homoclinic ($HOM_{2,3}$, blue stars) bifurcations and saddle-node bifurcation of periodic orbits (SNOP, orange purple square) are marked on the diagram. Numbered solutions are presented in panels (b-d). The orange curves trace the frequency of the oscillations. (b) Periodic orbits marked in panel (a), $L^0$ (red curve), fold curves $F_{1,2}$ (red points) and the critical surface $S^0$ (green surface) are projected on the ($v_0, v_2, v_3$)-space. Arrows indicate the flow direction and its time-scale (single-headed for super-slow, double-headed for slow dynamics. The homoclinic points $HOM_2$ and $HOM_3$ are marked by cyan and dark blue stars. (c) Periodic orbits marked in panel (a), $L^0$ (red surface), fold curves $F_{1,2}$ (black curves) and folded singular points $p_{1,2}$ (red dots) are projected on the ($y_7, v_2, v_0$)-space. Arrows indicate the corresponding time-scale. (d) Time series of the periodic orbits on panels (b, c) with respective color codes. Period is normalized to 1 ($\bar{t}/\bar{T} = 1$, where $\bar{T}$ represents period of a cycle)</figcaption></figure>

For a parameter set ensuring the type-II excitability ($C_3 = 80$, for instance), the fast spikes of the bursting oscillations disappear and the final oscillation turns out to be a relaxation type running in slow and super-slow timescales. These relaxation oscillations terminate via a complete canard-explosion near the singular Hopf bifurcation point $H_4$. This happens in a similar manner for all $C_5$ values under consideration. Figure 9 provides an example for $C_5 = 450$. As the large amplitude periodic solutions decrease in amplitude, they start to follow $L_m^0$ and take the shape of canard-with-head solutions (the second and third
<figure><img src="image_14.png" /><figcaption>**Figure 9** Example canard orbits near the type-II excitable regime. (a) Bifurcation diagram for $C_3 = 80$, $C_5 = 450$ and $B \in [16, 18]$ (see Fig. 5a for the whole diagram). Stable and unstable solutions are represented by continuous and dashed curves, respectively. Hopf bifurcation ($H_4$, red dot) is marked on the diagram. Numbered orbits on the lower branch of periodic solutions are presented in panels (c–d). The orange curve traces the frequency of the oscillations emerging at $H_4$. (b) Periodic orbits marked in panel (a), $L^0$ (red curve), fold curves $F_{1,2}$ (red points) and the critical surface $S^0$ (green surface) are projected on the ($v_0, v_2, v_3$)-space. Arrows indicate the flow direction and its time-scale (single-headed for super-slow, double-headed for slow dynamics). (c) Periodic orbits marked in panel (a), $L^0$ (red surface), fold curves $F_{1,2}$ (black curves) and folded singular points $p_{1,2}$ (red dots) are projected on the ($y_7, v_2, v_0$)-space. Arrows indicate the corresponding time-scale. (d) Time series of the periodic orbits on panels (b, c) are shown with respective color codes. Period is normalized to 1 ($\tilde{t}/\tilde{T} = 1$, where $\tilde{T}$ represents period of a cycle)</figcaption></figure>

orbits). The maximal canard of this canard family is the fourth orbit that stays along to the super-slow manifolds as long as possible. After the fourth orbit, we observe canard-without-head orbits (the 5th and the 6th orbits) that shrink to $p_1$. The frequency of the oscillations along the canard explosion ranges in 1.8–3.5 Hz. We also notice a region of bistability between large amplitude bursting oscillations and equilibrium points. Once again the canard solutions construct the boundary between them. For a parameter set giving relaxation oscillations in this region (e.g. $C_5 = \{80, 250, 350\}$ in Fig. 5), the relaxation oscillations shrink $H_4$ via a ‘classical’ canard explosion, similar to the one in the 2D van der Pol system, without having any fast component in $v_3$.
## 4 Local field potential in critical regimes

In the previous section we have shown two different regions in parameter space where canard solutions determine boundaries and organize transitions between different dynamical regimes. The narrow-band sinusoidal activity of 4.5–6 Hz emerging near $H_3$ and of 1.8–3.5 Hz emerging near $H_4$ are connected to large amplitude periodic multiple-timescale solutions through canard orbits. System (1a)–(1h) emits aperiodic large amplitude epileptic discharges under stochastic input ($p(t) = p + \xi$, with $\xi = N(0, 2^2)$) when it is initialized near the critical regions of $H_3$ and $H_4$ (Fig. 10–11). A parameter setting ensuring the type-I excitability without any canard solutions near $H_4$ gives a board band activity between the large amplitude spikes (Fig. 10a1–a3). On the other hand, taking the system to type-II excitability near $H_4$ introduces transient small amplitude oscillations of ≈ 3.5 Hz due to the presence of the canard cycles in this region (Fig. 10b1–b3). We observe transitions between large amplitude discharges and harmonic oscillations of ≈ 6 Hz when the system is initialized close to the Hopf bifurcation $H_3$ (Fig. 10c1–c3). Simulated PSPs at the level of the pyramidal cell subpopulation are given in Fig. 11.

Figures 12 and 13 show LFPs recorded by the SEEG electrodes in two different patients with drug-resistant focal epilepsy during presurgical evaluation (see Table 2 for the details). Multiple-contact depth electrodes were implanted according to the SEEG technique as a standard clinical procedure in the care of patients who consented the possible use of data for research purpose. The positioning of the electrodes is determined in each patient from hypotheses about the localization of the epileptogenic areas. Implantation accuracy

<figure><img src="image_15.png" /><figcaption>**Figure 10** LFP traces of system (1a)–(1h) near critical transitions under stochastic input. (a1) Transitions between multiple time-scale oscillations and background regime for a type-I setting at $B = 23$, $C_3 = 35$, $C_5 = 200$. Panel (a2) zooms between two large amplitude discharges (blue) and panel (a3) shows the normalized power spectral density of the signal. (b1) Transitions between multiple time-scale oscillations and background regime with slow oscillations of ≈ 3.5 Hz for a type-II $B = 17.8$, $C_3 = 80$, $C_5 = 200$. Panel (b2) zooms between two large amplitude discharges (red) and panel (b3) shows the normalized power spectral density of the signal. (c1) Transitions between multiple time-scale oscillations and sinusoidal oscillations for $B = 4.7$, $C_3 = 80$, $C_5 = 200$. Panel (c2) zooms between two large amplitude discharges (cyan) and panel (c3) shows the normalized power spectral density of the signal</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 11** Corresponding PSPs of the LFPs given in Fig. 10. (a1)–(a3) The EPSP, slow IPSP (IPSP<sub>s</sub>) and fast IPSP (IPSP<sub>f</sub>) for the type-I setting given in Fig. 10a. (b1)–(b3) The EPSP, slow IPSP (IPSP<sub>s</sub>) and fast IPSP (IPSP<sub>f</sub>) for the type-II setting given in Fig. 10b. (c1)–(c3) The EPSP, slow IPSP (IPSP<sub>s</sub>) and fast IPSP (IPSP<sub>f</sub>) for the setting given in Fig. 10c</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 12** SEEG signals recorded in a patient with epilepsy during the inter-ictal and ictal transition. (a) Transition from inter-ictal to ictal period in the first patient. Background activity (BKG) observed further away from seizure (panel b1) has a broad-band frequency distribution (normalized power density spectrum in panel b2). A sporadic spike is preceded by narrow band low amplitude *resembling* canard-mediated oscillations (CMO, marked in red, zoomed in panel (c1)) at ≈ 3.5 Hz (normalized power density spectrum in panel (c2))</figcaption></figure>

peri-operatively is verified by an X-ray CT scan. A post-operative CT scan without con-
trast product is then used to verify the precise 3D location of each electrode contact. After
SEEG exploration, intracerebral electrodes are removed. An MRI is performed on which
<figure><img src="image_18.png" /><figcaption>**Figure 13** SEEG signals recorded in a patient with epilepsy during the inter-ictal and ictal transition. (a) Transition from inter-ictal to ictal period in the second patient. Background activity (BKG) observed further away from seizure (panel **b1**) has a broad-band frequency distribution (normalized power density spectrum in panel **b2**). A sporadic spike is preceded by narrow band low amplitude resembling canard-mediated oscillations (CMO, marked in red, zoomed in panel (c1)) at ≈ 7 Hz (normalized power density spectrum in panel (c2))</figcaption></figure>

<table><caption>**Table 2** Summary of patients' features</caption><thead><tr><th>Feature</th><th>Patient 1 (Fig. 12)</th><th>Patient 2 (Fig. 13)</th></tr></thead><tbody><tr><td>Age at SEEG</td><td>16 y</td><td>14 y</td></tr><tr><td>Gender</td><td>Female</td><td>Female</td></tr><tr><td>MRI</td><td>Right occipito-temporal focal cortical dysplasia</td><td>Left hippocampal sclerosis</td></tr><tr><td>Syndrome</td><td>Temporal lobe epilepsy (temporal plus)</td><td>Mesial temporal lobe epilepsy</td></tr><tr><td>Recorded Cerebral Region</td><td>Medial part of the middle temporal gyrus (electrode contact: B4)</td><td>Internal temporal pole (electrode contact: Pt*1)</td></tr><tr><td>Surgical outcome (Engel Class)</td><td>II (cortectomy)</td><td>IA (anterior temporal lobectomy)</td></tr></tbody></table>

the trajectory of each electrode remains visible. Finally, a CT-scan/MRI data fusion is performed to anatomically locate each contact along each electrode trajectory. The patient had electrodes implanted in the temporal region. For this study, signals were selected as they exhibited clear transitions in electrophysiological patterns. In particular, we selected pre-ictal events followed by a fast discharge typical of the seizure onset, which is one of the markers of the imbalanced relation of excitation and inhibition [16, 71] that involves excitability variations.

In Fig. 12, a narrow band activity of theta-band of 3.5 Hz is followed by a large amplitude epileptic discharge between two sporadic discharges as we advance towards sustained pre-ictal discharges. Such narrow band activity may be a signature of canard-mediated regions where slowly varying system's parameters and/or remote interactions lead to transitions between small-amplitude-low-frequency oscillations and large amplitude discharges. In Fig. 13 a narrow band activity of about 7 Hz is followed by a large amplitude epileptogenic
discharge. We also notice that the form of the epileptic discharges in Fig. 12 and Fig. 13 are different, which may indicate that the systems would have different characteristics. Interestingly, we have identified parameter regions for the corresponding frequency bands in the model (see simulated LFPs and PSPs in Fig. 10 and Fig. 11, respectively). Hence, we think that the properties of transient narrow band oscillations may be related to the excitability properties and level of synaptic projections (scaled by the coupling coefficients in the model) of the epileptogenic zone.

## 5 Conclusion

In this article, we extended the multiple time-scale analyses previously initiated in [48]. Here, we both investigated canard transitions present in a neurophysiologically-relevant NMM and analyzed their consequences in terms of subsequent signatures in LFPs. In this three-time-scale model, the canard transitions occur in the 6-dimensional two-time-scale reduced system of slow and super-slow variables. They are associated with degenerate FSN II singularities and singular Hopf bifurcations. They organize initiation of relaxation/bursting oscillations from harmonic oscillations of 4.5–6 Hz or from equilibrium points, and determine the boundaries between them. We showed that the system switches between type-I and type-II excitability near the transitions between the equilibrium points and relaxation/bursting oscillations. We further noticed that the canard regimes of type-II excitability (and partially of type-I) yield low-frequency (near 3.5 Hz) oscillations in the LFP under stochastic input.

These model predictions motivated a close analysis of SEEG recordings performed in epileptic patients. In this paper, results illustrative of both signatures are reported only in two patients. Interestingly, in brain structures clearly involved in the transition from interictal to ictal activity, we observed a narrow band activity between sporadic discharges before the seizure initiation, which strongly differed from the preceding background activity. Although the parameter set used in this paper was not aimed for modeling these recordings specifically, it is striking to see such a matching between the mathematical analysis and the actual recordings.

It has been evidenced that impaired excitation–inhibition balance shapes the activity of neural networks and, therefore, causes the emergence of “pathological” electrophysiological patterns such as pre-ictal spikes and seizures in the context of epilepsy (see for a review [72]). Indeed, epileptogenic brain regions are typical example of such excitation–inhibition imbalanced networks [73]. We showed that the level of EPSP on the subpopulation of SOM+ interneurons determines the type of the excitability. In particular, the system is type-I excitable if the average number of synaptic contacts from the excitatory pyramidal cells to the GABAergic SOM+ interneurons is low, and type-II excitable if the average number of synaptic contacts from the excitatory pyramidal cells to the GABAergic SOM+ interneurons is high. It is then the decreasing GABAergic inhibition (modeled by decreasing inhibitory drive by the subpopulation of the SOM+ interneurons) that is responsible for transitions from background to epileptiform discharges. Interestingly, such model parameter variations are plausible and linked to the failure of inhibitory barrages observed in epileptic tissues [74] and generation of slow waves preceding the fast activity [75, 76]. Properties of emerging epileptic discharges (e.g. their shape and frequencies), and possible “silent” phases in between are strongly connected to the type of the excitability. In the context of epilepsy, transitory regimes between the background activity and epileptic discharges are crucial for understanding the underlying mechanisms [11, 77]. Epileptic
biomarkers during such regimes, such as high-frequency oscillations [78], shape features of epileptic spikes [79] or maybe frequency-specific oscillations reported here, are essen- tial for identification of epileptogenic networks and for further development of therapeutic procedures. Verification of the presence of such oscillations across different patients and accurate modeling of the inter-ictal activity are needed, of course, for suggesting them as biomarkers. This is the topic of future investigations.

As epilepsy can be considered as a dynamic disease [73, 80, 81], mathematical mod-
els of different cellular levels inherit multiple time-scale thinking [82–84]. We note a few
studies on the slow–fast transitions in NMMs. Desroches et al. extended NMMs [85] by
considering the synaptic gain of SOM+ interneurons as a slowly changing variable. They
showed that this configuration introduced regions of torus canards. Jafarian et al. [86] pro-
posed a NMM which incorporates slow variations in ionic currents leading to spontaneous
paroxysmal activity. Hebbink et al. [87] investigated response of the NMM of Wendling
et al. [16] to slowly varying inputs under which the systems yields bursting oscillations.
Weigenand et al. remarked the role of canard solutions in fast transitions in sleep wake
patterns of K-complexes in a NMM of sleep-wake patterns [43]. Our paper shows that
canard-mediated solutions are naturally present in the NMM of Wendling et al. [16]. Im-
portantly, as this model implements two main sub-types of interneurons (dendrite- and
soma-projecting), it is generic and can be considered for studying the dynamics of other
regions than hippocampus, such as neocortical areas, and in different contexts, such as
consciousness [47] and Alzheimer’s disease [20]. Furthermore, canard regimes reported
in this study are governed by the interactions between the pyramidal cell and SOM+ in-
terneuron subpopulations that follow a two-time-scale structure. It would be natural to
observe canard-mediated transitions in another generally used NMM of Jansen and Rit
[13] for modeling the brain activity. Hence, the canard-mediated fine structures we have
demonstrated here could be relevant for a number of situations and lead to markers of sub-
sequent critical transitions. The reported degenerate FSN II singularity leading to canard
trajectories is due to the general structure of NMMs, which are defined via second-order
differential equations. The dynamics associated with the degenerate FSN II singularity
merits further investigations and will be considered as a future work. Finally, organisation
of homoclinic canard orbits, possible codimension-two bifurcations and interactions with
the fold points will be studied in forthcoming works.

Acknowledgements

We would like to thank Prof. Fabrice Bartolomei (APHM, Timone Hospital, Clinical Neurophysiology, Marseille, France) for providing the clinical data, Mathieu Desroches (Inria Sophia Antipolis - Méditerranée, MathNeuro Team, Sophia Antipolis, France) for his helpful comments and referees for their useful suggestions.

Funding

EKE was supported by NIH (Application number: R01 NS092760-01A1; 18/02/2019 to 31/12/2020). She is member of the Galvani project (ERC-SyG 2019; grant agreement No 855109) since 01/01/2021.

Abbreviations

NMM, Neural Mass Model; EEG, Electroencephalography; LFP, Local Field Potential; GSPT, Geometric Singular Perturbation Theory; SOM+, Somatostatin positive; PV+, Parvanium positive; CA1, Cornu Ammonis 1; PSP, Post-Synaptic Potential; EPSP, Excitatory Post-Synaptic Potential; IPSP, Inhibitory Post-Synaptic Potential; SEEG, Stereoelectroencephalography; FSN II, Folded Saddle-Node type II; DSRS, Desingularized Slow Reduced System; LP, Limit Point; H, Hopf; HOM, Homoclinic; CP, Cusp Point; BT, Bogdanov–Takens; GH, Generalized Hopf; SNIC, Saddle-Node of Invariant Cycle; SNPO, Saddle-Node of Periodic Orbits.

Availability of data and materials

The codes used for numerical analysis are available from the GitHub database
(https://github.com/elifkoksal/NMM_BurstingDynamics).
**Ethics approval and consent to participate**

The SEEG recordings were carried out as part of normal clinical care of patients. Patients were informed that their data
may be used for research purposes.

**Competing interests**

The authors declare that they have no competing interests.

**Consent for publication**

Not applicable.

**Authors' contributions**

EKE developed the theoretical framework, conducted the mathematical analysis, interpreted the results and was a major contributor in writing the manuscript. FW contributed to interpreting the results and writing the manuscript. All authors read and approved the final manuscript.

**Authors' information**

Université de Rennes 1, INSERM, Laboratoire Traitement du Signal et de L'Image (LTSI) - U1099, Campus de Beaulieu - Batiment 22, 35042 Rennes, France.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Received: 16 September 2020 Accepted: 17 August 2021 Published online: 16 September 2021

**References**

1. Tukker JJ, Beed P, Schmitz D, Larkum ME, Sachdev R. Up and down states and memory consolidation across somatosensory, entorhinal, and hippocampal cortices. Front Syst Neurosci. 2020;14:22. https://doi.org/10.3389/fnsys.2020.00022.

2. Schulz DJ, Baines RA, Hempel CM, Li L, Liss B, Misonou H. Cellular excitability and the regulation of functional neuronal identity: from gene expression to neuromodulation. J Neurosci. 2006;26(41):10362-7. https://doi.org/10.1523/JNEUROSCI.3194-06.2006.

3. Ratté S, Hong S, De Schutter E, Prescott SA. Impact of neuronal properties on network coding: roles of spike initiation dynamics and robust synchrony transfer. Neuron. 2013;78(5):758-72. https://doi.org/10.1016/j.neuron.2013.05.030.

4. Brunel N, van Rossum MCW. Lapicque's 1907 paper: from frogs to integrate-and-fire. Biol Cybern. 2007;97(5):337-9. https://doi.org/10.1007/s00422-007-0190-0.

5. Lapique L. Recherches quantitatives sur l'excitation électrique des nerfs traitée comme une polarisation. J Physiol Pathol Gén. 1907;9:620-35.

6. Hodgkin AL. The local electric changes associated with repetitive action in a non-medullated axon. J Physiol. 1948;107(2):165-81. https://doi.org/10.1113/jphysiol.1948.sp004260.

7. Beraneck M, Idoux E. Reconsidering the role of neuronal intrinsic properties and neuromodulation in vestibular homeostasis. Front Neurol. 2012;3:25. https://doi.org/10.3389/fneur.2012.00025.

8. Sanabria ERG, Su H, Yaari Y. Initiation of network bursts by $Ca^{2+}$-dependent intrinsic bursting in the rat pilocarpine model of temporal lobe epilepsy. J Physiol. 2001;532(1):205-16. https://doi.org/10.1111/j.1469-7793.2001.0205g.x.

9. Jung S, Jones TD, Lugo JN, Sheerin AH, Miller JW, D'Ambrosio R, Anderson AE, Poolos NP. Progressive dendritic hcn channelopathy during epileptogenesis in the rat pilocarpine model of epilepsy. J Neurosci. 2007;27(47):13012-21. https://doi.org/10.1523/JNEUROSCI.3605-07.2007.

10. Shah MM, Anderson AE, Leung V, Lin X, Johnston D. Seizure-induced plasticity of h channels in entorhinal cortical layer III pyramidal neurons. Neuron. 2004;44(3):495-508. https://doi.org/10.1016/j.neuron.2004.10.011.

11. Freestone DR, Kuhlmann L, Grayden DB, Burkitt AN, Lai A et al. Electrical probing of cortical excitability in patients with epilepsy. Epilepsy Behav. 2011;22:110-8.

12. Chang W-C, Kudlacek J, Hlinka J, Chvojka J, Hadrava M, Kumpost V, Powell AD, Janca R, Murana MI, Karoly PJ, Freestone DR, Cook MJ, Palus M, Otahal J, Jefferys JGR, Jiruska P. Loss of neuronal network resilience precedes seizures and determines the ictogenic nature of interictal synaptic perturbations. Nat Neurosci. 2018;21:1742-52. https://doi.org/10.1038/s41593-018-0278-y.

13. Jansen BH, Rit VG. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol Cybern. 1995;73:357-66.

14. David O, Friston KJ. A neural mass model for MEG/EEG: coupling and neuronal dynamics. NeuroImage. 2003;20(3):1743-55. https://doi.org/10.1016/j.neuroimage.2003.07.015.

15. Ursino M, Cona F, Zavaglia M. The generation of rhythms within a cortical region: analysis of a neural mass model. NeuroImage. 2010;52(3):1080-94. https://doi.org/10.1016/j.neuroimage.2009.12.084.

16. Wendling F, Bartolomei F, Bellanger JJ, Chauvel P. Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition. Eur J Neurosci. 2002;15(9):1499-508. https://doi.org/10.1046/j.1460-9568.2002.01985.x.

17. Suffczynski P, Kalitzin S, Lopes Da Silva FH. Dynamics of non-convulsive epileptic phenomena modeled by a bistable neuronal network. Neuroscience. 2004;126(2):467-84. https://doi.org/10.1016/j.neuroscience.2004.03.014.

18. Molaee-Ardekani B, Benquet P, Bartolomei F, Wendling F. Computational modeling of high-frequency oscillations at the onset of neocortical partial seizures: from altered structure to dysfunction. NeuroImage. 2010;52(3):1109-22. https://doi.org/10.1016/j.neuroimage.2009.12.049.

19. Wendling F, Benquet P, Bartolomei F, Jirsa V. Computational models of epileptiform activity. J Neurosci Methods. 2016;260:233-51. https://doi.org/10.1016/j.jneumeth.2015.03.027.
20. Bhattacharya BS, Coyle D, Maguire LP. A thalamo-cortico-thalamic neural mass model to study alpha rhythms in Alzheimer's disease. Neural Netw. 2011;24(6):631-45. https://doi.org/10.1016/j.neunet.2011.02.009.

21. Liu F, Wang J, Liu C, Li H, Deng B, Fietkiewicz C, Loparo KA. A neural mass model of basal ganglia nuclei simulates pathological beta rhythm in Parkinson's disease. Chaos, Interdiscip J Nonlinear Sci. 2016;26(12):123113. https://doi.org/10.1063/1.4972200.

22. Liu C, Zhu Y, Liu F, Wang J, Li H, Deng B, Fietkiewicz C, Loparo KA. Neural mass models describing possible origin of the excessive beta oscillations correlated with Parkinsonian state. Neural Netw. 2017;88:65-73. https://doi.org/10.1016/j.neunet.2017.01.011.

23. Moehlis J. Canards for a reduction of the Hodgkin-Huxley equations. J Math Biol. 2006;52(2):141-53. https://doi.org/10.1007/s00285-005-0347-1.

24. Rubin J, Wechselberger M. Giant squid-hidden canard: the 3D geometry of the Hodgkin-Huxley model. Biol Cybern. 2007;97(1):5-32. https://doi.org/10.1007/s00422-007-0153-5.

25. Hasan CR, Krauskopf B, Osinga HM. Saddle slow manifolds and canard orbits in $R^4$ and application to the full Hodgkin-Huxley model. J Math Neurosci. 2018;8(1):5. https://doi.org/10.1186/s13408-018-0060-1.

26. Krupa M, Szmolyan P. Relaxation oscillation and canard explosion. J Differ Equ. 2001;174(2):312-68. https://doi.org/10.1006/jdeq.2000.3929.

27. Curtu R, Rubin J. Interaction of canard and singular Hopf mechanisms in a neural model. SIAM J Appl Dyn Syst. 2011;10(4):1443-79. https://doi.org/10.1137/110823171.

28. Köksal Ersöz E, Desroches M, Guillamon A, Rinzel J, Tabak J. Canard-induced complex oscillations in an excitatory network. J Math Biol. 2020. https://doi.org/10.1007/s00285-020-01490-1.

29. De Maesschalck P, Wechselberger M. Neural excitability and singular bifurcations. J Math Neurosci. 2015;5:16. https://doi.org/10.1186/s13408-015-0029-2.

30. De Maesschalck P, Desroches M. Numerical continuation techniques for planar slow-fast systems. SIAM J Appl Dyn Syst. 2013;12(3):1159-80. https://doi.org/10.1137/120877386.

31. Desroches M, Krupa M, Rodrigues S. Inflection, canards and excitability threshold in neuronal models. J Math Biol. 2013;67(4):989-1017. https://doi.org/10.1007/s00285-012-0576-z.

32. Desroches M, Freire E, Hogan SJ, Ponce E, Thota P. Canards in piecewise-linear systems: explosions and super-explosions. Proc R Soc A, Math Phys Eng Sci. 2013;469(2154):20120603. https://doi.org/10.1098/rspa.2012.0603.

33. Mitry J, McCarthy M, Kopell N, Wechselberger M. Excitable neurons, firing threshold manifolds and canards. J Math Neurosci. 2013;3(1):12. https://doi.org/10.1186/2190-8567-3-12.

34. Wechselberger M, Mitry J, Rinzel J. Canard theory and excitability. In: Kloeden PE, Pötzsche C, editors. Nonautonomous dynamical systems in the life sciences. vol. 2102. Cham: Springer; 2013. p. 89-132. https://doi.org/10.1007/978-3-319-03080-7_3.

35. Burke J, Desroches M, Barry AM, Kaper TJ, Kramer MA. A showcase of torus canards in neuronal bursters. J Math Neurosci. 2012;2(1):3. https://doi.org/10.1186/2190-8567-2-3.

36. Desroches M, Kaper TJ, Krupa M. Mixed-mode bursting oscillations: dynamics created by a slow passage through spike-adding canard explosion in a square-wave burster. Chaos, Interdiscip J Nonlinear Sci. 2013;23(4):046106. https://doi.org/10.1063/1.4827026.

37. Desroches M, Krupa M, Rodrigues S. Spike-adding in parabolic bursters: the role of folded-saddle canards. Phys D, Nonlinear Phenom. 2016;331:58-70. https://doi.org/10.1016/j.physd.2016.05.011.

38. Guckenheimer J, Kuehn C. Computing slow manifolds of saddle type. SIAM J Appl Dyn Syst. 2009;8(3):854-79. https://doi.org/10.1137/080741999.

39. Kramer MA, Traub RD, Kopell NJ. New dynamics in cerebellar Purkinje cells: torus canards. Phys Rev Lett. 2008;101(6):068103. https://doi.org/10.1103/PhysRevLett.101.068103.

40. Nowacki J, Osinga HM, Tsaneva-Atanasova K. Dynamical systems analysis of spike-adding mechanisms in transient bursts. J Math Neurosci. 2012;2(1):7. https://doi.org/10.1186/2190-8567-2-7.

41. Osinga HM, Tsaneva-Atanasova KT. Dynamics of plateau bursting depending on the location of its equilibrium. J Neuroendocrinol. 2010;22(12):1301-14. https://doi.org/10.1111/j.1365-2826.2010.02083.x.

42. Rubin J, Wechselberger M. The selection of mixed-mode oscillations in a Hodgkin-Huxley model with multiple timescales. Chaos, Interdiscip J Nonlinear Sci. 2008;18(1):015105. https://doi.org/10.1063/1.2789564.

43. Weigenand A, Costa MS, Ngo H-VV, Claussen JC, Martinetz T. Characterization of K-complexes and slow wave activity in a neural mass model. PLoS Comput Biol. 2014;10(11):1003923. https://doi.org/10.1371/journal.pcbi.1003923.

44. Krupa M, Vidal A, Desroches M, Clément F. Mixed-mode oscillations in a multiple time scale phantom bursting system. SIAM J Appl Dyn Syst. 2012;11:1458-98. https://doi.org/10.1137/110860136.

45. Vo T, Bertram R, Wechselberger M. Multiple geometric viewpoints of mixed mode dynamics associated with pseudo-plateau bursting. SIAM J Appl Dyn Syst. 2013;12(2):789-830. https://doi.org/10.1137/120892842.

46. Desroches M, Kirk V. Spike-adding in a canonical three-time-scale model: superslow explosion and folded-saddle canards. SIAM J Appl Dyn Syst. 2018;17(3):1989-2017. https://doi.org/10.1137/17M1143411.

47. Bensaid S, Modolo J, Merlet I, Wendling F, Benquet P. COALIA: a computational model of human EEG for consciousness research. Front Syst Neurosci. 2019;13:59. https://doi.org/10.3389/fnsys.2019.00059.

48. Köksal Ersöz E, Modolo J, Bartolomei F, Wendling F. Neural mass modeling of slow-fast dynamics of seizure initiation and abortion. PLoS Comput Biol. 2020;16(11):1008430. https://doi.org/10.1371/journal.pcbi.1008430.

49. Doedel EJ, Champneys A, Fairgrieve TF, Yu AB, Kuznetsov AP, Oldeman BE, Paffenroth RC, Sandstede B, Wang XJ, Zhang C. Auto-07p: continuation and bifurcation software for ordinary differential equations. 2007. http://cmvl.cs.concordia.ca/auto/.

50. Ermentrout B. Simulating, analyzing, and animating dynamical systems. Software, environments and tools. Philadelphia: SIAM; 2002. https://doi.org/10.1137/1.9780898718195.

51. Fenichel N. Geometric singular perturbation theory for ordinary differential equations. J Differ Equ. 1979;31(1):53-98. https://doi.org/10.1016/0022-0396(79)90152-9.

52. Cardin PT, Teixeira MA. Fenichel theory for multiple time scale singular perturbation problems. SIAM J Appl Dyn Syst. 2017;16(3):1425-52. https://doi.org/10.1137/16M1067202.

53. Szmolyan P, Wechselberger M. Canards in R3. J Differ Equ. 2001;177(2):419-53. https://doi.org/10.1006/jdeq.2001.4001.
54. Milik A, Szmolyan P. Multiple time scales and canards in a chemical oscillator. In: Jones CKRT, Khibnik AI, editors. Multiple-time-scale dynamical systems. The IMA volumes in mathematics and its applications. New York: Springer; 2001. p. 117–40. https://doi.org/10.1007/978-1-4613-0117-2_5.

55. Guckenheimer J. Return maps of folded nodes and folded saddle-nodes. Chaos, Interdiscip J Nonlinear Sci. 2008;18:015108. https://doi.org/10.1063/1.2790372.

56. Krupa M, Wechselberger M. Local analysis near a folded saddle-node singularity. J Differ Equ. 2010;248(12):2841–88. https://doi.org/10.1016/j.jde.2010.02.006.

57. Benoît E, Callot JF, Diener F, Diener M. Chasse au canard. Collect Math. 1981;31–32(1–3):37–119.

58. Dumortier F, Roussarie R. Canard cycles and center manifolds. Memories of the American mathematical society. vol. 121. Providence: Am. Math. Soc.; 1996.

59. Benoît E. Systèmes lents-rapides dans $\mathbb{R}^3$ et leur canards. Astérisque. 1983;2:109–10.

60. Wechselberger M. À propos de canards (Apropos canards). Trans Am Math Soc. 2012;364(6):3289–309. https://doi.org/10.1090/S0002-9947-2012-05575-9.

61. Wechselberger M. Existence and bifurcation of canards in $\mathbb{R}^3$ in the case of a folded node. SIAM J Appl Dyn Syst. 2005;4(1):101–39. https://doi.org/10.1137/030601995.

62. Krupa M, Popović N, Kopell N. Mixed-mode oscillations in three time-scale systems: a prototypical example. SIAM J Appl Dyn Syst. 2008;7(2):361–420. https://doi.org/10.1137/070688912.

63. Mitry J, Wechselberger M. Folded saddles and faux canards. SIAM J Appl Dyn Syst. 2017;16(1):546–96. https://doi.org/10.1137/15M1045065.

64. Albizuri JU, Desroches M, Krupa M, Rodrigues S. Inflection, canards and folded singularities in excitable systems: application to a 3D FitzHugh–Nagumo model. J Nonlinear Sci. 2020;30(6):3265–91. https://doi.org/10.1007/s00332-020-09650-9.

65. Wieczorek S, Ashwin P, Luke CM, Cox PM. Excitability in ramped systems: the compost-bomb instability. Proc R Soc A, Math Phys Eng Sci. 2011;467(2129):1243–69. https://doi.org/10.1098/rspa.2010.0485.

66. Shil'nikov L. On a new type of bifurcation of multidimensional dynamical systems. Sov Math Dokl. 1969;10.

67. Kuznetsov YA. Elements of applied bifurcation theory. 2nd ed. New York: Springer; 1998.

68. Rinzel J. Excitation dynamics: insights from simplified membrane models. Fed Proc. 1985;44(15):2944–6. https://doi.org/10.1186/s13408-015-0029-2.

69. Rinzel J, Ermentrout BG. Analysis of neural excitability and oscillations. In: Koch C, Segev I, editors. Methods in neuronal modeling. Cambridge: MIT Press; 1998.

70. Zhao Z, Huaguang G. Transitions between classes of neuronal excitability and bifurcations induced by autapse. Sci Rep. 2017;7:6760. https://doi.org/10.1038/s41598-017-07051-9.

71. Lagarde S, Buzori S, Trebuchon A, Carron R, Scavarda D, Milh M et al. The repertoire of seizure onset patterns in human focal epilepsies: determinants and prognostic values. Epilepsia. 2019;60(1):85–95. https://doi.org/10.1111/epi.14604.

72. Lopatina OL, Malinovskaya NA, Komleva YK, Gorina YV, Shuvaev AN, Olovyannikova RY, Belozor OS, Belova OA, Higashida H, Salmina AB. Excitation/inhibition imbalance and impaired neurogenesis in neurodevelopmental and neurodegenerative disorders. Rev Neurosci. 2019;30:807–20. https://doi.org/10.1515/revneuro-2019-0014.

73. Lopes da Silva FH, Blanes W, Kalitzin SN, Parra J, Suffczynski P, Velis DN. Dynamical diseases of brain systems: different routes to epileptic seizures. IEEE Trans Biomed Eng. 2003;50:540–8. https://doi.org/10.1109/TBME.2003.810703.

74. Trevelyan AJ, Sussillo D, Watson BO, Yuste R. Modular propagation of epileptiform activity: evidence for an inhibitory veto in neocortex. J Neurosci. 2006;26:12447–55. https://doi.org/10.1523/JNEUROSCI.2787-06.2006.

75. Timofeev I, Grenier F, Steriade M. The role of chloride-dependent inhibition and the activity of fast-spiking neurons during cortical spike wave electrographic seizures. Neuroscience. 2002;114:1115–32.

76. Timofeev I, Steriade M. Neocortical seizures: initiation, development and cessation. Neuroscience. 2004;123:299–336. https://doi.org/10.1016/j.neuroscience.2003.08.051.

77. Meisen C, Loddenkemper T. Seizure prediction and intervention. Neuropharmacology. 2020;172:107898.

78. Shamas M, Benquet P, Merlet I, Khalil M, El Falou W, Nica A, Wendling F. On the origin of epileptic high frequency oscillations observed on clinical electrodes. Clin Neurophysiol. 2018;129(4):829–41. https://doi.org/10.1016/j.clinph.2018.01.062.

79. Huneau C, Benquet P, Dieuset G, Biraben A, Martin B, Wendling F. Shape features of epileptic spikes are a marker of epileptogenesis in mice. Epilepsia. 2013;54(12):2219–27. https://doi.org/10.1111/epi.12406.

80. Lopes da Silva FH, Blanes W, Kalitzin N, Parra J, Suffczynski P, Velis DN. Epilepsies as dynamical diseases of brain systems: basic models of the transition between normal and epileptic activity. Epilepsia. 2003;44(12):72–83. https://doi.org/10.1111/j.0013-9580.2003.12005.x.

81. Meisel C, Kuehn C. Scaling effects and spatio-temporal multilevel dynamics in epileptic seizures. PLoS ONE. 2012;2:e30371. https://doi.org/10.1371/journal.pone.0030371.

82. Chizhov AV, Zefirov AV, Amakhin DV, Smirnova EY, Zaitsev AV. Minimal model of interictal and ictal discharges 'epileptor-2'. PLoS Comput Biol. 2018;14(5):1006186. https://doi.org/10.1371/journal.pcbi.1006186.

83. Ullah G, Wei Y, Dahlem MA, Wechselberger M, Schiff SJ. The role of cell volume in the dynamics of seizure, spreading depression, and anoxic depolarization. PLoS Comput Biol. 2015;11(8):1004414. https://doi.org/10.1371/journal.pcbi.1004414.

84. Hübel N, Dahlem MA. Dynamics from seconds to hours in Hodgkin-Huxley model with time-dependent ion concentrations and buffer reservoirs. PLoS Comput Biol. 2014;10(12):1003941. https://doi.org/10.1371/journal.pcbi.1003941.

85. Desroches M, Faugeras O, Krupa M. Slow-fast transitions to seizure states in the Wendling-Chauvel neural mass model. Opera Med Physiol. 2015;2(3–4):228–34.

86. Jafarian A, Freestone DR, Nešić D, Grayden DB. Slow-fast Duffing neural mass model. In: 41st annual international conference of the IEEE engineering in medicine and biology society (EMBC). 2019. p. 142–5.

87. Hebbink J, van Gils SA, Meijer HGE. On analysis of inputs triggering large nonlinear neural responses slow-fast dynamics in the Wendling neural mass model. Commun Nonlinear Sci Numer Simul. 2020;83:105103. https://doi.org/10.1016/j.cnsns.2019.105103.