## Abstract

Brain regions generating seizures in patients with refractory partial epilepsy are referred to as the epileptogenic zone (EZ). During a seizure, paroxysmal activity is not restricted to the EZ, but may recruit other brain regions and propagate activity through large brain networks, which comprise brain regions that are not necessarily epileptogenic. The identification of the EZ is crucial for candidates for neurosurgery and requires unambiguous criteria that evaluate the degree of epileptogenicity of brain regions. To obtain such criteria and investigate the mechanisms of seizure recruitment and propagation, we develop a mathematical framework of coupled neural populations, which can interact via signaling through a slow permittivity variable. The permittivity variable captures effects evolving on slow timescales, including extracellular ionic concentrations and energy metabolism, with time delays of up to seconds as observed clinically. Our analyses provide a set of indices quantifying the degree of epileptogenicity and predict conditions, under which seizures propagate to nonepileptogenic brain regions, explaining the responses to intracerebral electric stimulation in epileptogenic and nonepileptogenic areas. In conjunction, our results provide guidance in the presurgical evaluation of epileptogenicity based on electrographic signatures in intracerebral electroencephalograms.

## Introduction

Partial seizures arise from a localized area or network, called the epileptogenic zone or epileptogenic network (EZ; Fig. 1*A*; Talairach and Bancaud, 1966; Bartolomei et al., 2001; Spencer, 2002). Partial seizures can be asymptomatic (Cook et al., 2013), and their spread to downstream brain regions is often linked to the emergence of clinical symptoms (Bonini et al., 2014), including cognitive impairment and loss of consciousness (Arthuis et al., 2009; Bartolomei and Naccache, 2011). How brain areas are recruited during seizure propagation is not well understood. Intracranial depth or stereotactic EEGs are commonly used to delineate the EZ in drug-resistant patient candidates for neurosurgery (Bartolomei et al., 2002). In clinical practice, direct stimulation of brain regions with intracranial electrodes is used to localize epileptogenic regions (Talairach and Bancaud, 1973; Munari et al., 1993) and assess their degree of epileptogenicity (Wieser et al., 1979; Bernier et al., 1990; Chauvel et al., 1993; Halgren and Chauvel, 1993; Inoue et al., 1999; Ure et al., 2009). Time delays of seizure recruitment have also been considered to be indicative for the strength of epileptogenicity, but remain controversial as there is a large degree of spatial and temporal variation of propagation, even within the same patient (Bartolomei et al., 2005). Seizure propagation from the epileptogenic toward neighboring zones has been observed experimentally in animal models (Trevelyan et al., 2006) and in humans (Schevon et al., 2012). More common clinical observations report seizure propagation across large-scale networks (Bartolomei et al., 2013) evolving on considerably slower timescales (seconds compared with milliseconds for neuronal firing).

Neural population models, also called neural mass models, have provided important insights into seizure mechanisms, in particular the nature of transitions between ictal and interictal states such as bifurcations or effects of noise (Wendling et al., 2002; Kalitzin et al., 2010; Touboul et al., 2011; Taylor et al., 2013). Extracellular factors (ions, neuromediators, tissue oxygenation, etc.) constitute good candidates to account for seizure evolution on slow timescales (Heinemann et al., 1986; Suh et al., 2006; Zhao et al., 2011) and have been studied theoretically in the context of seizure-like activity and spreading depression (Kager et al., 2000; Somjen et al., 2008, 2009). Alternative approaches rely on first principles in the theory of nonlinear dynamic systems (Kramer and Truccolo, 2012; Wang et al., 2012; Jirsa et al., 2014) and allow analysis of the mathematical description of the origins of seizures, but abandon the proximity to the biophysics. Here we will exploit one of these phenomenological models, the Epileptor (Jirsa et al., 2014), and investigate seizure recruitment when Epileptors are coupled via signaling through slow permittivity variables. To this end, we will derive a phenomenological representation of a homeostatic coupling. Our analyses will provide a systematic characterization of coordinated behaviors between two brain regions through the development of parameter spaces. These parameter spaces provide us with means to understand the different behaviors of seizure recruitment and explain some of its variability, including responses to intracerebral electric stimulation as performed in the clinic.

## Materials and Methods

##### Model.

Acknowledging the timescale separation present in seizure evolution, Jirsa et al. (2014) used the theory of fast–slow systems in nonlinear dynamics to develop a taxonomy of seizures and identified from experimental seizure data of various species a predominant class of bifurcation pairs, integrated into a phenomenological dynamic model called Epileptor. The Epileptor comprises a slow permittivity variable accounting for, presumably predominantly, extracellular effects related to energy consumption and tissue oxygenation and captures details of the autonomous slow evolution of interictal and ictal phases, as well as various details of seizure evolution during each phase. Notably in the context of this article, the slow permittivity variable provides a definition linked to gradual degrees of epileptogenicity. The Epileptor is a five-dimensional model and comprises three different timescales accounting for various electrographic patterns (Fig. 2*A*): on the fastest timescale, two state variables (ensemble 1; variables *x*_{1} and *y*_{1} in Eq. 1) exhibit bistable dynamics between oscillatory activity modeling fast discharges and a stable node representing interictal activity. On the intermediate timescale, two state variables (ensemble 2; *x*_{2} and *y*_{2} in Eq. 1) model the spike and wave events (SWE) and form the second neuronal ensemble (state variable *x*_{2} in Fig. 2*B*). On the slowest timescale (order of tens of seconds), the evolution of a very slow permittivity variable (variable *z* in Eq. 1) guides the neural population through the seizures including seizure onset and offset (Fig. 2*B*, state variable *z*). The first ensemble is linearly inhibited by the second ensemble in order for fast discharges to occur only during the wave part of the SWE (via *f*_{1}(*x*_{1}, *x*_{2})), the second ensemble is excited by the first ensemble through a low-pass filter coupling to generate SWE and interictal spikes (via *g*(*x*_{1})). Both ensembles are coupled through the permittivity variable, whereas the first ensemble acts directly upon the permittivity variable (via *h*(*x*_{1})) and thus closes the interaction loop of the circuit. A separatrix (i.e., a boundary between two regions of qualitatively different behaviors) divides the state space of the first ensemble between ictal and interictal states and acts as a barrier. As the permittivity variable evolves over time, seizure onset occurs through a saddle-node bifurcation showing a DC shift at the transition between the interictal and the ictal state (Fig. 2*B*, state variable *x*_{1}). The DC shift has been recorded *in vivo* and *in vitro* (Ikeda et al., 1999; Vanhatalo et al., 2003; Jirsa et al., 2014). Seizure offset occurs through a homoclinic bifurcation showing the logarithmic scaling of interspike intervals when approaching seizure offset. The time course of the local field potential is related to the total activity, *x*_{1} + *x*_{2} (Fig. 2*C*).

The permittivity variable absorbs all intracellular and extracellular effects linked to homeostatic regulation of the cell such as extracellular levels of ions (Heinemann et al., 1986), tissue oxygenation (Suh et al., 2006), and metabolism (Zhao et al., 2011). Importantly the nature of these biological variables can vary, only their dynamics and respective relationship is enforced by the Epileptor. In particular, seizures can be triggered via multiple routes such as synaptic noise or stimulation (Jirsa et al., 2014).

The duration of the ictal and interictal state is approximately the same in the Epileptor, which is not the case in clinical situations, where the interictal period is typically longer. For the upcoming quantitative discussion of ictal and interictal periods in the context of the Epileptogenicity Index (EI), a better quantitative treatment necessitates a modification of the Epileptor equations, without altering the actual bifurcation structure. To this end, we use a slight modification of the Epileptor equations and replace the linear influence of *x*_{1} on the permittivity variable by the nonlinear function *h* causing a symmetry breaking between ictal and interictal period with an increase of the latter. We do not consider any further quantitative effects of this symmetry breaking, but wish to point out that realistic seizures in the human do not follow a simple periodic pattern as observed in some experimental preparations (Jirsa et al., 2014). In some cases interictal distributions have been shown to follow Poisson or gamma distributions for human seizures (Milton et al., 1987; Shayegh et al., 2014). Models of continuous random variables that show threshold behavior can generally generate Poisson distributions. For a window of sufficiently constant permittivity, these constraints are satisfied in the Epileptor model and seizures and interictal events will be Poisson distributed if other conditions such as bistability are fulfilled. The discharges of the second population *x*_{2} have a very small influence on the permittivity variable in the simulations and have been neglected here for simplicity.

The Epileptor equations with the here used modification read then as follows:
where
and the degree of epileptogenicity is represented through the value *x*_{0}. For *x*_{0} < 2.91, the Epileptor shows seizure activity autonomously and is said to be epileptogenic; for values of *x*_{0} > 2.91 the Epileptor is in its (healthy) equilibrium state. To simulate the system of differential equations above, we used a Euler–Maruyama integration scheme with an integration step of 0.05. For stochastic simulations we used additive white Gaussian noise in the variables *x*_{2} and *y*_{2} with mean 0 and variance 0.0025; 256 time steps are equivalent to 1 s of real time to obtain realistic frequencies of oscillation and seizure lengths and to match intracranial EEG sampling frequency. When specified in the Results section, we used a bandpass Butterworth filter of order 5, with cutoff frequencies of 0.16 and 97 Hz at −3 dB (Fig. 2*D*). The filter was chosen to be equivalent to the filter used to process intracranial EEG signals (Bartolomei et al. (2008) for processing methods used with intracranial EEG signals).

##### Coupling.

Arguments based on timescale separation in the theory of nonlinear dynamic systems (Haken, 1987) suggest that fast couplings through synapses or gap junctions will not qualitatively change the behavior on the slow timescales (see Discussion). For this reason multiscale seizure recruitment necessitates the inclusion of slower homeostatic mechanisms affecting the permittivity of cell tissue. Coupling through permittivity variables then has to model the extracellular/intracellular effects of local and distant discharges.

In the following we take inspiration from these approaches and introduce nonlocal interareal permittive coupling between Epileptors. The Epileptor, formally a phenomenological model of epileptiform neural activity, does not link directly to such described biophysical mechanisms. Instead we approximate the effect of fast neuronal discharges as a homeostatic perturbation of the slow permittive variable *z* of the Epileptor in the target region. On the slow timescale the local discharge influences upon the permittivity are captured through the expression *h* in the Epileptor. The distant discharges at remote locations *j* are transmitted via axons to target region *i* and perturb there the homeostatic equilibrium. This disturbance of ion homeostasis will depend on both local and distant neuronal discharges and can formally be expressed as a coupling term *C _{ij}*(x

_{1,j}, x

_{1,i}) influencing the rate of change of the permittivity variable

*z*

_{1}of region

*i*for

*i*,

*j*∈[1,

*N*] × [1,

*N*] where

*N*is the number of regions of interest. Figure 1

*B*shows this architecture where the permittivity is represented by the extracellular environment [which has been mostly modeled as extracellular potassium concentration so far (Fröhlich et al., 2008), but here it will not be limited to it]. The permittive coupling term

*C*(x

_{ij}_{1,j}, x

_{1,i}) subsumes all disturbances of ion regulation. We consider relevant working points

*x*

_{i}

^{0}and

*x*

_{j}

^{0}(corresponding locally to the mean of the ictal or interictal states) and their respective perturbations ξ

_{i(t)}, ξ

_{j(t)}. Then we express the coupling term as

*C*(

_{ij}*x*

_{i}

^{0}+ ξ

*,*

_{i}*x*

_{j}

^{0}+ ξ

*), with |ξ*

_{j}*|≪1 and |ξ*

_{i}*|≪1 indicating that the perturbations are small. We locally expand the permittive coupling around the working points in a Taylor series and consider its first order effects, i.e.,*

_{j}*C*(

_{ij}*x*

_{1,j},

*x*

_{1,i}) =

*C*(

_{ij}*x*

_{i}

^{0},

*x*

_{j}

^{0}) +

*and ξ*

_{i}*by*

_{j}*x*−

_{i}*x*

_{i}

^{0}and

*x*−

_{j}*x*

_{j}

^{0}, respectively, and reshuffle the subsequent terms such that the coupling is expressed by

*x*

_{⊥}=

*x*−

_{i}*x*, which now characterizes the coupling through the deviations orthogonal to the synchronization manifold

_{j}*x*=

_{i}*x*=

_{j}*x*(Dhamala et al., 2004) and hence is well suited to address seizure recruitment effects. All other modifications via linear influences

*x*

_{1}and constant parameters

*C*(x

_{ij}_{i}

^{0}, x

_{j}

^{0}),

*x*

_{i}

^{0}and

*x*

_{j}

^{0}are absorbed in the local permittivity term

*h*of the original Epileptor equation. Under these considerations, the final model equations for

*N*-coupled Epileptors with permittive coupling read as follows: with

*f*

_{1},

*f*

_{2}, g and

*h*as in Equations 2–5, respectively. There is no need to consider delays via signal transmission in the context of seizure recruitment, because the dynamics of

*z*is sufficiently slow to render these delays negligible.

##### 2D reduction of coupled Epileptors.

To study the recruitment effects on the slowest timescale of the permittivity, we apply averaging methods to the coupled Epileptor equations, which reduce the dimension of our model. Numerical analyses of the full system demonstrate that the effect of the second neuronal ensemble is negligible under this approximation, hence we neglect it in the analytical discussion also for simplicity, but verify all results against the full system numerically. Then the Epileptor equations (Eq. 6) become
with
The timescale difference is guaranteed by τ_{0}≫1. Then the fast variables, *x*_{1,i} and *y*_{1,i}, rapidly collapse on the slow manifold whose dynamics is governed by *z _{i}*. Since we are interested in the dynamics on this slow manifold, we can rescale the time by

*T*=

*t*/τ

_{0}and study the limit, which results, in τ

_{0}→ + ∞,

*y*

_{1,i}−

*f*

_{1}(

*x*

_{1,i}) −

*z*+

_{i}*I*

_{1,i}= 0 and

*y*

_{0}− 5

*x*

_{1,i}

^{2}−

*y*

_{1,i}= 0. Averaging over the fast oscillations, we keep

*y*

_{1,i}−

*f*

_{1}(

*x*

_{1,i}) −

*z*+

_{i}*I*

_{1,i}≠ 0, use

*y*

_{1,i}=

*y*

_{0}− 5

*x*

_{1,i}

^{2}, and obtain the following: with On the slow timescale, the dynamics of the system is nearly identical if we replace Equation 7 by the simpler expression

*f*

_{1}(

*x*

_{1,i}) =

*x*

_{1,i}

^{3}+ 2

*x*

_{1,i}

^{2}for all

*x*

_{1,i}∈

*R*. Then the neural mass model equations read as follows: with

*h*(

*x*

_{1,i}) =

*x*

_{0,i}+ 3/

_{0}= 2857 and

*I*

_{1,i}= 3.1.

##### Trigonometric reduction of slow Epileptor dynamics.

To get further insight into the bifurcations involved on the slow timescale, we derive a model equivalent to the uncoupled system in Equation 1 whose phase flow is constrained to the unit circle (the phase flow is shown in Fig. 3). We compute the term in *x*_{0,i} and the coupling term in *K* explicitly from Equation 6 by performing a circular approximation using *x _{i}* = cos(φ

*) and*

_{i}*z*= sin(φ

_{i}*) (the validity of these approximations is discussed in Roy et al., 2011). We obtain the following: with*

_{i}*j*the coupled Epileptor, τ

_{0}= 1000,

*x*

_{0,}

*− 2.1 taking values between 0 and 2 so*

_{i}*x*

_{0,}

*, the excitability parameter, takes values between 2.1 and 4.1 as in our model (Eq. 6), and*

_{i}*K*, the coupling parameter, taking values between 0 and 10. When the rate of change of a state variable is zero (i.e., its time derivative is zero), then we call this a fixed point. For reasons of simplicity, we introduce here a π/4 shift to move the fixed point between π and 5π/4. For

*x*

_{0,}

*= 3 and*

_{i}*K*= 0, a Saddle-Node on Invariant Circle (SNIC) bifurcation occurs: according to the value of the excitability parameter

*x*

_{0,}

*, the reduced model can thus be epileptogenic or not. In the oscillatory regime, the reduced model switches autonomously between interictal and ictal states as a result of the dynamics in Equation 9. All results are computationally validated against simulations of the complete system (Eq. 6).*

_{i}##### Analytical derivation of timing of seizures.

We can analytically calculate the length of the seizures and the length of the interictal periods by evaluating the integrals:
Under the assumption that cos(φ* _{j}*) − cos(φ

*) =*

_{i}*U*is constant, the strength of the coupling is constant over the intervals [0, π/4] and [π, 5π/4] compared with the DC shift, i.e., we assume cos(3π/4) − cos(0)≪cos(π) − cos(3π/4) and cos(5π/4) − cos(π)≪cos(2π) − cos(5π/4). Using the change of variable tan(ϕ/2), the integrals can then be analytically evaluated: with A = (

*x*

_{0,i}− 2)/

*K*·

*U*.

##### Analytical derivation of parameter space.

We summarize here the methods that entered into the derivation of the results in Figure 9. To determine the shape of the border between regime IV and III and between regime V and II, we used the 2D reduction (Eq. 8). On the border between regime IV and III, the second Epileptor, not epileptogenic (*x*_{0,2} < 2.91), becomes epileptogenic, equivalent to an Epileptor with *x*_{0,2} = 2.91, which we express as *x*_{0,2} − *K* · (x_{1,1} −x_{1,2}) = 2.91. *x*_{1,2} is the leftmost intersection of the nullclines (i.e., curves where the time derivative is zero for one state variable), solution of the equation − (*x*_{1}^{s})^{3} − 2(*x*_{1}^{s})^{2} + 4.1 − *x*_{0,2} = 0, using that when *x*→ − ∞, *z _{i}* =

*x*

_{0,i}on the sigmoid nullcline, computed numerically. As

*x*

_{1,1}varies in time, we take the temporal average for

*x*

_{1,1}. In the (

*K*,

*x*

_{0,2}) space, we find an almost linear relationship between

*x*

_{0,2}and

*K*(Fig. 9

*B*, red). Similarly, between regime V and II, the first Epileptor becomes nonepileptogenic if

*x*

_{0,1}− K · (x

_{1,2}−x

_{1,1}) = 2.91. On the border between regime V and II,

*x*

_{1,1}= −4/3 and

*x*

_{1,2}is the

*x*coordinate of the stable fixed point of the second Epileptor when its excitability is equal to

*x*

_{0,2}(solution of the equation − (

*x*

_{1}

^{s})

^{3}− 2(

*x*

_{1}

^{s})

^{2}+ 4.1 −

*x*

_{0,2}= 0, computed numerically; Fig. 9

*B*, cyan).

For the border between regime I and II, we used the trigonometric model; on the border, the time difference between intrinsic frequencies of both Epileptors is exactly equal to the time difference between the coupled and the uncoupled Epileptor:
*T*_{crise,1,K≠0} − *T*_{crise,1,K=0} ≈ 0, which gives
We numerically solve this equation using Equation 10 for different values of *K* and found the blue and green curve in Figure 9*B*.

##### Epileptogenicity Index (EI).

Bartolomei et al. (2008) suggested computing a ratio of spectral content and time delay of electrographic signatures as the EI. Specifically, they considered the energy ratio averaged over time to detect rapid discharges, divided by the time delay of seizure onset in the secondary region with respect to the primary seizure onset zone. We make two adaptations of the EI for the current Epileptor architecture: we simplify the definition of the energy ratio of high- to low-frequency discharges *E _{low}* = 1 for the Epileptor model to ease detection of faster frequencies and we identify the time delay between seizure onsets in the two regions via onset of their depolarization shifts.

##### Stimulation.

Stimulation was performed through the parameter *I*_{1} in Equation 1 (rectangular pulse of 1 unit strength, duration 300 ms) during the interictal period in. A seizure-like event was considered to be triggered if a seizure started <800 ms after the stimulation.

## Results

### Regimes of recruitment

To gain an understanding of the dynamic mechanisms underlying seizure recruitment, we focus on two coupled Epileptors and derive their characteristic behaviors. We systematically varied the parameters for the coupling strength *K* and the degree of excitability *x*_{0,}* _{i}* in the case of two Epileptors with

*i*∈ {1, 2} and simulated their coordinated dynamics. A recording of a typical empirical intracranial EEG time series is shown in Figure 4

*A*. A typical example of a simulated time series is shown in Figure 4

*B*. In a given patient, the spatial extent of a focal seizure can be variable (Fig. 4

*A*). To reproduce this feature, we explored systematically the activity displayed by the second Epileptor for different values of the coupling strength

*K*, excitability

*x*

_{0,1}of the first Epileptor and excitability

*x*

_{0,2}of the second Epileptor. If not mentioned otherwise, the first Epileptor was always chosen epileptogenic (

*x*

_{0,1}= 2.5). For each parameter value, time series were simulated over at least 10 seizures and the relevant regimes identified. Two seizures were considered coexistent if a seizure of one Epileptor started before the end of the seizure of the other Epileptor. This notion of seizure coexistence is distinct from the synchronization of discharges during a seizure, which is not considered here due to our focus on seizure recruitment.

We observed five different regimes of seizure coexistence (Fig. 5*A*,*B*). In regime I the time difference between seizure occurrence in the first and the second Epileptor is shifting (phase drifting). Regime I happens for low values of coupling, where each Epileptor is epileptogenic but with different values of excitability (*x*_{0,1} ≠ *x*_{0,2}), thus different intrinsic frequencies. With such low coupling strength, Epileptors do not affect each other and are independent. Regime II is an entrainment regime, where the leading (most epileptogenic) Epileptor rapidly recruits the following Epileptor (epileptogenic or not). Regime III shows a p:q phase-locking regime, i.e., p seizures of the first Epileptor correspond to q seizures of the second Epileptor. Conversely to regime II, the weak coupling term hinders the systematic recruitment of the second Epileptor before the first Epileptor ends its seizure. In regime IV, the first Epileptor still expresses epileptiform activity, but the second Epileptor (here not epileptogenic) is not recruited, because the coupling term is too weak for the second Epileptor to enter an ictal/interictal cycle. Finally, in regime V, no epileptiform activity is apparent: the second Epileptor inhibits the first Epileptor by exerting an “inhibitory coupling,” which blocks seizures. For higher values of *K*, the boundary between regime II and V remains horizontal; for higher values of *x*_{0,2}, regime V rapidly becomes the only regime available (data not shown).

We confirmed systematically that all regimes exist for other values of *x*_{0,1} (Fig. 5*C*). The boundaries between the different regimes change quantitatively according to the value of the parameter *x*_{0,1} but stay qualitatively the same speaking for the robustness of the topology of the parameter space. As *x*_{0,1} increases, regime V becomes larger and moves upward, because the first Epileptor is closer to the nonepileptogenic regime and easier to inhibit. The difference between intrinsic frequencies of both Epileptors will also change, increasing regime I and reducing regime III.

### Indices of epileptogenicity: delays of recruitment, EI, and seizure lengths

We define the delay of recruitment as the time between the onset of the seizure in a particular area and the onset of the seizure in another recruited area. We computed for several seizures the mean and SD of the delays of recruitment between seizure onset for each value of *x*_{0,2} and *K* (Fig. 6). In the entrainment regime (Fig. 5, regime II), the delay of recruitment is inversely proportional to the excitability of the second Epileptor and to the coupling strength (Fig. 6*B1*,*C1*). In the phase drift and the p:q coupling regime (Fig. 5, regime I and III, respectively), the SD takes substantial values because the delay changes from seizure to seizure (Fig. 6*B2*,*C2*). Then delays of recruitment can be either highly variable or stereotypical according to the values of excitability and coupling parameters.

As an additional measure we computed the EI for the second Epileptor according to the position in the parameter space (Fig. 6*D*). The EI decreases as the degree of excitability of the second Epileptor decreases and increases as the coupling strength increases. This result is consistent, as the EI increases when the second Epileptor is prone to be rapidly recruited in a seizure. A small region for very low values of coupling shows sparse high EI. This is an artifact of regime I: both Epileptors are independent from each other, resulting occasionally in very short delays between seizures, which artificially increases EI. Clinically, interictal periods are several orders of magnitude longer than in the model and such a situation is improbable.

We also quantified the length of the seizures for both Epileptors in the same parameter space (Fig. 6*E*). The length of the seizures of both Epileptors decreases as a function of decreasing the excitability of the second Epileptor, except for very low values of coupling where the first Epileptor is not affected by the change in the second Epileptor. The length of seizures appears to be high for high values of EI.

In the two previous sections, we showed that the excitability of downstream structures and the connectivity strength between upstream and downstream structures are key determinants of seizure recruitment. Since both parameters can vary considerably on a daily basis, coupled Epileptors provide a conceptual framework to study why and how seizures spread. Different values of the excitability and coupling parameters constrain the recruitment of downstream regions from seizure to seizure, which can be either stereotyped (e.g., in regime II) or variable (e.g., in regime III), both spatially (region recruited or not) and temporally (delay of recruitment).

### Stimulation

We explored the effect of stimulation in the Epileptor network for different values of coupling strength *K* and excitability *x*_{0,2}. As reported before (for review, see Jirsa et al. (2014)), there is a refractory period after a seizure, but a sufficient large stimulus can trigger a seizure (Fig. 7*A*). When Epileptors are coupled, stimulation of the first (Fig. 7*B*) or the second Epileptor (Fig. 7*C*) does not always trigger a seizure (if the stimulus is not strong enough to drive the system to the separatrix), and when a seizure is triggered, it does not necessarily recruit the second Epileptor, depending of the proximity of the second Epileptor to the separatrix, which is given by the value of the second permittivity variable *z*_{2}. We thus explored systematically the influence of the coupling strength *K* and current state of the permittivity variable *z*_{2} by stimulating the first Epileptor at random instants and found three regimes of different probability of recruitment: (1) never, (2) often, (3) always (Fig. 7*D*). The ability of a seizure triggered by stimulation to recruit a second Epileptor thus depends on the state of this second Epileptor.

### Analytical reduction

To get better analytic insights into the different regimes of the parameter space, we discuss the dynamics of the reduced system (Eq. 8) and validate the numerical simulations against the complete system (Eq. 6). Both systems show good correspondence (Fig. 8*A*), with slight differences in the intrinsic frequencies. The analysis of the trajectories of the reduced isolated system offers an intuitive understanding (Fig. 8*B*). We plot two nullclines that are defined as zero flux in either the *x* or *z* direction: the sigmoid nullcline for different values of *x*_{0} in red and the cubic nullcline in blue. The intersections of the nullclines identify the fixed points of the system. The interictal state and the ictal state correspond to the left and the right branches of the cubic nullcline, respectively. According to the value of the parameter *x*_{0}, the sigmoid nullcline in Figure 8*B* moves up and down, changing the number and stability of the fixed points. Two typical trajectories are possible. For *x*_{0} > 2.91, the trajectory is attracted to the stable fixed point in the interictal state on the left cubic nullcline and the Epileptor is said not epileptogenic, meaning not triggering epileptic seizure without external input. For *x*_{0} < 2.91, the two leftmost fixed points disappear through a SNIC bifurcation, and there is only one unstable fixed point left. The Epileptor enters an oscillatory regime and the Epileptor is said to be epileptogenic, triggering seizures autonomously; a trajectory is shown in green in Figure 8*B*. It follows a limit cycle alternating slow build up on the outer branches of the cubic nullcline and fast transitions between these nullclines, switching between the interictal and ictal state.

The reduced 2D representation in state space allows an intuitive understanding of the phenomena observed in the various parameter regimes. For instance, what is the effect of an epileptogenic Epileptor on a nonepileptogenic Epileptor? When a seizure occurs in the first Epileptor, it moves down the sigmoid nullcline of the second Epileptor whose fixed point disappears through a bifurcation (a Hopf bifurcation for this particular coupling function) and a seizure is triggered (Fig. 8*D*). If the second Epileptor is too far from the separatrix, several seizures of the first Epileptor can be necessary to trigger a seizure in the second Epileptor (Fig. 8*E*). If the second Epileptor is already epileptogenic, the occurrence of the seizure is accelerated (data not shown). Conversely, if one Epileptor ends its seizure while the other is still seizing, the former will exercise an inhibitory force on the latter. Increasing the coupling may also suppress seizures in the first Epileptor (regime V): due to the difference of their positions in the interictal state, as the coupling increases, the sigmoidal nullcline of the first (second) Epileptor goes down (up) respectively. Regime V will appear if the first Epileptor's nullcline crosses the bifurcation.

Note that in the oscillatory regime of the single Epileptor, cubic and sigmoid nullclines are close to the saddle-node bifurcation, resulting in a slow trajectory through a bottleneck, also called a ghost attractor (Fig. 8*D*; Deco and Jirsa (2012)). As we get closer to the SNIC bifurcation by moving the sigmoid nullcline up, the trajectories slow down in this bottleneck, resulting in an interictal period becoming longer than ictal periods. This bottleneck is due to the nonlinearity of the sigmoid nullcline (as modified from the Epileptor model in Jirsa et al. (2014) whose nullclines are shown in Fig. 8*C*). The effect of the proximity to the bifurcation can be quantified analytically by further reducing the 2D model to a 1D trigonometric model (Eq. 9). Using this trigonometric model, we compute the length of the seizures and the interictal periods analytically for *K* = 0 (Eq. 10) and compare the analytical results with those obtained from an Epileptor simulations using Equations 1–5 (Fig. 9*A*). We plot both results as a function of the epileptogenicity parameter *x*_{0} and find an excellent correspondence. In particular, as known from SNIC bifurcations (Jirsa et al. (2014)), the period length scales following a square root law.

Finally, we explain the shapes of the borders between the different regimes of the parameter space if we use the above reductions (see Materials and Methods). The analytically derived borders are shown in different colors in Figure 9*B* superimposed on the computational border (black) found previously.

In summary, two factors shape the topology of the parameter space, and thus the qualitative behaviors of coupled Epileptor networks: (1) the relaxation oscillator structure of the model and its two timescales and (2) the excitatory or inhibitory influence of the coupling function on the model stability.

## Discussion

Epileptogenic brain areas recruit other brain areas outside of the EZ during a seizure on a large range of time scales. We highlighted the role of permittivity coupling on a slow timescale in seizure recruitment and identified different regimes of seizure recruitment. We provided a set of indices including the delay of recruitment, EI, seizure lengths, and stimulation responses as indicators of tissue epileptogenicity. The ensemble of indices predicts conditions, under which seizures propagate to nonepileptogenic brain regions, explaining also the responses to intracerebral electric stimulation in epileptogenic and nonepileptogenic areas. Crucial to these findings is the timescale separation between fast seizure discharges and the slow evolution of permittivity. Approximations taking advantage of this timescale separation (Haken, 1987) allow us to lump system variables into subsets of slow and fast variables using adiabatic elimination. Prime candidates for slowly evolving variables involved in the recruitment of neuronal activity are mechanisms acting on the permittivity of neural tissue including extracellular ion concentrations, tissue oxygenation, and ATP. Pure synaptic and electric coupling are not sufficient to explain long delays of tens of seconds in seizure recruitment. The slow coupling via permittivity represents the mechanisms involved in the spread of the seizure and does not enable mechanisms for synchronization of discharges on faster timescales. Hence, in these frequency ranges, the coherence of seizing brain regions will be generally low during seizures in absence of synaptic and electric couplings. High synchronization (as seen in particular in absence seizures) may, however, be achieved by considering other additional coupling mechanisms, acting on faster timescales.

Various mechanisms have been proposed to explain seizure propagation and recruitment on timescales ranging from fast (1–100 ms) to very slow (several hundreds of seconds; Blume et al., 2001). Slice studies demonstrated that GABA antagonists reduce inhibition and allow excitatory activity to build up into spontaneous seizure-like activity (Sofronov and Golovko, 1992; Bausch and McNamara, 2000; Pinto et al., 2005). The lack of inhibition allows for seizure-like activity with fast propagation speeds of up to 100 mm/s (Chervin et al., 1988; Wadman and Gutnick, 1993). Hall and Kuhlmann (2013) proposed that slow propagation speeds (0.1–10 mm/s) are achieved in low [Mg^{2+}] solutions via presynaptic effects on excitatory (NMDA) receptors, which cause a decrease in the rate of presynaptic adaptation of excitatory synapses (DeBiasi et al., 1996), postsynaptic effects on excitatory (NMDA) receptors increasing postsynaptic excitatory conductance (Nowak et al., 1984), and a decrease in the rate of presynaptic adaptation of inhibitory (GABA) synapses. Several other biophysical parameters are known to change slowly during or preceding seizures, including extracellular levels of ions (Heinemann et al., 1986), tissue oxygenation (Suh et al., 2006), and metabolism (Zhao et al., 2011). Seizures recorded *in vivo* in cats (Moody et al., 1974) and awake baboons (Pumain et al., 1985) are characterized by the same time-dependent evolution of extracellular [K^{+}].

Several of these proposed mechanisms of seizure recruitment have been tested in computational models, mostly with a particular focus on either propagation through continuous sheets of neural populations (Kramer et al., 2005, 2007; Ursino and La Cara, 2006; Kim et al., 2009; Hall and Kuhlmann, 2013) or seizure recruitment of single cell activity within a population (Miles et al., 1988; Golomb and Amitai, 1997; Compte et al., 2003; Bazhenov et al., 2008; Fröhlich et al., 2008). For either approach, it is not evident how much we can infer about seizure recruitment of downstream brain regions that are distant, and often the biophysics of the slow processes is not well known enough to allow for explicit modeling or simply too complex (Marder and Taylor, 2011), which is why we have chosen here to systematically investigate seizure recruitment using a more abstract dynamic model and coupling. Though of a phenomenological and abstract nature, the network model allows highlighting the repertoire of recruitment scenarios observable in clinical settings allowing the delineation of the EZ, which is crucial for successful resective surgery. To characterize the EZ, clinical practice uses various indices including the high-frequency content of discharges, the presence of pre-ictal spikes, and time delays between the seizure onsets of two regions. Direct interrogation of the epileptogenicity of a particular region comprises its stimulation to trigger seizures. So far these procedures lack a theoretical foundation and their validity is challenged by the variability encountered in seizures (Chauvel et al., 1993; Halgren and Chauvel 1993), even within the same patient (Bartolomei et al., 2005). Neuro-electric patterns following stimulation may vary, whereas the most frequently observed pattern is a gradual increase of spikes and spikes and waves frequency, with or without occurrence of low-voltage fast activity (Munari et al., 1993). Seizures appearing morphologically similar may or may not be preceded by pre-ictal spikes, and may or may not propagate to secondary areas. Further, there is a large heterogeneity of excitability of brain tissue, such as the hippocampus and premotor regions being highly excitable and thus rendering it difficult to compare stimulation results across brain regions. For instance, stimulation of the amygdala may cause epileptiform after-discharges in the hippocampus, but not seize itself. In lateral temporal lobe epilepsy, stimulation of the mesial limbic structures is at least as efficient as stimulation of the seizure onset region to trigger seizures (Chauvel et al., 1993). Obviously our results cannot explain the whole range of phenomena, but may provide a first theoretical step to a better understanding of these interaction patterns. The systematic exploration of parameters explains at least part of the range of the variability encountered in patients. The understanding of the topology of behaviors in the parameter space spanned by these variables provides a chart of how to manipulate biophysical parameters to alter network behavior. As a map of navigation in the parameter space, it also offers a guide to estimate the degree of epileptogenicity. We discuss potential uses of this map in the following.

The use of stimulation provides less evident results regarding its use as an indicator of epileptogenicity. Jirsa et al. (2014) predicted that the ability of stimulation to trigger a seizure depends upon the strength of the threshold separating ictal and nonictal state. The strength is quantified by the distance to the separatrix and shows hysteresis. Similarly, when brain regions are coupled through permittivity, the distance to the separatrix together with the coupling strength are the main determinants of stimulation effects (Fig. 7*D*). Then, counterintuitive phenomena are encountered. For instance, a nonepileptogenic area (we refer here to the epileptogenic as defined parametrically in our model) can be stimulated and produce a seizure-like event, because the distance to the separatrix is small (Fig. 7*C*; where the second Epileptor is not epileptogenic), even though it does not spontaneously shows seizure-like activity, whereas a stimulation of an epileptogenic brain region just after a seizure is likely to stay without effects. Stimulation is an excellent tool to explore epileptogenicity, because it is the only way to interfere with the epileptogenic network; as such, it may be considered as indicative of causality, but not conclusive per se, because in each case the network embedding and the timing of the stimulation define a context that needs to be taken into account, especially in clinical situations.

Another counterintuitive phenomenon is encountered when multiple brain regions are mutually coupled through permittivity—some of these regions are epileptogenic—then the overall effect may render the network behavior nonepileptogenic (Fig. 5, regime V). For instance, despite well defined, preoperative electroclinical syndromes of temporal lobe epilepsy, many patients relapse unexpectedly, either immediately or remotely from the time of surgery (Hennessy et al., 2000). Often these relapses cannot be understood solely on the basis of local arguments (such as missed lesions identified in postoperative MRI) and necessitate other explanatory approaches as the one proposed here on the basis of network effects. A lesioning of the coupling between the two areas will effectively alter the network behavior and cause part of the network to show epileptiform activity. Figure 10, *c*, arrow, captures this form of surgical intervention as follows: originally the network is prepared in regime V. The surgical intervention reduces the coupling strength toward very small (or zero) values and moves the network into regime IV in the parameter space. As a consequence, Epileptor 1 starts displaying seizure-like activity, whereas Epileptor 2 shows rest activity. Conversely, strengthening the connection between epileptogenic and nonepileptogenic neural masses can stop the seizures and demonstrates the importance of the network effect in seizure propagation.

From these examples, it is evident that the precise evaluation of the degree of epileptogenicity and seizure recruitment necessitates an ensemble of diverse metrics and indices. In Figure 10 we show a summary of the ensemble of indices of epileptogenicity that we have investigated here. The set of indices may not be complete, but it is indicative for each region and has the potential to act as a biomarker for tissue epileptogenicity. Changes of parameters in this space allow tracing out the effects of medication or surgical intervention. This approach does not predict biophysically the functional pathways of medication, but allows on a more abstract level (via changes of coupling and excitability) to predict the interference of interventions. In conjunction, our results will aid in developing new dynamically based approaches in the presurgical evaluation of epileptogenicity based on electrographic signatures in intracerebral EEG.

## Footnotes

This work is supported by INSERM and the James S. McDonnell Foundation.

The authors declare no competing financial interests.

- Correspondence should be addressed to either Timothée Proix or Viktor K. Jirsa, Institut de Neurosciences des Systèmes, Aix Marseille Université, 27, Boulevard Jean Moulin, 13005 Marseille, France. timothee.proix{at}etu.univ-amu.fr or viktor.jirsa{at}univ-amu.fr