CC BY
83
Slow-Fast Transitions to Seizure States in the Wendling-Chauvel Neural Mass Model
Mathieu Desroches*1, Olivier Faugeras2, and Martin Krupa1
1 MathNeuro team, Inria Sophia-Antipolis - Mediterranee Research Centre, 06902 Sophia Antipolis Cedex, France;
2 MathNeuro and TOSCA teams, Inria Sophia-Antipolis - Mediterranee Research Centre, Sophia Antipolis Cedex, France. * Corresponding e-mail: [email protected]
Abstract. We revisit the Wendling-Chauvel neural mass model by reducing it to eight ODEs and adding a dierential equation that accounts for a dynamic evolution of the slow inhibitory synaptic gain. This allows to generate dynamic transitions in the resulting nine-dimensional model. The output of the extended model can be related to EEG patterns observed during epileptic seizure, in particular isolated pre-ictal spikes and low-voltage fast oscillations at seizure onset. We analyse the extended model using basic tools from slow-fast dynamical systems theory and relate the main transitions towards seizure states to torus canards, a type of solutions that has been shown to explain the spiking to bursting transition in many neural models. We nd that the original ten-dimensional Wendling-Chauvel model can be reduced to eight dimensions, two variables being scaled versions of two other variables of the model. We then obtain a model with four PSP blocks, which is consistent with the block-diagrams typically presented to describe this model. Instead of varying the slow inhibitory synaptic gain parameter B quasi-statically, or just performing numerical bifurcation analysis in B as the structure of the fast subsystem of an hypothetical extended system, we construct a true slow dynamics for B, depending sensitively on the main PSP output of the model, Y0. Near fold bifurcation of limit cycles of the original model, the solution to the extended model performs fast low-amplitude oscillations close to both attracting and repelling branches of limit cycles, which is the signature of a torus canard phenomenon.
Keywords: Neural mass models, Epileptic seizure, Inhibition, Slow-fast dynamics, Torus canards
Background
Neural Mass Models
Neural-mass models were initiated by Freeman (Freeman, 1972, 1975) and Lopes da Silva (Lopes Da Silva et al, 1974). These are population models with lumped parameters, aimed at capturing the dynamics of EEG recordings. Jansen and Rit (Jansen and Rit, 1995) have considered the coupling of two neural-mass models, corresponding to two neuronal populations, namely, pyramidal cells and interneurons. Their model aimed at capturing the transition to pathological EEG signals in the case of certain cases of epileptic seizure, in particular in temporal lobe epilepsy. Later on, Wendling and Chauvel (Wendling et al, 2002) extended the Jansen-Rit model to account for one more inhibitory population, having both \fast (somatic) inhibition" (FSI) and \slow (dendritic) inhibition" (SDI). This extension is based on experimental evidences and supposedly captures more features of the transition to seizure than the Jansen-Rit model, namely, small-amplitude fast oscillations at onset. In the present paper, we reduce the Wendling-Chauvel model to its minimal conguration, that is with eight state variables, and allow a key parameter, related to SDI, to vary slowly. In this way, we can reproduce time series from (Wendling et al, 2002; Wendling and Chauvel, 2008) corresponding to the dynamical transition to seizure via pre-ictal spikes and low-voltage fast oscillations. By doing so, we propose a mechanism that underpins these complex oscillations, based on the multiple timescales present in the model and related to the canard phenome-
non (Krupa and Szmolyan, 2001), well known to capture rapid changes between vastly dierent regimes of neural activity (Ermentrout and Wechselberger, 2009; Kramer et al, 2008; Krupa et al, 2008). This approach relies on analysing the bifurcation structure of the so-called fast subsystem, obtained when the slow variable's dynamics is frozen. That variable hence becomes a parameter that can then be varied statically to uncover an underlying bifurcation structure. In our case, the fast subsystem is the Wendling-Chauvel model and the parameter in question is the slow inhibitory synaptic gain. In the full model, one observes a dynamic transition through the fast subsystem's bifurcation structure, and such dynamic bifurcations are accompanied by "delayed" transitions that correspond to canard segments; see Section 4. Before presenting our extended model and analysing its slow-fast dynamics, we will rst review the main aspects of bifurcation analysis applied to epilepsy models.
Bifurcation Analysis and beyond
Bifurcation theory is an important mathematical tool to analyse the long-term behaviour of dynamical systems depending on one or several parameters. It provides a skeleton of models in terms of parameter dependence and, in particular, characterises qualitative changes in the form and associated dynamical regime of solutions that a system admits, such changes being referred to as bifurcations. For instance, one will often look for transitions from rest to periodic states (Hopf bifurcations) and for the possible appearance of additional frequencies in oscillatory solutions (period-doubling bifurcations, torus bifurcations).
An important feature of neuronal dynamics is that of being excitable and associated bifurcation scenarios generating excitability are important to characterise using bifurcation theory (e.g. subcritical Hopf bifurcation) ; see (Breakspear et al, 2006) for an example in the context of epilepsy. Changes of solution type, upon parameter variation, can also be observed when systems are operating close to instabilities (Milton et al, 2004; Zetterberg et al, 1978), that is, in a parameter regime close to, but not at, a bifurcation point. The idea that a disease such as epilepsy can occur because the system operates close to an instability has been put forth in, e.g., the two above references. Related to this concept is that of critical transitions, where a system is close to a bifurcation point but jumps to another activity regime without reaching the immediate vicinity of that particular bifurcation point, which typically can be reproduced in systems with noise. This scenario has also been invoked as a trigger for seizure appearance (Kramer et al, 2012).
In trying to characterize epilepsy using dynamical systems theory, the idea of dynamical disease has been introduced by Mackey and Milton in (Mackey and Milton, 1987). It was then applied to epilepsy in various works (Froohlich et al, 2006; Lopes Da Silva et al, 2003). What this concept means is simply that one should consider certain parameters in these models as slowly-varying and that the dynamic transitions observed in the resulting model are related to bifurcations undergone by the original model when varying these parameters. This naturally leads to the notion of dynamic bifurcation (Benoit, 1991), closely related to that of canard solutions already evoked in the previous section. Most importantly, these concepts allow to explain dynamic transitions between dierent regimes of activity in a given model, making the hypothesis that a slowly-varying quantity pushes dynamically the system from a resting state to a seizure state.
The Jansen-Rit and the Wendling-Chauvel Models
The Jansen-Rit (JR) Model
The model (Jansen et al, 1993; Jansen and Rit, 1995) is composed by two subsets of cells: pyramidal (PY) cells and interneurons. The main variables are PSPs, there
are three of them, which eventually gives a six-dimensional model. Each subset is characterised by two functions: the "pulse-to-wave" (P2W) function - which is linear and represented by a second-order linear ordinary dierential equation (ODE)|, and the \wave-to-pulse" (W2P) function, which is a static nonlinear function chosen to be a sigmoid. The P2W is organised as the output of a second-order linear dierential equation; it then represents a linear transfer function that changes presynaptic information into postsynaptic information, and it acts as a second-order lowpass lter. The W2P allows to convert the average level of membrane potential of neurons in a subset to an average pulse density of potentials red by these neurons. Therefore, the model's equations have the following form:
y - y3
y3=AaS(y1-y2)-2ay3-a2yo y. i - y4
y 4=Aa(p(t) + C2S(C1y0)) - 2ay4 - a2y1
y 2 - y5
y 5-Bb(CS (Cy)) - 2by5 - b2y2
(1)
where yff y1 and y2 are the PSPs, p(t) is the external input, the C's are constant representing average synaptic contacts between populations, and S is a sigmoid function. The main output of the model is y, it represents a summated average of PSPs on PY cells and reflects an EEG signal. The JR model can be viewed as a non-linear feedback system driven by a noise input p(t) that globally represents the average density of aerent action potentials from neighbouring or distant populations.
Wendling and Chauvel's Extension of the JR Model
Wendling, Chauvel and co-authors consider two populations of interneurons, namely slow (dendritic-projecting) interneurons (whose main associated parameter is B in the model) and (somatic-projecting) fast interneurons (whose main associated parameter is G in the model). Their model (Touboul et al, 2011; Wendling et al, 2002; Wendling and Chauvel, 2008) has the following form:
Figure 1. Minimal Wendling-Chauvel model. Block-diagram (left) and schematic diagram (right) of the minimal Wendling-Chauvel neural-mass model (3).
10 20 30 40 B 50 6 7 s v B
Figure 2. Fast oscillations in the minimal Wendling-Chauvel model. Torus canard (a) and spiking (b) solutions of system (3).
(2)
y - y5
- AaS(yfy2-y3)-2ay5-a2yo
Vi i - Ve
y e- Aa(p(t) + C2S(C1v0)) - 2ay6 - a2y1
V 2 - V7
y 7-BbC4S (Cyy) - 2by7 - b2y2
V 3 - V8
y.8- GgC7S (C5Vo- CeV4) - 2sy8 - g2Vy
V 4 - y9
y 9-BbS (Cyy) - 2by9 - b2V4
The main reason for the extension of system (1) by Wendling and Chauvel is the presence of two types of Gaba-ergic inhibition: fast, somatic inhibition (Gabaa; fast) and slow dendritic-projecting interneurons (Gabaa; fast), and a dierence in their behaviour prior to the beginning of the seizure. In the context of temporal lobe epilepsy (TLE) experiments, it was shown (Cossart et al, 2001) that dendritic inhibition decreases in the pre-ictal period. At the same time the fast inhibition increases, possibly in part due to the reduction in the inhibitory projection to the fast spiking interneurons.
Extension of a "Minimal" Wendling-Chauvel Model
Reduction of the Wendling-Chauvel Model
One can reduce the Wendling-Chauvel (2) by getting rid of the two equations for y4 and y9. Indeed, a simple rescal-ing shows that these two variables are proportional to variables y2 and y7, respectively. Therefore, the minimal three-population model extension of the JR model has the form
y - y4
y 4 - AaS(y1-y 2-y 3)-2ay 4-a2y„ y. i - ys
y 6- Aa(p(t) + C2S(C1y0)) - 2ay6 - a2y1
y 2 - y6
y 2 -BbCß (Cy) - 2by6 - b2y2
y 3 - y7
y7- gSC7S (Cyo- CJC4Ü2) - 2by7- b2y4
(3)
where the ratio C6/C4 corresponds to the aforementioned scaling factor. We show in Fig. 1 diagrams representing this minimal model. The model has the right dimensionality, that is, eight state variables, corresponding to four PSPs each modelled by a second-order system.
Extension of System (3)
The main dynamical capabilities of system (3) can be understood using the concept of dynamic bifurcation; see Section 4 below for details. That is, the behaviour of the model (upon a dynamic slow variation of parameter B) is characterised by a slow passage through the bifurcation diagram for the so-called fast subsystem, obtained when considering B as a parameter, which is obviously the case in (3). Hence, in (Wendling and Chauvel, 2008), Wendling and Chauvel show a brute-force bifurcation diagram (only stable branches are computed, by direct simulation with a suciently long time of integration so that any transients have disappeared) of their model, upon variation of parameter B; see figure 10 of (Wendling and Chauvel, 2008). Using numerical continuation, one can complement such a brute-force diagram and show the unstable branches, which are important since they correspond to boundaries between dierent activity regimes. We show the result in Fig. 2. Furthermore, one can make B a slowly-varying quantity and initiate a modeling process for its time evolution by creating a slow dierential equation for B. This is the purpose of our extension of system (3), which has the form
V »- V4
y4 - AaS(Vi~y 2-y 3)-2ay i-a2ytl V, i - Vs
y e- Aa(p(t) + C2S(C1yo)) - 2ay5 - a2y1 V2-Ve (4)
y 2 -BbCJS (C3yo) - 2by6 - b2y2
V 3 - V 7
y 7 - GgC7S (CsVo - Ce/C4V2) - 2gy7 - g2Vy
B --£(B - Bo + aSB(Vo)),
where SB is a sigmoid function. Therefore, we push the approach of static variation of B exposed in (Wendling and Chauvel, 2008) one step further by considering this variation of B as the result of a slow dynamical process. This slow process involves both a decrease towards a baseline level Bo and a subsequent increase nonlinearly triggered by the behaviour of the pyramidal cells. Note that this is a very rst attempt to model the variation of slow dendritic inhibition during a seizure event, which is not yet entirely grounded in experimental justications but rather a stepping stone towards a systematic modeling of it.
0.3
6 7 8 9 ß 10
Figure 3. Torus canard orbit in the minimal Wendling-Chauvel model with slowly-varing dendritic inhibition. The orbit is projected onto the ((Y0;B) space and superimposed onto the bifurcation diagram of the fast subsystem (using slow-fast dissection). This plot is a zoom of Fig. 2 (a) in the region of low-voltage fast oscillations at seizure onset. The point TR indicates the location of the torus bifurcation of the full system, which lies very close to the fold of cycles bifurcation points of the fast subsystem as expected in a torus canard phenomenon (Burke et al, 2012).
Dynamic Bifurcations and Canards
Since the early 1980s, the canard phenomenon has given rise to numerous work, both theoretical (Benoit et al, 1981;Benoit, 1991; Krupa and Szmolyan, 2001) and also in link with applications, in particular in neuroscience (Krupa et al, 2008; Moehlis, 2006). Canards are special solutions of dynamical systems with multiple timescales that follow for long time intervals repelling (locally invariant) slow manifolds. As a result, they are associated with rapid transitions - more precisely, exponentially small parameter variations in the timescale parameter 0 < £ << 1 | between dierent activity regimes, in particular from rest to spiking (Moehlis, 2006) in systems with at least one slow and at least one fast variables. This brutal transitions are typically referred to as explosions and, in the neural context, they can related to excitability properties of the cells under investigation (Desroches et al, 2013; Wechselberger et al, 2013). Three-dimensional slow-fast dynamical systems also feature many interesting canard phenomena in link with neural activity. In systems with two slow variables, the concept of mixed-mode oscillations (MMOs) (Desroches et al, 2012b) has been very successful to model the in-
terplay between subthreshold oscillations and spiking behaviour in various cells as well as the sensitivity of the oscillatory pattern on parameter changes (Krupa et al, 2008). On the other hand of the spectrum, three-dimensional systems with one slow variables naturally display dynamic bifurcations (Benoît, 1991), which correspond to slow passages through bifurcation points of the fast subsystem (where the slow variable is frozen and considered as a parameter). Dynamic bifurcations generate delays in the resulting full dynamics (e.g. dynamic Hopf bifurcation) and since the pioneering work of Rinzel (Rinzel, 1986), this concept has become a key tool to understand bursting oscillations from a timescale viewpoint. Indeed, the method of slow-fast dissection introduced by Rinzel (Rinzel, 1986) in the mid-1980s aims to compare the transitions between quiescence and burst phases of a bursting system, with passage through bifurcation points of the fast subsystem. This idea was further exploited by Izhikevich (Izhikevich, 2000) in the early 2000 to come up with a rather complete classication of bursting patterns.
The link between canards and bursting through spike-adding has been suggested by Terman in (Terman, 1991), then more recently revisited in (Guckenheimer and Kuehn, 2009) and, viewing it as another form of canard-explosive phenomenon, in (Desroches et al, 2012a). A related link between canard dynamics and bursting is through transition between bursting and spiking, which involve canard solutions within the burst, termed torus canards. This phenomenon has been observed in a number of neural models (Burke et al, 2012) and it is currently being analysed theoretically, in particular the transition regimes from MMOs to torus canards (Burke et al, 2016). Torus canards are solutions for which the burst stays close to both attracting and repelling families of limit cycles of the fast system. For a system to display torus canard dynamics, a fold bifurcation of cycles is required in the fast system. Described in an elementary way by Izhikevich in the context of elliptic bursters (2001), then unveiled more clearly by Kramer, Traub and Kopell in 2008 in the context of a Purkinje cell model (Kramer et al, 2008), it was then found systematically in many neural models whose fast subsystem displays a fold bifurcation of cycles (in link with a torus bifurcation in the full system) in (Burke et al, 2012). Beyond the mathematical interest for such complicated canard solutions that also include fast oscillatory dynamics, the concept of torus
Figure 4. Torus canard output time series. The time series for variable Y0 of the torus canards shown in Fig. 3 is shown.
10 15 21) 25 B 30 14 16 IS B 20
Figure 5. Torus canard in the Jansen-Rit model. By increasing the (xed) value of the external input p, one can obtain a bifurcation diagram for the fast subsybtem which has the correct structure to allow for torus canards in the extended systems where B is slowly evolving. Panel (a) shows the solution of the extended system together with the bifurcation diagram of the fast subsystem for p = 0:25; panel (b) shows a similar projection for p = 1.
canards is the natural extension to that of canard cycle to understand the boundaries of the bursting regimes in a large class of models. We now give numerical evidence that torus canards occur in the extended Wendling-Chauvel neural-mass model (4), in link with the transition to seizure states.
Torus Canards in the Extended Model
A simple numerical bifurcation analysis of the minimal WC model (4) reveals a structure that is suggestive of canard dynamics in the extended system where the dendritic inhibition varies slowly. This is a natural step forward from a quasi-static variations of these parameters as performed in (Touboul et al, 2011; Wendling et al, 2002) and which already revealed a bifurcation with both stationary and oscillatory branches of solutions, in particular multiple stable oscillatory branches separated by fold bifurcations of cycles. Therefore, when these parameters are slowly-varying, one can observe by direct simulation of the extended model dynamic transitions between dierent regimes of activity: quasi-stationary, quasi-periodic with small-amplitude oscillations, quasi-periodic with large-amplitude oscillations. This allows to reproduce times series featuring transitions to epileptic seizures, similar to the simulations from the 2002 paper by Wendling and Chauvel where they vary these quantities quasi-statically.
Going even one step beyond a slow drift on the dendritic inhibition, we propose a slow dynamics on this main parameter as the result of a dierential equations, which opens the way for a biophysical modeling of its evolution in a seizure regime. In particular, torus canards orbits organise the transition to solutions with epileptic bursts. In Figs. 2 and 3, we present the result of direct simulations of our extended model (4). The solution is plotted on top of the bifurcation diagram of the fast subsystem, as is customary when performing a slow-fast dissection of the problem, and we can observe a solution that evolves around the fold bifurcation of cycles of the fast subsystem, hence containing a torus canard segment; see Fig. 2 (a) and a zoomed view in Fig. 3. The same solution is presented in the time domain in Fig. 4. In Fig. 2 (b), we show a solution that converges towards a
fast-spiking regime, and which corresponds to a slightly dierent parameter set as in the previous case. This is a numerical evidence that the torus canard regime organises the transition between spiking and epileptic bursts in this extended model.
We also revisit the argument by Wendling and Chauvel according to which their model extends that of Jansen and Rit in that it allows to reproduce small-amplitude high-frequency oscillations at seizure onset. We show in Fig. 5 that this eect is due to the presence of a second inhibitory population. Indeed, the fold bifurcations of cycles responsible for the appearance of such small-amplitude burst (transition orbits being torus canards) disappear if the external input p increases past a certain value. However, this value of p increases substantially as the somatic inhibition parameter gets larger. That is to say, when the somatic inhibition parameter is zero (corresponding to the Jansen-Rit model), a small increase of p destroys the possibility for these small-amplitude bursts to exist, whereas when it is non-zero and becomes larger, p has to increase beyond reasonable values for that phenomenon to happen, even though it can happen. That is, the Jansen-Rit model can sustain torus canard solutions if the dendritic inhibition varies slowly, however the associated parameter regime is unrealistic. On the contrary, in the Wendling-Chauvel model, when the somatic inhibition parameter is suciently large, these small-amplitude bursts are always present at seizure onset and the adequate bifurcation structure of the fast subsystem is guaranteed. This is precisely why, in this modeling framework, the presence of a second inhibitory population is crucial to capture seizure-like states.
Conclusions and Perspectives
In this short article, we have shown that the Wendling-Chauvel neural-mass model can be reduced by two equations without losing any salient feature of its dynamics, simply eliminating two redundant equations. A bifurcation study of this model when considering the slow inhibitory synaptic gain as the main parameter, reveals a structure that allows for multiple oscillatory solutions to exist. Previous studies like (Grimbert and Faugeras, 2006; Touboul et al, 2011; Wendling et al, 2002)
have stopped at this fast subsystem bifurcation structure, then varying the main quasi-statically or driving a noisy external input in order to generate time series that resemble EEG data during seizure events. In this work, we push the slow-fast analysis further by proposing a rst dynamical modeling of the slow dendritic inhibition B within a neural-mass model, making it evolve slowly through the fast subsystem bifurcation diagram, in particular near bifurcation points. In this way, we can produce seizure-like time series and understand dynamically how they appear and disappear depending on the value of secondary parameters (in particular, the constant value of the somatic inhibition). Performing this dynamic bifurcation study, we can also understand why the second inhibitory population is important to obtain better synthetic time series for transitions to seizure.
Furthermore, we highlight the fact that canard theory, which has already been invoked to explain rapid transitions in neural rhythms, could also play an important role in modeling epileptic seizures. This complements the idea of modeling epilepsy as a dynamical disease, which then manifests itself through dynamic bifurcation that can be related to the presence of maximal canard solutions acting as boundaries between healthy and pathological states.
Future work will include making the slow evolution of dendritic (and possibly somatic) inhibition more precise, with the aim to gain more biological plausibility. We also want to compare neural-mass models to network models (e.g. that of Lopes da Silva et al. (Lopes Da Silva et al, 2003)), where slow-fast analysis can also be performed, canard solutions playing an equally important role, but where the variables are easier to control and to relate to biophysical elements.
Acknowledgements
The authors would like to thank Fabrice Wendling for useful discussions.
References
Benoît E (1991) Dynamic bifurcations: proceedings of a conference held in Luminy, France, March 5-10, 1990. 1493,Springer Verlag Benoît E, Callot JL, Diener F, Diener M (1981) Chasse au
canard. Collectanea Mathematica 32(1-2):37-119 Breakspear MJ, Roberts JA, Terry JR, Rodrigues S, Ma-hant N, Robinson PA (2006) A unifying explanation of primary generalized seizures through nonlinear brain modeling and bifurcation analysis. Cerebral Cortex 16(9):1296-1313 Burke J, Desroches M, Barry AM, Kaper TJ, Kramer MA (2012) A showcase of torus canards in neuronal bursters. The Journal of Mathematical Neuroscience 2(1):1 Burke J, Desroches M, Granados A, Kaper TJ, Krupa M, Vo T (2016) From canards of folded singularities to torus canards in a forced van der pol equation. Journal of Nonlinear Science 26(2):405-451 Cossart R, Dinocourt C, Hirsch JC, Merchan-Perez A,
De Felipe J, Ben-Ari Y, Esclapez M, Bernard C (2001) Dendritic but not somatic gabaergic inhibition is decreased in experimental epilepsy. Nature Neuroscience 4(1):52-62 Desroches M, Burke J, Kaper TJ, Kramer MA (2012a) Canards of mixed type in a neural burster. Physical Review E 85(2):021,920 Desroches M, Guckenheimer J, Krauskopf B, Kuehn C, Osinga HM, Wechselberger M (2012b) Mixed-mode oscillations with multiple time scales. SIAM Review 54(2):211-288
Desroches M, Krupa M, Rodrigues S (2013) Inflection, canards and excitability threshold in neuronal models. Journal of Mathematical Biology 67(4):989-1017 Ermentrout GB, Wechselberger M (2009) Canards, clusters, and synchronization in a weakly coupled inter-neuron model. SIAM Journal on Applied Dynamical Systems 8(1):253-278 Freeman WJ (1972) Waves, pulses, and the theory of neural masses. Progress in Theoretical Biology 2(1) Freeman WJ (1975) Mass action in the nervous system.
Academic Press, New York Froohlich F, Bazhenov M, Timofeev I, Steriade M, Se-jnowski TJ (2006) Slow state transitions of sustained neural oscillations by activity-dependent modulation of intrinsic excitability. The Journal of Neuroscience 26(23):6153-6162 Grimbert F, Faugeras O (2006) Bifurcation analysis of jansen's neural mass model. Neural Computation 18(12):3052-3068 Guckenheimer J, Kuehn C (2009) Computing slow manifolds of saddle type. SIAM Journal on Applied Dy-namicalSystems 8(3):854-879 Izhikevich EM (2000) Neural excitability, spiking and bursting. International Journal of Bifurcation and Chaos 10(06):1171-1266 Jansen BH, Rit VG (1995) Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biological Cybernetics 73(4):357-366 Jansen BH, Zouridakis G, Brandt ME (1993) A neuro-physiologically-based mathematical model of ash visual evoked potentials. Biological Cybernetics 68(3):275-283
Kramer MA, Traub RD, Kopell NJ (2008) New dynamics in cerebellar purkinje cells: torus canards. Physical Review Letters 101(6):068,103 Kramer MA, Truccolo W, Eden UT, Lepage KQ, Hochberg LR, Eskandar EN, Madsen JR, Lee JW, Mahesh-wari A, Halgren E, Chu CJ, Cash SS (2012) Human seizures self-terminate across spatial scales via a critical transition.
Proceedings of the National Academy of Sciences of the
USA 109(51):21,116-21,121 Krupa M, Szmolyan P (2001) Relaxation oscillation and canard explosion. Journal of Dierential Equations 174(2):312-368
Krupa M, Popovic N, Kopell NJ, Rotstein HG (2008) Mixed-mode oscillations in a three time-scale model for the dopaminergic neuron. Chaos: An Interdisci-
plinary Journal of Nonlinear Science 18(1):015,106 Lopes Da Silva F, Hoeks A, Smits H, Zetterberg LH (1974) Model of brain rhythmic activity. Kybernetik 15(1):27-37
Lopes Da Silva F, Blanes W, Kalitein SN, Parra J, Suc-zynski P, Velis DN (2003) Epilepsies as dynamical diseases of brain systems: basic models of the transition between normal and epileptic activity. Epilepsia 44(s12):72-83
Mackey MC, Milton JG (1987) Dynamical diseases. Annals
of the New York Academy of Sciences 504(1):16-32 Milton JG, Foss J, Hunter JD, Cabrera JL (2004) Controlling neurological disease at the edge of instability. In: Quantitative neuroscience, Springer, pp 117-143 Moehlis J (2006) Canards for a reduction of the Hodgkin-Huxley equations. Journal of Mathematical Biology 52(2):141-153
Rinzel J (1986) A formal classication of bursting mechanisms in excitable systems. In: Proceedings of the International Congress of Mathematicians, pp 1578-1593
Terman D (1991) Chaotic spikes arising from a model of bursting in excitable membranes. SIAM Journal on Applied Mathematics 51(5):1418-1450 Touboul J, Wendling F, Chauvel P, Faugeras O (2011) Neural mass activity, bifurcations, and epilepsy. Neural Computation 23(12):3232-3286 Wechselberger M, Mitry J, Rinzel J (2013) Canard theory and excitability. In: Nonautonomous dynamical systems in the life sciences, Springer, pp 89-132 Wendling F, Chauvel P (2008) Transition to ictal activity in temporal lobe epilepsy: insights from macroscopic models. In: Computational Neuroscience in Epilepsy, Academic Press, pp 356-386 Wendling F, Bartolomei F, Bellanger JJ, Chauvel P (2002) Epileptic fast activity can be explained by a model of impaired gabaergic dendritic inhibition. European Journal of Neuroscience 15(9):1499-1508 Zetterberg LH, Kristiansson L, Mossberg K (1978) Performance of a model for a local neuron population. Biological Cybernetics 31(1):15-26
