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

Abstract

We report a detailed analysis on the emergence of bursting in a recently developed neural mass model that includes short-term synaptic plasticity. Neural mass models can mimic the collective dynamics of large-scale neuronal populations in terms of a few macroscopic variables like mean membrane potential and firing rate. The present one is particularly important, as it represents an exact meanfield limit of synaptically coupled quadratic integrate and fire (QIF) neurons. Without synaptic dynamics, a periodic external current with slow frequency $\varepsilon$ can lead to burst-like dynamics. The firing patterns can be understood using singular perturbation theory, specifically slow–fast dissection. With synaptic dynamics, timescale separation leads to a variety of slow–fast phenomena and their role for bursting becomes inordinately more intricate. Canards are crucial to understand the route to bursting. They describe trajectories evolving nearby repelling locally invariant sets of the system and exist at the transition between subthreshold dynamics and bursting. Near the singular limit $\varepsilon = 0$, we report peculiar jump-on canards, which block a continuous transition to bursting. In the biologically more plausible $\varepsilon$-regime, this transition becomes continuous and bursts emerge via consecutive spike-adding transitions. The onset of bursting is complex and involves mixed-type-like torus canards, which form the very first spikes of the burst and follow fast-subsystem repelling limit cycles. We numerically evidence the same mechanisms to be responsible for bursting emergence in the QIF network with plastic synapses. The main conclusions apply for the network, owing to the exactness of the meanfield limit.

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]. 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 and Fire (QIF) neuron, which represents a canonical model for Saddle-Node on Invariant Circle (SNIC) bifurcations and so-called type I excitability. 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 led 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. Various applications of the MPR model have 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 be straightforwardly 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 [38].

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 quiescent 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.

Fig. 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. (Color figure online)

In order to understand how these bursts emerge, we have to encapsulate two main aspects. First, in the upcoming Sect. 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 Sect. 3. As a next step, this generic methodology for timescale separated problems is applied to the present model in Sect. 4. Dissection is crucial for understanding the results in Sect. 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 Sect. 6. Finally, in Sect. 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 [39]. 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 nonzero 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 [40, 41] and Morris–Lecar model [42].

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 Eqs. (1).

$$\begin{aligned} \dot{V}_i&=V_i^2 + \eta _i + J r(t) + I_1(t) \end{aligned}$$

(1a)

$$\begin{aligned}&\text {if }V_i>V_{\mathrm{thresh}}: V_i \leftarrow V_{\mathrm{r}}\nonumber \\ r(t)&=\frac{1}{N}\sum _{j=1}^{N} \sum _{k: t_j^{(k)}<t} \delta (t-t_j^{(k)}) \end{aligned}$$

(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_{\mathrm{thresh}}$ at which the reset rule applies leading to a reset of the potential to $V_{\mathrm{r}}$.

When performing the thermodynamic limit $N\rightarrow \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_{\mathrm{thresh}}=-V_{\mathrm{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 Eqs. (2).

$$\begin{aligned} \dot{r}&=\frac{\Delta }{\pi } + 2rv \end{aligned}$$

(2a)

$$\begin{aligned} \dot{v}&=v^2 - (\pi r)^2 + J r+ {\bar{\eta }} + I_1(t) \end{aligned}$$

(2b)

Despite having a rather simple state-space structure, the MPR model in absence of synaptic dynamics 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 Sect. 3.

2.2 Neural mass model with short-term plasticity

Our present study investigates bursting patterns in an extended version of Eqs. (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}\hbox { ms}$ (depression timescale).

Facilitation, as opposed to depression, leads to enhanced PSPs and is related to 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}\hbox { ms}$.

We will focus on an implementation of STP into the model on macroscopic level (m-STP) as suggested in [20], in order to 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-neuron 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 Eqs. (3). For mean-field approximations taking into account STP on single-neuron level we refer to [43].

$$\begin{aligned} \dot{V}_i&=V_i^2 + \eta _i + Jux r +I_1(t) \end{aligned}$$

(3a)

$$\begin{aligned} \dot{x}&= \frac{1-x}{\tau _\mathrm {d}}-uxr \end{aligned}$$

(3b)

$$\begin{aligned} \dot{u}&= \frac{U_0-u}{\tau _\mathrm {f}}+U_0(1-u)r \end{aligned}$$

(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 Eqs. (4), in the following referred to as neural mass with STP (NMSTP), represents an exact meanfield limit of the QIF network given in Eqs. (3). The state variables are the firing rate r, mean membrane potential v, amount of resources x and utilization factor u.

$$\begin{aligned} \dot{r}&=\frac{\Delta }{\pi } + 2rv \end{aligned}$$

(4a)

$$\begin{aligned} \dot{v}&=v^2 - (\pi r)^2 + Ju x r + {\bar{\eta }} + I_1(t) \end{aligned}$$

(4b)

$$\begin{aligned} \dot{x}&= \frac{1-x}{\tau _\mathrm {d}}-uxr \end{aligned}$$

(4c)

$$\begin{aligned} \dot{u}&= \frac{U_0-u}{\tau _\mathrm {f}}+U_0(1-u)r \end{aligned}$$

(4d)

Despite the fact that Eqs. (4) already possesses multiple timescales via $\tau _\mathrm {d}$ and $\tau _\mathrm {f}$, the system will evolve on the fastest timescale of the full problem in presence of a slow periodic drive. However, as we will discuss later, the 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 Sects. 5.2 and 6.3).

2.3 Dynamics under constant forcing

Most of the parameters values used for Eqs. (4) will remain fixed and, if not stated differently, given in Table 1. Note that the unit of time is given by the membrane time constant $\tau _\mathrm {m}$. For more details on the numerical methods, we refer to Appendix A.

Table 1 Parameters and their values, which are fixed throughout this work, if not stated differently

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. 2a the resulting bifurcation diagram is displayed.

The extended model Eqs. (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 Sect. 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 consider 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).

Fig. 2

Solution families of the system with constant forcing: a Bifurcation diagram r versus $I_1$ of Eqs. (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_{\mathrm{L}}$, lower orange dot), creating a branch of unstable limit cycles (LC, orange dashed line). Two saddle-node bifurcations $F_{\mathrm{L}}$ and $F_{\mathrm{U}}$ (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 > rsim 0.7$ at a supercritical Hopf-Bifurcation ($H_{\mathrm{U}}$, upper orange dot). The unstable LC stabilizes (solid orange line) via a saddle-node bifurcation of cycles (purple dot) and vanishes at the second Hopf Bifurcation. The orange line marks the maximum firing rates of the LC branch. An enlargement of the bifurcation structure nearby $F_{\mathrm{L}}$ and $F_{\mathrm{U}}$ is shown in the inset, together with the unstable LC branch emerging at the Hopf point $H_{\mathrm{L}}$. c–e Periodic solution (r(t), v(t), x(t), u(t)) versus 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 Eqs. (4), the gray ones of the network Eqs. (3). b Spike scatter plot for 20,000 representative neurons in the network out of $N=100{,}000$. (Color figure online)

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_{\mathrm{L}}$) followed by two saddle-node (fold) bifurcation at $F_{\mathrm{L}}$ and $F_{\mathrm{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. 2a the instability persists throughout the S-shaped curve up until the upper supercritical Hopf bifurcation $H_{\mathrm{U}}$. The lower Hopf bifurcation $H_{\mathrm{L}}$ generates a family of unstable limit cycles that undergoes a fold bifurcation of cycles, giving rise to stable periodic solutions.

One of these trajectories (r(t), v(t), x(t), u(t)) is presented in Fig. 2b–f as a function of time. It is superimposed onto the corresponding variables calculated by simulating a QIF network governed by Eqs. (3) and consisting of $N=100{,}000$ 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 population 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. 2b, 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. 2a 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. 2b–f justifies to perform the upcoming analysis using solely the NMSTP. We will return to the implication of NMSTP dynamics for the network in Sect. 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}\hbox { ms}/\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 Eqs. (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.

$$\begin{aligned} \dot{I_1}&=\varepsilon g_1(I_1,I_2)=\varepsilon \left[ I_1(a-I_1^2-I_2^2)+I_2\right] \end{aligned}$$

(5a)

$$\begin{aligned} \dot{I_2}&=\varepsilon g_2(I_1,I_2)=\varepsilon \left[ I_2(a-I_1^2-I_2^2)-I_1\right] \end{aligned}$$

(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 Eqs. (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. (5) and (4).

To understand the impact of this slow forcing, it is advantageous to superimpose solutions of the full problem on the r versus $I_1$ bifurcation diagram of the unforced system, this is at the core of the slow–fast dissection introduced by Rinzel [32,33,34]. An example of a purely slow trajectory is shown in Fig. 3$\hbox {a}_1$–$\hbox {c}_1$ and labeled $\varvec{\gamma }_0(t)$. The firing rate r(t), shown in panel ($\hbox {b}_1$), increases and decreases following the same pattern as the forcing $I_1(t)$ in panel ($\hbox {a}_1$). Moreover, in $r-I_1$ projection, shown in panel ($\hbox {c}_1$), 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 understood as a quasi-static motion, the more complex solutions $\varvec{\gamma }_1(t)$ to $\varvec{\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.

Fig. 3

Typical solutions $\varvec{\gamma }_0(t)$ to $\varvec{\gamma }_4(t)$ of the full system: a Periodic forcing current $I_1(t)$ and b firing rate r(t) versus time t. c Same trajectories superimposed on the bifurcation diagram of the unforced system in r–$I_1$. In ($\hbox {c}_1$) the blue trajectory evolves towards larger $I_1$ values in the vicinity of the stable equilibrium branch and turns back nearby the Hopf bifurcation $H_{\mathrm{L}}$. The trajectory $\varvec{\gamma }_3(t)$ in ($\hbox {c}_2$) moves towards larger $I_1$ values after reaching the upper unstable equilibrium branch up to a maximum value, given by the forcing amplitude, before following the branch down again towards smaller $I_1$ values. Note that in column (3) the firing rate axis scale reaches until larger values than in columns (1, 2). The parameter values are as follows: $\varepsilon =10^{-5}$; $\varvec{\gamma }_0$: $A\approx 0.2487$; $\varvec{\gamma }_1$ to $\varvec{\gamma }_3$: in increasing order exponentially close to $A\approx 0.2507$; $\varvec{\gamma }_4$: $\varepsilon =10^{-3}$, $A\approx 0.2553$. (Color figure online)

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 Eqs. (6) and here referred to as full system,

$$\begin{aligned} \dot{{\mathbf {X}}}_\mathrm {f}&=~{\mathbf {F}}({\mathbf {X}}_\mathrm {f},{\mathbf {X}}_\mathrm {s}) \end{aligned}$$

(6a)

$$\begin{aligned} \dot{{\mathbf {X}}}_\mathrm {s}&=\varepsilon {\mathbf {G}}({\mathbf {X}}_\mathrm {f},{\mathbf {X}}_\mathrm {s}) \end{aligned}$$

(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 \;\rightarrow \; {\mathbb {R}}^k$ and ${\mathbf {G}}({\mathbf {X}}_\mathrm {f},{\mathbf {X}}_\mathrm {s}): {\mathbb {R}}^k \times {\mathbb {R}}^l \;\rightarrow \; {\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 Eqs. (7) by parametrizing it in slow time $\tau :=\varepsilon t$,

$$\begin{aligned} \varepsilon {\mathbf {X}}^\prime _\mathrm {f}&={\mathbf {F}}({\mathbf {X}}_\mathrm {f},{\mathbf {X}}_\mathrm {s}) \end{aligned}$$

(7a)

$$\begin{aligned} {\mathbf {X}}^\prime _\mathrm {s}&={\mathbf {G}}({\mathbf {X}}_\mathrm {f}, {\mathbf {X}}_\mathrm {s}). \end{aligned}$$

(7b)

The derivative with respect to slow time $\tau$ is denoted . The two representations Eqs. (6) and Eqs. (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 Eqs. (8). This limit can be used to understand the dynamics of the full system for which ${\mathbf {X}}_\mathrm {f}$ evolves fast and results in the following $k+l$ ODEs, l of which are trivial:

$$\begin{aligned} \dot{{\mathbf {X}}}_\mathrm {f}&={\mathbf {F}}({\mathbf {X}}_\mathrm {f},{\mathbf {X}}_\mathrm {s}) \end{aligned}$$

(8a)

$$\begin{aligned} \dot{{\mathbf {X}}}_\mathrm {s}&={\mathbf {0}} \end{aligned}$$

(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 Eqs. (7), yields the slow subsystem, namely:

$$\begin{aligned} {\mathbf {0}}&={\mathbf {F}}({\mathbf {X}}_\mathrm {f},{\mathbf {X}}_\mathrm {s}) \end{aligned}$$

(9a)

$$\begin{aligned} {\mathbf {X}}^\prime _\mathrm {s}&={\mathbf {G}}({\mathbf {X}}_\mathrm {f},{\mathbf {X}}_\mathrm {s}). \end{aligned}$$

(9b)

Eqs. (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:

$$\begin{aligned} S_0=\{({\mathbf {X}}_\mathrm {f},{\mathbf {X}}_\mathrm {s})~\vert ~{\mathbf {F}}({\mathbf {X}}_\mathrm {f},{\mathbf {X}}_\mathrm {s})=0\} , \end{aligned}$$

(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 Eqs. (6) with both slow and fast segments are $\varepsilon$-perturbations of singular orbits and the prototypical slow–fast cycles are relaxation oscillations [44]. The way these cycles emerge in parameter space is rather peculiar and involve the famous canard solutions, which we briefly review in the context of the classical van der Pol (VdP) system below.

3.1 Classical canards in the van der Pol oscillator

The VdP system with $0<\varepsilon \ll 1$ displays various dynamical regimes, from a stable equilibrium branch, via Hopf cycles and canards, to relaxation oscillations. The subthreshold regime, i.e, the branch of stable equilibria, terminates at a supercritical Hopf bifurcation, at which stable small Hopf cycles emanate. These cycles are not yet of relaxation type and only exist in an $O(\varepsilon )$ distance from the Hopf point. This regime 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 [45] and the associated canards separate the small Hopf cycles from relaxation oscillations [46]. They evolve for some time near a repelling branch of $S_0$, before jumping to one of the attracting branches. Two types of canards can be distinguished, the ones with a head and those without a head. After the canard explosion a wide regime of relaxation oscillation follows and terminates at a second (supercritical) Hopf bifurcation.

3.1.1 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 [47], 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$ [48, 49]. They describe orbits following a stable family of fast subsystem cycles and switching to the unstable one past the fold. For an overview of different types of canards we refer to [50].

A hybrid of classical canards and torus canards, so-called mixed-type canards, were reported in [51]. 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.

In the following, in order to understand the emergence of canards and bursting in the NMSTP under slow forcing, i.e., Eqs. (4) and (5), we will treat the problem within the slow fast framework. Moreover, the introduction of an auxiliary system referred to as desingularized reduced system (DRS) reveals that the presence of a so-called folded-saddle singularity is responsible for the emergence of canards and has influence on shaping bursting trajectories.

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,

$$\begin{aligned} {\mathbf {F}}({\mathbf {X}}_\mathrm {f},{\mathbf {X}}_\mathrm {s})&= \begin{pmatrix} \frac{\Delta }{ \pi }+2rv \\ v^2 + J uxr - (\pi r)^2 +I_1 \\ (1-x)/\tau _\mathrm {d}-uxr \\ (U_0-u)/\tau _\mathrm {f}+U_0(1-u)r \end{pmatrix} \end{aligned}$$

(11)

$$\begin{aligned} {\mathbf {G}}({\mathbf {X}}_\mathrm {s})&= \varepsilon \begin{pmatrix} g_1(I_1,I_2) \\ g_2(I_1,I_2) \end{pmatrix} \nonumber \\&=\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}. \end{aligned}$$

(12)

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

$$\begin{aligned} {\mathbf {F}}({\mathbf {X}}_\mathrm {f},{\mathbf {X}}_\mathrm {s})&={\mathbf {0}}. \end{aligned}$$

(13)

We can therefore already infer the shape of $S_0$. It corresponds to the cartersian product $S^*\times \{I_2\vert I_2\in {\mathbb {R}}\}$, where $S^*$ denotes the S-shaped branch of equilibria of the fast subsystem and is given in Eq. (14),

$$\begin{aligned} S^*:=\{(r,v,x,u,I_1)\vert {\mathbf {F}}({\mathbf {X}}_\mathrm {f},{\mathbf {X}}_\mathrm {s})=0\}. \end{aligned}$$

(14)

Hence the associated fold set ${\mathcal {F}}$ has two 1D connected components, namely the two lines ${\mathcal {F}}_{\mathrm{U}}:=\{F_{\mathrm{U}}\} \times \{I_2\vert I_2\in {\mathbb {R}}\}$ and ${\mathcal {F}}_{\mathrm{L}}:=\{F_{\mathrm{L}}\} \times \{I_2\vert 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}}_{\mathrm{L}}$ and ${\mathcal {F}}_{\mathrm{U}}$, can be seen in Fig. 4, together with the critical manifold $S_0$ in $I_1$–$I_2$–r projection and the trajectories $\varvec{\gamma }_0(t)$ to $\varvec{\gamma }_4(t)$ of Fig. 3.

Fig. 4

Slow–fast dissection and critical manifold: Solutions $\varvec{\gamma }_0(t)$ to $\varvec{\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. (Color figure online)

4.1 Singular dynamics: fast subsystem

Essentially, since $S_0$ consists of fast subsystem equilibria, the bifurcation diagram in Fig. 2a is a projection of $S_0$ onto the $I_1$–r plane and we can associate local stability properties of the fast subsystem with the critical manifold. Stable (unstable) parts of the equilibrium branch in Fig. 2a 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 the fast dynamics dominates and (r, v, x, u) will be repelled and attracted accordingly. The stability changes along the set of Hopf bifurcation given by ${\mathcal {H}}_{\mathrm{L}}:=\{H_{\mathrm{L}}\} \times \{I_2\Vert I_2\in {\mathbb {R}}\}$ and ${\mathcal {H}}_{\mathrm{U}}:=\{H_{\mathrm{U}}\} \times \{I_2\Vert I_2\in {\mathbb {R}}\}$ . Attracting (repelling) parts of $S_0$ are marked as green (light green) surfaces throughout this chapter, as can be seen in Fig. 4.

We can already exploit this in order to examine the full system dynamics further. To start with, $\varvec{\gamma }_0(t)$ in Fig. 4$\hbox {a}_1$ evolves entirely in the vicinity of the lower attracting sheet of the critical manifold $S_0$. This means that the dynamics is purely slow and in approximation bounded to $S_0$. For $\varvec{\gamma }_1(t)$ we see the typical canard dynamics, since it follows the middle repelling part of $S_0$ for some time, before jumping to the bottom attracting sheet. Moreover, $\varvec{\gamma }_2(t)$ in Fig. 4$\hbox {a}_2$ resembles the behavior of a canard with a head. However, while in the VdP system the upper branch of the fast subsystem is attracting, we find the same behavior here, but with a repelling upper sheet of $S_0$. We will refer to this surprising solution of the full system as jump-on canard. Similarly, $\varvec{\gamma }_3(t)$ exhibits a jump to the upper sheet, but in absence of a canard segment along the middle sheet. Finally, the bursting solution $\varvec{\gamma }_4(t)$ in Fig. 4$\hbox {a}_3$ has segments along the bottom sheet of $S_0$ and displays a short canard segment along the repelling middle sheet, before it pierces through the critical manifold $S_0$ and starts to burst. In order to understand the emergence of canards, jump-on canards and bursting, we will proceed with the analysis of the slow subsystem in the following.

4.2 Singular dynamics: slow subsystem

4.2.1 Slow flow

We make further use of the dissection by studying the slow flow on the critical manifold $S_0$. In the slow subsystem Eqs. (9) the state space is reduced to $S_0$, described by four algebraic conditions in Eq. (13), 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. (12). 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. (15),

(15)

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

(16)

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. (11) and (12), respectively. Values of ${\mathbf {X}}_\mathrm {f}$ on $S_0$ are determined via the algebraic conditions Eq. (13), 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 Eqs. (17)

$$\begin{aligned} r^\prime&=g_1(I_1,I_2)Ar(1+ru\tau _\mathrm {d})(1+rU_0\tau _\mathrm {f})/D \end{aligned}$$

(17a)

$$\begin{aligned} I_1^\prime&=g_1(I_1,I_2) \end{aligned}$$

(17b)

$$\begin{aligned} I_2^\prime&=g_2(I_1,I_2) \end{aligned}$$

(17c)

We can see the relevance of the denominator D, given in Eq. (18), by noting that $D=0$ defines the fold set of $S_0$ at which the fast subsystem Jacobian $\partial {\mathbf {F}}/\partial {\mathbf {X}}_\mathrm {f}$ is singular:

$$\begin{aligned} D= 2&\left[ (\pi r)^2+v^2\right] (r \tau _\mathrm {d}u+1) (r \tau _\mathrm {f}U_0+1) \nonumber \\&-Jx r (r \tau _\mathrm {f}U_0+u). \end{aligned}$$

(18)

Saddle-node bifurcations of the fast subsystem are characterized by the identical condition $\det (\partial {\mathbf {F}}/\partial {\mathbf {X}}_\mathrm {f})=0$ and are equivalent to the fold curves ${\mathcal {F}}_{\mathrm{L}}$ and ${\mathcal {F}}_{\mathrm{U}}$ shown in Fig. 4. Therefore, the slow flow is undefined along these lines and the slow subsystem fails to describe the slow dynamics. This is for example relevant for the canard trajectories $\varvec{\gamma }_1(t)$ and $\varvec{\gamma }_2(t)$ in Figs. 3 and 4, which projected onto the $(r,I_2)$ plane cross ${\mathcal {F}}_{\mathrm{L}}$ transversely.

4.2.2 Desingularization

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

(19a)

(19b)

(19c)

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}}_{\mathrm{L}}$ and ${\mathcal {F}}_{\mathrm{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}}_{\mathrm{L}}$ and ${\mathcal {F}}_{\mathrm{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 Eqs. (17). 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 Eqs. (19) 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.

4.2.3 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 Eqs. (17). In the DRS (Eqs. (19)) 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. 5a 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.

Fig. 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 Eqs. (17)]. The black dots $p_0$ (unstable focus), $p_1$ (saddle), $p_2$ (center) denote equilibria of the DRS Eqs. (19). 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}}_{\mathrm{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 $\varvec{\gamma }_2(t)$ and $\varvec{\gamma }_3(t)$ are solutions of the full system. The objects $p_1,p_2$ and ${\mathcal {M}}_{\mathrm{FS}}$ depend on the choice of A, here they correspond to the value used to obtain $\varvec{\gamma }_3(t)$. Other parameters values are as in Figs. 3 and 4. (Color figure online)

Located on the bottom sheet of $S_0$, $p_0$ results from the Hopf form given in Eqs. (5) and is an unstable focus at ($I_1,I_2)=(0,0)$. On the other hand, $p_2$ lays on the upper fold line ${\mathcal {F}}_{\mathrm{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}}_{\mathrm{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}}_{\mathrm{FS}}={\mathcal {M}}_{\mathrm{FS}}^{\mathrm{L}}\cup {\mathcal {M}}_{\mathrm{FS}}^{\mathrm{U}}$ forms, consisting of two parts, which are connected via $p_1$; see the orange curve in Fig. 5a, b. The connection ${\mathcal {M}}_{\mathrm{FS}}^{\mathrm{L}}$ is located on the lower sheet, while ${\mathcal {M}}_{\mathrm{FS}}^{\mathrm{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}}_{\mathrm{L}}$ and ${\mathcal {F}}_{\mathrm{U}}$. Accordingly, the folded saddle $p_1$ has significant influence on the dynamics of the slow subsystem along ${\mathcal {M}}_{\mathrm{FS}}$, as described in the following.

• (i) First of all, trajectories in the DRS evolving on ${\mathcal {M}}_{\mathrm{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}}_{\mathrm{L}}$ trajectories are attracted to and above repelled from $p_1$.

• (iii) The double homoclinic connection ${\mathcal {M}}_{\mathrm{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 [52].

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

As a consequence of the previous properties (i)–(iv) and given a specific value of A a solution of the slow subsystem exists, that evolves along the folded homoclinic ${\mathcal {M}}_{\mathrm{FS}}$ below the bottom fold line ${\mathcal {F}}_{\mathrm{L}}$ and extends, while remaining on $S_0$, beyond the folded-saddle $p_1$ until the upper fold ${\mathcal {F}}_{\mathrm{U}}$. This solution of the slow subsystem represents a canard segment in the context of the full problem. More specifically it is associated with the maximal canard introduced below.

4.2.4 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}}_{\mathrm{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 $\varvec{\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_{\mathrm{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}}_{\mathrm{FS}}^{\mathrm{U}}$ is disconnected at the intersections with ${\mathcal {F}}_{\mathrm{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 $\varvec{\gamma }_0(t)$ to $\varvec{\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 [53]. 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 $(\partial {\mathbf {F}}/\partial {\mathbf {X}}_\mathrm {f})\vert _{S_0}$ to have non-zero real part [54]. The theorem can therefore not be applied on ${\mathcal {F}}_{\mathrm{L}}$ and ${\mathcal {F}}_{\mathrm{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, $\varvec{\gamma }_1(t)$ to $\varvec{\gamma }_3(t)$ exemplify such connection numerically: segments that evolve close to the bottom sheet of $S_0$ connect to segments close to 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 [46], matched asymptotics [55] and so-called blow-up techniques [56]. 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 $\varvec{\gamma }_1(t)$ in Figs. 3 and 4, where the forcing amplitude A is large enough to surpass the lower Hopf bifurcation $H_{\mathrm{L}}$ and the lower fold $F_{\mathrm{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 Sect. 6.

The trajectory $\varvec{\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. Moreover, the jump-on dynamics can occur after a regular canard segment, as for $\varvec{\gamma }_2$ or independent of that, like for $\varvec{\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 detail in Sect. 6.3.

The two solutions $\varvec{\gamma }_2(t)$ and $\varvec{\gamma }_3(t)$ are shown in Figs. 5b and 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 $\varvec{\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}}_{\mathrm{L}}$ and reaches the lower fold curve ${\mathcal {F}}_{\mathrm{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. 2a. Instead of entering a period of bursting, the cyan trajectories approaches the upper branch of unstable equilibria. As soon as it jumps onto $S_0$, the slow subsystem becomes a valid limit anew and the curve remains on $S_0$ until it reaches ${\mathcal {F}}_{\mathrm{U}}$, where it jumps down to the stable sheet.

In the second case, $\varvec{\gamma }_2$ in Fig. 5b and 6 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}}_{\mathrm{U}}$ and finally jumps down. The global motion is identical to that of a canard with a head in VdP.

Fig. 6

Folded-saddle and jump-on canards: $S_0$ in $I_1$–v–r space (green curve). The gray surface ${\mathcal {M}}_{rv}$ is defined by $r=-\frac{\Delta }{2\pi v}$. The curves $\varvec{\gamma }_{2,3}(t)$ are solutions of the full system and their jump-on points $\varvec{\gamma }_{2,3}^*$ on the upper repelling sheet of $S_0$ are marked by dots. Parameters values are as in Fig. 5. (Color figure online)

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 Eqs. (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}}_{\mathrm{U}}$, where the slow flow is undefined. Here the trajectory can be joined to a fast bit which connects from ${\mathcal {F}}_{\mathrm{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 $\varvec{\gamma }_3$ crosses ${\mathcal {F}}_{\mathrm{L}}$ far from $p_1$ and the slow segment on the attracting sheet stops; $\varvec{\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 $\varvec{\gamma }_2$ and $\varvec{\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 $\varvec{\gamma }^*_{2,3}=(r_{2,3}^*,v_{2,3}^*, x_{2,3}^*, u_{2,3}^*)$ of the jump-on canard solutions $\varvec{\gamma }_2$ and $\varvec{\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 $\varvec{\gamma }^*_{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. (20) and obtained by setting ${\dot{r}}=0$ in Eq. (4a). On one hand, for the slow pieces of the curve this observation is as expected, since by definition this conditions holds on $S_0$:

$$\begin{aligned} \dot{r}=0 \Rightarrow r=-\frac{\Delta }{2\pi v} \text {\quad for } v\ne 0. \end{aligned}$$

(20)

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 Eqs. (4) 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, resulting in Eqs. (21), and approximates the intermediate scale dynamics of jump-on canards, which takes place in the vicinity of ${\mathcal {M}}_{rv}$.

$$\begin{aligned} 0&=\frac{\Delta }{\pi } + 2rv \end{aligned}$$

(21a)

$$\begin{aligned} \dot{v}&=v^2 - (\pi r)^2 + Ju x r + {\bar{\eta }} + I_1\end{aligned}$$

(21b)

$$\begin{aligned} \dot{x}&= \frac{1-x}{\tau _\mathrm {d}}-uxr \end{aligned}$$

(21c)

$$\begin{aligned} \dot{u}&= \frac{U_0-u}{\tau _\mathrm {f}}+U_0(1-u)r \end{aligned}$$

(21d)

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^*, x^*, u^*)$ [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^*, x^*, u^*)$ are of saddle type with eigenvalues $\lambda _1 \gg \lambda _2\gg -\lambda _3 > 0$ and therefore have associated 1D stable manifolds ${\mathcal {M}}_{\mathrm{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$.

In the singular limit one can concatenate these solutions of Eqs. (21) with the regular canard segments that populate the middle sheet of $S_0$ and thereby construct singular jump-on canards. Overall, 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. For $0<\varepsilon \ll 1$ the singular jump-on canards persists as $\varepsilon$ perturbations, analogous to the classical 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 $\varvec{\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. 2a, 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 [57,58,59], bursting can effectively be initiated at that fold; see e.g. Fig. 7($\hbox {b}_4$). 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 cannot 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 [52]; 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 [60]. 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

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 $\varvec{\gamma }_0(t)$ in Figs. 3 and 4, and is continued towards larger amplitudes.

Fig. 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 versus 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$. The dashed vertical line is located at $A=A^*=0.25531851205$. b Solutions $(r(t),v(t),x(t),v(t),I_1(t),I_2(t))$ in $r-I_1$ 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 orange surface (wireframe) represents the family of stable (unstable) limit cycles of the fast subsystem [orange branch in column (b)]. Note that here the x-axis is inverted. In (b) the black dots denote ${\mathcal {F}}_{\mathrm{L}}$ and ${\mathcal {F}}_{\mathrm{U}}$, the orange dot ${\mathcal {H}}_{\mathrm{L}}$, while the black dots in (c) show the folded singularity $p_1$, assuming $A=A^*$. Spikes contributing to $n_\mathrm {s}$ are marked with blue dots. The dashed red curves show the solutions of the panels above. The A values in (b–d) are in increasing order from top to bottom and near $A^*$. Full system solutions obtained at $\varepsilon =10^{-3}$. (Color figure online)

As a solution measure the $L_2$-norm of this family is plotted vs. A in Fig. 7a. The first part until $A=A^*\approx 0.25531851205$ is in the subthreshold regime. Around $A^*$ 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^*$.

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^*$ 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_{\mathrm{L}}$. Unstable periodic solutions of this branch will be denoted by $\Gamma ^{\mathrm{r}}$, stable ones after the fold of cycles by $\Gamma ^{\mathrm{a}}$. Embedding these solutions $\Gamma ^{\mathrm{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}}^{\mathrm{a}}\cup {\mathcal {P}}^{\mathrm{r}}$, which consists of attracting and repelling parts ${\mathcal {P}}^{\mathrm{a,r}} = \{(\Gamma ^{\mathrm{a,r}}(I_1,t), I_1)~\vert ~t\in [0,T(I_1)]\} \times \{I_2~\vert ~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$.

6.1.1 Onset of fast oscillations $\varvec{(\mathrm{b}_1, \mathrm{c}_1,\mathrm{d}_1)}$

At the smallest of the chosen A values near $A^*$ 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 $(\mathrm{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 $\mathrm{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 ^{\mathrm{a, r}}$ into account, a remarkable feature of the dynamics can be seen in $(I_1, v, d)$ space (panel ($\hbox {d}_1$)). The spikes of the burst appear to follow the family of unstable limit cycles, thus evolving near the repelling surface ${\mathcal {P}}^{\mathrm{r}}$. Unexpectedly, after piercing through the critical manifold, the bursting solution stays in the proximity of ${\mathcal {P}}^{\mathrm{r}}$, instead of being repelled from it. As $I_1$ slowly drifts towards smaller values, the trajectory remains close to ${\mathcal {P}}^{\mathrm{r}}$. These windings around ${\mathcal {P}}^{\mathrm{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 ($\hbox {b}_1$) 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^*$.

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

Fig. 8

Emergence of bursting: Enlargement of Fig. 7($\mathrm{d}_3$). Critical manifold in $(I_1,v,x)$-space with attracting (repelling) sheets of $S_0$ visualized as a green solid (dotted) line; the orange surface (wireframe) represents the family of stable (unstable) limit cycles of the fast subsystem, emerging from the lower Hopf bifurcation $H_{\mathrm{L}}$ (orange dot). The blue and red curves show solutions of the full system near the canard explosion. Forcing amplitude for the blue trajectory ($A\approx 0.25531851417$) is larger than for red one ($A\approx 0.2553185136$) and both exponentially close to the canard explosion at $A^*$. Solutions obtained for $\varepsilon =10^{-3}$. (Color figure online)

6.1.2 Explosivity and spike-adding $\varvec{(\mathrm{b}_2, \mathrm{c}_2,\mathrm{d}_2)}$

At a slightly larger A value, close to $A^*$, 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 ($\mathrm{b}_2, \mathrm{c}_2,\mathrm{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 ($\mathrm{b}_2, \mathrm{c}_2,\mathrm{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}}^{\mathrm{r}}$, thus adding an additional layer of sensitivity.

Compared to the previous case in Fig. 7($\mathrm{b}_1, \mathrm{c}_1,\mathrm{d}_1$), the full system solution winds around ${\mathcal {P}}^{\mathrm{r}}$ more often until smaller values of $I_1$, before jumping back to $S_0$. The trajectory essentially remains close to ${\mathcal {P}}^{\mathrm{r}}$, but reaches up higher. This way, by passing through $A^*$, 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}}^{\mathrm{r}}$ increases, as it winds around it, indicating some extent of repulsion near the surface; see also panel ($\hbox {b}_2$). This way, the full dynamics starts to escape from ${\mathcal {P}}^{\mathrm{r}}$ and gets attracted to ${\mathcal {P}}^{\mathrm{a}}$.

6.1.3 Emergence of bursting $\varvec{(\mathrm{b}_3, \mathrm{c}_3,\mathrm{d}_3)}$

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

6.1.4 Bursting $\varvec{(\mathrm{b}_4, \mathrm{c}_4,\mathrm{d}_4)}$

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

Here, an additional spike-adding mechanism beyond the critical value $A^*$ acts and is related to the period $T(I_1)$ of ${\mathcal {P}}^{\mathrm{a}}$ as a function of $I_1$. When A is beyond the canard explosion at $A^*$ the full system dynamics approach ${\mathcal {P}}^{\mathrm{a}}$ right after surpassing the lower fold $F_{\mathrm{L}}$, without the excursion on ${\mathcal {P}}^{\mathrm{r}}$. As the amplitude A is increased, this attraction to ${\mathcal {P}}^{\mathrm{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}}^{\mathrm{a}}$. On top of this, the fraction of time for which the trajectory remains near ${\mathcal {P}}^{\mathrm{a}}$ increases. Both effects add more spikes to the burst and result in the spike-adding arches observed in Fig. 7a.

The presented results already show the complexity of how bursts are generated, that is, via a transition through the canard explosion at $A=A^*$, which rather surprisingly leads the canard segment to evolve around the repelling object ${\mathcal {P}}^{\mathrm{r}}$. This is followed by spike growth via repulsion from ${\mathcal {P}}^{\mathrm{r}}$ and attraction to ${\mathcal {P}}^{\mathrm{a}}$, until all oscillations evolve near ${\mathcal {P}}^{\mathrm{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}}^{\mathrm{r}}$. This clearly resembles mixed-type canards as described in [51]. Very similar MTTCs have been reported in [47], where they evidently arise in the singular limit $\varepsilon \rightarrow 0$. In particular Fig. 2 of [47] 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, they are only observed for intermediate timescale separations.

Fig. 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^*$. For all branches (br.) apart from br. 4, the $A^*$ 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 \times 10^{-4}$. (Color figure online)

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}}_{\mathrm{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}}^{\mathrm{r}}$; see Fig. 8.

Secondly, the solutions $\Gamma ^{r} \subset {\mathcal {P}}^{\mathrm{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 Sect. 5.2), form due to the intrinsic timescales of the fast subsystem. Consequently, this means that stable manifolds can be associated with ${\mathcal {P}}^{\mathrm{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}}^{\mathrm{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}}^{\mathrm{r}}$; the large $\varepsilon$ allows to bridge $S_0$ to ${\mathcal {P}}^{\mathrm{r}}$ and finally torus canard dynamics allow to follow ${\mathcal {P}}^{\mathrm{a}}$ closely, adding more and more spikes in the transition region around $A^*$.

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^*$. They also show how the subsequent spike-adding process, reflected by the arches in Fig. 7a, 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. Nevertheless bursting can be observed and possibly emerges in a discontinuous manner in parameter space.

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\times 10^{-3}$ to $\varepsilon = 1\times 10^{-5}$ and they are aligned to the canard explosion at $A=A^*$ (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^*$, take different paths. For $\varepsilon > rapprox 2.7\times 10^{-4}$ the branches display a continuous transition from subthreshold oscillations to bursting; see red curves in Fig. 9. For $\varepsilon \lessapprox 2.7\times 10^{-4}$ on the other hand the branches evolve into families of jump-on canards (cyan curves in Fig. 9). A value of $\varepsilon \approx 2.667521298 \times 10^{-4}$ separates the two $\varepsilon$ regimes.

However, it is clear by considering branches 4 and 5 with identical $\varepsilon$, 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\times 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\times 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 Sect. 5.2). They block a transition towards ${\mathcal {P}}$ and as a consequence, bursting remains absent for these solutions families.

For $\varepsilon > rapprox 2.7\times 10^{-4}$ however, jump-on canards cease to exist. In the amplitude regime where they would be expected, the system approaches ${\mathcal {P}}^{\mathrm{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 Eqs. (4) is an exact limit of the underlying QIF network Eqs. (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. 2b–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,7,8, 20, 61, 62]. 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. We suspect the additional state variables and their contributions to the non-linearity of the system to be responsible for this.

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 Eqs. (3) and investigate the QIF population dynamics nearby the canard explosion of the neural mass.

Fig. 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 versus the forcing amplitude $A-A^*$. 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 Eqs. (3) with $N=100{,}000$ neurons and in presence of sinusoidal forcing $I_1(t)=A\sin (\varepsilon t)$. In ($\hbox {b}_1$,$\hbox {c}_1$) the forcing amplitude is given by $A=0.265$; in ($\hbox {b}_2$,$\hbox {c}_2$) by $A=0.270$. Row (b) shows the time series of the solution in terms of x(t) versus 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 (orange wireframe and surface). The orange dot denotes ${\mathcal {H}}_{\mathrm{L}}$, the black dots ${\mathcal {F}}_{\mathrm{L}}$ and ${\mathcal {F}}_{\mathrm{U}}$. The x-axis is inverted in (c)

In Fig. 10a 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. 10b, 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^*, v^*,x^*,u^*)$ 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^*\quad \text {for all }i=1,\dots ,N$, $x=x^*$ and $u=u^*$. 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. 10a. 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($\hbox {b}_1$, $\hbox {c}_1$). They mainly evolve in the proximity of the invariant set ${\mathcal {P}}^{\mathrm{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^*$.

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_{\mathrm{L}}$ is evident. This is clear by looking at the black curve in Fig. 10($\hbox {c}_1$). More strikingly, a few low-amplitude oscillations are picked up in the vicinity of the repelling invariant set ${\mathcal {P}}^{\mathrm{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 Sect. 6.

At the larger amplitude $A=0.270$, shown in Fig. 10($\hbox {b}_2$, $\hbox {c}_2$), 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

We made 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.

Through the inclusion of synaptic dynamics, in form of STP, we employ a more biologically plausible model in this work than the original MPR model, while preserving the exactness of the meanfield limit with respect to the QIF network. This adds more details to the QIF network and neural mass, but at the same time complexifies 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.

Appendix A: Numerical methods

The results in Figs. 2a and 3–10 are obtained using the continuation software AUTO-07p [63], to numerically calculate solutions and solution branches of the fast subsystem Eq. (11) and full system Eqs. (11) and (12)

For Figs. 1, 2b–f, and 10 direct simulations of the QIF network with STP Eqs. (3) were performed, as well as simulations of the neural mass with STP Eq. (11). The latter was integrated using Python and the adaptive Dormand–Prince method [64]. The QIF network simulations use values of the Lorentzian distributed constant currents $\eta _i$, that are deterministically set as given below.

$$\begin{aligned} \eta _i = {\bar{\eta }}+\Delta \tan \left( \frac{\pi (2i-N-1)}{2(N+1)}\right) \end{aligned}$$

(A1)

Furthermore the numerical integration requires simplifications of the problem, in order to maintain reasonable computation times. For this we followed the method suggested in [5]. First of all, the network is integrated using the Euler scheme, with a timestep $dt=10^{-4}$. Second, the limit $V_{\mathrm{thresh}}\rightarrow \infty$ is not realizable numerically. However, one can approximate the time $T_{\infty }$ required for the membrane potential of a QIF neuron to evolve from a value $V_{\mathrm{max}}$ to infinity as $T_{\infty }\approx \frac{1}{V_{\mathrm{max}}}$, given that the total input I to that neuron fulfills $\sqrt{I}\ll V_{\mathrm{max}}$. The time from negative infinity to $V_i=-V_{\mathrm{max}}$ is given by the same expression. As soon as a neuron crosses the threshold $V_i>V_{\mathrm{max}}$ it is reset to $V_i=-V_{\mathrm{max}}$ and enters into a refractory period of $2T_{\infty }=2/V_{\mathrm{max}}$, for which its dynamics is deactivated. This is to account for the integration from $V_i=V_{\mathrm{max}}$ to infinity and from negative infinity to $V_i=-V_{\mathrm{max}}$.

A spike is registered at half of this refractory period, via an instantaneous change of all membrane potentials and the synaptic variables, as stated below.

$$\begin{aligned} V_i(t^+)&=\frac{J}{N} u(t^-) x(t^-) \end{aligned}$$

(A2)

$$\begin{aligned} x(t^+)&=x(t^-)-\frac{1}{N} u(t^-)x(t^-) \end{aligned}$$

(A3)

$$\begin{aligned} u(t^+)&=u(t^-)+\frac{1}{N} U_0(1-u(t^-)) \end{aligned}$$

(A4)

For clarity the time before the increment is denoted by $t^-$, the time after by $t^+$.