https://orcid.org/0000-0003-1402-1470
Figures
Abstract
The mechanisms underlying the emergence of seizures are one of the most important unresolved issues in epilepsy research. In this paper, we study how perturbations, exogenous or endogenous, may promote or delay seizure emergence. To this aim, due to the increasingly adopted view of epileptic dynamics in terms of slow-fast systems, we perform a theoretical analysis of the phase response of a generic relaxation oscillator. As relaxation oscillators are effectively bistable systems at the fast time scale, it is intuitive that perturbations of the non-seizing state with a suitable direction and amplitude may cause an immediate transition to seizure. By contrast, and perhaps less intuitively, smaller amplitude perturbations have been found to delay the spontaneous seizure initiation. By studying the isochrons of relaxation oscillators, we show that this is a generic phenomenon, with the size of such delay depending on the slow flow component. Therefore, depending on perturbation amplitudes, frequency and timing, a train of perturbations causes an occurrence increase, decrease or complete suppression of seizures. This dependence lends itself to analysis and mechanistic understanding through methods outlined in this paper. We illustrate this methodology by computing the isochrons, phase response curves and the response to perturbations in several epileptic models possessing different slow vector fields. While our theoretical results are applicable to any planar relaxation oscillator, in the motivating context of epilepsy they elucidate mechanisms of triggering and abating seizures, thus suggesting stimulation strategies with effects ranging from mere delaying to full suppression of seizures.
Author summary
Despite its simplicity, the modelling of epileptic dynamics as a slow-fast transition between low and high activity states mediated by some slow feedback variable is a relatively novel albeit fruitful approach. This study is the first, to our knowledge, characterizing the response of such slow-fast models of epileptic brain to perturbations by computing its isochrons. Besides its numerical computation, we theoretically determine which factors shape the geometry of isochrons for planar slow-fast oscillators. As a consequence, we introduce a theoretical approach providing a clear understanding of the response of perturbations of slow-fast oscillators. Within the epilepsy context, this elucidates the origin of the contradictory role of interictal epileptiform discharges in the transition to seizure, manifested by both pro-convulsive and anti-convulsive effect depending on the amplitude, frequency and timing. More generally, this paper provides theoretical framework highlighting the role of the slow flow component on the response to perturbations in relaxation oscillators, pointing to the general phenomena in such slow-fast oscillators ubiquitous in biological systems.
Citation: Pérez-Cervera A, Hlinka J (2021) Perturbations both trigger and delay seizures due to generic properties of slow-fast relaxation oscillators. PLoS Comput Biol 17(3): e1008521. https://doi.org/10.1371/journal.pcbi.1008521
Editor: Peter Neal Taylor, Newcastle University, UNITED KINGDOM
Received: November 19, 2020; Accepted: February 22, 2021; Published: March 29, 2021
Copyright: © 2021 Pérez-Cervera, Hlinka. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the manuscript.
Funding: AP and JH acknowledge the long-term strategic development financing of the Institute of Computer Science (RVO:67985807) of the Czech Academy of Sciences (www.avcr.cz) and Czech Academy of Sciences Praemium Academiae awarded to M. Paluš, JH acknowledges the support of project Nr. LO1611 with a financial support from the MEYS (www.msmt.cz) under the NPU I program. AP acknowledges the support from the Basic Research Program of the National Research University Higher School of Economics. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
The dynamics underlying complex processes usually involve many different time scales due to multiple constituents and their diverse interactions. When modelling such systems, the general distinction of at least two time-scales (fast and slow) is a useful and common conceptualization. Many examples of slow-fast dynamics can be found in cell modelling, ecosystems, climate or chemical reactions [1–4] and more recently in epilepsy [5], of particular interest for this paper.
Epilepsy is a chronic neurological disorder characterized by a marked shift in brain dynamics caused by an excessively active and synchronized neuronal population [6, 7]. Although several dynamical pathways have been proposed to explain the transition to seizure [8–12], in general, epileptic tissue is modelled as a system having two stable states: one corresponding to the low activity state and the other corresponding to high activity (that is to seizure) [13]. Besides external perturbations or noise, transitions between these two stable states can also be modelled considering the existence of a parameter evolving on some slow time scale. Whereas on the fast time scale the system can be seen as a bistable system, the variations of the slow parameter lead to bifurcations providing transitions between states [14].
During the last decade, there has been an increasing number of models approaching epilepsy through slow-fast time scales [15–21]. Recently, the slow-fast dynamics has been proposed to explain the role of the interictal epileptiform discharges (IEDs) in the generation of seizures [22]. The IEDs can be thought of as endogenous inputs affecting the target tissue. However, the effect of IEDs on the tissue activity is quite controversial: where some studies show that IEDs can prevent seizures [23, 24], other studies claim their seizure facilitating role [25, 26]. In the above mentioned work [22], the amplitude and frequency dependence of the effect of perturbations in a simple epilepsy model switching between seizure and non-seizure states due to a slow feedback variable, provided a conceptual reconciliation of the pro-convulsive and anti-convulsive effect of IEDs.
In this paper we elucidate this phenomenon in detail and provide theoretical foundations of this apparent perturbation effect paradox by studying the phase response of a generic relaxation oscillator. We perform this theoretical approach by means of the phase reduction [27]. In addition to simplifying the dynamics, the usage of phase reduction techniques allows the computation of its isochrons and phase response curves (PRCs), which clarify the dependence of the effect of perturbations of the oscillator on the perturbation timing, and also allows the study of possible synchronization regimes [28]. By studying a generic slow-fast system displaying relaxation oscillations we show, analytically, how the slow component of the vector field shapes their isochrons and PRCs, thus ultimately determining its response to perturbations. Therefore, our results, clarify the multifaceted effect of IEDs in epilepsy, and can be straightforwardly applied to understand the temporal dependency of perturbations over any model belonging to the wide family of models relying on slow-fast dynamics.
The paper is structured as follows. First, we present a general introduction to relaxation oscillators introducing the basic notation which will be used throughout the paper. Then, we describe the phenomenological epilepsy model and show how, through its phase analysis, we can unveil the mechanism integrating the contradictory role of IEDs in epilepsy. Next, we show, via a complete theoretical analysis, which factors determine the geometry of isochrons of planar relaxation oscillators and study the response of perturbations of relaxation oscillators. We support our theoretical findings studying the response of perturbations for a different reduced epileptor model and discuss our results in the context of epilepsy. We conclude the paper by explaining the computational techniques in the Materials and Methods section.
Results
Basic introduction to relaxation oscillations
The main aim of this Section is to introduce the reader to the basics of slow-fast systems and in particular to relaxation oscillations. For further details we refer the reader to [29–32]. We will consider systems in this form
(1)
the flow of which will be denoted as ϕt(x, y). Notice that ˙ indicates the derivative with respect to the time, t. As 0 ≤ ϵ ≪ 1, the variables x and y evolve on different time-scales, namely the fast time, t, and the slow time τ = ϵt. Next, we use this distinction between time-scales to illustrate how a system in the form (1) with the extra assumption of f(x, y) = 0 being a cubic manifold, generates a limit cycle (denoted as Γϵ) showing relaxation oscillations [33] (see also Fig 1). Consider a point p = (x, y). First, since ϵ ≪ 1, we can take the limit ϵ → 0 and approximate the dynamics of system (1) by the layer dynamics
(2)
The trajectory of p will initially (approximately) follow the layer dynamics in (2) so it will quickly converge to its set of equilibrium points, defined as the slow manifold
(3)
which in the limit ϵ → 0 corresponds to the nullcline (
) of the fast variable. As we considered the slow manifold
in (3) to be cubic, that is S-shaped, it will have two fold points (given by ∂x f(x, y) = 0), which we denote as
respectively, separating the repelling and attracting branches, denoted as
and
, respectively
(4)
Note that the attracting part of the slow manifold
in fact consists of a top and bottom branch
. Once the system has approached the slow manifold, its dynamics are given by the slow variable
(5)
where ′ denotes the derivative with respect to the slow time τ = ϵt. Furthermore, for points in
satisfying ∂x f(x, y) ≠ 0 we know from the implicit function theorem that we can write a function x = m(y) from f(x, y) = 0, so we can express (5) as
(6)
Therefore, once the trajectory has converged to the slow manifold,
, the y variable evolves following the dynamics in (6), while the x variable is given by x = m(y). So trajectories slowly move along
until reaching the fold points
. There, they become governed by the fast dynamics, leading to an almost instantaneous transition to the other stable branch. Indeed, as Fig 1 shows, this is the mechanism underlying the generation of a stable periodic orbit Γϵ showing relaxation oscillations, that is, the motion over Γϵ consists of the alternation of long intervals of quasi-static behaviour (corresponding to the stable branches
of
) and almost instantaneous transitions between the branches [31].
The slow manifold, , is a S-shaped curve having two stable branches
(solid red line) and one repelling
(dashed red line) (see Eq (4)). Stable and unstable branches of
are separated by the fold points
. A given point, p, (see A) will quickly converge to the attracting branch of the slow manifold
(see B). Then, it evolves along
following (6) until reaching the fold point
(see C) where it traverses fast to the other branch
(see D). Then, following again the slow dynamics, the trajectory approaches
(see E) where it goes back to
(see F). Therefore, the system (1) in the singular limit (ϵ → 0) admits a singular periodic orbit Γϵ (in blue) generating relaxation oscillations.
Phenomenological epilepsy model
As we discussed in the introduction, the mechanism of relaxation oscillations (see Fig 1) has been recently used in [22] to explain the apparent contradictory role of IEDs in epilepsy. In this work, the authors propose the following simple phenomenological epilepsy model, further referred to simply as the phenomenor:
(7)
where v and a represent the firing rate and the excitability of a neuronal population, respectively. The dynamical changes in the excitability depend on the difference between v and h through the function f(x) = (tanh(cx) − a0), that is, an hyperbolic tangent whose slope is given by c. When v values are below h, the excitability increases, whereas when v values exceed h, the excitability decreases. Hence, h can be thought of as a threshold. For this study, h = hm a − hn. We keep fixed the particular set of parameters
(8)
for which the system (7) displays a limit cycle denoted as Γpheno with a period of T ≈ 508.42; although the qualitative behaviour of the model stays the same for a wide range of parameters. Indeed, as τa ≪ τx and the fast nullcline
—which corresponds with the slow manifold
in (3)—describes a cubic curve, dynamics over Γpheno consists of a periodic switching between the states of low and high activity within relaxation oscillations.
The following Fig 2 illustrates the mechanism proposed in [22] by which the phenomenor (7) reconciles the antagonistic role of IEDs. Consider the IEDs as a random train of pulses whose inter pulse interval distribution, ts, follows a normal distribution with mean value, Ts, and standard deviation, σ: . Whether or not a given perturbation causes an immediate transition to seizure depends on whether the perturbation manages to make the trajectory cross from the lower branch above the middle branch of the v-nullcline. If this happens, the trajectory rapidly converges to the upper branch, i.e. transitioning to the seizure regime. However, the response of the system dramatically changes depending on the amplitude, A, and mean inter pulse interval, Ts, of IEDs (see Fig 2A and 2B).
Panels A and B show, in red, the v-nullcline whose stable branches correspond to the stable low and high activity states of the system. The unstable part of the v-nullcline (dashed red line) separates the basin of attraction of both branches. As was illustrated in [22], Figure 5, whether the pulses make the system cross the unstable part of the v-nullcline determines the opposite nature of IEDs. For a random train with amplitude A = 0.25 and the system goes to seizure (panel A). By contrast, for a random train with amplitude A = 0.5 and
the system avoids the seizure state (panel B). By plotting the change in seizure rate Δ as a function of both the amplitude, A, and the mean inter-perturbation interval, Ts (panel C), we can distinguish between pro-convulsive regimes (yellow and white areas) in which the transition is potentiated, and preventive regimes (red and black areas) in which the transition is delayed or completely suppressed. We refer the reader to Materials and Methods section for the specific details about the computation of panel C.
The effect of a single pulse applied to the system, while on the lower branch, is either to keep the trajectory on the lower branch or to cause a transition to the upper branch. Therefore, the total effect of a train of pulses depends on the proportion of pulses causing transitions. Indeed, this dependence can be seen by plotting the change in the seizure rate Δ as a function of both the amplitude, A, and the mean inter-perturbation interval, Ts (Panel C).
Phase dynamics
Oscillations correspond to attracting limit cycles whose dynamics can be described by a single variable: the phase. As we now expose, the study of the dynamics on a limit cycle by means of the phase variable provides a more intuitive and simplified view of its synchronization properties. Consider an autonomous system of ODEs
(9)
whose flow is denoted by ϕt(z). Assume that Z is an analytic vector field and that system (9) has a T-periodic hyperbolic attracting limit cycle, Γ. This T-periodic limit cycle, Γ, can be parametrized by the phase variable θ = t/T as
(10)
so that it has period 1, that is, γ(θ) = γ(θ + 1). While originally defined only on the limit cycle, the phase can be extended to the whole basin of attraction of Γ (which we will denote by
). Indeed, as we consider attracting limit cycles, any point in
converges towards Γ as time tends to infinity. Therefore, we will say that two points
have the same asymptotic phase if
(11)
We further define the isochron
as the set of points having the same asymptotic phase θ [34], that is,
(12)
Let us now consider the effect of an instantaneous delta-like pulse over the T-periodic limit cycle Γ,
(13)
It is clear that the perturbation will just change the trajectory from one point z to another point
. As we illustrate in Fig 3, since the isochrons foliate the whole basin of attraction
of Γ, we can say that the perturbation moved the trajectory from one isochron
to another isochron
, thus causing a phase shift
. However, the phase shift will depend on the amplitude of the pulse and on the phase at which it was applied. This dependency is captured by the Phase Response Curves (PRCs). They are calculated by applying the same pulse to the limit cycle at different phases and registering how much the phase is advanced (or delayed). Let z = γ(θ) be a point on the limit cycle Γ. If we consider an instantaneous pulse as (13), it is clear that it will move z to
+ Δz. Thus, the PRC is defined as
(14)
As Fig 3A shows, the isochrons of Γpheno portray the distribution of phases along the basin of attraction
. Whereas the isochrons for the upper branch of the cycle are almost vertical, the isochrons for the lower branch of the cycle show a more interesting geometry: they start vertical until crossing the a-nullcline, when they all bend. The shape of the PRCs as the amplitude, A, of the pulse increases is determined by this particular geometry of the isochrons. Since there is an almost constant distance of 0.1 between the lower branch of Γpheno and the slow nullcline
, we can distinguish between two cases. For perturbations of A < 0.1, the perturbed trajectories only reach the part of isochrons consisting in almost vertical lines. Therefore, the corresponding phase shift Δθ for perturbations on the upper and lower branches of Γpheno is almost negligible. Hence, the PRC for these phases will be close to zero. Indeed, only in the vicinity of the fold points
the PRCs will show larger values (see zoom window in Fig 3C).
Panel (A) shows the limit cycle Γpheno, 16 equispaced isochrons, the v and a nullclines (dashed black and green curves, respectively) and the fixed point, P, at their intersection. The distribution of isochrons clarifies the time dependency of perturbations: as panel (B) shows, a pulse of amplitude A applied at a time t1 (t2) causes a negative (positive) phase shift, delaying (promoting) the transition to seizure. This time dependency can be directly inferred from panel (A): a pulse of amplitude A applied at a point on the cycle z1 = γ(θ1) = γ(t1/T) (z2 = γ(θ2) = γ(t2/T)) displaces the trajectory to a point (
). Since
(
) the perturbation causes a phase shift
(
) delaying (advancing) the phase of the oscillator. The panel (C) shows the PRCs for the phenomenor for positive voltage pulses of different amplitudes summarising the timing (phasic) effect of a given perturbation.
By contrast, for perturbations of A > 0.1, the change on the geometry of isochrons for points on the lower branch remarkably changes the shape of PRCs. Perturbations on the lower branch will have a delaying effect unless they bring trajectories above the middle branch of the v-nullcline—which corresponds with in (4)—so they advance phase (see points z1 and z2 at Fig 3A and 3B illustrating delaying and advancing effects, respectively). The delaying or advancing effect of a given pulse of amplitude, A, is delimited across a discontinuity for its corresponding PRC at the exact phase θ* for which
(see Fig 3C).
The isochrons and PRCs computed for the phenomenor provide insight about how the combination of both the amplitude, A, and the mean inter pulse interval, Ts, generate the different seizure propensity regimes in Fig 2C. As isochrons in Fig 3 show, unless they bring trajectories above , the effect of positive voltage pulses at a point, z = γ(θ), on the lower branch is to cause a delay Δθ < 0. However, for large enough mean inter-pulse intervals Ts, although perturbations delay the system, they are not frequent enough to stop it from eventually transitioning to seizure (see Fig 2A). Moreover, larger pulses are able to cause the trajectory to cross the v-nullcline earlier through the cycle (way before the fold point). Thus, the larger the amplitude of the pulse, the more common are these transitions. By contrast, for small enough inter-pulse intervals, Ts, the transition to seizure can be delayed or even stopped across the accumulation of the delays caused by each single pulse (see Fig 2B). Thus we can conclude that the mechanism underlying the description of the phenomenor of the role of IEDs, relies on the one hand on its cubic v-nullcline structure, allowing for relaxation oscillations and on the other hand on the prevalence of delays for positive perturbations at points on the lower branch not crossing the middle branch of the v-nullcline.
Phase analysis of relaxation oscillators
As explained in the previous Section, the accurate description of the role of IEDs provided by the phenomenor is based on the prevalence of delays for perturbations in the ‘non-epileptic’ state, i.e. on the bottom branch of the cycle Γpheno. Since this determining feature of the model—the prevalence of delays—is based on the bending in a particular direction of the isochrons, we aim to identify which elements in the model are key to cause this particular isochron geometry. As we show next, we perform this identification by taking advantage of the dynamical properties underlying any relaxation oscillator.
The 
geometry of isochrons 
.
Next, we discuss some generalities shaping the isochrons of planar relaxation oscillators. To begin, it is worth recalling that if two points and z = γ(θ) belong to the same isochron,
, they have to meet at the same point of the cycle after a large enough time, t (see Eq 12). For this reason, the determination of the shape of isochrons requires to study the converging dynamics towards Γϵ which we recall consist of trajectories covering
distances in the fast direction and
distances in the slow direction.
Since we aim to study the isochrons for relaxation oscillators, we can take advantage of the time-scale separation to be more precise concerning this convergence. Consider a point . In a first approximation one can assume that the convergence of
is achieved simply following the layer dynamics (2). If that was the case, since the layer dynamics consider the variable y as frozen, the isochrons will always be lines of y constant that we denote as
. However, for correctly determining the shape of isochrons, we have to take into account that neither the convergence towards the limit cycle Γϵ is instantaneous nor the dynamics on y during convergence are negligible. As a result, the isochrons are expected to be
corrections of
. Indeed, it is worth to note that generalities determining the sign of those
corrections will explain the prevalence of delaying (or advancing) effects of delta-like pulses in the fast direction.
Regarding the time needed for solutions to converge to the limit cycle, although the convergence towards a normally hyperbolic attracting limit cycle is ensured [35], for the case of slow-fast dynamics we can give even more details about this convergence by means of Tikhonov’s theorem [36] (see also [37, 38]). Roughly speaking, Tikhonov’s theorem states that after a time , all orbits starting in a neighbourhood
of the slow manifold
will have reached a neighbourhood of
of
.
Once we know the time, , needed to converge, we can compute the motion of the converging point
in the slow direction. The travelled distance in the y direction by
to approach a
neighbourhood of Γϵ is given by
(15)
where
refers to the solution of the layer system (2). During the time, th, needed to converge, the point z = γ(θ) on the limit cycle has travelled a distance Δy given by
(16)
where in the second equality we utilize the fact that for points on the slow manifold we can use Eq (6).
Now, let us illustrate how the difference between the distances travelled by the base point and the converging point, , will determine the sign of the
correction for isochron
at the point
. If we write the expression for
:
(17)
we can see that the difference
is directly determined by the difference between the speeds in the slow direction for the base z and converging
points during the time
needed for approaching Γϵ. Basically, since both points z and
have to meet at the same point after the same time, the one travelling slower, needs to travel less distance. The difference,
, corresponds exactly to the
correction to
(see Fig 4).
In the limit ϵ → 0 isochrons are lines of y constant denoted by . However, since ϵ ≠ 0 but small, the isochrons are
perturbations of
. As we show in the right panel, the sign of the
corrections depends on the difference of speeds between the converging point
and the base point z during the convergence time th. In this case, to approach Γϵ,
has to cross layers of x whose values are smaller than the ones surrounding Γϵ. For this reason
travels slower than z. Since
and z have to meet after a time th at the same point on Γϵ, but
travels slower than z, then
needs to travel a short distance. This determines the sign of the
correction. Furthermore, if the slow vector field is monotonous along the fast direction, the farther the point
, the slower (faster) it travels, so the slope of the isochrons will have the same sign for all the points
satisfying fast convergence, thus determining the effect of perturbations in the fast direction.
However, at the moment we have a local argument just justifying the shape of isochrons for a given point . Nevertheless, we can globalize this argument by assuming some conditions for g(x, y). In particular, as we show now, if the slow vector field g(x, y) is monotonous in the fast direction, then the slope of the isochrons will have the same sign for all the points
satisfying fast convergence.
The asymptotic phase defined in (11) allows to assign a phase to any point , by defining the following function Θ(z)
(18)
whose level curves indeed correspond to the isochrons. Let us assume we can invert Θ(x, y), so we can define the following function
(19)
The slope of isochron
, which we denote by
, is then given by
(20)
as we also have
(21)
we can write the following expression for the slope
(22)
As the term
can be written in integral form (see Eq (17)), the slope
, can be evaluated as the derivative of the difference of two sums (integrals)
(23)
As we see, assuming that the vector field g(x, y) is strictly increasing (decreasing) function with x it is easy to discuss the sign of
. If the trajectory followed by the approaching point, satisfies
for 0 < δt ≤ th, then, the second integral will be smaller than the first one. Since this difference will increase with
, then
. Furthermore, the larger the changes in g(x, y), the larger the slope. We remark that in the case
, we can argue identically to obtain
.
In conclusion, we have illustrated the relationship between geometry of isochrons for relaxation oscillations and the slow vector field. First, we have shown how the tilt of the isochron at a given point
depends on the difference of speeds between
and the base point z during convergence. Furthermore, we showed that if the monotonicity of the vector field does not change, the tilt of the isochrons does not change sign as well.
We can illustrate these theoretical results by revisiting the isochrons for the phenomenor. As Fig 5 shows, the parameters in (8) were chosen so that the tanh in (7) acts almost as a step function. As a result, the speed in the slow direction dramatically changes when crossing the slow nullcline. Since there are almost no differences between speeds for points below the slow nullcline, the isochrons are almost vertical. By contrast, this large difference of speeds once the slow nullcline is crossed, results in a remarkable bending of the isochrons for points on the lower branch of Γpheno.
For the phenomenological epilepsy model (7) with the set of parameters in (8) the figure shows: (A) Limit cycle Γpheno and its isochrons
(left). (B) Values of the slow vector field (corresponding to
in (7)) for points
.
PRCs.
Since the shape of PRCs is determined by the geometry of isochrons, next we discuss the extensions of our previous analysis of isochrons to PRCs. First, we can consider the limit ϵ → 0. In this case the isochrons would be vertical lines. Therefore, for points in the lower branch, unless the pulse brings trajectories above the mid-branch of the slow nullcline, its corresponding phase shift will be zero. For those points going to the other branch, the phase shift will be proportional to the skipped segment of the cycle, thus generating the characteristic shape of PRCs for relaxation oscillators [39] (see black curve in Fig 6 right). However, our knowledge of the geometry of isochrons can extend this result. Without loss of generality we discuss the case
. In this case, we know that perturbations acting over points on the lower branch not crossing
will delay the system (see Fig 4). As a consequence, the PRC will have negative values for all the phases θ in the lower branch such that θ < θ* where
. Although the particular shape of the delaying segment of the PRC will depend on the particular slow vector field chosen, in general, we expect the crossing of the slow and fast nullclines to generate a single unstable fixed point (denoted by P) inside Γϵ. It is worth to mention that since isochrons will approach P through
[40], we expect the bending of a particular isochron to increase as it approaches
. As a consequence, we expect the maximal delay values of a PRC to concentrate near the critical phase θ*. Finally, if we consider perturbations over points in the upper branch, arguing similarly as in Fig 4, we can conclude that the effect of pulses of positive amplitude is to advance trajectories.
Next we sketch the PRCs for pulses of amplitude A > 0 for the case . For phases θ < θ*, where θ* satisfies
, due to the slope of isochrons the effect of the pulses will be to delay trajectories. Since isochrons approach the unstable point P through
, the closer the phase θ to θ*, the larger the bending of the isochrons and thus the larger the corresponding delay value. For phases
, there is an advancement proportional to the fraction of cycle skipped. This prevalence of advancements is also seen for points in the upper branch. For phases in a neighbourhood of the fold point
, we expect a transition between advancement and delays not drawn because our analysis is only valid for normally hyperbolic points.
Phase locking.
So far we have theoretically identified the factors shaping the isochrons for relaxation oscillators. Furthermore, we have discussed how the particular geometry of the isochrons for relaxation oscillators is reflected in the corresponding PRCs. Next, we aim to continue extending our theoretical approach to determine generalities underlying the mechanism by which external perturbations suppress the original oscillatory dynamics. We recall that, in the particular case of epilepsy, we are studying the suppression of the original oscillation through the perturbation-triggered delays which cause the system to remain in the lower activity state and thus to prevent the transition to seizure; one mechanism for this is the existence of stable phase-locked solutions of the perturbed system.
A delta-like pulse of amplitude, A, reaching the cycle at a phase, θ, will map it to a new phase fA(θ) = θnew, where the map fA(θ) writes as
(24)
If the perturbation was a train of periodic pulses with an inter stimulus interval given by Ts, we can describe the phase dynamics of the system by
(25)
where θ0 = θ. The fixed points of the above map (25), which are given by
(26)
correspond to the phase locking states of the system.
In the following, we shall use Eq (26) to find the phase-locked solutions that prevent seizures. Note that not all the phase-locked states predicted by Eq (26) prevent seizures. For example, consider a train of pulses whose inter stimulus interval is Ts ≈ T. Since it is forcing the system with an almost resonant frequency, the system entrains to it. However, that particular phase-locked state will not prevent seizures since the system will essentially perform one full oscillation before the next perturbation occurs.
The particular phase-locked solution that does prevent seizures is sketched in Fig 7. Consider a pulse displacing a point z = γ(θ) to . If we denote by th the time that
needs to approach Γϵ, we need
with
. That is, we need the perturbed trajectory to reach the cycle at a previous phase. Assuming fast convergence, we can write
(27)
Since we need
as a necessary condition for phase locking, then, if we assume without loss of generality that the motion over Γϵ is counter-clockwise, the above integral has to be negative. For that to happen, the perturbation has to necessarily send trajectories above the slow nullcline. Indeed, if we denote by t* the time needed to cross the slow nullcline, then, the particular class of locking we are interested in has to satisfy
(28)
Since the first integral is negative and the second is positive, above Eq (28) shows that the appearance of phase locking requires the perturbed trajectories to be sent to a point such that the distance travelled during convergence in the negative direction overcomes the distance travelled in the positive direction, so the total displacement is negative. Then there is a time Ts = −TΔθ > 0 (with
) for which the next pulse will kick the system at the same initial point z = γ(θ) (see Fig 7). The repetition of this process keeps the trajectory on the lower branch, and prevents the seizure emergence by suppressing the original oscillatory dynamics. Importantly, we highlight the strong influence of the slow vector field on the appearance of this locking mechanism. Indeed, the smaller the distance between the slow nullcline and the lower branch, the smaller the amplitude of perturbations needed for locking the system.
To suppress the original oscillation and keep the system in the lower branch of Γϵ the amplitude A of the pulse has to be large enough so besides causing a delay Δθ, it displaces trajectories above enough the slow-nullcline so the distance travelled in the negative direction overcomes the distance travelled in the positive direction, thus causing a negative net displacement. The locking appears by repeating this mechanism after Ts = TΔθ intervals so the new pulse always hits the system at the same initial point.
We can check the validity of this result by revisiting the results for the phenomenor. Fig 8A shows the relative seizure rate increase Δ for a Ts periodic train of pulses. We can see how the locking preventing the transition to seizure starts for values , which is the approximate distance between the lower branch and the slow-nullcline. Furthermore, for a fixed amplitude
, if we consider the maximum delay value (denoted by Δθ*) of the corresponding PRC and compute the inter pulse interval value given by
, it is clear that for inter-pulse intervals
, the system is likely to transition to seizure because the delays are not large enough to stop the system. Therefore, we expect the pair
to delimit the locking regime. By computing the PRCs for all the amplitude values satisfying
, we can calculate the corresponding
values and thus generate a curve in the (A, Ts) space—which indeed corresponds with the bifurcation curve of the map (25)—showing a nice agreement with the boundaries of the locking area (see purple line in Fig 8A).
We plot the change in the seizure rate Δ for a random train of pulses following a Gaussian distribution of mean time Ts and standard deviation σ denoted as . Panels (A) and (B) correspond to the deterministic periodic case
and to the random case
, respectively. For panel (A) we plot a purple solid line corresponding to the bifurcation of the phase map (25). We plot the same curve as a dashed purple curve in panel (B) illustrating the resilience of the deterministic phase-locked states to noise. For the specific details about the numerical computation of this Figure, we refer the reader to Materials and Methods section.
Results for the random case in Fig 8B can be interpreted by means of the periodic case. The random dynamics can be computed as well by using a similar map to (25) but substituting Ts for ts values in the distribution . In this case, the system does not ‘lock’ in the same way the deterministic system does, that is through fixed points in Eq (26). However, one might try to interpret random dynamics by means of the periodic case. As we explained for the deterministic case, the maximum delay value Δθ* of the PRC allows to compute a characteristic value of
, such that perturbations
will lead to seizures. Therefore, the robustness of the deterministic locking states to noise decreases, as the probability of ts being above
increases, which happens with increasing the width σ of the inter pulse distribution.
Epileptor model
Apart from long standing detailed epilepsy models [41], a recent successful phenomenological model of epileptic dynamics is the epileptor model [17]. This model consists of 5 differential equations (4 fast and 1 slow) so it can not only display a wide range of dynamical regimes explaining many different pathways to seizure [42], but importantly it also contains a phenomenological yet explicit deterministic mechanism for spontaneous switching between seizing and non-seizing regime. In order to show the generality of the results derived from our theoretical approach and to demonstrate their consequences in models of epilepsy, we will study the following 2D reduction of the epileptor model [43]:
(29)
where v and z represent the firing rate and the permittivity of a neuronal population, respectively. For this model we will work with the sets of parameters
and
in Table 1.
For the set of parameters , the system will display a limit cycle Γi of period Ti.
Identically as the phenomenor, since the time constant for the z variable is small τz ≪ 1, and describes a cubic curve, the three sets of parameters
and
lead to relaxation oscillators denoted as Γ+, Γ0 and Γ− respectively. The three different sets of parameters
and
were chosen to illustrate the influence of the slow vector field on the response of perturbations of the system. Indeed, we denoted the parameters as
and
because they set the nullcline to have positive, horizontal and negative slope, respectively. Fig 9 shows the isochrons and PRCs for the three sets of parameters
and
. Since the slow vector field of the reduced epileptor is monotonic in the fast variable v, the slope of the isochrons does not change sign for any of the considered cases, and again, it causes a prevalence of delays for perturbations of positive amplitude over points on the lower branch which is captured by the PRCs (see Fig 9). We remark the similarity between the computed PRCs in Fig 9 and the ones sketched in Fig 6.
For the sets of parameters and
in Table 1 we show: Limit cycle Γ+, Γ0 and Γ− and its isochrons
(left). The phase response curves for pulses in the v direction for different values of A (right). For the three cases we plot 16 equispaced isochrons. Consistently with our previous analysis, since the monotonicity of the slow vector field does not change, the slope of isochrons does not change sign. For numerical details about the computation of both the isochrons and the PRCs see Materials and methods section.
Response to perturbations
Next, we show how while the unperturbed behaviour of the cycles Γ+, Γ0 and Γ− remains qualitatively identical, that is, they show relaxation oscillations, their response to the same train of pulses will be completely different. As we will argue, these remarkable differences can be explained by the different sets of parameters and
causing different slow vector fields for each cycle. Identically as in the phenomenor case, we consider a random train of pulses whose inter pulse interval follows a normal distribution of mean Ts and standard deviation σ, denoted by
and compute the change of the seizure rate Δ for a train with
and
.
The simulation results are summarized in Fig 10. Consistently with the theoretical results, we can see a direct correspondence between the mean distance between the lower branch and the slow nullcline and the appearance of areas suppressing the oscillation. For this reason, Γ− locks for smaller amplitude values than for Γ0 and Γ+. Furthermore, although the bending of the isochrons is small and so are the corresponding delays Δθ, because of its large T value (see Table 1), the range of Ts = −TΔθ values for which Γ− shows locking is even larger than for Γ0 and Γ+. We also remark the good agreement between the bifurcation curves of map (25) and the areas suppressing the transition to seizure.
We show the change in the seizure rate Δ for a random train of pulses whose mean inter impulse interval follows a normal distribution with a mean time Ts and standard deviation σ. Panels A, B, C correspond to the sets of parameters
and
in Table 1. Left Figures correspond to the periodic case
and right Figures to the random case
. Consistently with our theoretical analysis there is a direct correspondence between the mean distance between the lower branch and the slow nullcline and the minimal pulse amplitude A for which perturbations may lead to lock the system. Purple solid lines, bounding locking regimes, correspond to the bifurcations of the map (25). By drawing the same curve for the random case, we illustrate the resilience of locking states to noise. For the specific details about the numerical computation of this Figure, we refer the reader to Materials and Methods section.
Regarding the interpretation of the random perturbation train scenario, we can interpret results approximately by means of the results for the periodic perturbation scenario. Similarly as we argued in the phenomenor case (see Fig 8), the robustness of a given locking state to noise will depend on whether the critical value of , (where Δθ* corresponds with the maximal delay value of the PRC) is or not within the width σ of the distribution
. The higher the probability of occurrence of
values, the likely is the system to switch to the upper branch. The differences in the resilience of the deterministic locking areas for Γ+, Γ0 and Γ− in Fig 10, can be explained by the different values of the period for the 3 cycles (see Table 1). Despite the PRCs for the three cycles show a similar range of values for the delays Δθ, the differences come when these delays are transformed in inter impulse intervals through Ts = TΔθ. The shorter the period T, the smaller the critical
value. Since in the three cases the ts distributions have the same width, the smaller the critical
value, the higher the probability of occurrence of
values. As a consequence, the resilience of locking states for Γ+ and Γ0 is weaker than for Γ− in which the distribution Ts = −TΔθ is larger because of its larger period.
Comparison between the phenomenor and the reduced epileptor
Although both the phenomenor in Eq (7) and the reduced epileptor in Eq (29) model seizure dynamics through relaxation oscillations, it is worth to mention the different role of the slow variable in the models. In the phenomenor the variable a describes the excitability of the tissue (the higher excitability, the more likely the spontaneous seizure initiation), whereas in the (both original and reduced) epileptor the z variable (dubbed as permittivity) has the opposite polarity: for its low values, the system switches to seizure as its only stable state. As a consequence, although the dynamical mechanism of the two models generate is virtually identical, the monotonicity of the slow vector field and the rotation direction over the cycle is flipped (see Fig 11). However, in both models, the motion and the tilt of the isochrons are related in such a way that the prevalent effect of positive voltage perturbations over the lower branch of the cycle is to slow-down the oscillations, or in particular to delay the seizures.
Each cycle is depicted in purple, the v-nullcline in black and the slow nullcline in green. Notice that the direction of the slow variable in both models is flipped, and thus is also the motion over the cycles and the sign of the derivative of the slow vector field in the fast direction .
To further compare both models from a general mathematical perspective, we come back to the notation for a generic planar slow-fast system we defined in (1), where g(x, y) corresponds to the slow vector field and x, y to the fast and slow variables, respectively. The main differences between both models rely on their different time constant τ values and the specific slow vector field functions g(x, y). Because of the correspondence between τ ≪ 1 and ϵ, we expect the isochrons to be bounded in domains (see Fig 4). However, from our analysis it also follows that the bending of the isochrons, although being contained in
domains, will be also determined by dependence of g(x, y) on the fast variable x between the perturbed and the base trajectories (see Eq (23)). To illustrate these role of τ and
, let us compare Γpheno with Γ−. In both cases, the slow nullcline was near the lower branch, so we have a (qualitatively) similar geometry for both phase spaces. For this reason the response to perturbations was qualitatively similar in both cases (compare Figs 8 and 10C). However, the larger range of Ts values for which perturbations over Γpheno avoid seizure can be explained by both the larger τ and the strong change in the monotonicity of g(x, y) for the PE. The combination of both effects causes a larger bending of the isochrons and thus larger delays. Therefore, for (qualitatively) similar geometries, the differences in both the time constants values and in the strength of variations in the fast component of the slow-vector field have a substantial effect on the amplitude of the phase response of the system to inputs.
Discussion
In this paper we applied a phase approach to analyse planar relaxation oscillators, motivated by models of epileptic dynamics. Indeed, the study of neural oscillators by means of the phase reduction has been extensively utilized in neuroscience from the level of single neurons to the network scale [28, 44–46]. In this work, the computation of isochrons and PRCs of the phenomenological seizure dynamics model introduced in [22] fully clarified the mechanism integrating the antagonistic potential effects of IEDs. Furthermore, the theoretical analysis of the phase response of a generic planar relaxation oscillator manifested the crucial role of the slow vector field on the geometry of their isochrons. Due to the direct link between isochrons and PRCs, we have been able to study the relationship between the slow vector field and the different response behaviour a planar oscillator can display depending on the amplitude and frequency of perturbations. For the cases considered, whereas the distance between the slow nullcline and the bottom branch of the cycle indicated the minimum value of amplitude values suppressing the original oscillation, the minimum value of PRCs (that is, the maximum delay) was related to the maximum interpulse intervals for which this locking mechanism holds. Furthermore, besides confirming our results, the study of variants of the reduced epileptor model [43] showed how vastly different responses to perturbations can be exhibited by models differing only in the slow-nullcline position, but possessing almost identical unperturbed behaviour, i.e. equivalent limit cycle oscillations, thus demonstrating the key role of the slow vector field in the response of perturbations for planar relaxation oscillators.
We acknowledge that due to the motivation by models of epilepsy, we showcased the theory only on a small set of example dynamical systems previously used for modelling the cyclical transition between an ictal and interictal state, which showed quite similar dynamics, including having one linear and one cubic nullcline, and a monotonous slow component of the flow field. A quick glance at other slow-fast relaxation oscillator models however suggests, that these properties are far from uncommon in many other models. Moreover, careful consideration of the theoretical arguments however shows, that the specific linear or cubic shape is indeed not crucial for the general observations to hold. Also, careful consideration of the theoretical arguments shows, that the monotonicity of slow vector field is firstly quite natural (the function needs to change from positive to negative values between the two stable branches of the stable manifold; it may likely do so just monotonically); and moreover not necessarily needed—if the change is not monotonic, the dependence of the PRC on the size of the perturbation just becomes more complicated, however the (sign of) the PRC is still given by the integral of the slow component along the recovery trajectory.
Another apparent limitation is that we focused on the effect of positive pulses acting on the bottom branch of the cycle. However, the approach straightforwardly extends to planar oscillators having more complex slow vector fields and to pulses of different sign applied either to the lower or higher branch. Indeed, we suggest that for a given slow vector field the applied geometrical approach is instrumental in providing an intuitive insight concerning the isochrons and therefore the PRCs. In that sense, our analysis extends previous results on PRCs and isochrons of planar relaxation oscillators beyond the weak and singular limit [39, 47]. Theoretically more interesting, while also more demanding, is the generalization to higher dimensional oscillators, providing richer geometrical structure of the flow, perturbations and trajectories. However, previous simulation-based results on the full Epileptor model [22] suggest that the potential dual effect of perturbations on oscillatory behaviour is preserved even in higher dimensions, although richer behaviour might show for other models or perturbation scenarios.
Regarding epilepsy, our results indicate the key influence of the slow vector field on the propensity for seizure emergence. We acknowledge our analysis relied on reduced planar models. However, we plan to make advantage of recent methodologies computing isochrons of high dimensional systems [48] to extend our approach to different high dimensional models as [17, 19–21, 49]. In general, the high dimensionality of these models permits to describe more accurately seizures initiation and termination [14, 50]. We believe the continuation of this line of research may provide an alternative vision to the questions these models approach. Furthermore, because of the usage of the phase variable and the determination of PRCs, we think this approach can also help to determine more accurately coupling functions for studies approaching epilepsy from the coupling of different oscillatory units [51].
Importantly, the quest for deeper and intuitive understanding of the effect of perturbation on epileptic network dynamics is not just an intriguing mathematical exercise, but an indispensable part of an important while difficult journey to understand the mechanisms of seizure initiation, and the possible ways to preclude this initiation by therapeutic stimulation interventions [52]. Of course, while the general conceptual insights are on their own relevant for general understanding the possible dynamical phenomena in response to perturbations, the observed role of the slow component of the field and in particular the nullcline suggests that any computational models of epilepsy dynamics should also attempt to reasonably approximate these aspects (and not only the unperturbed behaviour), if aspiring for providing relevant predictions concerning treatment protocols or just outcomes of endogenous perturbations and inter-regional interactions. This opens also the question of how to practically estimate these properties from experimental data, be it through stimulation protocols or purely observation data; this seems to be a natural avenue for obtaining more realistic models of epileptic dynamics.
In conclusion, we have outlined and carried out phase response analysis of planar relaxation oscillator models of epileptic dynamics that opens not only a path in epilepsy research with many interesting analytical, computational, experimental and potentially clinical implications, but also provides a framework applicable to gain insight in the plethora of other computational biology problems in which slow-fast relaxation oscillator models are pertinent.
Materials and methods
This section contains some technical details concerning the numerical implementation of computations used to provide the presented results. Integration of ordinary differential equations was done using a 8th-order Runge-Kutta Fehlberg method (rk78) with a tolerance of 10−14.
Counting of seizures
In this Section we explain how we generate the diagrams in Figs 2C, 8, and 10 showing the effect of perturbations on the transition to seizure for both the phenomenor (PE) and the reduced epileptor (RE). As we explain along the manuscript, both models describe epileptic dynamics through a relaxation oscillator of period T whose dynamics on the upper stable branch correspond to seizures. As for both models, the upper branch of the cycle terminates at the upper fold point (which in both cases corresponds to v = 0), we have proceeded this way: for each case we integrate the corresponding system for a time t ≫ T with a time step Δt and apply a pulse of amplitude A in the v direction at equispaced Ts intervals. Each time the following condition is satisfied: v(t − Δt) < 0 < v(t), we consider the system transitions to seizure. Finally, the change in seizure rate Δ is computed by dividing the number of seizures in the perturbed case by the number of seizures for the unperturbed case (which is 1 seizure per period). For the random case we proceed the same way just perturbing the system at intervals
.
From the adopted criteria for counting seizures it follows that very large perturbations might cause or constitute a seizure per se independently of the Ts value. For that reason, since we were interested in the relationship of both the amplitude A and the inter pulse interval Ts we limited the simulated amplitude in the above mentioned Figures to maximum of A = 0.8 for the PE and A = 1.8 for the RE (note that perturbations of A ≈ 1 for the PE and A ≈ 2 for the RE will cause seizures independently of the Ts value).
Computation of isochrons
To compute isochrons of slow-fast systems, we assume we have an analytic vector field having a T-periodic hyperbolic attracting limit cycle Γ which we parametrise by γ(θ) (see Eqs (9) and (10)). To find γ(θ), we construct a Poincaré section and use a Newton method to find a fixed point of the corresponding Poincaré map. By doing this, we obtain a point z0 ∈ Γ and the period T. Then, we integrate the system (9) with initial condition z(0) = z0 for a time T to obtain
.
Next, we need to compute the linearisation N(θ) of the isochrons around Γ. To that aim, typically one solves a variational-like equation [53]. However, in slow-fast systems the cycle is strongly attracting (indeed, its characteristic multiplier is with k > 0) [54]. For this reason, obtaining N(θ) via numerical integration requires to deal with very small numbers, so one needs high precision algorithms and large number of decimals.
As an alternative to numerical integration, we took advantage of the fact that ∇Θ(z) is perpendicular to the level curves of Θ(z), which indeed correspond to the isochrons. Therefore, we can use the infinitesimal PRC (iPRC), that is ∇Θ(γ(θ)), to compute N(θ) through the following equation [53]:
(30)
where v⊥ refers to a perpendicular vector to v and <⋅, ⋅> to the usual dot product. Instead of computing the iPRC ∇Θ(γ(t/T)) by integrating the adjoint equations (which also display numerical instabilities) we compute it by means of the procedure described in next subsection.
Finally, we globalise the isochrons via the backwards integration of N(θ) (we refer the reader to [53] for more details about the globalisation procedure).
Computation of PRCs
The PRCs in this paper were computed using a continuation method. The computation of PRCs by direct integration of the perturbed trajectories, usually measures the phase shift over the maximum of a certain variable. That is, they require to integrate a relaxation time Trel large enough so the perturbed trajectories reach the maximal values over the cycle. By contrast, as we now show, continuation methods just require the perturbed trajectories to reach a point on the cycle. Therefore, one needs to integrate a shorter time Trel. Specifically in slow-fast systems, in which the periods of the system are large, the usage of continuation methods saves a lot of computational effort. To compute PRCs, we have used the continuation method introduced in [55], which we now briefly review for the sake of completeness.
A pulse acting on a point z = γ(θ) ∈ Γ will displace the trajectory to . Then, after a time Trel large enough, the trajectory will be again on the limit cycle but with another phase
. Mathematically
(31)
where
, and
. Then the PRC is PRC(θ, A) = fA(θ) − (θ + Trel/T).
The idea of the method is to obtain fA(θ) by solving Eq (31). To that aim, one can use the following algorithm which computes the PRC for a perturbation of amplitude A by means of a Newton method. The computation of PRCs via continuation is achieved using the computed PRC as an initial seed for computing the PRC for a new amplitude A′ = A + ΔA. Given the parameterization of the limit cycle γ(θ), and fA(θ) an approximate solution of Eq (31), we perform the following operations:
- Compute E(θ) = FA(γ(θ)) − γ(fA(θ)).
- Compute ∂θ γ(fA(θ)) = TZ(γ(fA(θ))).
- Compute
.
- Set fA(θ) ← fA(θ) + ΔfA(θ).
- Repeat steps 1-4 until the error E is smaller than the established tolerance. Then PRC(A, θ) = fA(θ) − (θ + Trel/T).
We refer the reader to [55] for the implementation of this methodology for not pulsatile perturbations. To compute the iPRC by means of this algorithm one has to consider perturbations of A small and fA(θ) = θ + Trel/T as initial seed. Then, ∇Θ(γ(θ)) = PRC(A, θ)/A.
References
- 1.
Rubin JE, Terman D. Geometric singular perturbation analysis of neuronal dynamics. In: Handbook of dynamical systems. vol. 2. Elsevier; 2002. p. 93–146.
- 2. Rinaldi S, Scheffer M. Geometric analysis of ecological models with slow and fast processes. Ecosystems. 2000;3(6):507–521.
- 3. Ashwin P, Ditlevsen P. The middle Pleistocene transition as a generic bifurcation on a slow manifold. Climate dynamics. 2015;45(9-10):2683–2695.
- 4. Moehlis J. Canards in a Surface Oxidation Reaction. Journal of Nonlinear Science. 2002;12(4).
- 5. Coombes S, Terry JR. The dynamics of neurological disease: integrating computational, experimental and clinical neuroscience. European Journal of Neuroscience. 2012;36(2):2118–2120.
- 6. Fisher RS, Acevedo C, Arzimanoglou A, Bogacz A, Cross JH, Elger CE, et al. ILAE official report: a practical clinical definition of epilepsy. Epilepsia. 2014;55(4):475–482. pmid:24730690
- 7. Jiruska P, De Curtis M, Jefferys JG, Schevon CA, Schiff SJ, Schindler K. Synchronization and desynchronization in epilepsy: controversies and hypotheses. The Journal of physiology. 2013;591(4):787–797.
- 8.
Milton J, Jung P. Epilepsy as a dynamic disease. Springer Science & Business Media; 2013.
- 9. Breakspear M, Roberts J, Terry JR, Rodrigues S, Mahant N, Robinson P. A unifying explanation of primary generalized seizures through nonlinear brain modeling and bifurcation analysis. Cerebral Cortex. 2006;16(9):1296–1313.
- 10. Da Silva FL, Blanes W, Kalitzin SN, Parra J, Suffczynski P, Velis DN. Epilepsies as dynamical diseases of brain systems: basic models of the transition between normal and epileptic activity. Epilepsia. 2003;44:72–83.
- 11. Baier G, Goodfellow M, Taylor PN, Wang Y, Garry DJ. The importance of modeling epileptic seizure dynamics as spatio-temporal patterns. Frontiers in physiology. 2012;3:281.
- 12. Goodfellow M, Schindler K, Baier G. Intermittent spike–wave dynamics in a heterogeneous, spatially extended neural mass model. Neuroimage. 2011;55(3):920–932.
- 13. Suffczynski P, Kalitzin S, Da Silva FL. Dynamics of non-convulsive epileptic phenomena modeled by a bistable neuronal network. Neuroscience. 2004;126(2):467–484.
- 14. Saggio ML, Spiegler A, Bernard C, Jirsa VK. Fast–slow bursters in the unfolding of a high codimension singularity and the ultra-slow transitions of classes. The Journal of Mathematical Neuroscience. 2017;7(1):7.
- 15. Wang Y, Goodfellow M, Taylor PN, Baier G. Phase space approach for modeling of epileptic dynamics. Physical Review E. 2012;85(6):061918.
- 16. Taylor PN, Baier G. A spatially extended model for macroscopic spike-wave discharges. Journal of computational neuroscience. 2011;31(3):679–684.
- 17. Jirsa VK, Stacey WC, Quilichini PP, Ivanov AI, Bernard C. On the nature of seizure dynamics. Brain. 2014;137(8):2210–2230.
- 18. Desroches M, Faugeras O, Krupa M. Slow-fast transitions to seizure states in the Wendling-Chauvel neural mass model. Opera Medica et Physiologica. 2016;(3-4).
- 19. Chizhov AV, Zefirov AV, Amakhin DV, Smirnova EY, Zaitsev AV. Minimal model of interictal and ictal discharges “Epileptor-2”. PLoS computational biology. 2018;14(5):e1006186.
- 20. Baier G, Taylor PN, Wang Y. Understanding epileptiform after-discharges as rhythmic oscillatory transients. Frontiers in computational neuroscience. 2017;11:25.
- 21. Ersöz EK, Modolo J, Bartolomei F, Wendling F. Neural mass modeling of slow-fast dynamics of seizure initiation and abortion. PLoS computational biology. 2020;16(11):e1008430.
- 22. Chang WC, Kudlacek J, Hlinka J, Chvojka J, Hadrava M, Kumpost V, et al. Loss of neuronal network resilience precedes seizures and determines the ictogenic nature of interictal synaptic perturbations. Nature neuroscience. 2018;21(12):1742. pmid:30482946
- 23. Barbarosie M, Avoli M. CA3-driven hippocampal-entorhinal loop controls rather than sustains in vitro limbic seizures. Journal of Neuroscience. 1997;17(23):9308–9314.
- 24. Karoly PJ, Freestone DR, Boston R, Grayden DB, Himes D, Leyde K, et al. Interictal spikes and epileptic seizures: their relationship and underlying rhythmicity. Brain. 2016;139(4):1066–1078. pmid:26912639
- 25. Avoli M, de Curtis M. GABAergic synchronization in the limbic system and its role in the generation of epileptiform activity. Progress in neurobiology. 2011;95(2):104–132.
- 26. Huberfeld G, de La Prida LM, Pallud J, Cohen I, Le Van Quyen M, Adam C, et al. Glutamatergic pre-ictal discharges emerge at the transition to seizure in human epilepsy. Nature neuroscience. 2011;14(5):627–634. pmid:21460834
- 27.
Hoppensteadt FC, Izhikevich EM. Weakly connected neural networks. vol. 126. Springer Science & Business Media; 2012.
- 28.
Izhikevich EM. Dynamical systems in neuroscience. MIT press; 2007.
- 29. Desroches M, Guckenheimer J, Krauskopf B, Kuehn C, Osinga HM, Wechselberger M. Mixed-mode oscillations with multiple time scales. Siam Review. 2012;54(2):211–288.
- 30.
Berglund N, Gentz B. Noise-induced phenomena in slow-fast dynamical systems: a sample-paths approach. Springer Science & Business Media; 2006.
- 31.
Grasman J. Asymptotic methods for relaxation oscillations and applications. vol. 63. Springer Science & Business Media; 2012.
- 32.
Mishchenko E. Differential equations with small parameters and relaxation oscillations. vol. 13. Springer Science & Business Media; 2013.
- 33. Krupa M, Szmolyan P. Extending geometric singular perturbation theory to nonhyperbolic points—fold and canard points in two dimensions. SIAM journal on mathematical analysis. 2001;33(2):286–314.
- 34. Guckenheimer J. Isochrons and phaseless sets. Journal of Mathematical Biology. 1975;1(3):259–273.
- 35. Fenichel N. Persistence and smoothness of invariant manifolds for flows. Indiana Univ Math J. 1971/1972;21:193–226.
- 36. Tikhonov AN. Systems of differential equations containing small parameters in the derivatives. Matematicheskii sbornik. 1952;73(3):575–586.
- 37. Gradshtein I. Application of AM Lyapunov’s theory of stability to the theory of differential equations with small coefficients in the derivatives. Matematicheskii Sbornik. 1953;74(2):263–286.
- 38.
Pontryagin LS, Rodygin L. Periodic solution of a system of ordinary differential equations with a small parameter in the terms containing derivatives. In: Doklady Akademii Nauk. vol. 132. Russian Academy of Sciences; 1960. p. 537–540.
- 39.
Sacré P, Franci A. Singularly perturbed phase response curves for relaxation oscillators. In: 2016 IEEE 55th Conference on Decision and Control (CDC). IEEE; 2016. p. 4680–4685.
- 40. Langfield P, Krauskopf B, Osinga HM. Solving Winfree’s puzzle: The isochrons in the FitzHugh-Nagumo model. Chaos: An Interdisciplinary Journal of Nonlinear Science. 2014;24(1):013131.
- 41. Wendling F, Bartolomei F, Bellanger J, Chauvel P. Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition. European Journal of Neuroscience. 2002;15(9):1499–1508.
- 42. El Houssaini K, Bernard C, Jirsa VK. The Epileptor model: a systematic mathematical analysis linked to the dynamics of seizures, refractory status epilepticus and depolarization block. Eneuro. 2020;.
- 43. Proix T, Bartolomei F, Guye M, Jirsa VK. Individual brain structure and modelling predict seizure propagation. Brain. 2017;140(3):641–654.
- 44. Ashwin P, Coombes S, Nicks R. Mathematical frameworks for oscillatory network dynamics in neuroscience. The Journal of Mathematical Neuroscience. 2016;6(1):2.
- 45.
Schultheiss NW, Prinz AA, Butera RJ. Phase response curves in neuroscience: theory, experiment, and analysis. Springer Science & Business Media; 2011.
- 46. Rinzel J, Ermentrout GB. Analysis of neural excitability and oscillations. Methods in neuronal modeling. 1998;2:251–292.
- 47. Izhikevich EM. Phase equations for relaxation oscillators. SIAM Journal on Applied Mathematics. 2000;60(5):1789–1804.
- 48. Pérez-Cervera A, M-Seara T, Huguet G. Global phase-amplitude description of oscillatory dynamics via the parameterization method. Chaos: An Interdisciplinary Journal of Nonlinear Science. 2020;30(8):083117.
- 49. Gast R, Schmidt H, Knösche TR. A mean-field description of bursting dynamics in spiking neural networks with short-term adaptation. Neural Computation. 2020;32(9):1615–1634.
- 50. El Houssaini K, Ivanov AI, Bernard C, Jirsa VK. Seizures, refractory status epilepticus, and depolarization block as endogenous brain activities. Physical Review E. 2015;91(1):010701.
- 51. Proix T, Bartolomei F, Chauvel P, Bernard C, Jirsa VK. Permittivity coupling across brain regions determines seizure recruitment in partial epilepsy. Journal of Neuroscience. 2014;34(45):15009–15021.
- 52. Fisher RS, Velasco AL. Electrical brain stimulation for epilepsy. Nature Reviews Neurology. 2014;10(5):261–270.
- 53. Guillamon A, Huguet G. A computational and geometric approach to phase resetting curves and surfaces. SIAM Journal on Applied Dynamical Systems. 2009;8(3):1005–1042.
- 54. Fenichel N. Asymptotic stability with rate conditions. Indiana University Mathematics Journal. 1974;23(12):1109–1137.
- 55. Pérez-Cervera A, Tere M, Huguet G, et al. A geometric approach to Phase Response Curves and its numerical computation through the parameterization method. Journal of Nonlinear Science. 2019;29(6):2877–2910.