## Abstract

One characteristic of epilepsy is the variety of mechanisms leading to the epileptic state, which are still largely unknown. Refractory status epilepticus (RSE) and depolarization block (DB) are other pathological brain activities linked to epilepsy, whose patterns are different and whose mechanisms remain poorly understood. In epileptogenic network modeling, the Epileptor is a generic phenomenological model that has been recently developed to describe the dynamics of seizures. Here, we performed a detailed qualitative analysis of the Epileptor model based on dynamical systems theory and bifurcation analysis, and investigate the dynamic evolution of “normal” activity toward seizures and to the pathological RSE and DB states. The mechanisms of the transition between states are called bifurcations. Our detailed analysis demonstrates that the generic model undergoes different bifurcation types at seizure offset, when varying some selected parameters. We show that the pathological and normal activities can coexist within the same model under some conditions, and demonstrate that there are many pathways leading to and away from these activities. We here archive systematically all behaviors and dynamic regimes of the Epileptor model to serve as a resource in the development of patient-specific brain network models, and more generally in epilepsy research.

- bifurcation analysis
- depolarization block
- dynamical systems theory
- epilepsy
- neural mass model
- refractory status epilepticus

## Significance Statement

Epilepsy is characterized by patient-specific electrophysiological discharges. The range of mechanisms and pathways leading to the same type of seizure, however, is large. The Epileptor model has found many applications in epilepsy research and clinical applications, because it allows the classification and dynamic modeling of seizure types independent of the knowledge of its underlying biophysical mechanisms. It is based purely on the dynamic features of the seizure. We provide here a complete functional atlas of all possible behaviors of the Epileptor model, which serves as a useful resource in modeling brain networks in epilepsy. More, we explore the contribution of the Epileptor model to better understand the dynamics of the refractory status epilepticus and depolarization block phenomena, which are linked to epilepsy.

## Introduction

Epilepsy is a condition of the nervous system in which neuronal populations manifest as repeated epileptic seizures lasting a few minutes. These seizures are spontaneous and commonly accompanied by fast abnormal discharges (10 ms timescale), after which the brain activity slowly returns to normal. Epileptic seizures may present in different forms, as well as transitions from and to these pathological states. Some seizures are controlled by medication, particularly antiepileptic drugs (AEDs). However, there are seizures that last >1 h without returning to baseline, and do not respond to AEDs (31-43%), resulting then in so-called refractory status epilepticus (RSE) (Mayer et al., 2002; Holtkamp et al., 2005; Rossetti et al., 2005). Patients with RSE are at an increased risk of death. Indeed, the short-term mortality rates of RSE are estimated to be between 16% and 39% (Mayer et al., 2002; Holtkamp et al., 2005; Rossetti et al., 2005; Novy et al., 2010).

Other neuroelectric phenomena are linked to epilepsy, among which spreading depression (SD) is the most prominent. SD is characterized by a slowly propagating depolarization wave [or depolarization block (DB)] in neuronal networks, followed by a shutdown of brain activity (Pietrobon and Moskowitz, 2014). DB is a state in which the neuronal membrane is depolarized, but neurons stop firing. Spreading depression may occur during migraine and some seizures (Rogawski, 2008; Pietrobon and Moskowitz, 2014), and was first described by Leo (1944), who observed a depression of electroencephalographic (EEG) activity that moved across the cortex.

The mechanisms underlying the genesis of epilepsy, DB, and RSE are still largely unknown. In this article, we take an integrative approach toward the understanding of these neuroelectric phenomena. Our thinking is inspired by the wish to understand the underlying dynamic underpinnings rather than the biophysiological basis.

In this context, a large number of experimental and computational models have been proposed to clarify the basic mechanisms of seizures and DB. Most computational models rely on biophysically realistic parameters, trying to reproduce experimental data (Kager et al., 2000; Traub et al., 2001; Destexhe, 2008; Cressman et al., 2009; Ullah et al., 2009). Although these models provide important advances in seizures and DB research, they rarely produce general rules. Furthermore, there has been no attempt at evaluating whether or not seizures, RSE, DB, and normal brain activities can coexist, and, if so, under what conditions.

The previously mentioned models are rooted in physiological mechanisms generating a fairly limited range of behaviors. The physiological foundation is critical when the intended therapeutic intervention shall make use of the physiology. Examples include the identification of signaling pathways leading to the control of neurotransmitters linked to excitability such as lamotrigine or topiramate acting on calcium channels. Another type of intervention acts on the brain as a network and harnesses the capacity to modulate networks, such as stimulation, resection, or disconnection. In this case, physiological realism of a mechanism is to be replaced by dynamic realism as the network communication depends more on the type of signal rather than how it is generated. These types of models are phenomenological. A neural mass model of partial seizures called Epileptor was previously developed, which has found many applications in brain network modeling of epilepsy patients and is the network node model used in the European clinical trial EPINOV (www.epinov.com; Jirsa et al., 2014). Given the wide application of this model and its relevance for clinical research, we here archive in detail its dynamic repertoire. The Epileptor comprises one susbsystem (called subsystem 1) with two state variables responsible for generating fast discharges, another subsystem (called subsystem 2) with two state variables generating sharp-wave events (SWEs). The subsystem 1 (fast) and subsystem 2 (slow) are linked to a state variable, *z*, evolving on a very slow timescale called the permittivity variable (Jirsa et al., 2014). Interestingly, the transition from and to the pathologic states can occur autonomously, under the slow *z* evolution. Fast discharges and sharp-wave events are pathological features commonly associated with seizures despite their different forms. The goal of this article is, first, to perform a systematic mathematical analysis of the Epileptor; and, second, to determine the range of behaviors present in the Epileptor, with the further reaching goal to ask whether “normal” brain activities, seizures, RSE, and DB can coexist within the same model. To this aim, we present a qualitative analysis of the Epileptor model based on dynamical systems theory and bifurcation analysis.

## Materials and Methods

We provide here a detailed bifurcation analysis of the Epileptor model, which is a neural mass model of partial seizures, to analyze seizure dynamics (Jirsa et al., 2014). The Epileptor model consists of a system of coupled nonlinear differential equations with five state variables. It comprises two 2D subsystems and one slow variable *z*. Subsystems 1 and 2 are responsible for fast discharges and SWEs, respectively (Jirsa et al., 2014). Analysis of the separated subsystems 1 and 2 was performed to provide more in-depth information on the Epileptor dynamics.

### Analysis of the Epileptor

#### Epileptor equations

The Epileptor equations generate seizure-like events (SLEs) written in the following form: (1) (2) (3) (4) (5)where (6) (7) (8)

In what follows, unless otherwise stated, the values of the parameters are as follows:

Subsystem 1: *a* = 1; *b* = 3; *c *= 1; *d *=* *5;
; *m *=* * 0

Subsystem 2:
;
;
; *γ* = 0.01

Slow *z* dynamics: *r * = * *0.00,035; *s * = * *4; *x*_{0} = −1.6.

Equations 1 and 2 describe the subsystem 1 dynamics and equations; Equations 4 and 5 describe subsystem 2. The state variables *x*_{1} and *y*_{1} comprise the subsystem 1 responsible for fast discharges, and *x*_{2} and *y*_{2} the subsystem 2 involved in spike-wave events (Jirsa et al., 2014). The third differential equation represents a slow adaptation variable. The Epileptor model evolves under a slow timescale, *r* (small). Equations 1–3 represent a fast-slow subsystem.

#### Numerical integration

All stochastic simulations are performed with linear additive Gaussian white noise with a zero mean and a variance of 0.0025 using the Euler–Maruyama method. Deterministic simulations are performed using the Runge-Kutta method with a maximal time step of 0.01.

#### Finding equilibrium points and their stability

In order to understand the Epileptor dynamics, we identify the equilibrium points and investigate their stability. This analysis is performed on the deterministic Epileptor model. We find the equilibrium points
by solving
. We obtain a system of algebraic equations that we solve using the Matlab “solve” symbolic solver. The variables to solve for are the state variables of the Epileptor model, which are
. The inputs of the solve function are the parameter values of the Epileptor model. Here, we analyze the existence of the equilibrium points when varying the *m* and *x*_{0} parameters. The solutions are the vector equilibrium points. We first determine the equilibrium points of the whole system when the parameter *x*_{0} varies, for two values of the parameter *m* (*m * = * *0 and *m * = * *0.5). We determine the stability of the equilibrium points *E* by evaluating the eigenvalues of the Jacobian matrix *J* at the equilibrium point *E*, which is stable if all the real parts of the eigenvalues of *J* are negative (Izhikevich, 2007).

### Parameter space of equilibrium points

We explore the stability of the coexisting equilibrium points for a range of *m* and *x*_{0} values. For each values combination, we determine the equilibrium points using the solve function and their stability using the Jacobian matrix. The complete results are drawn within a (*m*, *x*_{0}) parameter space for equilibrium points, which consists of various areas, each of which is characterized by coexisting equilibrium points with different stability. Moreover, when the parameter
is varied, the Epileptor behavior changes. We therefore explore two parameter spaces of equilibrium points, for
(default Epileptor value) and
.

#### Bifurcation diagram

One important step in the analysis of a mathematical model is the geometrical analysis of bifurcations, which here was performed on the deterministic Epileptor model. First, we draw a (*z*, *x*_{1}) bifurcation diagram of the Epileptor for constant *m*, here *m * = * *0.5. The (*z*, *x*_{1}) bifurcation diagram comprises the equilibrium points of the Epileptor, where *x*_{1} is the first coordinate of the vector of equilibrium points that we determined with respect to *z*. All equilibrium point solutions were obtained using the Matlab solve symbolic solver, where only the (
) equations were considered for constant *z*. We discretized the space for 801 points in the *z*-dimension and performed numerical continuation computing the linear stability using the Jacobian matrix. The (*z*, *x*_{1}) curve is divided into different branches, of which each is numerically continued separately representing a different stability type of the equilibrium points solutions. Bifurcations were identified at stability changes for given *z*-values and are indicated by dots in the corresponding visualizations. As the main parameters *m* and
control the dynamics of the Epileptor, we typically present bifurcation diagrams for varying values of *m* and
. All bifurcation diagrams were verified using XPPAUT and explicit numerical simulations of the model system (typically 2000–4000 units of time, of which the initial transients are removed (Figs. 1, 2), for which representative trajectories were plotted in the corresponding diagrams.

#### Finding a fast-slow limit cycle

We found a large attractor in the (*y*_{1}, *ψ*, *z*) phase space below the SLE attractor shown in Figure 1 The *z* equation of the Epileptor model (Jirsa et al., 2014) was introduced as follows:
(9)

When using this *z* equation, the system evolves toward the attractor and then diverges with time, which we show in a (*z*, *x*_{1}) bifurcation diagram. Thus, we modify the *z* equation to stabilize the final state. We interpret graphically the divergence of the Epileptor by introducing the averaging method, which is used to locate the periodic orbits in the phase space. The periodic orbits are intersections of the *z*-nullcline and the
-curve. The *z*-nullcline is given by
, and the
-curve is the average value of *x*_{1} for each *z* constant written in the following form:
(10)where
(11)

We also determine the periodic orbits by using the Pontryagin’s averaging technique (Shilnikov and Kolomiets, 2008). To illustrate this technique, we introduce a slow averaged nullcline written as follows: (12)

where periodic orbits are the zeros of the
. We determine the stability of periodic orbits using the following derivative:
(13)which represents the dynamics of the averaged equation. A periodic orbit is stable when Equation 13 is negative, and graphically if the graph of
decreases at the given zero. First, we find the periodic orbits as the parameter *x*_{0} varies in a (*x*_{0}, *z*) bifurcation diagram of periodic orbits. We determine the bifurcation that results in the appearance or disappearance of the periodic orbits. Second, we explore two parameter spaces of periodic orbits for
and
.

#### Stabilizing equilibrium points in the Epileptor model

The parameter *x*_{0} can change the stability of equilibrium points and then can control the Epileptor behavior. Using the fact that a trajectory moves around the equilibrium point when it is stable, we use the bifurcation analysis to find the values of *x*_{0} for which an equilibrium point is stable. There is a qualitative distinction between the stable equilibrium points, which we discuss in more detail later.

#### Coexisting attractors in the Epileptor model

The Epileptor behavior depends on the stability of its equilibrium points. More, there is a pre-existing limit cycle (LC) attractor, which is a stable periodic orbit. In this case, the trajectory can have a behavior that is controlled by the stability of the equilibrium points or can jump to the limit cycle attractor. Thus, there is a coexistence of two attractors. We identify all coexisting attractors according to the values of the parameters *m* and *x*_{0}.

### Analysis of subsystem 1

We provide a detailed analysis of the subsystem 1 dynamics without coupling (i.e., subsystem 1 is not coupled to subsystem 2).

#### Subsystem 1 equations

The equations of subsystem 1 are given by:
(14)
(15)where
(16)
*z* is constant (
).

#### Subsystem 1 equilibrium points and stability

We analytically find the equilibrium points (*x*_{1}, *y*_{1}) by solving the following equations:
(17)

We determine the stability of the equilibrium points by evaluating the eigenvalues of the Jacobian matrix *J*. An equilibrium point is stable if all the real parts of the eigenvalues of *J* are negative (Izhikevich, 2007).

We graphically find the equilibrium points by intersecting the *x*_{1}- and *y*_{1}-nullclines. We determine the nullclines of subsystem 1 and show how to find the equilibrium points and their stability in a phase plane.

#### Subsystem 1 bifurcation diagram

When finding the equilibrium points of the subsystem 1, we observe a qualitative change of the phase plane, which is interpreted as a bifurcation. We identify the different types of bifurcations that exist according to *z*, which is considered as a parameter control, and the parameter *m*.

Using the bifurcation diagrams, we analyze these bifurcations and the qualitative behavior of the susbsystem 1 with respect to *z*. The bifurcation analysis is performed on the deterministic subsystem 1. First, we plot the whole bifurcation diagram of the subsystem 1 for *m * = * *0, which consists of two curves. The curves consist of the subsystem 1 equilibrium points for each value of *z*. When the parameter *m* changes, the shape of one of the curves changes on an interval of *z*. As a result, we then plot only this curve in bifurcation diagrams for *m * = * *0 and *m * = * *2, and describe the trajectories behavior for each value of *m*. Below, we summarize the types of bifurcations that exist with respect to *m*.

#### Fast-slow subsystem equations

The equilibrium points of subsystem 1 depend on *z*. When *m * = * *0 and *z * = * *3.1, three equilibrium points coexist: a saddle, a stable node and a stable focus (see *m * = * *0; see Fig. 39*A*). When *m * = * *1.5 and *z * = * *3.1, then equilibrium points coexist: a saddle, a stable node, and an unstable focus (see *m * = * *1.5; see Fig. 39*A*). A stable node and a stable focus are equilibrium points of a resting state and a nonoscillatory state, respectively. When the resting and nonoscillatory states coexist, the trajectories converge to one of them, depending on the initial conditions (see *m * = * *0 and *m * = * *1; see Fig. 39*A*).

To ensure that trajectories switch between two stable states, we introduce the following equation:
(18)to the subsystem 1
(19)
(20)
*r* and *s* are positive constant parameters (
).*z* changes the input
of the subsystem 1 to
. Indeed, when *z* decreases, then
is increased, and only a stable focus exists (see Fig. 39*B*). If *m* increases, the stable focus becomes unstable, surrounded by a stable limit cycle. Hence, only a stable limit cycle exists (see Fig. 39*B*). When *z* increases, then
is reduced, and the stable limit cycle coexists with a stable node (see Fig. 39*A*). The saddle acts as a separatrix (*S*) between them. When *z* further increases, the stable limit cycle disappears through a homoclinic bifurcation, HB (see Fig. 39*C*), and hence only a stable node exists. Thus, *z* mimics a slow adaptation of the subsystem 1 to produce resting and oscillatory states, and ways to switch between them. The dynamics of subsystem 1 is fast and *z* is a slow variable, hence Equations 18–20 represent a fast-slow subsystem.

#### Finding equilibrium points of the fast-slow subsystem

We analytically find the equilibrium points by solving Equations 19 and 20, where *z* is a solution of
(Eq. 18). We determine the stability of the equilibrium points by analyzing the Jacobian matrix *J*. The equilibrium point is stable if all the real parts of the eigenvalues of *J* are negative (Izhikevich, 2007). We graphically find the equilibrium points of the fast-slow subsystem by using the (*z*, *x*_{1}) bifurcation diagram and the *z*-nullcline, which is related to the parameter *x*_{0}. We show how the *z*-nullcline moves in the bifurcation diagram when varying *x*_{0}, and use the bifurcation diagram to find the equilibrium points of the fast-slow subsystem for each value of *x*_{0}.

#### Finding periodic orbits of the fast-slow subsystem

As follows from the Epileptor analysis, we determine the periodic orbits of the fast-slow subsystem as well as their stability by using the averaging method and the Pontryagin’s averaging technique. The slow averaged nullcline is given by the following: (21)where (22)

is a solution of the subsystem 1 with *z* constant.

Here, we plot the graph of
to explain how we find the periodic orbits and their stability. The periodic orbits are the zeros of the graph of
. It is stable if the graph of
decreases at the given zero. More, we find the periodic orbits as *x*_{0} varies and determine the bifurcation leading to their appearance or disappearance.

#### Stabilizing equilibrium points in the fast-slow subsystem

Our approach to stabilizing the equilibrium points of the fast-slow subsystem is the same as that of the Epileptor model: analyze bifurcation diagrams of the fast-slow subsystem and then find the values of *x*_{0} for which an equilibrium point is stable. We show how stable equilibrium points can correspond to different states.

#### Coexisting attractors in the fast-slow subsystem

The fast-slow subsystem exhibits a bistability of two attractors as the Epileptor model. The first attractor is the fast-slow limit cycle. The behavior of the second attractor depends on the parameters *m* and *x*_{0}, which control the stability of the equilibrium point. We discuss the coexisting attractors when *m* and *x*_{0} have different values.

### Analysis of subsystem 2

In this section, we analyze the dynamics of subsystem 2 (without coupling), which generates SWEs. We explore the equilibrium points, evaluating their stability and determining the bifurcation types.

#### Subsystem 2 equations

The subsystem 2 equations are given by the following: (23) (24)where (25)

#### Subsystem 2 equilibrium points and stability

We analytically find the equilibrium points (*x*_{2}, *y*_{2}) by solving the following equations:
(26)
(27)

We determine the stability of the equilibrium points by evaluating the eigenvalues of the Jacobian matrix *J*_{2}. *J*_{2} is defined at the equilibrium point *x*_{2}, which is stable if all the real parts of the eigenvalues of *J*_{2} are negative (Izhikevich, 2007).

The equilibrium points lie graphically at the intersection of the *x*_{2}- and *y*_{2}-nullclines. We determine the nullclines of the subsystem 2 and show how the equilibrium points and their stability changes as the parameter
varies.

#### Subsystem 2 bifurcation diagram

As the phase plane changes qualitatively when varying the parameter , we discuss this change in a bifurcation diagram where is the parameter control. The analysis is performed on the deterministic subsystem 2.

## Results

To get a better understanding of the dynamics of the generation and termination of epileptic seizures, we adopted a computational model that reproduces epileptic activity and perform a detailed analysis using a mathematical approach.

### The Epileptor

#### Epileptor dynamics

##### Epileptor model behavior

The Epileptor equations generate SLEs, which are characterized by an onset and offset (Jirsa et al., 2014). We plot time series of the Epileptor system *ψ*, subsystem 1 *ψ*_{1}, and subsystem 2 *ψ*_{2} in Figure 1*A*. We find the major elements of an SLE: onset, timescale, offset of SLEs, and their recurrence. Figure 1*A* shows that during the ictal phase, fast discharges decrease in frequency with time. Ictal states are separated by a period of normal brain activity (non-ictal state). State variables *x*_{1} and *y*_{1} are responsible for generating fast discharges in the ictal state with a fast timescale (Fig. 1*A*; see *ψ*_{1}). State variables *x*_{2} and *y*_{2} are responsible for generating the SWEs with an intermediate timescale (Fig. 1*A*; see *ψ*_{2}). We plot the Epileptor trajectory in a (*y*_{1}, *ψ*, *z*) phase space (Fig. 1*B*). Ictal and normal states (NSs) coexist and their coexistence necessitates a separation in the state space so-called “separatrix.” The ictal onset occurs when the trajectory collides with the separatrix after a transient normal state. The seizure offset occurs when the trajectory collides with the separatrix after a transient ictal state. The separatrix acts as a barrier between ictal and normal states. The slow state variable *z* is responsible for the alternation of both states under a slow timescale (
).

##### Equilibrium points

Using the Matlab solve symbolic solver, we localize the equilibrium points *E* as *x*_{0} varies in a (*x*_{0}, *z**) diagram for *m * = * *0.5 (Fig. 2*A*) and *m * = * *0 (Fig. 2*B*). *z*^{*} is the third component of the equilibrium points *E*. For *m * = * *0.5 and *m * = * *0, a unique equilibrium point exists when *x*_{0} = −1.6, while three equilibrium points coexist when *x*_{0} = −0.9 (Fig. 2*A*,*B*). The Jacobian matrix *J* of the Epileptor at the equilibrium point E is given by the following:
where
(28)
(29)
(30)
(31)
(32)

The stability of the Epileptor equilibrium points depends on *x*_{0} (Fig. 2*A*,*B*). For *x*_{0} = −1.6 and *m * = * *0 or *m * = * *0.5, the equilibrium point is a saddle. The stable manifold of the saddle equilibrium point corresponds to a separatrix between ictal and normal states (Fig. 1*B*). The trajectory behavior in Figure 1*A* corresponds to a recurrent alternation between ictal and normal states, which characterizes SLEs.

When

*x*_{0}= −2.5 and*m*= 0 or*m*= 0.5, the equilibrium point is a stable node. The trajectory behavior corresponds to a normal activity (Fig. 2*C*).When

*x*_{0}= −0.9 and*m*= 0, the equilibrium points are: an unstable focus, a saddle, and a stable focus (Fig. 2*B*). The Epileptor remains in a nonoscillatory state (Fig. 2*D*).When

*x*_{0}= −0.9 and*m*= 0.5, three equilibrium points coexist. The equilibrium points are as follows: one saddle and two unstable foci (Fig. 2*A*). The trajectory behavior corresponds to a periodic solution (Fig. 2*E*).

##### Parameter space of equilibrium points

The stability of equilibrium points changes when we vary *m* and *x*_{0} (Fig. 2*A*,*B*). Trajectories can exhibit: (1) a normal activity when an equilibrium point is a stable node, (2) a nonoscillatory state when an equilibrium point is a stable focus, (3) SLEs when a unique saddle equilibrium point exists, and (4) a periodic solution when equilibrium points are unstable. To explore the various coexisting equilibrium points, we plot a (*m*, *x*_{0}) parameter space of equilibrium points in Figure 3*A*.

The parameter space comprises nine areas. The equilibrium point is an unstable focus in area 1, a stable node in area 5, and a saddle in area 4. The Epileptor model has three equilibrium points in areas 2, 3, 6, 7, 8, and 9. In area 2, a stable node, a saddle and an unstable focus coexist. In area 3, a stable focus, a saddle, and an unstable focus coexist. In area 6, three saddles coexist. In area 7, two saddles and one unstable focus coexist. In area 9, two unstable foci and one saddle coexist. In area 8, an unstable focus, an unstable node, and a saddle coexist.

We plot a (*m*, *x*_{0}) parameter space of equilibrium points for
in Figure 3*B*. The parameter space comprises seven areas. The equilibrium point is an unstable focus in area 1, a stable node in area 5, a stable focus in area 6, and a saddle in area 7. The Epileptor model has three equilibrium points in areas 2, 3, and 4. In area 2, a stable node, a saddle, and an unstable focus coexist. In area 3, a stable focus, a saddle, and an unstable focus coexist. In area 4, two unstable foci and one saddle coexist.

##### Main parameters of the Epileptor model

*m* and *x*_{0} play roles of particular importance on the Epileptor dynamics. We showed above that the equilibrium points of the Epileptor change as *x*_{0} varies, thereby allowing the obtaining of different behaviors of the system, including normal activity, nonoscillatory state (Fig. 2), and the alternation between the normal and ictal states (Fig. 1).

*m* can be considered as a parameter that controls the Epileptor dynamic during the ictal period. The equilibrium points of the Epileptor model are graphically defined by the intersection of the nullclines. The *x*_{1}-nullcline (
) is a straight line for
. The sign of the slope of the straight line varies the direction of movement of the state variable *x*_{1} in the phase portrait (
). Then, the stability of the equilibrium points (
) changes when varying the slope, which depends on *m*, and hence the dynamic of the Epileptor can pass from fast discharges to nonoscillatory state during the ictal period as *m* varies.

Below, we theoretically demonstrate how the system can switch between normal and epileptic activities by varying *x*_{0} and using bifurcation diagrams. Moreover, we explore the significant role of *m* in controlling the frequency of discharges and in generating a DB.

##### Bifurcation diagram of the Epileptor for *m *=* *0.5

The Epileptor can alternate between ictal and normal states (see *ψ*; Fig. 1) when the equilibrium point is a saddle (see *x*_{0} = −1.6; Fig. 2). The alternation is interpreted mathematically as a bifurcation, which is a qualitative change of the system behavior (Izhikevich, 2007). We identify two bifurcation types in an SLE. The first bifurcation corresponds to the transition from normal to ictal activity. The second one corresponds to the transition from ictal to normal activity. We identify bifurcation types and find the equilibrium points in a bifurcation diagram. Using the Matlab solve symbolic solver, we plot a (*z*, *x*_{1}) bifurcation diagram of the Epileptor model in Figure 4 for *m * = * *0.5. *x*_{1} is the first component of the equilibrium points *E*. *z* acts as a control parameter (*r * = * *0). Figure 4 first shows a *Z*-shaped curve (*Z*-curve). *Z*-lower and *Z*-upper branches consist of stable nodes and unstable foci, respectively. The *Z*-middle branch separates *Z*-lower and *Z*-upper branches and consists of saddles. The *Z*-middle and *Z*-upper branches collide as *z* increases in a saddle-node bifurcation SN_{3} (Fig. 4, inset). The *Z*-middle and *Z*-lower branches collide as *z* decreases in a saddle-node bifurcation SN_{1}.

Moreover, there are two branches above the *Z*-shaped curve (Fig. 4, inset). The lower branch consists of saddles and terminates as *z* decreases in a saddle-node bifurcation SN_{4}. The upper branch comprises two sub-branches: one sub-branch (dashed) consists of unstable foci and another (dash-dotted) consists of saddles. The two sub-branches collide as *z* increases in a saddle-node bifurcation SN_{2}. The equilibrium points *E* of the Epileptor model lie at the intersection of the *z*-nullcline and the curve of equilibrium points (Fig. 4). The *z*-nullcline (
) is given by the following:
(33)and depends on *x*_{0}. The *z*-nullcline moves downward (*x*_{0} decreases) or upward (*x*_{0} increases) in the bifurcation diagram; hence, it intersects the curve of equilibrium points at different sites. When *x*_{0} = −1.6, then the *z*-nullcline intersects the *Z*-middle branch, which consists of saddles (Fig. 4). Consistent with Figure 2*A*, the equilibrium point is a saddle. The trajectory behavior corresponds to SLEs, which is illustrated in Figure 4 (right). When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation SNs_{1}. Then the trajectory switches to the *Z*-upper branch, which consists of unstable foci surrounded by stable periodic orbits. Hence the trajectory exhibits an oscillatory solution on the *Z*-upper branch, which approaches the *Z*-middle branch as *z* increases. The oscillatory solution is homoclinic to one of the saddles along the *Z*-middle branch, then it is destroyed as *z* is increased. The trajectory terminates in a homoclinic bifurcation HB and switches to the *Z*-lower branch. *Z*-lower and *Z*-upper branches correspond to normal and ictal states, respectively. The *Z*-middle branch acts as a separatrix between ictal and normal states. The transition from normal to ictal states (first bifurcation) occurs through a saddle-node bifurcation, SN_{1}. The transition from ictal to normal states (second bifurcation) occurs though a homoclinic bifurcation, HB. The homoclinic bifurcation *z*(HB) corresponds to an offset threshold and the saddle-node bifurcation
corresponds to an onset threshold. Therefore, the transitions between ictal and normal states occur through a fold/homoclinic bifurcation. The system is bistable on [SN_{1}, HB].

##### Bifurcation diagram when *m* decreases

To explore the role of *m* in the dynamics of the Epileptor, we plot a (*z*, *x*_{1}) bifurcation diagram, when *m* decreases. Let *m * = * *0 (Fig. 5), the plot first shows a *Z*-shaped curve. *Z*-Lower, *Z*-middle, and *Z*-upper branches consist of stable nodes, saddles, and unstable foci, respectively.

The *Z*-upper branch is surrounded by stable periodic orbits (see
and
curves), which terminate as *z* increases in a saddle-node on invariant circle (SNIC) bifurcation with
. Above the *Z*-curve, the lower branch consists of saddles and the upper branch comprises two sub-branches (Fig. 5, inset). The first sub-branch (solid) consists of stable foci and terminates as *z* decreases in a SNIC bifurcation. The second sub-branch (dash-dotted) consists of saddles. The two sub-branches collide as *z* increases in a saddle-node bifurcation SN_{2}. We plot a trajectory in Figure 5 (right), which corresponds to SLEs. When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation SN_{1} and the trajectory switches to the *Z*-upper branch exhibiting an oscillatory solution, which terminates as *z* increases in a SNIC bifurcation (see inset). Then the trajectory is at the upper branch above the *Z*-curve, which consists of stable foci. The trajectory exhibits a nonoscillatory solution, which terminates as *z* increases in a saddle-node bifurcation SN_{2} and the trajectory switches to the *Z*-lower branch. The *Z*-lower branch corresponds to the normal state. The *Z*-upper branch and the upper branch above the *Z*-curve correspond to the ictal state (see inset). The transition from normal to ictal states (first bifurcation) occurs through a saddle-node bifurcation, SN_{1}. The transition from ictal to normal states (second bifurcation) occurs through a saddle-node bifurcation, SN_{2}. Fast discharges characterizing the ictal state correspond to the oscillatory solution, which disappears through a saddle-node on invariant circle bifurcation SNIC. Therefore, the transitions between ictal and normal states are said to be of a fold/circle bifurcation. The system is bistable on [SN_{1}, SN_{2}].

We conclude that the bifurcation diagrams of the Epileptor model is a *Z*-shaped curve. Above, there are the following three branches: two branches consisting of saddles, and one branch consisting of stable foci for *m *≤* *0 and unstable foci for *m *>* *0. The Epileptor model undergoes a SNIC bifurcation for *m *≤* *0 and an HB bifurcation for *m *>* *0. The transitions between ictal and normal states occur through a fold/circle bifurcation for *m *≤* *0 and through a fold/homoclinic bifurcation for *m * >* * 0. The system is bistable on [SN_{1}, SN_{2}] for *m * ≤ * *0 and on [SN_{1}, HB] for *m *>* *0.

##### Bifurcation diagram for *m *=* *0 and

We now analyze the bifurcation diagram when
decreases. Let
, we plot a (*z*, *x*_{1}) bifurcation diagram for *m * = * *0 in Figure 6, which shows a *Z*-shaped curve. *Z*-lower, *Z*-middle, and *Z*-upper branches consist of stable nodes, saddles, and unstable foci, respectively. Decreasing *z*, *Z*-lower and *Z*-middle branches collide in a saddle-node bifurcation, SN_{1}. Increasing *z*, *Z*-upper and *Z*-middle branches collide in another saddle-node bifurcation, SN_{2}. The *Z*-upper branch terminates as *z* decreases in a saddle-node bifurcation SN_{4}. Moreover, two branches appear as *z* decreases. The lower branch (dashed) consists of unstable foci, and the upper branch (dash-dotted) consists of saddles. Increasing *z*, lower and upper branches collide in a saddle-node bifurcation SN_{3}. Decreasing *z*, the upper branch terminates in a saddle-node bifurcation SN_{4}.

Let *x*_{0} = −1.6, the *z*-nullcline is at the *Z*-middle branch, then the equilibrium point is a saddle. We plot a trajectory in Figure 6, which shows transitions between *Z*-lower and *Z*-upper branches. When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation, SN_{1}, and the trajectory switches to the *Z*-upper branch, which consists of unstable foci. An unstable focus is surrounded by a stable periodic orbit (see
and
curves). Hence, when the trajectory is at the *Z*-upper branch, it exhibits an oscillatory solution, which terminates as *z* increases in a homoclinic bifurcation HB. Then the trajectory switches to the *Z*-lower branch. *Z*-Lower (solid) and *Z*-upper (dashed) branches correspond to normal and ictal states, respectively. The transition from normal to ictal states occurs through a saddle-node bifurcation SN_{1}. The transition from ictal to normal states occurs through an HB. Thus, the transitions between ictal and normal states occur through a fold/homoclinic bifurcation. The system is bistable on [SN_{1}, HB].

##### Bifurcation diagram when *m* decreases, for

Let *m* = −0.5, we plot a (*z*, *x*_{1}) bifurcation diagram in Figure 7. The *Z*-upper branch is divided into two sub-branches separated by a Hopf (H) bifurcation (Fig. 7). The first sub-branch (solid) consists of stable foci and terminates as *z* increases in a saddle-node bifurcation SN_{2}. The second sub-branch (dashed) consists of unstable foci and terminates as *z* decreases in a saddle-node bifurcation SN_{4}. An unstable focus is surrounded by a stable periodic orbit (see
and
curves), which terminates as *z* increases in an H bifurcation. When *x*_{0} = −1.6, then the *z*-nullcline is at the *Z*-middle branch, hence the equilibrium point is a saddle. We plot a trajectory in Figure 7, which shows transitions between *Z*-lower and *Z*-upper branches. When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation, SN_{1}, and the trajectory switches to the *Z*-upper sub-branch (dashed) exhibiting an oscillatory solution, which terminates as *z* increases in an H bifurcation. Then, the trajectory is at the *Z*-upper (solid) sub-branch exhibiting a nonoscillatory solution, which terminates as *z* increases in a saddle-node bifurcation SN_{2}. *Z*-lower (solid) and *Z*-upper branches correspond to normal and ictal states, respectively. The oscillatory solution corresponds to fast discharges. The transition from normal to ictal states occurs through a saddle-node bifurcation, SN_{1}. The transition from ictal to normal states occurs through a saddle-node bifurcation SN_{2}. The oscillatory solution (fast discharges) terminates in an H bifurcation. The transitions between ictal and normal states are said to be of a fold/Hopf bifurcation. The system is bistable on [SN_{1}, SN_{2}].

Let *m* = −1, we plot a (*z*, *x*_{1}) bifurcation diagram in Figure 8, which has a *Z*-shaped curve. The Hopf bifurcation point H disappears, then the *Z*-upper branch consists of stable foci only, which is the equilibrium point of nonoscillatory state. Increasing *z*, the *Z*-upper (solid) branch terminates in a SNIC bifurcation (Fig. 8). Below the *Z*-upper branch, two branches appear as *z* decreases: the upper (dash-dotted) branch consists of saddles, and the lower (dashed) one consists of unstable foci. An unstable focus is surrounded by a stable periodic orbit (see
and
curves), which terminates as *z* increases in a SNIC bifurcation. Let *x*_{0} = −1.6, the *z*-nullcline is at the *Z*-middle branch, then the equilibrium point is a saddle. We plot a trajectory in Figure 8, which shows transitions between *Z*-lower and *Z*-upper branches. When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation, SN_{1}, and the trajectory switches to the *Z*-upper (solid) branch, which consists of stable foci, thereby exhibiting a nonoscillatory solution. Then, the stable focus disappears as *z* increases through a saddle-node bifurcation, SN_{2}, and the trajectory switches to the *Z*-lower branch. The transition from *Z*-lower to *Z*-upper branches occurs through a saddle-node bifurcation, SN_{1}, and the transition from *Z*-upper to *Z*-lower branches occurs through a saddle-node bifurcation, SN_{2}. The transitions between *Z*-lower and *Z*-upper branches occur through a fold/fold bifurcation. The system is bistable on [SN_{1}, SN_{2}]. Here, the SLE attractor reduces to a periodic switch between a nonoscillatory state and a NS.

We conclude that when
, the transitions between *Z*-lower and *Z*-upper branches occur through a fold/homoclinic bifurcation,
; a fold/Hopf bifurcation for *m* = −0.5; and a fold/fold bifurcation,
. The system is bistable on [SN_{1}, SN_{2}] for *m *<* * 0 and on [SN_{1}, HB] for *m * ≥* * 0. Moreover, the transitions between ictal and normal states occur through a fold/circle bifurcation when
, and do not when
.

#### Finding LC: a stable limit cycle

##### LC in the Epileptor model

We found a large geometrical object in the (*y*_{1}, *ψ*, *z*) phase space below the SLE attractor shown in Figure 1. When using the *z* equation (Eq. 34), the Epileptor diverges with time (Fig. 9*C*). The initial conditions *i*_{2} are on the left of the separatrix (*S*), which means that *i*_{2} are below the separatrix (S) and the SLE attractor in the phase space. To explain this divergence, we plot the *z*-nullcline and the
-curve in a (*z*, *x*_{1}) bifurcation diagram for *m * = * *0.5 (Fig. 9*A*). The *z*-nullcline corresponds to
(34)

When *x*_{0} = −1.6, then the *z*-nullcline intersects the
-curve at one periodic orbit labeled as S (Fig. 9*A*). We consider two different initial conditions *i*_{1} and *i*_{2} such as
. Only the trajectory *i*_{1} is plotted in the bifurcation diagram, which exhibits transitions between *Z*-lower and *Z*-upper branches through a fold/homoclinic bifurcation. We plot a time series for the initial condition *i*_{2} in Figure 9*C*, which shows that the trajectory diverges with time.

To stabilize the final state that the Epileptor evolves toward, we modified the *z*-equation by introducing the following equation:
(35)

We plot the corresponding *z*-nullcline (
; Eq. 35) in a (
) bifurcation diagram for *m * = * *0.5 (Fig. 9*B*). When *x*_{0} = −1.6, then the *z*-nullcline intersects the
-curve at two periodic orbits: S and LC (Fig. 9*B*). We consider the initial conditions *i*_{1} and *i*_{2} such as
. The trajectory *i*_{1} exhibits transitions between *Z*-lower and *Z*-upper branches, which occur through a fold/homoclinic bifurcation. We plot time series for the initial condition *i*_{2} in Figure 9*D*, which shows a stable LC. We plot the trajectory as its transients toward LC in Figure 4 and the corresponding time series in Figure 9*D*. LC is characterized by a large amplitude and a fast-slow invariant cycle (Fig. 9*D*).

We consider Equation 35 and plot the *z*-nullcline when
in a (
) bifurcation diagram (Fig. 10). The *z*-nullcline intersects three branches of the (
) curve, as follows: lower (dash-dotted) and upper (dashed) branches above the *Z*-curve, and the *Z*-upper (dashed) branch. Then three equilibrium points coexist: two unstable foci and one saddle. The *z*-nullcline intersects the
-curve at two periodic orbits, which are LC and S. Then the trajectory only exhibits an oscillatory solution, which stabilizes at the point LC.

##### Evolution of periodic orbits

We plot the
- and
-curves to limit the oscillatory solutions (Figs. 5–8). Decreasing *z*, trajectories converge to LC and the amplitude of periodic orbits increases, leading to LC with a large amplitude. The intersection of the *z*-nullcline and the
-curve gives rise to S and LC, which are saddle and stable periodic orbits, respectively.

Let *x*_{0} = −1.6, two periodic orbits exist for *m * = * *0.5 (Fig. 4) and *m * = * *0 (Fig. 5). We show only one periodic orbit in Figures 4 and 5, which corresponds to S. The second periodic orbit is shown in Figure 9*B*, which corresponds to LC. LC and S are stable and saddle periodic orbits, respectively. Trajectories plotted in Figures 4 and 5 exhibit SLEs (right) or converge to LC (left), depending on the initial conditions. SLEs and LC are considered as two attractors. The saddle periodic orbit (S) acts as a separatrix between the basin of attraction of SLEs and the basin of attraction of LC.

The *z*-nullcline moves downward (*x*_{0} decreases) or upward (*x*_{0} increases) in the bifurcation diagram, hence the existence of periodic orbits depends on *x*_{0}. We find S and LC in a (*x*_{0}, *z*^{*}) bifurcation diagram of periodic orbits (Fig. 11*A*). *x*_{0} acts as a control parameter. Finding S helps us to limit the basin of attraction of SLEs and the basin of attraction of LC. Finding LC helps us to determine the *z*-value at which LC stabilizes. Red (+) markers and black dots correspond to LC and S, respectively. Only LC exists for large *x*_{0}. Decreasing *x*_{0}, LC and S coexist, and collide as *x*_{0} is further decreased in a saddle-node of periodic orbits bifurcation (SNPO), then fades.

Figure 11*B* (*m* = −16) shows that for some *x*_{0}, a third periodic orbit appears, which is labeled as blue top (+) markers. This periodic orbit is stable, and is localized above LC and S (Fig. 11*B*). Since the amplitude of periodic orbits decreases as *z* increases, the third periodic orbit is a stable limit cycle with a small amplitude, which we denoted as small limit cycle (SLC). SLC and LC are two stable limit cycles that coexist for some *x*_{0}, and are separated by S. SLC exists after a SNPO bifurcation occurs when *m* = −16 and disappears when *x*_{0} is further decreased.

##### Parameter space of periodic orbits

We explore two (*m*, *x*_{0}) parameter spaces of periodic orbits: for
in Figure 12A and for
in Figure 12*B*. There are five areas separated by a bold line (SNPO bifurcation). LC exists above, but not below. LC exists in area I and coexists with S in area II. LC and S coexist with SLC in area III. Only SLC exists in area IV, which is localized below a SNPO bifurcation (bold line). Periodic orbits disappear in area V.

#### SLC dynamics

We used the averaging method to find periodic orbits, which lie at the intersection of the *z*-nullcline and the
-curve. When
, area III shows the coexistence of LC, SLC, and S (Fig. 12*B*).

The stability of the equilibrium points depends on *m* and *x*_{0} in area III. We analyze the coexistence of LC and SLC for each stability type of equilibrium points. We plot a (*z*, *x*_{1}) bifurcation diagram for *m *=* *1 in Figure 13*A*. When
, then the *z*-nullcline is at the middle branch and the equilibrium point is a saddle. The *z*-nullcline intersects the
-curve at three periodic orbits: LC (left), S (middle), and SLC (right). LC and SLC are two stable limit cycles with large and small amplitudes, respectively. S acts as a separatrix between LC and SLC. When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation SN_{1} and the trajectory switches to the *Z*-upper branch, which consists of unstable foci. Then the trajectory exhibits an oscillatory solution on the *Z*-upper branch and continues to SLC at which *z* stabilizes (Fig. 13*A*, inset). The transition to the *Z*-lower branch does not occur even if the equilibrium point is a saddle. In fact, *z* stabilizes at SLC before the homoclinic bifurcation HB. The trajectory behavior corresponds to a periodic solution.

We plot a (*z*, *x*_{1}) bifurcation diagram for *m* = −2 in Figure 13*B*. Let *x*_{0} = 0, the equilibrium point is an unstable focus. The *z*-nullcline intersects the
-curve at three periodic orbits (Fig. 13*B*): LC (left), S (middle), and SLC (right). S acts as a separatrix between LC and SLC. When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation SN_{1}, and the trajectory switches to the *Z*-upper (solid) branch, which consists of stable foci. Increasing *z*, the trajectory exhibits a nonoscillatory solution on the *Z*-upper branch (solid), which terminates in a SNIC bifurcation (Fig. 13*B*, inset). Then, the trajectory is at the dashed branch exhibiting an oscillatory solution and continues to SLC at which *z* stabilizes. The transition to the *Z*-lower branch does not occur and the trajectory behavior corresponds to a periodic solution.

We conclude that when
, then LC and SLC coexist in area III and S acts as a separatrix between them. The equilibrium point of the Epileptor model is either a saddle or an unstable focus. Time series of the Epileptor system *ψ*, subsystem 1 *ψ*_{1}, and subsystem 2 *ψ*_{2} are plotted in Figure 13*C* when the equilibrium point is a saddle, and in Figure 13*D* when it is an unstable focus. The SLC patterns depend on the stability of the equilibrium point. The subsystem 2 generates a resting state when the equilibrium point is a saddle (see *ψ*_{2}; Fig. 13*C*) and a spiking state when it is an unstable focus (see *ψ*_{2}; Fig. 13*D*).

#### Transition from epileptiform fast discharges to a DB

The Epileptor exhibits epileptiform fast discharges in the ictal state, which corresponds to an oscillatory solution at the *Z*-upper branch of the (*z*, *x*_{1}) curve. Let
, the *Z*-upper branch consists of unstable foci when *m * = * *0 (Fig. 6), and consists of stable foci when *m* = −1 (Fig. 8). Recall that the SLE reduces to a periodic switch between a nonoscillatory state and NS when *m* = −1.

Further decreasing *m*, the imaginary part of the complex–conjugate eigenvalues goes to zero, and then the *Z*-upper branch only consists of stable nodes. We plot a (*z*, *x*_{1}) bifurcation diagram for *m* = −8 (
) in Figure 14*A*, which shows a *Z*-shaped curve. *Z*-lower (solid), *Z*-middle (dash-dotted), and *Z*-upper (solid) branches consist of stable nodes, saddles, and stable nodes, respectively. The *Z*-upper branch terminates as *z* decreases in a SNIC bifurcation (Fig. 14*A*, inset). Below the *Z*-upper branch (Fig. 14*A*, inset), lower (dashed) and upper (dash-dotted) branches appear, which approach as *m* decreases to the *Z*-upper (solid) branch. The lower branch consists of unstable foci, and the upper branch consists of saddles. Increasing *z*, the lower and upper branches collide in a saddle-node bifurcation, SN_{3}. Note that
. Let *x*_{0} = −1.6, the equilibrium point is a saddle. When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation, SN_{1}, and the trajectory switches to the *Z*-upper branch, which consists of stable nodes. The trajectory exhibits a silent activity on the *Z*-upper branch, which terminates as *z* increases in a saddle-node bifurcation SN_{2}. Then, the trajectory switches to the *Z*-lower branch. The transition between *Z*-lower and *Z*-upper branches occur through a fold/fold bifurcation. The system is bistable on [SN_{1}, SN_{2}]. Time series of the Epileptor model *ψ*, subsystem 1 *ψ*_{1}, and subsystem 2 *ψ*_{2} are plotted in Figure 14*B*. The fast discharges during the ictal period disappear with decreasing *m* underlying a slow wave. This scenario corresponds to a depolarization block (or excitation block; Izhikevich, 2007), which is indicated by segment number 1 (see *ψ*; Fig. 14*B*). The NS is indicated by segment number 2. Therefore, the transitions between ictal states and NSs reduce to a periodic transition between DB and NS, which are both slow manifolds.

The DB occurs when
but not when
. To see this, we plot a (*z*, *x*_{1}) bifurcation diagram for *m* = −8 (
) in Figure 15*A*, which is with a *Z*-shaped curve. *Z*-lower (solid), *Z*-middle (dash-dotted), and *Z*-upper (dashed) branches consist of stable nodes, saddles, and unstable foci, respectively. Above the *Z*-curve, lower branch consists of saddles and upper branch consists of stable nodes (Fig. 15*A*, inset II). Decreasing *z*, lower and upper branches collide in a SNIC bifurcation. Let *x*_{0} = −1.6, the equilibrium point is a saddle. When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation SN_{1} and the trajectory switches to the *Z*-upper branch exhibiting an oscillatory solution (Fig. 15*A*, inset I), which terminates as *z* increases in a SNIC bifurcation (Fig. 15*A*, inset II). Then, the trajectory exhibits a nonoscillatory solution (Fig. 15*A*, inset II), which terminates as *z* increases in a saddle-node bifurcation, SN_{2} (Fig. 15*A*), and the Epileptor switches to the *Z*-lower branch. The transitions between *Z*-lower and *Z*-upper branches occur through a fold/circle bifurcation. The system is bistable on [SN_{1}, SN_{2}]. Time series of the Epileptor model *ψ*, subsystem 1 *ψ*_{1}, and subsystem 2 *ψ*_{2} are plotted in Figure 15*B*. The Epileptor exhibits fast discharges in the ictal state (see *ψ*), and hence does not enter into a DB.

To understand why a DB occurs for
and not for
, we analyze the [SN_{1}, SN_{2}] interval, on which the Epileptor model is bistable. We determine the branches of the (*z*, *x*_{1}) curve existing on this interval (
) for
and
. Distinct branches coexist on [SN_{1}, SN_{2}] depending on *m*:

When

*m*= 0, then three branches coexist on [SN_{1}, SN_{2}], consisting of stable nodes, saddles, and unstable foci for (Fig. 5), and they consist of saddles and unstable foci for (Fig. 6).When

*m*= −1, then three branches coexist on [SN_{1}, SN_{2}], which consist of stable foci, saddles, and unstable foci for (figure not shown) and (Fig. 8).When

*m*= −8, then three branches coexist for , which consist of stable nodes, saddles, and unstable foci (Fig. 15*A*). For , one branch exists which consists of stable nodes (Fig. 14*A*).

Therefore, a DB occurs when *m* = −8 and
because one branch exists on [SN_{1}, SN_{2}] (
), which consists of stable nodes. The SLE attractor is characterized by transitions between *Z*-lower and *Z*-upper branches. The *Z*-lower branch consists of stable nodes, which are the equilibrium points of the NS. When *m* = −8 and
, the *Z*-upper branch consists of stable nodes, which are the equilibrium points of the DB. Hence, the SLE attractor reduces as *m* further decreases to a periodic switch between DB and NS.

#### Stabilizing equilibrium points in the Epileptor model

The transitions between ictal and normal states of the SLE attractor occur when the equilibrium point is a saddle (*x*_{0} = −1.6). When *x*_{0} increases, the Epileptor stabilizes in the ictal (nonoscillatory) state, and when *x*_{0} decreases, the Epileptor model stabilizes in the normal state.

##### Stabilizing the equilibrium point of the nonoscillatory state

During the ictal state, epileptiform fast discharges (oscillatory state) disappear through three bifurcation types depending on *m* and
: a homoclinic bifurcation (Figs. 4, 6), a SNIC bifurcation (Fig. 5), or a Hopf bifurcation (Fig. 7). When the equilibrium point is a saddle (*x*_{0} = −1.6), the Epileptor switches differently to the normal state according to the bifurcation type: homoclinic, SNIC, or Hopf. The Epileptor switches to the normal state just after a homoclinic bifurcation, and remains in the nonoscillatory state after SNIC and Hopf bifurcations, thereby exhibiting a nonoscillatory solution. The latter disappears as *z* increases through a saddle-node bifurcation, SN_{2}, and then the Epileptor switches to the normal state. When increasing *x*_{0}, the stability of the equilibrium point changes, and the Epileptor remains in the nonoscillatory state. We illustrate this scenario in Figure 16*A–C*.

Let *m * = * *0 and *x*_{0} = −0.9 (Fig. 16*A*), the *z*-nullcline intersects three branches of the (*z*, *x*_{1}) curve: lower and upper (solid) branches above the *Z*-curve, and the *Z*-upper (dashed) branch. The *z*-nullcline intersects the upper branch above *Z* at *C*. Then three equilibrium points coexist: an unstable focus, a saddle, and a stable focus (C). When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation, SN_{1}, and the trajectory switches to the *Z*-upper branch exhibiting an oscillatory solution (fast discharges), which terminates as *z* increases in a SNIC bifurcation. Then, the trajectory exhibits a nonoscillatory solution at the upper (solid) branch above the *Z*-curve, and continues to *C* at which *z* stabilizes. Therefore, the transition to the *Z*-lower branch does not occur and the Epileptor remains in the nonoscillatory state. We plot this scenario in a (*Y*, *X*, *z*) phase space (see Fig. 28*A*, top trajectory). Corresponding time series are plotted in Figure 28*B*. We add numbers to identify trajectory segments. We indicate the normal state by (1), fast discharges by (2), and the final state by (3), which corresponds to the equilibrium point of the nonoscillatory state. Here, the NS exists but not the equilibrium point. Hence, after a transient normal state, the Epileptor exhibits fast discharges which disappear in a SNIC bifurcation, and then remains in the nonoscillatory state.

Let *m* = −0.5 and *x*_{0} = −0.8 (Fig. 16*B*), the *z*-nullcline intersects the *Z*-upper (solid) sub-branch at *C* which is a stable focus. When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation SN_{1} and the trajectory switches to the *Z*-upper (dashed) sub-branch exhibiting an oscillatory solution (fast discharges), which terminates as *z* increases in a Hopf bifurcation, *H*. Then, the trajectory exhibits a nonoscillatory solution at the *Z*-upper (solid) sub-branch, which consists of stable foci, and continues to *C* at which *z* stabilizes. The transition to the *Z*-lower branch does not occur, and the Epileptor remains in the nonoscillatory state. We plot this scenario in a (*Y*, *X*, *z*) phase space (see Fig. 29A, top trajectory). Corresponding time series are plotted in Figure 29*B*. We indicate the normal state by (1), fast discharges by (2), and the final state by (3), which corresponds to the equilibrium point of the nonoscillatory state. After a transient normal state, the Epileptor exhibits fast discharges which disappear in a Hopf bifurcation, and then remains in the nonoscillatory state.

Let *m* = −1, the SLE attractor reduces to transitions between nonoscillatory and normal states when *x*_{0} = −1.6, which occur through a fold/fold bifurcation (Fig. 8). When increasing *x*_{0}, the Epileptor model stabilizes its stable focus, which corresponds to the equilibrium point of the nonoscillatory state. We illustrate this case in Figure 16*C*. Let *x*_{0} = −0.8, the *z*-nullcline intersects the *Z*-upper branch at *C*, which is a stable focus. When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation, SN_{1}, and the trajectory switches to the *Z*-upper branch exhibiting a nonoscillatory solution. Then, the trajectory continues to *C* at which *z* stabilizes. The transition to the *Z*-lower branch does not occur and the Epileptor remains in the nonoscillatory state. We plot this scenario in a (*Y*, *X*, *z*) phase space (see Fig. 30*A*, top trajectory). Corresponding time series are plotted in Figure 30*B*. We indicate the normal state by (1) and the final state by (3), which corresponds to the equilibrium point of the nonoscillatory state. The imaginary part of the eigenvalues of *C* is responsible for the oscillations (2) around the stable focus. After a transient normal state, the Epileptor spirals into the equilibrium point of the nonoscillatory state exhibiting damped oscillations, and does not re-enter the normal state.

##### Stabilizing the equilibrium point of the DB

Figure 14 shows that the SLE attractor reduces to a periodic switch between DB and NS when further decreasing *m* (
). The equilibrium point is a saddle (*x*_{0} = −1.6). When increasing *x*_{0}, the Epileptor model stabilizes its stable node, which corresponds to the equilibrium point of DB. To see this, we plot a (*z*, *x*_{1}) bifurcation diagram when *x*_{0} = −0.8 in Figure 16*D*. The *z*-nullcline intersects the *Z*-upper (solid) branch at *C*, which is a stable node. When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation SN_{1} and the trajectory switches to the *Z*-upper branch and enters into a DB. The trajectory continues as *z* increases to *C* at which *z* stabilizes. Thus, the transition to the *Z*-lower branch does not occur and the Epileptor remains in DB. We plot this scenario in a (*Y*, *X*, *z*) phase space (see Fig. 27*A*, top trajectory). Corresponding time series are plotted in Figure 27*B*. We indicate the normal state by (2) and the final state by (4), which corresponds to the equilibrium point of DB. Then, after a transient NS, the Epileptor enters into DB and then remains in it.

##### Stabilizing the equilibrium point of the NS

When the equilibrium point is a saddle (*x*_{0} = −1.6), the normal state of the SLE attractor disappears through a saddle-node bifurcation SN_{1}. When decreasing *x*_{0}, the Epileptor stabilizes its stable node, which is the equilibrium point of the normal state. To see this, we plot a (*z*, *x*_{1}) bifurcation diagram when *x*_{0} = −2.5 in Figure 17. The *z*-nullcline intersects the *Z*-lower branch at *C*, which corresponds to the equilibrium point of the normal state. The ictal state (*Z*-upper branch) exists but not the equilibrium point, hence the Epileptor remains in the normal state, at *C* which is a stable node.

When
, the transition from ictal to normal states occurs according to two scenarios, depending on *m*.

*m*= 0.5 (Fig. 17*A*): When a trajectory is at the*Z*-upper branch (dashed), the oscillatory solution terminates as*z*increases in an HB bifurcation and the trajectory switches to the*Z*-lower branch.*m*= 0 (Fig. 17*B*): When a trajectory is at the*Z*-upper branch (dashed), the oscillatory solution terminates as*z*increases in a SNIC bifurcation, then the trajectory exhibits a nonoscillatory solution on the upper branch above the*Z*-curve, which terminates as*z*increases in a saddle-node bifurcation SN_{2}. Then the trajectory switches to the*Z*-lower branch.

When
, the transition from ictal to normal states occurs according to four scenarios, depending on *m*.

For

*m*≥ 0 (Fig. 17*C*): The trajectory exhibits an oscillatory solution on the*Z*-upper (dashed) branch, which terminates as*z*increases in a homoclinic bifurcation HB, and then it switches to the*Z*-lower branch.For

*m*= −0.5 (Fig. 17*D*): The trajectory exhibits an oscillatory solution on the*Z*-upper sub-(dashed) branch, which terminates as*z*increases in a Hopf bifurcation*H*. Then the trajectory exhibits a nonoscillatory solution on the*Z*-upper (solid) sub-branch, which terminates as*z*increases in a saddle-node bifurcation SN_{2}, and the trajectory switches to the*Z*-lower branch.For

*m*≤ −1, when a trajectory starts at*i*with (third scenario), then it exhibits a nonoscillatory solution on the*Z*-upper branch, which terminates as*z*increases in a saddle-node bifurcation, SN_{2}. Then the trajectory switches to the*Z*-lower branch (Fig. 17*E*). When a trajectory starts at*i*with (fourth scenario), then it exhibits an oscillatory solution, which terminates as*z*increases in a SNIC bifurcation. Then, the trajectory is at the*Z*-upper branch exhibiting a nonoscillatory solution, which terminates as*z*increases in a saddle-node bifurcation SN_{2}. Then, the trajectory switches to the*Z*-lower branch (Fig. 17*F*).

For all *m*, when the trajectory is at the *Z*-lower branch, which consists of stable nodes, it continues as *z* decreases to *C*, at which *z* stabilizes. Therefore, the transition to the *Z*-upper branch does not occur and the Epileptor remains in the NS.

#### Coexistence in the Epileptor model

Figure 11, *A* and *B*, shows the evolution of periodic orbits as *x*_{0} varies. The periodic orbits LC and S disappear as *x*_{0} decreases through a SNPO bifurcation. The dynamics of SLC when it exists, depends on the equilibrium points stability. Figure 4 shows that the *z*-nullcline and the
-curve intersect at two points: LC and S. Depending on initial conditions, trajectories either converge to LC or exhibit SLEs. We identify seven types of coexisting attractors in the Epileptor model.

##### Coexistence of LC and seizures (SLEs)

Depending on *m* and
, transitions between ictal and normal states, which characterize the SLE attractor, occur through a fold/homoclinic bifurcation, a fold/circle bifurcation, or a fold/Hopf bifurcation. We describe for each bifurcation type the coexistence of LC and the SLE attractor:

*A fold/homoclinic bifurcation:* We plot two trajectories in a (*z*, *x*_{1}) bifurcation diagram for *m * = * *0.5 (Fig. 4). The equilibrium point is a saddle (*x*_{0} = −1.6). The behavior of the first trajectory (right) corresponds to SLEs, which occur through a fold/homoclinic bifurcation. The *Z*-middle branch acts as a separatrix between ictal and normal states. The second trajectory (left) converges to LC. The saddle periodic orbit (S) separates LC and the SLE attractor. We plot LC, SLEs, and S in a (*Y*, *X*, *z*) phase space (Fig. 18*A*). SLEs and LC coexist and the separatrix between them corresponds to the saddle periodic orbit (S). SLEs lock in LC when *z* declines to below baseline shift (Fig. 18*A*). Time series are plotted in Figure 18*B* for SLEs and in Figure 18*C* for LC.

Figure 6 shows that SLEs occur through a fold/homoclinic bifurcation when
. Parameter settings in Figure 6 correspond to area II in Figure 12*B*. Hence, SLEs and LC coexist and are separated by S. Only the SLE attractor is plotted in Figure 6. We plot LC, SLEs, and S in a (*Y*, *X*, *z*) phase space (Fig. 19*A*). SLEs and LC coexist and are separated by a separatrix, which corresponds to the saddle periodic orbit (S). Time series are plotted in Figure 19*B* for SLEs and in Figure 19*C* for LC.

*A fold/circle bifurcation:* We plot two trajectories in a (*z*, *x*_{1}) bifurcation diagram (Fig. 5). The first trajectory (right) exhibits SLEs, and the second one (left) converges to LC. The equilibrium point is a saddle (*x*_{0} = −1.6). The transitions between ictal and normal states occur through a fold/circle bifurcation. Then SLEs and LC coexist and are separated by a saddle periodic orbit (S). We plot LC, SLEs, and S in a (*Y*, *X*, *z*) phase space (Fig. 20*A*). SLEs and LC coexist and are separated by a saddle periodic orbit (S), which is the separatrix. Time series are plotted in Figure 20*B* for SLEs and in Figure 20*C* for LC.

*A fold/Hopf bifurcation:* Figure 7 shows that SLEs occur through a fold/Hopf bifurcation when
(
). The equilibrium point is a saddle (*x*_{0} = −1.6). Parameter settings in Figure 7 correspond to area II in Figure 12*B*. Then SLEs and LC coexist and are separated by S. Figure 7 only shows the SLE attractor. We plot LC, SLEs, and S in a (*Y*, *X*, *z*) phase space (Fig. 21*A*). SLEs and LC coexist and the separatrix between them corresponds to a saddle periodic orbit (S). Time series are plotted in Figure 21*B* for SLEs and in Figure 21*C* for LC.

##### Coexistence of LC and a periodic switch between nonoscillatory state and NS

Figure 8 shows that the SLE attractor reduces to a periodic switch between nonoscillatory state and NS, which occurs through a fold/fold bifurcation when *m* = −1 (
). The equilibrium point is a saddle (*x*_{0} = −1.6). Parameter settings in Figure 8 correspond to area II in Figure 12*B*. Then LC and the periodic switch between nonoscillatory state and NS coexist, and are separated by S. Only the periodic switch between nonoscillatory state and NS is plotted in Figure 8. We plot LC, this periodic switch, and S in a (*Y*, *X*, *z*) phase space (Fig. 22*A*). LC and the periodic switch between nonoscillatory state and NS coexist and the separatrix between them corresponds to a saddle periodic orbit (S). Time series are plotted in Figure 22*B* for this periodic switch and in Figure 22*C* for LC.

##### Coexistence of two periodic orbits: LC and SLC

SLC is a stable periodic orbit with a small amplitude. When
, then SLC, S, and LC coexist in area III (Fig. 12*B*). Figure 13 shows that the *z*-nullcline intersects the
-curve at three points: LC, S, and SLC. Here, we plot the SLC behavior only, which depends on the stability of the equilibrium point. SLC and LC coexist in both cases: when the equilibrium point is a saddle (Fig. 13*A*) and when it is an unstable focus (Fig. 13*B*). We plot SLC, LC, and S in a (*Y*, *X*, *Z*) phase space (Figs. 23, 24). LC and SLC coexist and are separated by a saddle periodic orbit (S), which corresponds to the separatrix. When the equilibrium point is a saddle, the coexistence of LC and SLC is illustrated in Figure 23*A*. Time series are plotted in Figure 23*B* for SLC and in Figure 23*C* for LC. When the equilibrium point is an unstable focus, the coexistence of LC and SLC is illustrated in Figure 24*A*. Time series are plotted in Figure 24*B* for SLC and in Figure 24*C* for LC.

##### Coexistence of LC and a periodic switch between DB and NS

A (*z*, *x*_{1}) bifurcation diagram is plotted in Figure 14*A* for *m* = −8 (
). The trajectory characterizes the periodic switch between DB and NS. Parameter settings in Figure 14 correspond to area II in Figure 12*B*. Thus, LC and the periodic switch between DB and NS coexist, and are separated by a saddle periodic orbit (S). We plot LC, S, and the periodic switch between DB and NS in a (*Y*, *X*, *z*) phase space (Fig. 25*A*). We observe that DB occurs when the fast manifold of the SLE collapses to a point that traces out a line under the slow *z*-evolution. The ictal state thus reduces to a silent activity. LC and the periodic switch between DB and NS coexist in the phase space and the separatrix acts as a barrier between them. Time series are plotted in Figure 25*B* for the periodic switch between DB and NS, and in Figure 25*C* for LC. DB and NS are indicated by segment numbers 1 and 3, respectively.

##### Coexistence of LC and a NS

Figure 17 shows that when decreasing *x*_{0}, the *z*-nullcline is at the *Z*-lower branch, which consists of equilibrium points of the normal state (stable nodes). Let *m * = * *0, the equilibrium point is a stable node when *x*_{0} is below −2.1 (Fig. 3, area 5). LC and S coexist when *x*_{0} = −2.1 (Fig. 12*A*, area II). We plot LC, NS, and S in a (*Y*, *X*, *z*) phase space (Fig. 26*A*). LC and NS coexist in the phase space, and the separatrix acts as a barrier between them. Fast discharges and NS are indicated by segment numbers 1 and 4, respectively. After the transient seizure-like fast discharges, the Epileptor remains in the NS. Time series are plotted in Figure 26*B* for NS and in Figure 26*C* for LC.

##### Coexistence of LC and a DB

Figure 16*D* shows that when increasing *x*_{0}, the *z*-nullcline is at the *Z*-upper branch, which consists of equilibrium points of DB (stable nodes). Parameter settings in Figure 16*D* correspond to area II in Figure 12*B*, where LC and S coexist. We plot LC, DB, and S in a (*Y*, *X*, *z*) phase space (Fig. 27*A*). The transient NS and DB are indicated by segment numbers 2 and 4, respectively. After the transient normal state, the Epileptor model remains in the depolarization block. Time series are plotted in Figure 27*B* for DB and in Figure 27*C* for LC.

##### Coexistence of LC and a nonoscillatory state

Figure 16*A–C* shows that when increasing *x*_{0}, the equilibrium point C is at the *Z*-upper branch when
and at the upper branch above the *Z*-curve when
. *C* (stable focus) is the equilibrium point of the nonoscillatory state. Figure 16, *A* and *B*, shows that after the transient normal state, and the transient seizure-like fast discharges, the Epileptor remains in the nonoscillatory state. Figure 16*C* shows that after the transient normal state, the Epileptor spirals into the equilibrium point of the nonoscillatory state exhibiting damped oscillations and then stabilizes at C. Parameter settings in Figure 16*A* correspond to area II in Figure 12*A*, and parameter settings in Figure 16, *B* and *C*, correspond to area II in Figure 12*B*, and then LC and S coexist. Only the nonoscillatory state is plotted in Figure 16*A–C*. We plot LC, the nonoscillatory state, and S in a (*Y*, *X*, *z*) phase space (Figs. 28
*A*, 29
*A*, 30
*A*). LC (below) and the nonoscillatory state (above) coexist. The nonoscillatory state shown in Figures 28
*A*, 29
*A*, and 30
*A* is the same as that shown in Figure 16*A–C*, respectively. NS and the nonoscillatory state are indicated by segment numbers 1 and 3, respectively. Time series for each case are plotted in Figures 28*B*, 29*B*, and 30*B* for the nonoscillatory states, and in Figures 28*C*, 29*C*, and 30*C* for LC.

#### Parameter space of equilibrium points and periodic orbits

We identified seven types of coexisting attractors: a coexistence of SLE and LC (Figs. 18–21), a coexistence of LC and periodic switch between nonoscillatory state and NS (Fig. 22), a coexistence of SLC and LC (Figs. 23, 24), a coexistence of LC and periodic switch between DB and NS (Fig. 25), a coexistence of NS and LC (Fig. 26), a coexistence of DB and LC (Fig. 27), and a coexistence of nonoscillatory state and LC (Figs. 28–30). SLEs are characterized by transitions between ictal and normal states. Depending on *m*, *x*_{0} and
, these transitions occur through a fold/homoclinic bifurcation (Figs. 18, 19), a fold/circle bifurcation (Fig. 20), or a fold/Hopf bifurcation (Fig. 21).

To characterize the different types of the coexisting attractors, we explore a (*m*, *x*_{0}) parameter space of periodic orbits and equilibrium points in Figures 31 and 32 for
.

##### Parameter space of the Epileptor model for

The parameter space is divided into two parts which are separated by a bold line (SNPO bifurcation; Fig. 31). Above the bold line, LC exists but does not exist below it. Twelve areas exist. LC exists in area I and coexists with SLC in area II. Only SLC exists in area III. The nonoscillatory state exists in area IV and coexists with LC in area V. Only SLE exists in area VI and coexists with LC in areas VII and VIII. An SLE occurs through a fold/circle bifurcation in areas VI and VII. An SLE occurs through a fold/homoclinic bifurcation in area VIII. NS exists in area IX and coexists with LC in area X. LC coexists with a chaotic state in areas XI and XII.

For each area, the Epileptor model has different equilibrium points, depending on *m* and *x*_{0}. To explore this further, we determine the equilibrium points and the periodic orbits of each area (I-XII) by using 32, which shows 18 (1–18) areas:

*(I) LC:* Area I in Figure 31 is composed of areas 1, 7, 8, 9, 10, and 14 in Figure 32. The *z*-nullcline intersects the
-curve at one periodic orbit, which corresponds to LC. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at different equilibrium points. The equilibrium point is an unstable focus in area 1 and a saddle in area 14. The equilibrium points are one saddle and two unstable foci in area 7. The equilibrium points are one unstable focus and two saddles in area 9. The equilibrium points are three saddles in area 10. The equilibrium points are a saddle, an unstable focus, and an unstable node in area 8. Trajectories exhibit only a fast-slow cyclic behavior (LC), which is plotted in Figures 9*D* and 10.

*(II) LC and SLC:* Area II in Figure 31 corresponds to area 2 in Figure 32. The *z*-nullcline intersects the
-curve at three periodic orbits in area 2, which correspond to LC, S, and SLC, respectively. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at one equilibrium point, which is an unstable focus. LC and SLC coexist. The stable manifold of the saddle periodic orbit (S) separates the basin of attraction of LC and the basin of attraction of SLC. Trajectories exhibit two periodic solutions depending on initial conditions. The first solution corresponds to a fast-slow cyclic behavior (LC). The second solution corresponds to a periodic solution with a small amplitude (SLC).

*(III) SLC:* Area III in Figure 31 corresponds to area 3 in Figure 32. The *z*-nullcline intersects the
-curve at one periodic orbit in area 3, which corresponds to SLC. Area three is below the SNPO bifurcation (bold line), which means that SLC exists after a SNPO bifurcation. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at one equilibrium point, which is an unstable focus. The trajectories exhibit only a periodic solution with a small amplitude (SLC).

*(IV) Nonoscillatory state:* Area IV in Figure 31 corresponds to area 4 in Figure 32. Here, the *z*-nullcline and the
-curve do not intersect, hence periodic orbits LC, S, and SLC do not exist. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at different equilibrium points. The equilibrium points are a stable node, an unstable focus, and a saddle in area 4. Trajectories remain in the nonoscillatory state.

*(V) LC and a nonoscillatory state:* Area V in Figure 31 is composed of areas 5 and 6 in Figure 32. The *z*-nullcline intersects the
-curve at two periodic orbits in areas 5 and 6, which correspond to LC and S. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at different equilibrium points. The equilibrium points are a stable node, an unstable focus, and a saddle in area 5. The equilibrium points are a stable focus, an unstable focus, and a saddle in area 6. Trajectories either exhibit a fast-slow cyclic behavior (LC) or remain in the nonoscillatory state. The coexistence of LC and the nonoscillatory state is plotted in Figure 28*A*. The equilibrium points belong to area 6 in Figure 32. Time series are plotted in Figures 28*B* for the nonoscillatory state and in Figure 28*C* for LC.

*(VI) SLE with a fold/circle bifurcation:* Area VI in Figure 31 corresponds to area 11 in Figure 32. Here, the *z*-nullcline and the
-curve do not intersect, hence periodic orbits LC, S, and SLC do not exist. The *z*-nullcline intersects the *Z*-middle branch, which consists of saddles. Thus, a unique saddle equilibrium point exists. Since *m *≤* * 0, an SLE occurs through a fold/circle bifurcation.

*(VII) LC and SLE with a fold/circle bifurcation:* Area VII in Figure 31 corresponds to area 12 in Figure 32. Here, the equilibrium point is a saddle, and an SLE occurs through a fold/circle bifurcation. The *z*-nullcline intersects the
-curve at two periodic orbits, which correspond to LC and S. Then LC and an SLE with a fold/circle bifurcation coexist and are separated by a saddle periodic orbit (S). This coexistence is plotted in Figure 20A. Time series are plotted in Figure 20*B* for SLE and in Figure 20*C* for LC.

*(VIII) LC and SLE with a fold/homoclinic bifurcation:* Area VIII in Figure 31 corresponds to area 13 in Figure 32. Here, the equilibrium point is a saddle, and an SLE occurs through a fold/homoclinic bifurcation. The *z*-nullcline intersects the
-curve at two periodic orbits, which correspond to LC and S. Then LC and SLE with a fold/homoclinic bifurcation coexist and are separated by a saddle periodic orbit (S). This coexistence is plotted in Figure 18*A*. Time series are plotted in Figure 18*B* for SLE and in Figure 18*C* for LC.

*(IX) NS:* Area IX in Figure 31 corresponds to area 17 in Figure 32. The *z*-nullcline and the
-curve do not intersect, hence periodic orbits LC, S, and SLC do not exist. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at the *Z*-lower branch, and then the equilibrium point is a stable node in area 17. The Epileptor remains in NS after a transient period, which is shown in Figure 17

*(X) LC and NS:* Area X in Figure 31 corresponds to area 18 in Figure 32. The *z*-nullcline intersects the
-curve at two periodic orbits, which correspond to LC and S. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at the *Z*-lower branch, which is a stable node in area 18. Then LC and NS coexist and are separated by a saddle periodic orbit (S). This coexistence is plotted in Figure 26*A*. Time series are plotted in Figure 26*B* for NS and in Figure 26*C* for LC.

*(XI*, *XII) LC and a chaotic state:* Areas XI and XII in Figure 31 correspond to areas 15 and 16 in Figure 32. The *z*-nullcline intersects the
-curve at three periodic orbits in area 15. Then LC, S, and SLC coexist. S acts as a separatrix between the basin of attraction of LC and the basin of attraction of SLC. Then trajectories exhibit either LC or SLC in area 15. The *z*-nullcline intersects the
-curve at two periodic orbits in area 16, which correspond to LC and S. The *z*-nullcline intersects the *Z*-middle branch which consists of saddles. An SLE occurs through a fold/homoclinic bifurcation in area 16. Then LC and SLE coexist in area 16. S acts as a separatrix between the basin of attraction of LC and the basin of attraction of SLE in area 16.

We calculate the times between two successive spikes, which are referred to as the interspike intervals (ISIs) with respect to *x*_{0} in Figure 33*A* for *m * = * *1 and in Figure 33*B* for *m * = * *1.5. The system is chaotic when the ISI is irregular. Therefore, LC and the chaotic state coexist in areas 15 and 16 and are separated by a saddle periodic orbit (S).

The parameter space of the Epileptor model includes attractors as LC, SLE, and normal state, but neither DB nor SLE with a fold/Hopf bifurcation. In fact, using bifurcation analysis, we demonstrated that DB and Hopf bifurcation appeared when decreasing
. To characterize these attractors, we explore a (*m*, *x*_{0}) parameter space of periodic orbits and equilibrium points in Figures 34 and 35 for
.

##### Parameter space of the Epileptor model for

The parameter space is divided into two parts, which are separated by a bold line (SNPO bifurcation) (Fig. 34). Above the line, LC exists and does not exist below it. Thirteen areas exist. LC exists in area I and coexists with SLC in area II. Only SLC exists in area III. A DB exists in area IV and coexists with LC in area V. A periodic switch between a DB and a NS exists in area VII and coexists with LC in area VIII. A nonoscillatory state coexists with LC in area VI, and a periodic switch between a nonoscillatory state and a NS coexists with LC in area IX. LC and SLE coexist in areas X and XI. an SLE occurs through a fold/Hopf bifurcation in area X and through a fold/homoclinic bifurcation in area XI. A normal state exists in area XII and coexists with LC in area XIII.

For each area, the Epileptor model has different equilibrium points depending on *m* and *x*_{0}. To analyze this further, we determine the equilibrium points and the periodic orbits of each area (I-XIII) by using 35, which shows 20 (1–20) areas:

*(I) LC:* Area (I) in Figure 34 is composed of areas 1, 7, and 17 in Figure 35. The *z*-nullcline intersects the
-curve at one periodic orbit, which is LC. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at different equilibrium points. The equilibrium point is an unstable focus in area 1 and a saddle in area 17. The equilibrium points are one saddle and two unstable foci in area 7. Trajectories exhibit only a fast-slow cyclic behavior (LC).

*(II) LC and SLC:* Area (II) in Figure 34 is composed of areas 2 and 18 in Figure 35. The *z*-nullcline intersects the
-curve at three periodic orbits in areas 2 and 18, which correspond to LC, S, and SLC, respectively. Then LC and SLC coexist in areas 2 and 18. The stable manifold of a saddle periodic orbit (S) separates the basin of attraction of LC and the basin of attraction of SLC. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at one equilibrium point, which is an unstable focus in area 2 and a saddle in area 18. The trajectories exhibit two periodic solutions depending on the initial conditions. The first solution corresponds to a fast-slow cyclic behavior (LC). The second solution corresponds to a periodic solution with a small amplitude (SLC). Moreover, the SLC behavior depends on the stability of the equilibrium point. The SLC behavior is plotted in Figure 13*A* for area 18 and in Figure 13*B* for area 2. The coexistence of LC and SLC is plotted in a phase space, in Figure 23*A* for area 18, and in Figure 24*A* for area 2. Time series are plotted in Figure 23*B* for SLC and in Figure 23*C* for LC, for area 18. Time series are plotted in Figure 24*B* for SLC and in Figure 24*C* for LC, for area 2.

*(III) SLC:* Area III in Figure 34 corresponds to area 3 in Figure 35. The *z*-nullcline intersects the
-curve at one periodic orbit in area 3, which corresponds to SLC. Area 3 is below the SNPO bifurcation (bold line), which means that SLC exists after an SNPO bifurcation occurs. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at one equilibrium point, which is an unstable focus. Trajectories exhibit only a periodic solution with a small amplitude (SLC). The SLC behavior is plotted in Figures 13*B* and 24*B*.

*(IV) DB:* Area IV in Figure 34 is composed of areas 4 and 8 in Figure 35. Here, the *z*-nullcline and the
-curve do not intersect, hence periodic orbits LC, S, and SLC do not exist. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at different equilibrium points. The equilibrium point is a stable node in area 8. The equilibrium points are a stable node, an unstable focus, and a saddle in area 4. The Epileptor remains in DB after a transient period, which is plotted in Figures 16*D* and 27*B*.

*(V) LC and a DB:* Area V in Figure 34 is composed of areas 5 and 9 in Figure 35. The *z*-nullcline intersects the
-curve at two periodic orbits, which correspond to LC and S. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at different equilibrium points. The equilibrium points are a stable node, an unstable focus, and a saddle in area 5. The equilibrium point is a stable node in area 9. Trajectories either exhibit a fast-slow cyclic behavior (LC) or enter into a DB. The coexistence of LC and DB is plotted in Figure 27*A*. Time series are plotted in Figure 27*B* for DB and in Figure 27*C* for LC.

*(VI) LC and a nonoscillatory state:* Area VI in Figure 34 is composed of areas 6, 10, and 11 in Figure 35. The *z*-nullcline intersects the
-curve at two periodic orbits, which correspond to LC and S. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at different equilibrium points. The equilibrium points are a stable focus, an unstable focus, and a saddle in area 6. The equilibrium point is a stable focus in area 10 and a stable node in area 11. Trajectories either exhibit a fast-slow cyclic behavior (LC) or remain in a nonoscillatory state. The coexistence of LC and the nonoscillatory state is plotted in Figures 29*A* and 30*A*, which corresponds to area 10 in Figure 35. Time series are plotted in Figures 29*B* and 30*B* for the nonoscillatory state, and in Figures 29*C* and 30*C* for LC.

*(VII) Periodic switch between DB and NS:* Area VII in Figure 34 corresponds to area 12 in Figure 35. Here, the *z*-nullcline and the
-curve do not intersect, hence LC, S, and SLC do not exist. The *z*-nullcline intersects the *Z*-middle branch which consists of saddles. Then only a saddle equilibrium point exists in area 12. The ictal state of SLE reduces to DB with decreasing *m*. Periodic switch between DB and NS occurs through a fold/fold bifurcation, and is plotted in Figures 14 and 25*B*.

*(VIII) LC and a periodic switch between DB and NS:* Area VIII in Figure 34 corresponds to area 13 in Figure 35. The *z*-nullcline intersects the
curve at two periodic orbits, which correspond to LC and S. The equilibrium point is a saddle, and the ictal state of an SLE is reduced to DB, with decreasing *m*. Periodic switch between DB and NS occurs through a fold/fold bifurcation, and coexists with LC. Both are separated by a saddle periodic orbit (S). The coexistence of LC and periodic switch between DB and NS is plotted in Figure 25*A*. Time series are plotted in Figure 25*B* for periodic switch between DB and NS, and in Figure 25*C* for LC.

*(IX) LC and a periodic switch between nonoscillatory state and NS:* Area IX in Figure 34 corresponds to areas 14 in Figure 35. The *z*-nullcline intersects the
curve at two periodic orbits, which correspond to LC and S. Here, there is a unique saddle equilibrium point. The SLE attractor reduces to a periodic switch between nonoscillatory state and normal state, which occurs through a fold/fold bifurcation in area 14, and coexists with LC. Both are separated by a saddle periodic orbit (S). This coexistence is plotted in Figure 22*A*. Time series are plotted in Figure 22*B* for the periodic switch between nonoscillatory state and NS, and in Figure 22*C* for LC.

*(X) LC and SLE with a fold/Hopf bifurcation:* Area X in Figure 34 corresponds to area 15 in Figure 35. The *z*-nullcline intersects the
-curve at two periodic orbits, which correspond to LC and S. There is a unique saddle equilibrium point, and an SLE occurs through a fold/Hopf bifurcation (Fig. 7). Then LC and SLE with a fold/Hopf bifurcation coexist and are separated by a saddle periodic orbit (S). This coexistence is plotted in Figure 22*A*. Time series are plotted in Figure 21*B* for SLE and in Figure 21*C* for LC.

*(XI) LC and SLE with a fold/homoclinic bifurcation:* Area XI in Figure 34 corresponds to area 16 in Figure 35. The *z*-nullcline intersects the
-curve at two periodic orbits, which correspond to LC and S. A unique saddle equilibrium point exists, and an SLE occurs through a fold/homoclinic bifurcation in area 16 (Fig. 6). Then LC and SLE coexist and are separated by a saddle periodic orbit (S). This coexistence is plotted in Figure 19*A*. Time series are plotted in Figure 19*B* for SLE and in Figure 19*C* for LC.

*(XII) A NS:* Area XII in Figure 34 corresponds to area 19 in Figure 35. Here, the *z*-nullcline and the
-curve do not intersect, hence periodic orbits LC, S, and SLC do not exist. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at the *Z*-lower branch, and then the equilibrium point is a stable node in area 19. The Epileptor remains in NS after a transient period, which is shown in Figure 17.

*(XIII) LC and a NS:* Area XIII in Figure 34 corresponds to area 20 in Figure 35. The *z*-nullcline intersects the
-curve at two periodic orbits, which correspond to LC and S. The *z*-nullcline intersects the (*z*, *x*_{1}) curve at the *Z*-lower branch, and hence the equilibrium point is a stable node in area 20. Then LC and NS coexist and are separated by a saddle periodic orbit (S). This coexistence is plotted in Figure 26*A*. Time series are plotted in Figure 26*B* for NS and in Figure 26*C* for LC.

#### Epileptor behaviors as a function of *m* and *x*_{0}

Different areas exist in the parameter space of the Epileptor model as *m* and *x*_{0} vary (Figs. 31, 34). We show the behavior of the SLE attractor in Figure 36*A*, *a*, *b*, *b*_{1}, and *c*_{1}; the periodic switch between DB and NS in Figure 36 *A*, *a*_{1}; the periodic switch between nonoscillatory state and NS in Figure 36 *A*, *c*; the nonoscillatory state in Figure 36*B*, *a*; the normal state in Figure 36*B*, *b*; *s* and the chaotic state in Figure 37.

*SLE*: it occurs through:

A fold/circle bifurcation in Figure 36

*A*,*a*and*b*, and corresponds to areas VI and VII, respectively, in Figure 31.A fold/homoclinic bifurcation in Figure 36

*A*,*b*_{1}, and corresponds to area XI in Figure 34.A fold/Hopf bifurcation in Figure 36

*A*,*c*_{1}, and corresponds to area X in Figure 34.

*Periodic switch between DB and NS*: The switch occurs through a fold/fold bifurcation in Figure 36*A*, *a*_{1}, and corresponds to area VII in Figure 34, which shows a periodic switch between DB and NS.

*Periodic switch between a nonoscillatory state and NS:* The switch occurs through a fold/fold bifurcation in Figure 36*A*, *c*, and corresponds to area IX in Figure 34, which shows a periodic switch between nonoscillatory state and NS.

*Nonoscillatory state*: The Epileptor remains in the nonoscillatory state after a transient seizure-like fast discharges (Fig. 36*A*). Parameter settings correspond to area V in Figure 31 and to area 6 in Figure 32.

*Depolarization block*: The Epileptor remains in the DB after a transient NS (Fig. 27*A*). Parameter settings correspond to area V in Figure 34 and to area 5 in Figure 32.

*Normal state*: The Epileptor remains in the NS after a transient seizure-like fast discharges (Fig. 36*B*, *b*). Parameter settings correspond to area IX in Figure 31 and to area 17 in Figure 32.

*Chaotic state*: The chaotic state exists in areas XI and XII (
; Fig. 31). Time series of the Epileptor model *ψ*, subsystem 1 *ψ*_{1}, and subsystem 2 *ψ*_{2} are plotted in Figure 37. Figure 37*A* correspond to area XI in Figure 31, and Figure 37*B* correspond to area XII in Figure 31.

### Subsystem 1

Here, we present results on the analysis of the subsystem 1 dynamics without coupling. We analytically determine the equilibrium points and use tables to classify them, according to the trace and the determinant of the Jacobian matrix. We present the different tables as a function of the parameter *m*.

#### Analysis of subsystem 1

##### Subsystem 1 equilibrium points and stability

We analytically find the equilibrium points (*x*_{1}, *y*_{1}) by solving the following equations:
(36)where *x*_{1} is a solution of:
(37)
(38)

Let be the discriminant of Equation 37 and be the discriminant of Equation 38.

If , Equation 37 has one solution.

If , Equation 37 has one solution.

If , Equation 37 has three solutions.

The solutions are equilibrium points if .

If , Equation 38 has two solutions, and are written as: (39)

Let *X * = * z*−4, then
and *X* is a solution of
(40)where
(41)

We find
and
as *m* varies in Figure 38. Solutions of Equation 38 are equilibrium points if
and
. Then an equilibrium point *x*_{1} (*x*_{1} ≥ 0) exists if:
(42)where *X* is a solution of Equations 40 and 41.

To determine the equilibrium points stability, we analyze the eigenvalues of the following Jacobian matrix *J*:
where
(43)

*J* is defined at the equilibrium point (*x*_{1}, *y*_{1}). To find the eigenvalues of *J*, we solve the characteristic equation:
(44)where *Tr*(*J*) and *Det*(*J*) are the trace and the determinant of the matrix *J*, respectively.

The trace and the determinant of *J* are given by:
(45)
(46)

An equilibrium point is stable if all the real parts of the eigenvalues of J are negative (Izhikevich, 2007). We can classify the equilibrium points according to the trace and the determinant of *J* (Izhikevich, 2007). We determine the intervals on which the *Tr*(*J*) and *Det*(*J*) are negative and positive.

The roots of *Tr*(*J*) are as follows:
(47)

The roots of *Det*(*J*) are as follows:
(48)where *X* is a solution of Equations 40 and 41.

We conclude that an equilibrium point *x*_{1} (*x*_{1} ≥ 0) exists (see Eq. 42) if:
(49)

We find the stability of equilibrium points
and
in Tables 2–Tables 5, depending on *m*.

We graphically determine the equilibrium points in a phase plane. The equilibrium points lie at the intersection of the *x*_{1}- and *y*_{1}-nullclines. The *x*_{1}-nullcline (
) corresponds to a cubic curve for *x*_{1} < 0 and a straight line for *x*_{1} ≥ 0.

The cubic curve is written in the form: (50)and the straight line is written as follows: (51)

The *y*_{1}-nullcline (
) corresponds to a parabola given by:
(52)

We plot the phase plane of the subsystem 1 in Figure 39*A* for *z * = * *3.1 and in Figure 39*B* for *z * = * *0. We consider only the negative *x*_{1}-axis for the cubic curve, and plot a straight line of the *x*_{1}-nullcline (
) for different *m*.

When *z * = * *3.1, then
(Fig. 39*A*). We plot a straight line of the *x*_{1}-nullcline for *m * = * *0, *m *=* * 1, and *m * =* * 1.5. Then the *x*_{1}- and *y*_{1}-nullclines intersect at three points for each *m* value. We plot four trajectories with different initial conditions (*i*_{1}, *i*_{2}, *i*_{3}, and *i*_{4}). Trajectories *i*_{1} and *i*_{2} converge to one equilibrium point, which is a stable node. Trajectories *i*_{3} and *i*_{4} converge to one equilibrium point, which is a stable focus when *m * = * *0 (Fig. 39*A*, *m * = * *0). By increasing (*m * = * *1 and *m * = * *1.5), the stable focus loses its stability, and the equilibrium point is an unstable focus. Trajectories *i*_{3} and *i*_{4} converge to a stable limit cycle, which surrounds the unstable focus. The radius of the limit cycle increases as *m* increases (Fig. 39*A*, *m * = * *1 and *m * = * *1.5).

We conclude that when *z * = * *3.1, three equilibrium points coexist. The first one is a stable node and the second one is a stable or an unstable focus, depending on *m*. The third equilibrium point is a saddle. The stable manifold of the saddle corresponds to a separatrix between the first and second equilibrium point. The subsystem 1 undergoes a supercritical Andronov–Hopf bifurcation, *H*, as *m* increases at *m*(*H*)* * = * *0.514 [i.e., solution of *Tr*(*J*)* * = * *0; Eq. 45].

When *z * = * *0, then
decreases and the *x*_{1}-nullcline moves downward in the phase plane (*m* = −10; Fig. 39*B*). We plot a straight line of the *x*_{1}-nullcline (
) for *m * = * *0 and *m* = −10. Then the *x*_{1}- and *y*_{1}-nullclines intersect at one equilibrium point, which is a stable focus for *m* = −10 and an unstable focus for *m * = * *0 (Fig. 39*B*). The stable node and saddle disappear. We plot trajectories in the phase plane, which converge to a stable focus for *m* = −10 and to a stable limit cycle for *m * = * *0. The limit cycle surrounds the unstable focus. The subsystem 1 undergoes a supercritical Andronov–Hopf bifurcation, *H*, at *m*(*H*) = −8.6 [i.e., solution of *Tr*(*J*)* * = * *0; Eq. 45].

We now fix *m* on two values to identify the equilibrium points as *z* varies.

Let

*m*= 0, three equilibrium points coexist when*z*= 3.1: a stable node, a saddle, and a stable focus (Fig. 39*A*). When*z*= 0, then only an unstable focus exists (Fig. 39*B*). Hence, as*z*decreases, a stable focus becomes unstable, and the subsystem 1 undergoes a supercritical Andronov–Hopf bifurcation,*H*. Moreover, the stable node and saddle disappear when*z*= 0 (Fig. 39*B*). In fact, the stable node and saddle approach each other as*z*decreases and coalesce at*z*≈ 2.9 (figure not shown), which corresponds to (see Eq. 48). Decreasing*z*further, the stable node and saddle disappear through a saddle-node bifurcation.Let

*m*= 1.5, three equilibrium points coexist when*z*= 3.1: a stable node, a saddle, and an unstable focus (Fig. 39*A*). When*z*decreases to*z*≈ 2.9, then a saddle-node bifurcation occurs (figure not shown). When*z*= 2.7, the stable node and saddle disappear (Fig. 39*C*). In contrast, when*z*= 4 (≥ 3.1), three equilibrium points coexist: a stable node, a saddle, and an unstable focus (Fig. 39*C*). A stable limit cycle surrounds the unstable focus when*z*= 3.1 (Fig. 39*A*), and it does not when*z*= 4 (Fig. 39*C*). In fact, the stable limit cycle approaches as*z*increases to the saddle, touches the saddle, and then terminates. Hence, when*z*= 4, the stable limit cycle disappears through a homoclinic bifurcation, and the trajectories converge to a stable node.

We conclude that the stability of equilibrium points depends on *z* and *m*. In addition, we find three bifurcation types: a saddle-node bifurcation, a homoclinic bifurcation, and a supercritical Hopf bifurcation.

##### Subsystem 1 bifurcation diagram

We find the equilibrium points as *z* varies in a (*z*, *x*_{1}) bifurcation diagram plotted in Figure 40 for *m * = * *0. The (*z*, *x*_{1}) curve comprises two geometrical shapes. The right one comprises two branches: the lower (dash-dotted) branch consists of saddles; and the upper (dashed) branch consists of unstable foci. When decreasing *z*, the lower and upper branches collide in a saddle-node bifurcation, SN. The left geometrical shape corresponds to a *Z*-curve known from fast-slow subsystems (Ermentrout and Terman, 2010). Left and right shapes are separated by SN_{2} and SN, such that
and
(Eq. 42). The discriminant
of Equation 38 is negative between SN_{2} and SN. Lower and upper branches (right shape) consist of saddles and unstable nodes, respectively
. Hence, in what follows, we plot only the *Z*-curve (left shape) in the next bifurcation diagrams.

We plot a (*z*, *x*_{1}) bifurcation diagram for *m * = * *0 in Figure 41*A*. The geometrical shape of the bifurcation diagram corresponds to a *Z*-curve, which comprises four branches. *Z*-lower (solid) and *Z*-middle (dash-dotted) branches consist of stable nodes and saddles, respectively. Decreasing *z*, *Z*-lower and *Z*-middle branches collide in a saddle-node bifurcation point, SN_{1}. The *Z*-middle branch acts as a separatrix between *Z*-lower and *Z*-upper branches.

The *Z*-upper branch comprises two sub-branches separated by a Hopf bifurcation, *H*: one sub-branch (solid) consists of stable foci and another (dashed) consists of unstable foci (Fig. 41*A*). The first sub-branch corresponds to a nonoscillatory state, and the second sub-branch corresponds to an oscillatory state. The subsystem 1 then switches as *z* decreases from nonoscillatory to oscillatory states.

An unstable focus is surrounded by a stable limit cycle, which is limited by red (
) and green (
) curves. The stable limit cycle radius increases as *z* is decreased. The stable limit cycle terminates in a Hopf bifurcation, *H*, written in the form:
(53)which corresponds to
(see Eq. 47),
). The Hopf bifurcation point *z*(*H*) exists for *m * ≤* * 1 moving rightward (leftward) on the *Z*-upper curve as *m* increases (decreases), and does not exist for *m * >* * 1.

We plot a (*z*, *x*_{1}) bifurcation diagram for *m * = * *2 in Figure 41*B*. *Z*-lower (solid) and *Z*-middle (dash-dotted) branches consist of stable nodes and saddles, respectively. *Z*-upper (dashed) branch consists only of unstable foci, which are surrounded by stable limit cycles. The stable limit cycle amplitude is limited by red (
) and green (
) curves. The stable limit cycle terminates as *z* increases in a homoclinic bifurcation HB. Decreasing *z*, *Z*-lower and *Z*-middle branches collide in a saddle-node bifurcation SN_{1}. Increasing *z*, *Z*-upper and *Z*-middle branches collide in a saddle-node bifurcation SN_{2} (Fig. 41*B*, inset).

SN_{1} and SN_{2} correspond to saddle-node bifurcation points, and HB corresponds to a homoclinic bifurcation point. SN_{1} and SN_{2} points are given by the following:
(54)which corresponds to
(see Eq. 48),
),
(55)where *X* is a solution of Equations 40 and 41. We find saddle-node bifurcation points SN_{1} and SN_{2}, and a Hopf bifurcation point, H, in Figure 42 as *m* varies. Figure 42 shows the following three cases:

For

*m*≤ 0.29, then , and the subsystem 1 is bistable on [SN_{1}, SN_{2}].For , then , and the subsystem 1 is bistable on [SN

_{1}, SN_{2}].For

*m*> 1, then H does not exist, and the subsystem 1 is bistable on [SN_{1}, HB].

##### Fast-slow subsystem equilibrium points

We analytically find the equilibrium points by solving the Equations 37 and 38, where *z* is a solution of
(Eq. 18). We determine the stability of the equilibrium points by analyzing the Jacobian matrix *J* written as:
where
(56)
(57)
(58)

*J* is defined at the equilibrium point (*x*_{1}, *y*_{1}, *z*). We find the eigenvalues of *J* by solving the characteristic equation:
(59)

The trace, *Tr*(*J*), and the determinant, *Det*(*J*), of *J* are given by:
(60)
(61)

We graphically determine the equilibrium points using the nullclines. The *x*_{1}-nullcline corresponds to the Equations 50 and 51, and the *y*_{1}-nullcline corresponds to Equation 52.

The *z*-nullcline (
, Eq. 18) corresponds to a straight line for *z * ≥ * *0 written in the form:
(62)and to a curve for *z *<* *0 written as:
(63)

The equilibrium points lie at the intersection of the *z*-, *x*_{1}-, and *y*_{1}-nullclines. We graphically find the equilibrium points using the bifurcation diagram of the subsystem 1 (Fig. 41). We plot the *Z*-curve using the intersection of the *x*_{1}- and *y*_{1}-nullclines for each *z*. Therefore, the equilibrium points of the fast-slow subsystem lie at the intersection of the *z*-nullcline and the *Z*-curve (Fig. 43). A bifurcation diagram of the subsystem 1 is plotted in Figure 43*A* for *m * = * *0 and in Figure 43*B* for *m * = * *2. When *x*_{0} = −1.6, then the *z*-nullcline is at the *Z*-middle branch for *m * =* * 0 (Fig. 43*A*) and *m * = * *2 (Fig. 43*B*), and then only a saddle exists. When *x*_{0} = −1.9, then the *z*-nullcline is at the *Z*-middle branch for *m * =* * 2, and the equilibrium point is a saddle (Fig. 43*B*).

The *z*-nullcline moves downward (*x*_{0} decreases) or upward (*x*_{0} increases) in the bifurcation diagram, and intersects the (*z*, *x*_{1}) curve at different points. To determine these points as *x*_{0} varies, we plot a (*x*_{0}, *x*_{1}) diagram in Figure 44*A* for *m * = * *0 and in Figure 44*B* for *m * = * *2. The curve consists of equilibrium points of the fast-slow subsystem. Blue (bottom solid), black (plus sign markers), red (dashed), and green (top solid) branches consist of stable nodes, saddles, unstable foci, and stable foci, respectively. The branch of stable foci exists for *m * = * *0 but does not for *m * = * *2. SN_{1} and SN_{2} correspond to saddle-node bifurcation points, and *H* to a Hopf bifurcation point (Fig. 44).

##### Fast-slow subsystem behavior

The stability of equilibrium points depends on *m* and *x*_{0}. Different behavioral patterns are discovered as *m* and *x*_{0} vary: a resting state, a fold/fold bifurcation, a fold/homoclinic bifurcation, a fold/Hopf bifurcation and a periodic solution (Fig. 45). To observe this dynamical change, we integrated the fast-slow subsystem equations for different *m* and *x*_{0}. We plot time series of the fast-slow subsystem in Figure 45, and trajectories in Figure 46, as *m* and *x*_{0} vary.

*NS:* We plot a (*z*, *x*_{1}) bifurcation diagram for *m * = * *0 in Figure 46*D*. Let *x*_{0} = −2.6, the *z*-nullcline is at the *Z*-lower branch. The equilibrium point is a stable node. When a trajectory is at the *Z*-upper (dashed) sub-branch, then it exhibits an oscillatory solution, which terminates as *z* increases in a Hopf bifurcation *H*. Then, the trajectory exhibits a nonoscillatory solution on the *Z*-upper (solid) sub-branch, which terminates as *z* increases in a saddle-node bifurcation point, SN_{2}, and the trajectory switches to the *Z*-lower branch. Decreasing *z*, the trajectory continues to *C* at which *z* stabilizes. Time series are plotted in Figure 45, *E* and *E1*.

If *m * = * *2, then the trajectory exhibits an oscillatory solution on the *Z*-upper branch, which terminates as *z* increases in a homoclinic bifurcation and the trajectory switches to the *Z*-lower branch. The trajectory continues as *z* decreases to *C*, at which *z* stabilizes (figure not shown).

*Fold/Hopf bifurcation:* We plot a (*z*, *x*_{1}) bifurcation diagram for *m * = * *0.5 in Figure 46*A*. When *m * = * *0.5, then
(Fig. 42). Let *x*_{0} = −1.6, the *z*-nullcline is at the *Z*-middle branch. The equilibrium point is a saddle. When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation, SN_{1}. Then, the trajectory switches to the *Z*-upper (dashed) sub-branch exhibiting an oscillatory solution, which terminates as *z* increases in a Hopf bifurcation, *H*. The trajectory exhibits a nonoscillatory solution on the *Z*-upper (solid) sub-branch, which terminates as *z* increases in a saddle-node bifurcation, SN_{2}, and then switches to the *Z*-lower branch. Therefore, the transitions between *Z*-lower and *Z*-upper branches occur through a fold/Hopf bifurcation. Time series are plotted in Figure 45*F*.

*Fold/fold bifurcation:* We plot a (*z*, *x*_{1}) bifurcation diagram for *m * = * *0 in Figure 43*A*. When *m * = * *0, then
(Fig. 42). Let *x*_{0} = −1.6, the *z*-nullcline is at the *Z*-middle branch (Fig. 43*A*). The equilibrium point is a saddle. If a trajectory is at the *Z*-upper (dashed) sub-branch, then it exhibits an oscillatory solution, which terminates as *z* increases in a Hopf bifurcation *H*. The trajectory exhibits a nonoscillatory solution on the *Z*-upper (solid) sub-branch, which terminates as *z* increases in a saddle-node bifurcation SN_{2}, and then switches to the *Z*-lower branch. Decreasing *z*, the stable node disappears through a saddle-node bifurcation SN_{1} and the trajectory switches to the *Z*-upper (solid) branch, which consists of stable foci. Therefore, the transitions between *Z*-lower and *Z*-upper branches occur through a fold/fold bifurcation. Time series are plotted in Figure 45*A*.

*Fold/homoclinic bifurcation:* We plot a (*z*, *x*_{1}) bifurcation diagram for *m * = * *2 in Figures 43*B* and 46*B*. Here, a Hopf bifurcation points does not exist (Fig. 42). Let *x*_{0} = −1.9, the *z*-nullcline is at the *Z*-middle branch. The equilibrium point is a saddle. When a trajectory is at the *Z*-upper (dashed) branch, it exhibits an oscillatory solution, which terminates as *z* increases in a homoclinic bifurcation HB, and then switches to the *Z*-lower branch. Decreasing *z*, the stable node disappears through a saddle-node bifurcation, SN_{1}, and the trajectory switches to the *Z*-upper branch. Therefore, the transitions between *Z*-lower and *Z*-upper branches occur through a fold/homoclinic bifurcation. Time series are plotted in Figure 45*D*.

*Periodic solution (SLC _{1}):* We plot a (

*z*,

*x*

_{1}) bifurcation diagram for

*m*=

*2 in Figure 46*

*C*. Let

*x*

_{0}= –1.7, the

*z*-nullcline is at the

*Z*-middle branch. The equilibrium point is a saddle. When a trajectory is at the

*Z*-lower branch, the stable node disappears as

*z*decreases through a saddle-node bifurcation, SN

_{1}. Then, the trajectory switches to the

*Z*-upper (dashed) branch exhibiting as

*z*increases an oscillatory solution which stabilizes at

*C*. The

*z*-nullcline intersects the -curve at

*C*, hence the system solution is periodic and corresponds to a stable limit cycle, with a small amplitude denoted by SLC

_{1}. Time series are plotted in Figure 45

*C*.

Figure 45*B* shows the existence of a periodic solution with a large amplitude when *m * = * *2 and *x*_{0} = −1.6. Here, the equilibrium point is a saddle (Fig. 44*B*). The solution is a stable limit cycle with a fast-slow cyclic behavior, which corresponds to LC.

#### Finding LC

LC and SLC_{1} are stable limit cycles with large and small amplitudes, respectively. Here we show when periodic orbits exist and how they evolve, determining their stability.

##### Sensitivity to initial conditions

The Epileptor model is sensitive to initial conditions. In fact, integrating the fast-slow subsystem in Equations 18–20 from an initial condition results in different behaviors as *m* and *x*_{0} are varied (Fig. 45*A–F*, left). We considered another initial condition and used numerical integration techniques to solve the fast-slow subsystem equations. Figure 45*A1–F1*, (right), shows two different solutions as *m* and *x*_{0} vary. The first solution corresponds to a stable LC and appears in Figure 45, *A1–D1*, and *F1*. The second one corresponds to a normal state, which exists when *x*_{0} = −2.6 (Fig. 45*E*_{1}).

We now fix *m* and *x*_{0}, and compare the solutions of the fast-slow subsystem as the initial conditions vary:

When

*m*= 0 and*x*_{0}= −1.6, two solutions coexist: a fold/fold bifurcation (Fig. 45*A*) and LC (Fig. 45*A*)._{1}When

*m*= 0.5 and*x*_{0}= −1.6, two solutions coexist: a fold/Hopf bifurcation (Fig. 45*F*) and LC (Fig. 45*F*)._{1}When

*m*= 2 and*x*_{0}= −1.9, two solutions coexist: a fold/homoclinic bifurcation (Fig. 45*D*) and LC (Fig. 45*D*)._{1}When

*m*= 2 and*x*_{0}= −1.7, then two periodic solutions coexist: SLC_{1}(Fig. 45*C*) and LC (Fig. 45*C*)._{1}

Therefore, the fast-slow subsystem is sensitive to initial conditions and exhibits a bistability of LC and SLC_{1}, a fold/fold bifurcation, a fold/Hopf bifurcation or a fold/homoclinic bifurcation.

##### Finding periodic orbits

For *x*_{0} = −2.6, the normal state exists and LC does not (Fig. 45*E*,*E _{1}
*). However, when

*m*=

*2 and*

*x*

_{0}= −1.6, only the LC exists (Fig. 45

*B*,

*B*). We explain this using the averaging method. We introduce , which is the average value of

_{1}*x*

_{1}-coordinate associated with the subsystem 1, given by: (64)where .

We plot
in a (*z*, *x*_{1}) bifurcation diagram (Fig. 43). Periodic orbits lie at the intersection of the *z*-nullcline and the
-curve. A periodic orbit branch is referred to the *Z*-upper branch, which is surrounded by periodic orbits. This branch and the
-curve terminate as *z* increases in a Hopf bifurcation for *m * ≤* * 1 (Fig. 43A) and in a homoclinic bifurcation for *m * >* * 1 (Fig. 43*B*).

Let *x*_{0} = −1.6, Figure 43 shows that the *z*-nullcline intersects the
-curve at two points: LC and S when *m * = * *0, and only at LC when *m * = * *2. However, LC and S coexist when *m * = * *2 and *x*_{0} = −1.9 (Fig. 43*B*).

We graphically determine the periodic orbit by using the Pontryagin’s averaging technique (Shilnikov and Kolomiets, 2008). We introduce the following slow averaged nullcline: (65)where the periodic orbits correspond to the zeros of the .

We determine the slow averaged nullcline by substituting the *x*_{1}-coordinates by their averaged values in the *z* dynamics Equation 18 and setting the right-hand side of the equation equal to zero (Eq. 65). The slow averaged nullcline depends on *x*_{0} and *z*. We determine the periodic orbit stability by using the following derivative:
(66)that represents the dynamics of the averaged equation. A periodic orbit is stable when Equation 66 is negative and is unstable when Equation 66 is positive.

We plot the
-curve in Figure 47*A* for *m * = * *0 (Figure 47*A*,*a*), and *m * = * *2 (Figure 47*A,b*). *z* is a control parameter. The maximal value of *z* corresponds to a Hopf bifurcation point H for
and to a homoclinic bifurcation point HB for *m *> 1. The periodic orbits are the zeros of the
. A periodic orbit is stable when the
-curve decreases over *z* at the given zero and is unstable when the
curve increases. Stable and saddle periodic orbits are labeled as black circles and squares respectively.

Let

*m*= 0 (Fig. 47*A*,*a*), has a simple zero when*x*_{0}= 1, which is stable. Decreasing*x*_{0}to −0.25, another periodic orbit is born, which is unstable (saddle). When*x*_{0}= −2.45, stable and saddle periodic orbits coalesce forming a saddle-node periodic orbit, which then fades . Stable and saddle periodic orbits disappear through a SNPO bifurcation (Shilnikov and Kolomiets, 2008).Let

*m*= 2 (Fig. 47*A*,*b*), has a simple zero when*x*_{0}= −1.6, which is stable. When has two zeros: one stable and another saddle. When has three zeros: two stable periodic orbits separated by a saddle periodic orbit. By decreasing*x*_{0}, SNPO bifurcation occurs.

We conclude that LC and S are stable and saddle periodic orbits, respectively (Fig. 43). When *m * = * *2 and *x*_{0} = −1.7, then two zeros of
are LC and S. The third zero is a stable periodic orbit, which we denote by SLC_{1}.

##### Evolution of periodic orbits

LC and S are stable and saddle periodic orbits, respectively. SLC_{1} is a stable periodic orbit, which exists for *m * = * *2. We now analyze how the periodic orbits evolve when *x*_{0} varies slowly. We plot a (*x*_{0}, *z*^{*}) bifurcation diagram of periodic orbits in Figure 47*B* for *m * = * *0 (Figure 47*B*,*a*), and *m * = * *2 (Figure 47*B*, *b*). We localize the periodic orbits at *z*^{*}. *x*_{0} is a control parameter.

Let

*m*= 0, only LC exists [red bottom (+) markers] for large*x*_{0}. Decreasing*x*_{0}, S (black dotted squares) appears. Stable (LC) and saddle (S) periodic orbits approach each other as*x*_{0}decreases, collide in a saddle-node periodic orbit bifurcation, and then fades. LC and S disappear through a SNPO bifurcation (Fig. 47*B*,*a*).Let

*m*= 2, only LC exists for large*x*_{0}. Decreasing*x*_{0}, S and SLC_{1}[blue top (+) markers] appear. Hence, LC and SLC_{1}coexist and are separated by S. Decreasing*x*_{0}, SLC_{1}disappears. Further decreasing*x*_{0}, LC and S disappear through an SNPO bifurcation (Fig. 47*B*,*b*).

We can now explain the solutions found in Figure 45. We plot in the (*x*_{0}, *z*^{*}) bifurcation diagram the equilibrium points, labeled as dots (Fig. 47*B*). Blue (segment 1), black (segment 2), green (segment 3), and purple (segment 4) dots correspond to stable nodes, saddles, stable foci, and unstable foci, respectively. For instance, when *m * = * *0 and *x*_{0} = −2.6, a stable node exists and LC does not exist, whence trajectories exhibit a normal state in Figure 45, *E* and *E1*. We discuss all solutions found in Figure 45 in the sections: “Monostability in the fast-slow subsystem” and “Coexistence in the fast-slow subsystem”.

#### Periodic switch between a DB and a NS in the fast-slow subsystem

When *m * ≤* * 1, a Hopf bifurcation *H* partitions the *Z*-upper branch into two sub-branches: one consists of an unstable focus (
) and another consists of a stable focus (
; Fig. 41*A*). When decreasing *m*, the imaginary part of the complex-conjugate eigenvalues corresponding to the stable focus goes to zero, and then the stable focus becomes a stable node. We plot a (*z*, *x*_{1}) bifurcation diagram when *m* = −8 in Figure 48*A*. The *Z*-lower branch consists of stable nodes, which are the equilibrium points of the NS. The *Z*-upper branch consists of stable nodes (
) which are responsible for the silent activity, similar to a DB. The *z*-nullcline is at the *Z*-middle branch which consists of saddles (*x*_{0} = −1.6). When a trajectory is at the *Z*-lower branch, the stable node disappears as *z* decreases through a saddle-node bifurcation SN_{1} and then the trajectory enters into a depolarization block, which terminates as *z* increases in a saddle-node bifurcation, SN_{2}. Thus, the transitions between the *Z*-lower and *Z*-upper branches are reduced to a periodic switch between a NS and a DB. Time series are plotted in Figure 49*C*.

#### Stabilizing equilibrium points in the fast-slow subsystem

The transitions between *Z*-lower and *Z*-upper branches occur when the equilibrium point is a saddle (*x*_{0} = −1.6). When *x*_{0} is increased, the fast-slow subsystem stabilizes the equilibrium point of the *Z*-upper branch, and when *x*_{0} is decreased, the fast-slow subsystem stabilizes the equilibrium point of the *Z*-lower branch (normal state).

##### Stabilizing the equilibrium point of the nonoscillatory state

Let *x*_{0} = −0.5. We plot a (*z*, *x*_{1}) bifurcation diagram when *m * = * *0 in Figure 48*C*. The *z*-nullcline is at the *Z*-upper branch which consists of a stable focus. Decreasing z, the stable node disappears through a saddle-node bifurcation SN_{1} and the fast-slow subsystem switches to the *Z*-upper branch (dashed) exhibiting an oscillatory solution, which terminates in a Hopf bifurcation, H. The fast-slow subsystem exhibits a nonoscillatory solution at the *Z*-upper branch (solid) and continues to *C*, at which *z* stabilizes. The imaginary part of the eigenvalues of *C* is responsible for the oscillations around the stable focus, which is the equilibrium point of the nonoscillatory state. Time series are plotted in Figure 49*A* for the nonoscillatory state.

##### Stabilizing the equilibrium point of the DB

Let *x*_{0} = −0.5. We plot a (*z*, *x*_{1}) bifurcation diagram when *m* = −8 in Figure 48*D* The *z*-nullcline is at the *Z*-upper branch which consists of stable nodes. Decreasing *z*, the stable node disappears through a saddle-node bifurcation, SN_{1}, then the fast-slow subsystem switches to the *Z*-upper branch (solid) exhibiting DB and continues to *C*, at which *z* stabilizes. The imaginary part of the eigenvalues of *C* goes to zero, and therefore the oscillations around *C* disappear. Here, the stable node (*C*) is the equilibrium point of DB. Time series are plotted in Figure 49*D* for the DB.

##### Stabilizing the equilibrium point of the NS

Let *x*_{0} = −2.4. We plot a (*z*, *x*_{1}) bifurcation diagram when *m * = * *0 in Figure 48*B*. The *z*-nullcline is at the *Z*-lower branch which consists of stable nodes. When a trajectory is at the *Z*-upper branch (dashed), it exhibits an oscillatory solution, which terminates as *z* increases in a Hopf bifurcation, *H*, and then a nonoscillatory solution on the *Z*-upper branch (solid), which terminates as *z* increases in a saddle-node bifurcation SN_{2}. The fast-slow subsystem continues to *C*, at which *z* stabilizes. The stable node (*C*) is the equilibrium point of the NS. Time series are plotted in Figure 49*B* for the NS.

#### Monostability in the fast-slow subsystem

Figure 45 shows that the fast-slow subsystem presents a monostability in two cases:

*Case 1:* Only LC exists.
has a simple zero for *m * = * *2 and *x*_{0} = −1.6, which is LC (Fig. 47*A*,*b*). Then only LC exists, to which all trajectories converge (Fig. 45,*B* and *B _{1}
*).

*Case 2:* Only normal state exists. The *z*-nullcline and
-curve do not intersect when *m * = * *0 and *x*_{0} = −2.6 (Fig. 46*D*). Then LC disappears through an SNPO bifurcation (Fig. 47*B*,*a*). As a consequence, only a stable node exists to which all trajectories converge (Fig. 45,*E* and *E _{1}
*).

#### Coexistence in the fast-slow subsystem

Figure 45 shows that LC coexists with SLE or SLC_{1}. Using bifurcation analysis, we identify all coexisting attractors in the fast-slow subsystem.

##### Coexistence of LC and SLE

A coexistence occurs when
has two zeros, which are LC and S. LC coexists with SLE, which occurs through a fold/Hopf bifurcation (Fig. 45,*F* and *F1*) or a fold/homoclinic bifurcation (Fig. 45,*D* and *D1*). The equilibrium point is a saddle. A saddle periodic orbit (*S*) acts as a separatrix between LC and SLE (Fig. 47*B*).

##### Coexistence of LC and SLC_{1}

has three zeros, which are LC, S, and SLC_{1} (Fig. 47*A*,*b*). LC coexists with a periodic solution (SLC_{1}; Fig. 45*C* and *C _{1}
*). The equilibrium point is a saddle. A saddle periodic orbit (S) acts as a separatrix between LC and SLC

_{1}(Fig. 47

*B*,

*b*).

##### Coexistence of LC and a NS

Figure 47*B* shows the coexistence of LC and a stable node (blue dots), which is the equilibrium point of the normal state. Therefore, LC and the normal state coexist and are separated by S. Let *m * = * *0 and *x*_{0} = −2.4, time series are plotted in Figure 49*D* for the normal state and in Figure 49*D*_{1} for LC. The bifurcation diagram is plotted in Figure 48*B*. Only the normal state is shown.

##### Coexistence of LC and a nonoscillatory state

(Figure 47*B*,*a*), shows the coexistence of LC and a stable focus (green dots), which is the equilibrium point of the nonoscillatory state. Then, LC and the nonoscillatory state coexist and are separated by S. Time series are plotted in Figure 49*A* for the nonoscillatory state and in Figure 49*A*_{1} for LC. The bifurcation diagram is plotted in Figure 48*C*. Only the nonoscillatory state is shown.

##### Coexistence of LC and a periodic switch between nonoscillatory state and NS

When *m* is decreased and the equilibrium point is a saddle, the fast-slow subsystem switches between the nonoscillatory state and the NS through a fold/fold bifurcation (Fig. 45*A*). (Figure 47*B*,*a*), shows the coexistence of LC and a saddle (black dots), which is the equilibrium point of this periodic switch. Therefore, LC and the periodic switch between nonoscillatory state and NS coexist and are separated by S. Time series are plotted in Figure 45*A* for the periodic switch between nonoscillatory state and NS, and in Figure 45*A*_{1} for LC.

##### Coexistence of LC and a periodic switch between DB and NS

When *m* is further decreased and the equilibrium point is a saddle, the fast-slow subsystem switches between DB and NS (Fig. 48*A*). Time series are plotted in Figure 49*C* for the periodic switch between DB and NS, and in Figure 49*C1* for LC. Parameters *m* and *x*_{0} are the same, only initial conditions change. Then, LC and the periodic switch between DB and NS coexist.

##### Coexistence of LC and a DB

When *m* is decreased, the fast-slow subsystem enters into DB (Fig. 49*C*). When *x*_{0} is increased, the fast-slow subsystem remains in DB (Fig. 48*D*). Time series are plotted in Figure 49*D* for DB, and in Figure 49*D*_{1} for LC. Hence, LC and DB coexist.

#### Parameter space of equilibrium points and periodic orbits

To characterize the coexisting attractors, we explore a (*m*, *x*_{0}) parameter space of the fast-slow subsystem in Figure 50, using numerical techniques. The parameter space is divided in two parts separated by a bold line (SNPO bifurcation), above which LC exists, while it does exist below it. There are 10 areas. For large values of *m* and *x*_{0} (area I) only LC exists. The adjacent area II shows bistability of LC and a stable focus. Both attractors are separated by a saddle periodic orbit. In the middle, LC coexists with the SLE attractor separated by a saddle periodic orbit. SLE occurs via a saddle-node/saddle-node bifurcation (area VI) and a saddle-node/homoclinic bifurcation (area VII). In area VIII, the homoclinic bifurcation HB is not completed and gives rise to another (coexisting) stable periodic orbit with a small amplitude, SLC_{1}. In area V, DB coexists with LC. In area IX, normal brain activity and LC coexist. In area X, only normal activity exists.

*(I) LC:* The *z*-nullcline intersects the
-curve at one point which is LC, and intersects the (*z*, *x*_{1}) curve at different equilibrium points. The equilibrium point can be a stable focus, an unstable focus, or a saddle. All trajectories converge to LC. Time series are plotted in Figure 45, *b* and *b _{1}
*, for LC.

*(II) Coexistence of LC and a stable focus:* The *z*-nullcline intersects the (*z*, *x*_{1}) curve at one equilibrium point, which is a stable focus, and intersects the
-curve at two points: LC and S. The trajectories converge to LC or a stable focus, depending on the initial conditions. A saddle periodic orbit (S) limits the basin of attraction of LC and the basin of attraction of a stable focus. The stable focus is the equilibrium point of the nonoscillatory state. Then, LC and nonoscillatory state coexist and are separated by S.

If a trajectory is in the basin of attraction of the stable focus, then it exhibits a nonoscillatory solution following two scenarios. When the trajectory is at the *Z*-lower branch, the stable node disappears through a saddle-node bifurcation SN_{1} and the trajectory switches to the *Z*-upper branch, which comprises a Hopf bifurcation point, *H*. If
(first scenario), then the trajectory exhibits an oscillatory solution, which terminates as *z* decreases in a Hopf bifurcation, *H*. Then, the trajectory continues to the stable focus at which *z* stabilizes. If
(second scenario), then the trajectory exhibits a nonoscillatory solution, which continues as *z* decreases to the stable focus at which *z* stabilizes. The stable focus is the equilibrium point of the nonoscillatory state. The system is bistable on [SN_{1}, SN_{2}].

Using analytic techniques, we observe that the stable focus reduces as *m* decreases to a stable node (i.e., the imaginary part of the complex-conjugate eigenvalues corresponding to the stable focus goes to zero). Here, the stable node is the equilibrium point of the DB, and then DB is present in area II. LC and DB coexist as *m* decreases.

*(V) Coexistence of LC and a periodic switch between DB and NS:* The *z*-nullcline intersects the (*z*, *x*_{1}) curve at one equilibrium point, which is a saddle, and intersects the
-curve at two points: LC and S. A Hopf bifurcation point, *H*, exists for *m *≤* *1 (Fig. 42). Using analytic techniques (Fig. 42), we can deduce that
. The transitions between *Z*-upper and *Z*-lower branches occur through a fold/fold bifurcation. For small *m*, the *Z*-upper branch consists of stable nodes, which is the equilibrium point of a DB. The system is bistable on [SN_{1}, SN_{2}]. The transitions between *Z*-upper and *Z*-lower branches reduce to a periodic switch between DB and NS. The LC and periodic switch between DB and NS coexist and are separated by S.

*(VI) Coexistence of LC and SLE with a fold/fold bifurcation:* Using numerical techniques, there is a saddle-node bifurcation at both offset and onset seizures. We can deduce that SLE occurs with a fold/fold bifurcation in area VI. The *z*-nullcline intersects the
-curve at two points: LC and S. As a consequence, LC and SLE with a fold/fold bifurcation coexist and are separated by S. In contrast to area V, the *Z*-upper branch consists of a stable focus, which is the equilibrium point of the nonoscillatory state. The SLE attractor reduces then to a periodic switch between nonoscillatory state and NS, which occurs through a fold/fold bifurcation, coexisting with LC. Time series are plotted in Figure 45*A* for a periodic switch between nonoscillatory state and NS, and in Figure 45*A _{1}
* for LC.

Using analytic techniques (Fig. 42), we can deduce that as *m* increases,
. This means that although there is a saddle-node bifurcation at both offset and onset seizures, there is a Hopf bifurcation occurring during the ictal period before offset seizure. Thus, an SLE occurs as *m* increases through a fold/Hopf bifurcation in area VI, coexisting with LC. Time series are plotted in Figure 45*F* for SLE with a fold/Hopf bifurcation, and in Figure 45
*F1*for LC.

*(VII) Coexistence of LC and SLE with a fold/homoclinic bifurcation:* The *z*-nullcline intersects the (*z*, *x*_{1}) curve at one equilibrium point, which is a saddle, and intersects the
-curve at two points: LC and S. Since *m *>* *1, then the transitions between *Z*-upper and *Z*-lower branches occur through a fold/homoclinic bifurcation. For *m * = * *1, a Hopf bifurcation point, *H*, exists such that
(Fig. 42). When a trajectory is at the *Z*-upper branch, then it exhibits an oscillatory solution, which touches the *Z*-middle branch before the Hopf bifurcation, *H*. Then the transitions between *Z*-upper and *Z*-lower branches occur through a fold/homoclinic bifurcation for *m * = * *1. LC and a fold/homoclinic bifurcation coexist and are separated by S. Time series are plotted in Figure 45*D* for fold/homoclinic bifurcation and in Figure 45*D _{1}
* for LC.

*(VIII) Coexistence of LC and SLC _{1}:* The

*z*-nullcline intersects the (

*z*,

*x*

_{1}) curve at one equilibrium point, which is a saddle, and intersects the -curve at three points: LC, SLC

_{1}, and S. Then LC and SLC

_{1}coexist and are separated by S. If a trajectory is in the basin of attraction of SLC

_{1}, then it exhibits a periodic solution with a small amplitude. When the trajectory is at the

*Z*-lower branch, the stable node disappears through a saddle-node bifurcation, SN

_{1}, and the trajectory switches to the

*Z*-upper branch. Then, the trajectory exhibits, as

*z*increases, an oscillatory solution, which stabilizes at SLC

_{1}. Time series are plotted in Figure 45

*C*for SLC

_{1}and in Figure 45

*C*for LC.

_{1}
*(IX) Coexistence of LC and a NS:* The *z*-nullcline intersects the (*z*, *x*_{1}) curve at one equilibrium point which is a stable node, and intersects the
-curve at two points: LC and S. The stable node is the equilibrium point of the normal state. Then, LC and the normal state coexist and are separated by S. Time series are plotted in Figure 49*B*) for the normal state and in Figure 49*B _{1}
* for LC. The fast-slow subsystem switches to NS according to three scenarios, depending on

*m*. We discuss these scenarios below in area X.

In the following areas, LC and S disappear through a SNPO bifurcation.

*(III) DB:* The *z*-nullcline intersects the (*z*, *x*_{1}) curve at one equilibrium point, which is a stable node (Fig. 43*A*). The range of *m* indicates that a Hopf bifurcation point *H* exists, with
(Fig. 42). When a trajectory is at the *Z*-lower branch, the stable node disappears through a saddle-node bifurcation, SN_{1}, and the trajectory switches to the *Z*-upper branch exhibiting DB and continues to the equilibrium point, at which *z* stabilizes. Therefore, the fast-slow subsystem remains in the depolarization block.

*(IV) Periodic switch between DB and NS:* This area exhibits a periodic switch between DB and NS, like area V. In contrast, LC and S disappear through a SNPO bifurcation.

Therefore, the fast-slow subsystem only switches between DB and NS.

*(X) NS:* The *z*-nullcline intersects the (*z*, *x*_{1}) curve at one equilibrium point, which is a stable node. Time series of this solution are plotted in Figure 45, *E* and *E1*. Trajectories switch to the normal state according to three scenarios, depending on *m*.

For

*m*< 1 (first scenario), a Hopf bifurcation point,*H*, exists. If a trajectory is at the*Z*-upper sub-branch, which consists of unstable foci, it exhibits an oscillatory solution that terminates as*z*increases in a Hopf bifurcation,*H*. Then, the trajectory exhibits a nonoscillatory solution on the*Z*-upper sub-branch, which consists of stable foci. The nonoscillatory solution terminates as*z*decreases in a saddle-node bifurcation SN_{2}and the trajectory switches to the*Z*-lower branch.For

*m*= 1 (second scenario), a Hopf bifurcation point*H*exists. If a trajectory is at the*Z*-upper sub-branch, which consists of unstable foci, it exhibits an oscillatory solution, which touches the*Z*-middle branch before the Hopf bifurcation,*H*, as*z*increases. Then the trajectory switches to the*Z*-lower branch through an HB.For

*m*> 1 (third scenario), a Hopf bifurcation point,*H*, does not exist. If a trajectory is at the*Z*-upper branch, which consists of unstable foci, then it exhibits an oscillatory solution, which terminates as*z*increases in a homoclinic bifurcation, and the trajectory switches to the*Z*-lower branch.

The *Z*-lower branch consists of stable nodes, which are the equilibrium points of the normal state. At the *Z*-lower branch, the trajectory continues to a stable node, at which *z* stabilizes. Therefore, the fast-slow subsystem only remains in the normal state. Time series are plotted in Figure 45, *E* and *E1*, for the normal state.

Note that the SNPO curve (bold line) decreases as *m* increases.

#### Special cases

The fast-slow subsystem has one equilibrium point, which is a stable focus, an unstable focus, a saddle, or a stable node, depending on *m* and *x*_{0}. When LC coexists with an attractor and both are separated by S, the trajectories converge to LC or the attractor, depending on initial conditions. When S does not exist (*x*_{0} = −1.6; Fig. 47*B*,*b*), then all trajectories converge to LC (Fig. 45*B*,*B1*).

Figure 47*B*,*a*, shows that a stable focus and LC coexist, and S does not for *x*_{0} = −0.1 and *x*_{0} = −0.05. We plot trajectories in a (*z*, *x*_{1}) bifurcation diagram for *x*_{0} = −0.1 (Fig. 51*A*) and *x*_{0} = −0.05 (Fig. 51*B*). Even if S does not exist, trajectories converge to either LC or a stable focus, depending on the initial conditions.

### Subsystem 2

Here we present results of our analysis on the dynamics of subsystem 2 without coupling. We partition the results into two main parts: the equilibrium points and their stability, and bifurcations.

#### Subsystem 2 equilibrium points

We analytically determine the equilibrium points (*x*_{2}, *y*_{2}) by solving the following equations:
(67)

and (68)

We determine the stability of equilibrium points by evaluating the eigenvalues of the Jacobian matrix, *J*_{2}, as follows:
where
(69)

The trace
and the determinant
of *J*_{2} are given by:
(70)
(71)

The roots of are : (72)

The roots of are : (73)

We find the stability of equilibrium points in Table 6.

We graphically determine the equilibrium points (*x*_{2}, *y*_{2}) by using the *x*_{2}- and *y*_{2}-nullclines, where the *x*_{2}-nullcline corresponds to a cubic curve given by:
(74)and the *y*_{2}-nullcline corresponds to two straight lines given by:
(75)

We plot a phase plane of the subsystem 2 in Figure 52, for
(Figure 52*A*,*a*);
(Figure 52*A*,*b*); and
(Figure 52*A*,*c*).

Let
, the *x*_{2}- and *y*_{2}-nullclines collide in one point, which is a stable node (Fig. 52*A*,*a*). Increasing
, the *x*_{2}- and *y*_{2}-nullclines collide in three points: a stable node, a saddle, and an unstable focus (Fig. 52*A*,*b*). The stable node and saddle approach each other as
increases and coalesce at
forming an invariant circle (figure not shown). When
, the stable node and saddle disappear, and the invariant circle becomes a limit cycle through a SNIC bifurcation. Then, trajectories converge to the stable limit cycle (Fig. 52*A*,*c*).

We analytically determine the SNIC bifurcation point by solving the following equation: (76)

where *x*_{2} is a solution of
(77)

#### Subsystem 2 bifurcation diagram

The dynamics of subsystem 2 changes as
varies (Fig. 52*A*). We plot a (
, *x*_{2}) bifurcation diagram of the subsystem 2 in Figure 52*B*.
is a control parameter. The (
, *x*_{2}) curve is an *S*-shaped curve, which comprises four branches. *S*-lower and *S*-upper branches consist of stable nodes and unstable foci, respectively. The *S*-middle branch consists of saddles and acts as a separatrix between the *S*-lower and *S*-upper branches. The *S*-upper branch is divided into two sub-branches: a stable limit cycle surrounds one sub-branch (red) but does not surround the other sub-branch (green). We call the red sub-branch a limit cycle branch, which is limited by
- and
-curves. Increasing
, *S*-middle and *S*-upper branches collide in a saddle-node bifurcation, SN. Increasing further
, *S*-middle and *S*-lower branches collide in a SNIC bifurcation (
). The limit cycle branch terminates as
decreases in a SNIC bifurcation.

We identify three intervals on which the equilibrium points change:

]-

*∞*, SN[: only a stable node exists.]SN, SNIC[: a stable node, a saddle, and an unstable focus coexist. Here, the unstable focus is not surrounded by a stable limit cycle.

]SNIC,+

*∞*[: only an unstable focus exists, which is surrounded by a stable limit cycle.

The subsystem 2 is stable in the three intervals. Since the stable limit cycle emerges after a SNIC bifurcation, then there is no bistability in the subsystem 2.

We now link the dynamics of subsystem 2 to its behavior (see *ψ*_{2}; Fig. 1, main file) in the Epileptor model, which corresponds to transitions between an oscillatory state and a resting state. To explain these transitions, we plot time series *ψ*_{2} and equilibrium points of the subsystem 2 in Figure 53. The equilibrium points are labeled as dots. First, we observe that the bifurcation diagram of the subsystem 2 (Fig. 52*B*) is repeated over time. Blue, black, and green dots correspond to stable nodes, saddles, and unstable foci, respectively. Red dots, which are surrounded by stable limit cycles, correspond to unstable focus. A SNIC bifurcation occurs at the intersection of blue (stable nodes) and black (saddles) dots. Second, transitions from and to the limit cycle occur through a SNIC bifurcation, and under the evolution of
. In fact, the coupling function *g*(*x*_{1}) and *z* change the input
of the subsystem 2 to
(see Epileptor equations). We plot time series of *z* in Figure 53. Let
, the subsystem 2 is at the limit cycle branch exhibiting an oscillatory solution (Fig. 52*B*). When *z* increases, then
is decreased, and the oscillatory solution terminates in a SNIC bifurcation. Hence, the subsystem 2 exhibits a resting state. When *z* decreases, then
is increased, and the oscillatory solution emerges through a SNIC bifurcation. We conclude that the transitions between the oscillatory state and the resting state (see *ψ*_{2}; Fig. 1, main file) occur under the evolution of *g*(*x*_{1}) and *z*.

## Discussion

Epileptic seizures take many forms, as well as transitions from and to these pathological states. The mechanisms involved in epileptic dynamics are assumed to be different according to seizure patterns. We have used a neural mass model of partial seizures the “Epileptor” to analyze seizure dynamics. Using a mathematical approach, we showed that the model can generate different dynamic behaviors including SLEs, as the value of the parameters change. Most computational models tried to reproduce experimental data with biophysical realistic parameters (Kager et al., 2000; Traub et al., 2001; Destexhe, 2008; Cressman et al., 2009; Ullah et al., 2009). In contrast, parameters of the Epileptor model do not directly relate to biological variables. Their meaning can only be inferred based on dynamics. The SLE attractor consists of two states, an oscillatory state (limit cycle) interpreted as ictal activity, and a normal state interpreted as interictal activity. The ictal and normal states coexist and occupy two separated regions in the phase space. The SLE attractor corresponds to a periodic transition between the ictal and normal states. Transitions from ictal to normal states, and vice versa, can autonomously occur under the slow evolution of the variable *z* and without changing the parameters. The transitions from ictal to normal states occur through a bifurcation of the limit cycle, and the transitions from normal to ictal states occur through a bifurcation of the equilibrium point. The combination of two such bifurcations results in classes of (point-cycle) bursters as proposed by Izhikevich (2000). The Epileptor switches from normal to ictal states through a saddle-node bifurcation (Jirsa et al., 2014). In contrast, to return to the normal state, the Epileptor may undergo different bifurcation types, depending on the value of some selected parameters. There are two main parameters having effects on the Epileptor dynamics. The first, *m*, controls the Epileptor dynamics during the ictal period, and the second, *x*_{0}, controls the equilibrium points of the Epileptor. Bifurcations that define how the Epileptor qualitatively changes its behavior were identified using the scaling behavior of frequency and amplitude during seizure onset and offset. A taxonomy of seizures (SLEs) assembles the bifurcations into 16 different classes (Jirsa et al., 2014). SLE offsets in the Epileptor show the logarithmic scaling of homoclinic bifurcations that occurred more often (83%) in drug-resistant epileptic patients As a consequence, the fold/homoclinic bifurcation was identified as the predominant class of SLEs (Jirsa et al., 2014). Through a bifurcation diagram, we showed that the Epileptor undergoes a homoclinic bifurcation at seizure offset.

On the other hand, the behavior of SLEs in the 17% of patients with a non fold/homoclinic bifurcation was potentially consistent with two classes: fold/Hopf and fold/fold cycle (Jirsa et al., 2014). We evaluated whether other bifurcations may occur at seizure offset, when changing some parameters of the model. Upon change of the first main parameter, we indeed found that the Epileptor may undergo two different types of bifurcation at seizure offset: Hopf and SNIC bifurcations. Hence, the Epileptor can generate two classes: fold/Hopf and fold/circle. Our analysis demonstrated the existence of another burster, defined by Izhikevich (2000) as a point-point type (i.e., both bifurcations are of equilibrium points), called fold/fold bifurcation. This burster does not pertain to the 16 classes (point-cycle) defined by Izhikevich (2000), and then it does not describe any behavior of SLEs.

We proposed in this study a modification of the slow *z* variable equation (Eq. 34). This modification resolves the divergences observed in the model behavior for some initial conditions, in particular when the initial *z* value is negative. More, the modification unraveled the existence of a stable LC whose behavior results from the fast-slow structure of the Epileptor model. To explore this issue theoretically, we explained in this study how the modification of the *z* variable equation leads the model to generate a fast-slow cyclic behavior, in particular for negative initial *z* values, using bifurcation analysis and the averaging method. We demonstrated that LC does not destabilize, but disappears through a SNPO bifurcation. Interestingly, LC resembles the RSE and does not exist for certain initial conditions. As a consequence, if we consider LC as a prime candidate for RSE, our analysis of the LC dynamics would contribute to the understanding of the mechanisms underlying RSE and to the prediction of how to treat RSE.

In addition to SLEs and RSE, our analysis unraveled another pathological activity, called DB. Using bifurcation analysis, we explained how the Epileptor enters into DB. Indeed, we showed that reducing the first main parameter, the frequency of oscillations during the ictal period decreases until they disappear and DB occurs. In this case, the SLE attractor reduces to a periodic switch between DB and NS. The Epileptor model can therefore reproduce SLEs, RSE, and a periodic switch between DB and NS. Interestingly, it can remain in DB if we increase the value of the second main parameter *m*, and in NS if we decrease *m*.

We predicted that epileptic seizures, refractory status epilepticus, and “normal” brain activity can coexist in the brain, under some conditions. In fact, we demonstrated that the normal and ictal states of the SLE attractor coexist in the phase space, and, for certain conditions, LC also exists below. The transitions between them require a large change of *z*. The “fold/homoclinic bifurcation” was proposed as the predominant class of SLEs (Jirsa et al., 2014). We demonstrated that the coexistence of LC and SLEs occurs for the predominant class of the SLE attractor, as well as two other classes of the SLE attractor: fold/Hopf and fold/circle bifurcations. More, we demonstrated that only the SLEs with fold/circle type can exist alone in the phase space, for certain conditions.

The Epileptor model presents another type of bistability, which consists in the coexistence of LC and a periodic switch between DB and NS. Under certain conditions, this periodic switch between DB and NS can exist in the phase space without the presence of LC. Therefore, we could deduce that DB, normal brain activity and RSE coexist in the brain. However, epileptic seizures and DB do not coexist in the brain. In fact, the SLEs and periodic switch between DB and NS do not coexist in the phase space. DB was developed during the ictal period when reducing the first parameter value *m*.

We explored the parameter space of, first, the fast-slow subsystem and extend to the whole system Epileptor. For the Epileptor, we explore two parameter spaces, as the parameter control of the subsystem 2 is varied. In one of them, DB appears, when reducing the parameter value. The parameter space can provide pathways to switch among SLEs, RSE, DB, and NS. In fact, there is a region in the parameter space representing the coexistence of LC and DB in the phase space, another where only DB exists, and a third where only LC exists. Based on these parameters, we can then propose how to switch from SLEs to DB to RSE, and how to escape the RSE and DB. Interestingly, there is a region representing the coexistence of LC and NS in the phase space, and another where only the NS exists. Therefore, we could propose how to stay in the NS, and thus how to return to the normal brain activity.

In addition to the stable limit cycle LC, the Epileptor model generates another stable periodic orbit, which we denoted by SLC. This limit cycle has higher frequency and smaller amplitude than LC, and occupies the same region as the SLE attractor. SLC is different from LC; both coexist in the phase space. The transitions between them requires a large change of the *z* variable. Furthermore, there are two different patterns of SLC depending on the equilibrium point stability of the Epileptor, which can be a saddle as for the SLE attractor, or an unstable focus. SLC mimics the dynamics of epileptic seizures during the ictal period of the SLE attractor. For the SLE attractor, the Epileptor returns to the normal state after a transient epileptic seizures, while for SLC, epileptic seizures do not stop. Hence, we believe that SLC could characterize the status epilepticus, or might have another possible clinical and physiological explanation.

In conclusion, our study provides a better understanding of the Epileptor dynamics, exploring the different behaviors that the model can generate. Three classes can be generated by the Epileptor model, which are: fold/homoclinic, fold/Hopf and fold/circle bifurcations. We predict the coexistence of epileptic seizures (SLEs), RSE, DB, and normal brain activity (NS) in the brain, under a variety of conditions. Furthermore, the model predicts the existence of different paths between these neuroelectric activities which may help to explain the mechanisms underlying their genesis, making progress in clinical and brain research. (78) (79) (80)

## Footnotes

The authors declare no competing financial interests.

This research received financial support from the following agencies: Fondation pour la Recherche Médicale (Grant DIC20161236442 to Viktor K. Jirsa), European Commission’s Human Brain Project (Grant H2020-720270), and the SATT Sud-Est (TVB-Epilepsy). The work has been carried out within the FHU (Fédération Hospitalo-Universitaire) EPINEXT (Epilepsy and Disorders of Neuronal Excitability) with the support of the A*MIDEX project (ANR-11-IDEX-0001-02) funded by the “Investissements d'Avenir” French Governement program managed by the French National Research Agency (ANR).

- Received December 10, 2018.
- Revision received October 10, 2019.
- Accepted October 28, 2019.

- Copyright © 2020 El Houssaini et al.

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.