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

Fig 1

Phase space for relaxation oscillators.

The slow manifold, $$ \mathcal{S}$$, is a S-shaped curve having two stable branches $$ {\mathcal{S}}_{\pm }^a$$ (solid red line) and one repelling $$ {\mathcal{S}}^r$$ (dashed red line) (see Eq (4)). Stable and unstable branches of $$ \mathcal{S}$$ are separated by the fold points $$ {\mathcal{S}}_{\pm }^f$$. A given point, p, (see A) will quickly converge to the attracting branch of the slow manifold $$ {\mathcal{S}}_{-}^a$$ (see B). Then, it evolves along $$ {\mathcal{S}}_{-}^a$$ following (6) until reaching the fold point $$ {\mathcal{S}}_{-}^f$$ (see C) where it traverses fast to the other branch $$ {\mathcal{S}}_{+}^a$$ (see D). Then, following again the slow dynamics, the trajectory approaches $$ {\mathcal{S}}_{+}^f$$ (see E) where it goes back to $$ {\mathcal{S}}_{-}^a$$ (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:

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, t_s, follows a normal distribution with mean value, T_s, and standard deviation, σ: $$ t_s\sim \mathcal{N}(T_s,{\sigma }^2)$$. 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, T_s, of IEDs (see Fig 2A and 2B).

Fig 2

The antagonistic effect of IEDs on the transition to seizure.

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 $$ t_s\sim \mathcal{N}(50,3^2)$$ the system goes to seizure (panel A). By contrast, for a random train with amplitude A = 0.5 and $$ t_s\sim \mathcal{N}(30,2^2)$$ 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, T_s (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, T_s (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

Let us now consider the effect of an instantaneous delta-like pulse over the T-periodic limit cycle Γ,

As Fig 3A shows, the isochrons $$ {\mathcal{I}}_{\theta }$$ of Γ_pheno portray the distribution of phases along the basin of attraction $$ \mathcal{W}$$. 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 $$ \dot{a}=0$$, 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 $$ {\mathcal{S}}_{\pm }^f$$ the PRCs will show larger values (see zoom window in Fig 3C).

Fig 3

Isochrons and PRCs for the phenomenor model (7).

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 t₁ (t₂) 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 z₁ = γ(θ₁) = γ(t₁/T) (z₂ = γ(θ₂) = γ(t₂/T)) displaces the trajectory to a point $$ {\overline{z}}_1\in {\mathcal{I}}_{{\overline{\theta }}_1}$$ ($$ {\overline{z}}_2\in {\mathcal{I}}_{{\overline{\theta }}_2}$$). Since $$ {\overline{\theta }}_1<{\theta }_1$$ ($$ {\overline{\theta }}_2>{\theta }_2$$) the perturbation causes a phase shift $$ \Delta {\theta }_1={\overline{\theta }}_1-{\theta }_1<0$$ ($$ \Delta {\theta }_2={\overline{\theta }}_2-{\theta }_2>0$$) 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 $$ {\mathcal{S}}^r$$ in (4)—so they advance phase (see points z₁ and z₂ 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 $$ \gamma \left({\theta }^{*}\right)+A\in {\mathcal{S}}^r$$ (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, T_s, generate the different seizure propensity regimes in Fig 2C. As isochrons in Fig 3 show, unless they bring trajectories above $$ {\mathcal{S}}^r$$, 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 T_s, 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, T_s, 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 $$ \mathcal{O}\left(1\right)$$ geometry of isochrons $$ {\mathcal{I}}_{\theta }$$

Next, we discuss some generalities shaping the isochrons of planar relaxation oscillators. To begin, it is worth recalling that if two points $$ \overline{z}\in \mathcal{W}$$ and z = γ(θ) belong to the same isochron, $$ {\mathcal{I}}_{\theta }$$, 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 $$ \mathcal{O}\left(1\right)$$ distances in the fast direction and $$ \mathcal{O}\left(ϵ\right)$$ 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 $$ \overline{z}\in \mathcal{W}$$. In a first approximation one can assume that the convergence of $$ \overline{z}$$ 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 $$ {\mathcal{F}}_y$$. 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 $$ \mathcal{O}\left(ϵ\right)$$ corrections of $$ {\mathcal{F}}_y$$. Indeed, it is worth to note that generalities determining the sign of those $$ \mathcal{O}\left(ϵ\right)$$ 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 $$ t_h=\mathcal{O}\left(ϵ\right|logϵ\left|\right)$$, all orbits starting in a neighbourhood $$ \mathcal{O}\left(1\right)$$ of the slow manifold $$ \mathcal{S}$$ will have reached a neighbourhood of $$ \mathcal{O}\left(ϵ\right)$$ of $$ \mathcal{S}$$.

Once we know the time, $$ t_h=\mathcal{O}\left(ϵ\right|logϵ\left|\right)$$, needed to converge, we can compute the motion of the converging point $$ \overline{z}\in \mathcal{W}$$ in the slow direction. The travelled distance in the y direction by $$ \overline{z}$$ to approach a $$ \mathcal{O}\left(ϵ\right)$$ neighbourhood of Γ_ϵ is given by

where $$ {\phi }_t\left(\overline{x}\right)$$ refers to the solution of the layer system (2). During the time, t_h, needed to converge, the point z = γ(θ) on the limit cycle has travelled a distance Δy given by

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, $$ \Delta y-\Delta \overline{y}$$, will determine the sign of the $$ \mathcal{O}\left(ϵ\right)$$ correction for isochron $$ {\mathcal{I}}_{\theta }$$ at the point $$ \overline{z}$$. If we write the expression for $$ \Delta y-\Delta \overline{y}$$:

we can see that the difference $$ \Delta y-\Delta \overline{y}$$ is directly determined by the difference between the speeds in the slow direction for the base z and converging $$ \overline{z}$$ points during the time $$ t_h=\mathcal{O}\left(ϵ\right|logϵ\left|\right)$$ needed for approaching Γ_ϵ. Basically, since both points z and $$ \overline{z}$$ have to meet at the same point after the same time, the one travelling slower, needs to travel less distance. The difference, $$ \Delta y-\Delta \overline{y}$$, corresponds exactly to the $$ \mathcal{O}\left(ϵ\right)$$ correction to $$ {\mathcal{F}}_y$$ (see Fig 4).

Fig 4

The slow vector field shapes the isochrons for relaxation oscillators.

In the limit ϵ → 0 isochrons are lines of y constant denoted by $$ {\mathcal{F}}_y$$. However, since ϵ ≠ 0 but small, the isochrons are $$ \mathcal{O}\left(ϵ\right)$$ perturbations of $$ {\mathcal{F}}_y$$. As we show in the right panel, the sign of the $$ \mathcal{O}\left(ϵ\right)$$ corrections depends on the difference of speeds between the converging point $$ \overline{z}$$ and the base point z during the convergence time t_h. In this case, to approach Γ_ϵ, $$ \overline{z}$$ has to cross layers of x whose values are smaller than the ones surrounding Γ_ϵ. For this reason $$ \overline{z}$$ travels slower than z. Since $$ \overline{z}$$ and z have to meet after a time t_h at the same point on Γ_ϵ, but $$ \overline{z}$$ travels slower than z, then $$ \overline{z}$$ needs to travel a short distance. This determines the sign of the $$ \mathcal{O}\left(ϵ\right)$$ correction. Furthermore, if the slow vector field is monotonous along the fast direction, the farther the point $$ \overline{z}$$, the slower (faster) it travels, so the slope of the isochrons will have the same sign for all the points $$ \overline{z}\in {\mathcal{I}}_{\theta }$$ 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 $$ \overline{z}\in \mathcal{W}$$. 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 $$ \overline{z}\in {\mathcal{I}}_{\theta }$$ satisfying fast convergence.

The asymptotic phase defined in (11) allows to assign a phase to any point $$ z\in \mathcal{W}$$, by defining the following function Θ(z)

whose level curves indeed correspond to the isochrons. Let us assume we can invert Θ(x, y), so we can define the following function

The slope of isochron $$ {\mathcal{I}}_{\theta }$$, which we denote by $$ {\mathcal{K}}_{\theta }$$, is then given by

as we also have

we can write the following expression for the slope $$ {\mathcal{K}}_{\theta }$$

As the term $$ \Delta y-\Delta \overline{y}$$ can be written in integral form (see Eq (17)), the slope $$ {\mathcal{K}}_{\theta }$$, can be evaluated as the derivative of the difference of two sums (integrals)

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 $$ {\mathcal{K}}_{\theta }$$. If the trajectory followed by the approaching point, satisfies $$ g(\overline{x}\left(t_0\right),\overline{y}\left(t_0\right))<g(\overline{x}(t_0+\delta t),\overline{y}(t_0+\delta t))$$ for 0 < δt ≤ t_h, then, the second integral will be smaller than the first one. Since this difference will increase with $$ \overline{x}$$, then $$ {\mathcal{K}}_{\theta }>0\forall \overline{x}>{\gamma }^x(t/T)$$. Furthermore, the larger the changes in g(x, y), the larger the slope. We remark that in the case $$ g(\overline{x}\left(t_0\right),\overline{y}\left(t_0\right))>g(\overline{x}(t_0+\delta t),\overline{y}(t_0+\delta t))$$, we can argue identically to obtain $$ {\mathcal{K}}_{\theta }<0\forall \overline{x}>{\gamma }^x(t/T)$$.

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 $$ {\mathcal{I}}_{\theta }$$ at a given point $$ \overline{z}$$ depends on the difference of speeds between $$ \overline{z}$$ 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 $$ {\mathcal{P}}_{pheno}$$ 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.

Fig 5

Relationship between the curvature of isochrons and the values for the slow vector field.

For the phenomenological epilepsy model (7) with the set of parameters $$ {\mathcal{P}}_{pheno}$$ in (8) the figure shows: (A) Limit cycle Γ_pheno and its isochrons $$ {\mathcal{I}}_{\theta }$$ (left). (B) Values of the slow vector field (corresponding to $$ \dot{a}$$ in (7)) for points $$ z\in \mathcal{W}$$.

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 $$ {\mathcal{S}}^r$$ 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 $$ g_x^{\prime }(x,y)<0$$. In this case, we know that perturbations acting over points on the lower branch not crossing $$ {\mathcal{S}}^r$$ 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 $$ \gamma \left({\theta }^{*}\right)+A\in {\mathcal{S}}^r$$. 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 $$ {\mathcal{S}}^r$$ [40], we expect the bending of a particular isochron to increase as it approaches $$ {\mathcal{S}}^r$$. 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.

Fig 6

PRC of pulses A > 0 for relaxation oscillators.

Next we sketch the PRCs for pulses of amplitude A > 0 for the case $$ g_x^{\prime }(x,y)<0$$. For phases θ < θ*, where θ* satisfies $$ \gamma \left({\theta }^{*}\right)+A\in {\mathcal{S}}^r$$, due to the slope of isochrons the effect of the pulses will be to delay trajectories. Since isochrons approach the unstable point P through $$ {\mathcal{S}}^r$$, the closer the phase θ to θ*, the larger the bending of the isochrons and thus the larger the corresponding delay value. For phases $$ {\theta }^{*}<\theta <{\theta }_{-}^f$$, 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 $$ {\theta }_{-}^f$$, 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 f_A(θ) = θ_new, where the map f_A(θ) writes as

If the perturbation was a train of periodic pulses with an inter stimulus interval given by T_s, we can describe the phase dynamics of the system by

where θ₀ = θ. The fixed points of the above map (25), which are given by

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 T_s ≈ 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 $$ \overline{z}$$. If we denote by t_h the time that $$ \overline{z}$$ needs to approach Γ_ϵ, we need $$ {\varphi }_{t_h}\left(\overline{z}\right)=\gamma \left(\overline{\theta }\right)$$ with $$ \overline{\theta }<\theta $$. That is, we need the perturbed trajectory to reach the cycle at a previous phase. Assuming fast convergence, we can write

Since we need $$ \overline{\theta }<\theta $$ 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

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 T_s = −TΔθ > 0 (with $$ \Delta \theta =\overline{\theta }-\theta $$) 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.

Fig 7

Mechanism preventing the emergence of seizures.

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 T_s = 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 T_s periodic train of pulses. We can see how the locking preventing the transition to seizure starts for values $$ \overline{A}\approx 0.1$$, which is the approximate distance between the lower branch and the slow-nullcline. Furthermore, for a fixed amplitude $$ A>\overline{A}$$, if we consider the maximum delay value (denoted by Δθ*) of the corresponding PRC and compute the inter pulse interval value given by $$ T_s^{*}=T\Delta {\theta }^{*}$$, it is clear that for inter-pulse intervals $$ T_s>T_s^{*}$$, the system is likely to transition to seizure because the delays are not large enough to stop the system. Therefore, we expect the pair $$ (A,T_s^{*})$$ to delimit the locking regime. By computing the PRCs for all the amplitude values satisfying $$ A>\overline{A}$$, we can calculate the corresponding $$ T_s^{*}$$ values and thus generate a curve in the (A, T_s) 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).

Fig 8

Response of perturbations for the phenomenor.

We plot the change in the seizure rate Δ for a random train of pulses following a Gaussian distribution of mean time T_s and standard deviation σ denoted as $$ \mathcal{N}(T_s,{\sigma }^2)$$. Panels (A) and (B) correspond to the deterministic periodic case $$ \mathcal{N}(T_s,0)$$ and to the random case $$ \mathcal{N}(T_s,0.05T_s)$$, 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 T_s for t_s values in the distribution $$ t_s\sim \mathcal{N}(T_s,{\sigma }^2)$$. 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 $$ T_s^{*}$$, such that perturbations $$ t_s>T_s^{*}$$ will lead to seizures. Therefore, the robustness of the deterministic locking states to noise decreases, as the probability of t_s being above $$ T_s^{*}$$ 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]:

Table 1

For the set of parameters $$ {\mathcal{P}}_i$$, the system will display a limit cycle Γ_i of period T_i.

Different parameters for the reduced 2D Epileptor model in (29).

	τz	v0	Iapp	c	s	Lim. Cycle	Period
$$	1/2857	-2	3.1	-4	-1	Γ+	T+ ≈ 2181.6
$$	1/2857	-1.5	3.1	-16	-1	Γ0	T0 ≈ 695.7
$$	1/2857	-0.1	3.1	2.4	1	Γ−	T− ≈ 7333.3

Identically as the phenomenor, since the time constant for the z variable is small τ_z ≪ 1, and $$ \dot{v}=0$$ describes a cubic curve, the three sets of parameters $$ {\mathcal{P}}_{+},{\mathcal{P}}_0$$ and $$ {\mathcal{P}}_{-}$$ lead to relaxation oscillators denoted as Γ₊, Γ₀ and Γ₋ respectively. The three different sets of parameters $$ {\mathcal{P}}_{+},{\mathcal{P}}_0$$ and $$ {\mathcal{P}}_{-}$$ 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 $$ {\mathcal{P}}_{+},{\mathcal{P}}_0$$ and $$ {\mathcal{P}}_{-}$$ 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 $$ {\mathcal{P}}_{+},{\mathcal{P}}_0$$ and $$ {\mathcal{P}}_{-}$$. 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.

Fig 9

Isochrons and PRCs for the reduced epileptor.

For the sets of parameters $$ {\mathcal{P}}_{+},{\mathcal{P}}_0$$ and $$ {\mathcal{P}}_{-}$$ in Table 1 we show: Limit cycle Γ₊, Γ₀ and Γ₋ and its isochrons $$ {\mathcal{I}}_{\theta }$$ (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 Γ₊, Γ₀ 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 $$ {\mathcal{P}}_{+},{\mathcal{P}}_0$$ and $$ {\mathcal{P}}_{-}$$ 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 T_s and standard deviation σ, denoted by $$ \mathcal{N}(T_s,{\sigma }^2)$$ and compute the change of the seizure rate Δ for a train with $$ \mathcal{N}(T_s,0)$$ and $$ \mathcal{N}(T_s,0.05T_s)$$.

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 Γ₀ 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 T_s = −TΔθ values for which Γ₋ shows locking is even larger than for Γ₀ and Γ₊. We also remark the good agreement between the bifurcation curves of map (25) and the areas suppressing the transition to seizure.

Fig 10

Response of perturbations for the reduced epileptor (29).

We show the change in the seizure rate Δ for a random train of pulses whose mean inter impulse interval follows a normal distribution $$ \mathcal{N}(T_s,{\sigma }^2)$$ with a mean time T_s and standard deviation σ. Panels A, B, C correspond to the sets of parameters $$ {\mathcal{P}}_{+},{\mathcal{P}}_0$$ and $$ {\mathcal{P}}_{-}$$ in Table 1. Left Figures correspond to the periodic case $$ \mathcal{N}(T_s,0)$$ and right Figures to the random case $$ \mathcal{N}(T_s,0.05T_s)$$. 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 $$ T_s^{*}=T\Delta {\theta }^{*}$$, (where Δθ* corresponds with the maximal delay value of the PRC) is or not within the width σ of the distribution $$ t_s\sim \mathcal{N}(T_s,{\sigma }^2)$$. The higher the probability of occurrence of $$ t_s>T_s^{*}$$ values, the likely is the system to switch to the upper branch. The differences in the resilience of the deterministic locking areas for Γ₊, Γ₀ 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 T_s = TΔθ. The shorter the period T, the smaller the critical $$ T_s^{*}=-T\Delta {\theta }^{*}$$ value. Since in the three cases the t_s distributions have the same width, the smaller the critical $$ T_s^{*}$$ value, the higher the probability of occurrence of $$ t_s>T_s^{*}$$ values. As a consequence, the resilience of locking states for Γ₊ and Γ₀ is weaker than for Γ₋ in which the distribution T_s = −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.

Fig 11

Slow vector field for the reduced epileptor and the phenomenor.

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 $$ g_x^{\prime }(x,y)$$.

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 $$ \mathcal{O}\left(\tau \right)$$ (see Fig 4). However, from our analysis it also follows that the bending of the isochrons, although being contained in $$ \mathcal{O}\left(\tau \right)$$ 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 $$ g_x^{\prime }(x,y)$$, 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 T_s 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.

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 $$ S_{+}^f$$ (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 T_s 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 $$ t_s\sim \mathcal{N}(T_s,{\sigma }^2)$$.

From the adopted criteria for counting seizures it follows that very large perturbations might cause or constitute a seizure per se independently of the T_s value. For that reason, since we were interested in the relationship of both the amplitude A and the inter pulse interval T_s 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 T_s value).

Computation of isochrons

To compute isochrons of slow-fast systems, we assume we have an analytic vector field $$ \dot{z}=Z\left(z\right)$$ 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 z₀ ∈ Γ and the period T. Then, we integrate the system (9) with initial condition z(0) = z₀ for a time T to obtain $$ z(\theta T)≕\gamma (\theta )\text{for}\theta \in [0,1)$$.

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 $$ \mathcal{O}\left(e^{\mathrm{-}k/ϵ}\right)$$ 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]:

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 T_rel 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 T_rel. 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 $$ \overline{z}=z+A$$. Then, after a time T_rel large enough, the trajectory will be again on the limit cycle but with another phase $$ \overline{\theta }$$. Mathematically

The idea of the method is to obtain f_A(θ) 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 f_A(θ) an approximate solution of Eq (31), we perform the following operations:

Compute E(θ) = F_A(γ(θ)) − γ(f_A(θ)).

Compute ∂_θ γ(f_A(θ)) = TZ(γ(f_A(θ))).

Compute $$ \Delta f_A=\frac{<{\partial }_{\theta }\gamma \left(f_A\left(\theta \right)\right),E\left(\theta \right)>}{<{\partial }_{\theta }\gamma \left(f_A\left(\theta \right)\right),{\partial }_{\theta }\left(f_A\left(\theta \right)\right)>}$$.

Set f_A(θ) ← f_A(θ) + Δf_A(θ).

Repeat steps 1-4 until the error E is smaller than the established tolerance. Then PRC(A, θ) = f_A(θ) − (θ + T_rel/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 f_A(θ) = θ + T_rel/T as initial seed. Then, ∇Θ(γ(θ)) = PRC(A, θ)/A.