Active probing to highlight approaching transitions to ictal states in coupled neural mass models

Abstract

Altered auditory steady-state responses (ASSRs) have been found in patients with schizophrenia [17], bipolar disorder [18], and temporal lobe epilepsy [19]. Abnormal responses to TMS [26], single pulse electrical stimulation (SPES) and/or cortico-cortical evoked potentials (CCEP) have been used to explore cerebral functional connectivity and assist in delineating epileptogenic cortex for resection [27–29]. In both cases, the choice of stimulation parameters is a challenging task; widely different effects may be elicited by varying location, intensity, timing, frequency, polarity, and waveform of the stimulation component [59,64]. [...]computational neuronal models allow for the evaluation of a wider range of stimulation parameters that would be cumbersome, if frankly infeasible, to do experimentally. Perturbation stimuli could highlight intrinsic excitability and/or synchrony changes that lead up to seizures, but this effect is limited to the neuronal population being stimulated.

Introduction

The identification of robust electrophysiological features indicating approaching transitions from interictal to ictal states is one of the main challenges in epilepsy research [1,2]. The extraction of markers heralding the imminence of a seizure would be an important step for elucidating the mechanisms of epilepsy and for the development of seizure-prediction methods. The latter would alleviate one of the greatest burdens afflicting epilepsy patients, which is the unpredictability of seizures [3]. Achieving this would be a breakthrough to improve patients’ quality of life, whether by serving as seizure advisory systems [4] or by closing the loop to enable time-specific interventions to prevent upcoming ictal events [5].

Although several articles have reported promising seizure-prediction methods since the 90s, several of these have been questioned after the impossibility of reproducing the original positive results [6–8]. Such discrepancies occurred primarily in methods that were originally optimized and tested in specific datasets but failed to achieve similar performances when confronted with a greater variety of long-term EEG recordings. Achieving a sufficiently robust and reliable seizure forecasting method has thus been shown to be a task that is far from trivial [9] and has been described as “the long and winding road” [6].

In addition to being robust to the high number of brain states (sleep, attention, cognition, etc.), the development of seizure-prediction methods should consider that ictogenesis mechanisms may be distinct for different brain structures and different pathologies [6,10]. This motivates the use of individualized methods [11–13]; algorithms should be trained or tuned for each patient in order to account for the innumerous sources of inter-subject variability, which include different epilepsy and seizure types. Other challenges to overcome in developing seizure prediction methods are the scarcity of publicly available long-term EEG databases (although this is changing with the increased use of implantable devices) and the difficulty in extracting features that consistently reflect increased seizure susceptibility–in general, there is a high level of abstraction from signal features in relation to the underlying mechanisms leading up to seizures [6,14]. A possible approach to overcome some of these problems involves the use of active probing or active observation paradigms [7,15]; by applying stimuli and evaluating responses, more information could be obtained about the state of a system, which would otherwise be inaccessible with its passive observation. Such perturbational approaches have been shown to provide diagnostic value for different conditions and have served as valuable tools throughout the history of neurology and neuroscience research [16]. Altered auditory steady-state responses (ASSRs) have been found in patients with schizophrenia [17], bipolar disorder [18], and temporal lobe epilepsy [19]. The complexity of responses to transcranial magnetic stimulation (TMS) can be used to discriminate levels of consciousness [20] and has shown to have significant diagnostic potential in genetic generalized epilepsy [21]. Visual stimulation, through a binocular rivalry paradigm, reveals slower responses in autism [22]. Tracking cortical excitability in patients with epilepsy has been done with TMS [23,24] and electrical stimuli [25]. Abnormal responses to TMS [26], single pulse electrical stimulation (SPES) and/or cortico-cortical evoked potentials (CCEP) have been used to explore cerebral functional connectivity and assist in delineating epileptogenic cortex for resection [27–29]. However, given the wide variety of epilepsy types and seizure models available, the use of similar interventional or perturbational approaches for seizure forecasting is still limited.

The first use of stimuli to “probe” transitions to ictal states was done by Kalitzin et al. [10], with the use of intermittent photic stimulation (IPS) to test the hypothesis that response changes to these stimuli would precede the onset of seizure activity in photosensitive patients. An increase of the relative phase clustering index (rPCI) was found, which would reflect the hyperexcitability of the underlying dynamical system and would indicate the presence of nonlinear dynamics. Auditory stimuli were used to reveal the pathophysiology of mTLE, where altered neural synchronization induced by monaural pure-tone stimulation [30] or by ASSR [19] may provide useful diagnostic information. ASSRs were also used to analyze ictal and interictal state dynamics in an audiogenic seizure model, demonstrating that responses in this model are enhanced and hypersynchronous [31] in addition to the impairment of auditory processing integration [32]. Electrical stimuli in the form of either cortical or DBS (Deep Brain Stimulation) have been used to monitor neuronal excitability and predict the occurrence of seizures in patients with epilepsy [7,25,33]. DBS for pre-ictal probing was also used in canines [34] and in a Pentilenotetrazole (PTZ) model in rats [35]. In the latter, probing effects were potentiated by previously pairing stimulation with ictal activity.

Seizures in isolated CA1 sections are preceded by the progressive increase of amplitude and duration of responses evoked by stimulation of Schaffer’s collaterals [36]. This increased sensitivity and delayed recovery to perturbations characterizes the loss of a systems’ stability and may be related to an approaching critical transition. It has been hypothesized that the onset of seizures would be ruled by this type of dynamics, sharing similarities with other types of complex phenomena and showing common early-warning signs such as critical slowing down, flickering, increased correlation, variability, and lag-1 autocorrelation (which is the estimate of the correlation of a signal with a 1-sample delayed version of itself) [37,38]. Contrasting evidence has been found both against [39–41] and in favor of this hypothesis [42–45]. Some questions regarding this might be elucidated with the use of perturbational approaches, which have been used to evaluate the existence and proximity to such tipping points in living dynamical systems [46] and could also be useful tools to explore the underlying dynamics of seizure onset in different models.

Neuronal computational models have been used extensively in the last decades, such as in trying to uncover the mechanisms responsible for brain rhythms [47], neuronal dynamics [48], coding [49], and pathologies like Epilepsy [50–54]. By defining sets of rules and equations, models can be built to express either detailed or emerging general neuronal dynamics, offering insights about underlying mechanisms and enabling the testing of hypotheses that would be otherwise difficult or impossible to do experimentally. Individualized whole-brain models can be tailored to provide patient-specific interventions, such as in presurgical evaluation [55], or to predict and control seizure propagation [56,57]. Simulated microstimulation is proposed in [58] as a tool for tracking surround excitability levels to enable the operation of targeted intervention devices, as well as predicting surgical outcomes.

Established computational models provide a useful way of evaluating the effects of stimulation in neuronal systems [59,60]. Specifically, seizure-generating models permit exploration of how stimulation may be used for controlling seizures [61] or for the detection of transitions to ictal states with the aforementioned probing approach [9,62,63]. In both cases, the choice of stimulation parameters is a challenging task; widely different effects may be elicited by varying location, intensity, timing, frequency, polarity, and waveform of the stimulation component [59,64]. Thus, computational neuronal models allow for the evaluation of a wider range of stimulation parameters that would be cumbersome, if frankly infeasible, to do experimentally.

Despite the potential of this type of approach, few studies using probing to assess interictal to ictal transitions in neuronal computational models are found in the literature. Suffczynski et al. [62] used the model developed by Wendling et al. [65] to show that changes in excitability preceding the transition to a seizure are not observed in the ongoing EEG activity, but can be inferred by analysis combined with stimulation paradigms. Simulations also predicted the effects of electrical extracellular, local, bipolar stimulation (ELBS), showing that appropriately tuned stimulation activated preferentially GABAergic interneurons [59]. Furthermore, responses to ELBS in the computational model and patients with focal epilepsy were shown in this study to reveal hyperexcitable brain regions. The application of an active probing paradigm enabled a better observation of the underlying synchrony of a system of coupled Kuramoto oscillators [9], but this was only possible with suitably dense electrode recordings. That is, when electrodes are too spaced in-between, compared to the synchrony area, probing does not have the intended effect. The use of perturbation stimuli was also used with the Epileptor model to highlight feature changes signaling parameter changes towards ictal activity [63].

Using two coupled populations of Wendling’s neural mass model, this work aims to explore different conditions in which an active probing approach is useful for highlighting the dynamical system changes occurring before seizure onsets. We expect that, as model parameters are shifted from normal towards ictal state, neuronal excitability and synchrony changes would be observable with the use of low frequency probing stimuli. Different model settings of coupled neural masses were evaluated, varying stimulation target and amplitude, as well as which model parameter is shifted to bring the system towards ictal states.

Results in this work corroborate previous literature evidence from computational and experimental data indicating that the impending risk of seizure states can be assessed using active probing approaches. Perturbation stimuli could highlight intrinsic excitability and/or synchrony changes that lead up to seizures, but this effect is limited to the neuronal population being stimulated. The observed feature changes such as increased variance and lag-1 autocorrelation are common hallmarks of the phenomenon of critical slowing down, which may be related to the dynamics underlying transitions to seizure states and can serve as early-warning signs of the imminence of such tipping points. And while stimuli did not elicit sustained discharges of spikes (thus indicating the lack of proictal effects), interictal-like discharges were elicited just before the transition to ictal activity under different model settings, indicating these could be robust markers for forecasting the imminence of seizures.

Materials and methods

CA1 neural mass model

Some neuronal computational models aim to describe in detail membrane potential and ion channel dynamics of specific neurons, such as the classic Hodgkin-Huxley model [66]. Networks of such detailed neurons interacting through synapses and gap junctions can be built, but their size and complexity are limited by both computational costs and dimension of their parameter space. Mesoscopic models may overcome some of these limitations by describing the mean activity of neuronal populations. Instead of individual variables for membrane potentials, action potentials (APs), or post-synaptic potentials (PSP) for each neuron, lumped descriptions of these attributes can be represented for specific subpopulations included in larger-scale networks [67], which are so far unfeasible to represent with detailed microscopic models. Pioneering works on such mean-field (or lumped-parameter) models include those from Wilson & Cowan [68], Freeman [69], Lopes da Silva et al. [70], and Wright & Liley [71]. Other models of this kind have been used to study spike-wave (SW) seizures [72] and their abatement using single-pulse stimulation [73].

This work is built upon the neural mass model developed by Wendling et al. [65], which consists of four coupled subpopulations and based on the global cellular organization of the CA1 subfield of the hippocampus; excitatory main cells (pyramidal cells), excitatory interneurons, inhibitory neurons with slow kinetics (O-LM neurons, with IPSCs mediated by dendritic synapses) and inhibitory neurons with fast kinetics (soma-projecting basket cells). In this work, an instance of this model is referred to as a neuronal population (or as a model subset), while specific neural masses (representing excitatory, slow inhibitory or fast inhibitory) are referred to as subpopulations. Interactions between these neuronal masses are represented in Fig 1, with excitatory interneurons being implicitly represented as the self-feedback loop from Pyramidal Cells. The average activity of each subpopulation is described by two main input-output relations, referred as “wave to pulse” and “pulse to wave” functions (i) The first function (wave to pulse) is an asymmetric sigmoid which outputs the average efferent AP density of a population as a function of its average PSP or membrane potential. (ii) The second is a 2nd order linear transfer function that relates the average PSP to the average density of incoming APs and can be either excitatory (EPSP) or inhibitory (IPSP). Impulse responses are represented by he(t) = Aae-at for excitatory synapses (EXC), hi(t) = Bbe-bt for slow dendritic inhibition (SDI) and hg(t) = Gge-gt for fast somatic inhibition (FSI).

[Figure omitted. See PDF.]

Fig 1. Modified neural mass model.

(a) Representation of each coupled neuronal population and corresponding subpopulations. Excitatory interneurons are implicitly represented as the self-feedback loop from main pyramidal cells (EXCi to itself). Each model settings consists of either (I) unilateral (s2) or (II) simultaneous bilateral low-frequency probing stimulation (s1 and s2), while slowly changing one parameter towards ictal activity–A1, B1 or K.

https://doi.org/10.1371/journal.pcbi.1008377.g001

The model receives an input signal pi(t), a normally distributed random variable and represents the density of afferent pulses or APs to the main excitatory (pyramidal) cells of neuronal population (or model subset) i. The model output is usually chosen as the average PSP of the main pyramidal cells, which presumably constitute the main generators of field potentials that make up the EEG signal.

Model parameters and modifications

The main adjustable parameters of Wendling’s neural mass model are usually A, B, and G. These refer to synaptic gains and influence the amplitude of the average excitatory (for A) or inhibitory (for B and G) PSPs. Changing these gives rise to at least six different activity types [65]; (1) normal background, (2) sporadic spikes, (3) sustained discharges of spikes, (4) slow rhythmic activity, (5) low-voltage fast activity and (6) slow quasi-sinusoidal activity.

In this work, two main modifications were made to the original model. The first is the bidirectional coupling of two neuronal populations (or model subsets) through excitatory main cells, as shown in Fig 1. This aims to represent inter-hemispheric hippocampal coupling, where each model subset or population represents the CA1 region of one cerebral hemisphere. This interaction is influenced by a coupling factor K (varied across simulations) and a delay coefficient τd, set to 10 ms [74].

The second modification to the original model is the addition of low-frequency periodic stimuli. These probing stimuli are represented by si(t) and summed to the model’s Gaussian noise input pi(t), which is the incoming density of action potentials to the main excitatory cells of the neuronal population i. That is, the simplified modeling of the probing stimuli consists of an additional input function that represents incoming AP bursts due to the electrical stimulation of afferent regions such as the dentate gyrus or Schaffer collaterals. In the model, this is defined as periodic pulses with duration of 10 ms and amplitude Estim, delivered every 2 seconds. This differs from the perturbation approach used previously by Suffczynski et al. [62], which simulated the placement of extracellular bipolar stimulation electrodes to change the membrane voltage of all subpopulations with pulses applied at a frequency of 10 Hz. We opted to use a slower stimulation frequency to ensure that the neuronal circuits are not entrained, preventing eventual pro-ictal effects due to the stimulus [75] (for a review on the pro or anticonvulsant effects of various stimulation parameters, see Ref. [76]).

The following set of 10 differential equations rules the activity of each neuronal population i.where i ∈ [1,2] and refers to the number of the current neuronal population, si is the probing stimulus applied to population i, and Ki,j is the coupling factor between populations i and j. In this work, bidirectional equal coupling was used, with K1,2 = K2,1 = K.

Values and descriptions of model parameters are given in Table 1.

[Figure omitted. See PDF.]

Table 1. Parameter descriptions and values [65].

A, B, and K values are either fixed or varied linearly throughout the simulations.

https://doi.org/10.1371/journal.pcbi.1008377.t001

In order to evaluate the effectiveness of probing stimuli to highlight preictal states, simulations are run with parameters set to normal activity, but linearly shifting one parameter towards the sustained discharges of spikes state throughout the simulation. This can be brought by simply increasing excitation (A) or decreasing inhibition (B) [65], but also by the time-varying interplay of interactions between excitatory and inhibitory synaptic gains [77]. In the modified model with coupled population subsets, the coupling factor K can also be used for seizure initiation. Thus, simulations were realized with default parameters as given in Table 1, but linearly varying A1, B1, or K in order to approach but not cross the threshold between normal and sustained spiking discharge activity.

The Euler-Maruyama method was used to solve the model’s stochastic differential equations (SDE’s), using a step size of 1.95 ms (sampling frequency of 512 Hz). The final time for each simulation is 2000 seconds. The model was implemented in Python 3.7, using the NumPy [78], Joblib, scikit-learn [79], statsmodels, pandas [80], Matplotlib [81] and SciPy [82] packages. All scripts are available at https://github.com/vrcarva/WNMM_probing.

Model settings

Different model settings were used to evaluate the effects of stimulation on predicting pre-ictal states, as illustrated in Fig 1. Each setting is characterized by the stimulation targets (one or both neuronal populations) and the parameter used to shift population activity towards ictal states–by increasing A1 (EXC gain of population 1), decreasing B1 (SDI gain of population 1), or increasing K (inter-population coupling gain). Detailed parameters of each setting are given next. All other parameters are set according to Table 1.

1. Probing neuronal population 2

1. A. Increasing A1: A1 = 2.5 to 4.6, A2 = 4.0, B = 40, K = 0.3

2. B. Decreasing B1: A = 4.0, B1 = 45 to 30, B2 = 40, K = 0.3

3. K. Increasing K: A = 4.0, B = 40, K = 0 to 0.5

2. Probing neuronal populations 1 and 2

1. A. Increasing A1: A1 = 2.5 to 4.6, A2 = 4.0, B = 40, K = 0.3

2. B. Decreasing B1: A = 4.0, B1 = 45 to 30, B2 = 40, K = 0.3

3. K. Increasing K: A = 4.0, B = 40, K = 0 to 0.5

Thus, model setting I consists of probing population 2 while set 1 parameters are shifted towards ictal-like activity. This aims to evaluate if applying stimuli to different regions that are not ictogenic (but coupled to ictogenic ones) allows the identification of impending ictal transitions. Unilateral probing of the ictogenic population 1 is evaluated in additional model settings, included in S1 and S2 Figs. The use of bilateral probing is evaluated in model setting II, where both populations receive stimuli simultaneously. Fig 1 illustrates the interactions between neuronal populations and subpopulations of the proposed model settings.

Evaluating probing effectiveness

After each simulation, different features were extracted from 400 ms post-stimuli epochs of each neuronal population’s output. The set of features is composed of variance, skewness, kurtosis and lag-1 autocorrelation (lag-1 AC) of each simulated LFP, and mutual information between the activity of both coupled populations to assess synchronization [83,84]. A moving average filter of order 20 was used for smoothing the resulting feature time series.

In each simulation, the dynamics of the extracted features may indicate if the system is approaching the transition to an ictal state. In addition to visual inspection, where increasing or decreasing feature trends could indicate increased seizure susceptibility, it is important to quantify how good predictor a feature is and compare this quantitative measure with and without probing stimulation. This can be done by assessing the relationship between the shifted parameter and the selected feature. The Spearman’s rank correlation coefficient [85] was used for this in order to account for non-linear correlations. It is defined by:where k is the paired score, x and y are the ranks of extracted features (variance, skewness, kurtosis, Lag-1 AC, and mutual information) and shifted parameters (A1, B1, or K), respectively.

To evaluate the sensitivity of probing efficiency as a function of stimulus intensity, input pulse amplitude is gradually increased from 0 (i.e., passive paradigm) to 200 APs/s, in steps of 20. In each test, fifteen realizations were run, yielding distributions for each extracted feature time-series, from which correlations measures were calculated. This enables statistical comparison between the effects of increasing probing stimuli amplitude–this was done with Tukey's honestly significant difference (HSD) method, accounting for pairwise comparisons of probing efficiency between increasing stimulation amplitudes and the passive observation case (stimulus amplitude = 0).

A crucial assessment that must be done is whether the stimuli have pro-ictal effects. That is, probing stimuli should help to predict impending seizures, not elicit them. Thus, additional simulations were done with the same model settings and varying degrees of stimulus amplitude, but shifting parameters to elicit ictal activity in the form of sustained discharges of spikes. For each simulation, seizure onset was defined as the first time in the simulation where spike discharges occurred in succession (more specifically, five spikes with an inter-spike interval of at least 4 seconds occur in succession). Spikes were detected if a peak was found in the simulated signal with amplitude greater than 5 mV. In order to exclude stimuli responses, peaks found up to 200 ms after perturbation onset were discarded. Once again, multiple comparisons testing with Tukey’s HSD were made to test if stimulation of increasing amplitudes resulted in significantly different seizure onset times in relation to the passive case with no stimuli.

Results

Simulated signals are shown in Fig 2, with and without stimulation of both neuronal populations (120 APs/sec). First, the parameter A1 was varied linearly from 4.3 to 4.95, eliciting ictal activity at the end of the simulation. As a seizure approaches, the shape of the evoked responses changes gradually, up to a point in which the response is drastically increased. These evoked responses can be considered as interictal like spikes (ILS) [86] evoked by the stimuli when the system is approaching an ictal-like state. They are akin to afterdischarges encountered during electrical stimulation. There was no evidence that a sustained series of spontaneous discharges (e.g. an ictal run) was generated. Only unitary events were observed.

[Figure omitted. See PDF.]

Fig 2. Neuronal mass model outputs with and without active probing stimuli.

Simulation changing parameter A1 from 4.30 to 4.95 (B = 40, G = 20) to elicit ictal activity in the end, with (a) passive observation or (b) active probing of both populations. Outputs are shown for population 1 (a.1 and b.1) and population 2 (a.2 and b.2), with tick marks showing stimuli onsets for (b). Ictal-like activity in the form of sustained discharges of spikes is observed for population 1 at the end of both simulations. ILS appear before this transition and are elicited even further in advance by stimulation. (c) Bottom plots show overlapped and mean response waveforms for increasing values of A1 for (c.1) population 1 and (c.2) population 2. Gradual response changes are observed until the appearance of ILS for A1 = 4.5 and ictal activity at A1 = 5.0.

https://doi.org/10.1371/journal.pcbi.1008377.g002

Resulting features and probing effectiveness

Probing was then evaluated with model setting I, in which a population is made ictogenic (but without crossing the threshold to ictal activity) by shifting a parameter linearly throughout each simulation–increasing A1 (I-A), decreasing B1 (I-B), or increasing K (I-K)–while perturbation stimuli are applied to the input of population 2. Fig 3 shows the resulting feature series for I-A for simulations without probing, followed by the Spearman's rank correlation measure between A1 and each extracted feature series for various degrees of stimulus intensity.

[Figure omitted. See PDF.]

Fig 3. Feature series and probing efficiency for model setting I-A.

For population 1, excitability increase is correlated with variance even with passive observation. However, increased amplitude probing responses from population 2 reflect gradual changes of A1 towards ictogenesis. (a) shows a simplified illustration of the model setting. For (b) variance, (c) skewness, (d) kurtosis and (e) lag-1 AC features, top rows show feature values extracted from the output of each population, as the excitability gain parameter (A1) is increased. (f) shows the inter-population synchronization as measured by the mutual information feature. Bottom row shows Spearman’s rank correlation between each feature series and the shifted parameter (A1), as stimulus amplitude is increased. Black dots indicate significant differences from passive observation (stim. Amplitude = 0) according to Tukey’s HSD.

https://doi.org/10.1371/journal.pcbi.1008377.g003

For model setting I-A, the approaching transition to ictal activity of the ictogenic population 1 can be observed even without probing, due to the increase in signal variance. However, increasing stimulus amplitude (greater than 140 APs/s) affects other features for the “normal” population 2 –slight decreases in skewness and kurtosis can be observed with the approaching ictal transition, only when probing stimuli are used, suggesting a predictive effect of the perturbational approach in this setting. High stimulation amplitudes also highlight the increased mutual information between both model subsets.

For model setting I-B, Fig 4 shows slight changes in Lag-1 AC and variance (and to a lesser extent, in kurtosis as well) of population 1, regardless of stimulation. That is, as the SDI gain is decreased, no predictive effects of active probing are observed in the activity of the neuronal populations, and neither in the mutual information between both simulated signals.

[Figure omitted. See PDF.]

Fig 4. Feature series and probing efficiency for model setting I-B.

Decreased slow inhibition of population 1 is not highlighted by the extracted features if only the “normal” neural mass (set 2) is stimulated. (a) shows a simplified illustration of the model setting. For (b) variance, (c) skewness, (d) kurtosis and (e) lag-1 AC features, top rows show feature values extracted from the output of each population, as the slow inhibition gain parameter (B1) is decreased. (f) shows the inter-population synchronization as measured by the mutual information. Bottom row shows Spearman’s rank correlation between each feature series and the shifted parameter (B1), as stimulus amplitude is increased.

https://doi.org/10.1371/journal.pcbi.1008377.g004

For I-K, all features from the stimulated population correlate with the increased coupling factor K–and consequently with the time to ictal activity onset. However, as shown in Fig 5, these changes are not observed with the un-perturbed population 1. Altogether, the results for setting I indicate that the effectiveness of active probing to highlight preictal changes is limited to the model subset being perturbed. However, if unilateral stimulation is applied to the ictogenic population, as shown in S1 and S2 Figs, increased feature changes can be detected in both populations.

[Figure omitted. See PDF.]

Fig 5. Feature series and probing efficiency for model setting I-K.

Coupling gain increase is highlighted by features extracted from the stimulated population and by the mutual information between activities of both model subsets. (a) shows a simplified illustration of the model setting. For (b) variance, (c) skewness, (d) kurtosis and (e) lag-1 AC features, top rows show feature values extracted from the output of each population, as the coupling gain gain parameter (K) is increased. (f) shows the inter-population synchronization as measured by the mutual information feature. Bottom row shows Spearman’s rank correlation between each feature series and the shifted parameter (K), as stimulus amplitude is increased.

https://doi.org/10.1371/journal.pcbi.1008377.g005

Next, results are shown for model settings II, where both populations are simultaneously probed while one parameter is shifted towards ictal activity. This is done with the increase of A1 (EXC), the decrease of B1 (SDI), or the increase of coupling both populations (K). Fig 6 shows the resulting feature series and Spearman's rank correlation for model setting II-A.

[Figure omitted. See PDF.]

Fig 6. Feature series and probing efficiency for model setting II-A.

Active probing effects are visible for population 2, but only for increased stimulation amplitudes. (a) shows a simplified illustration of the model setting. For (b) variance, (c) skewness, (d) kurtosis and (e) lag-1 AC features, top rows show feature values extracted from the output of each population, as the excitability gain parameter (A1) is increased. (f) shows the inter-population synchronization as measured by the mutual information feature. Bottom row shows Spearman’s rank correlation between each feature series and the shifted parameter (A1), as stimulus amplitude is increased.

https://doi.org/10.1371/journal.pcbi.1008377.g006

For higher stimulation amplitudes in model setting II-A, feature changes are evident as the transition to ictal activity approaches. This can be seen in both populations, albeit to a lesser degree in the “normal” model subset, which parameters are held constant. Higher stimulus amplitudes result in sudden increases in variance near the threshold from normal to ictal activity, reflecting the appearance of increased sharp responses in the form of ILS as shown in Fig 2. Furthermore, increased stimuli amplitudes highlight increased mutual information between the simulated activities of both populations.

For the ictogenic population 1 in setting II-B, Fig 7 shows that active probing induces significant feature changes as the transition to ictal state approaches. Once again, the sharp increase in signal variance for higher stimulation amplitudes reflects enhanced responses or ILS shown in Fig 2. Preictal changes are much less visible in the “normal” simulated LFP from population 2, appearing only for skewness and kurtosis with stimulation amplitudes greater than 140 APs/s. On the other hand, consistent increases in mutual information between the activities of both populations are revealed by bilateral stimuli greater than 80 APs/s.

[Figure omitted. See PDF.]

Fig 7. Feature series and probing efficiency for model setting II-B.

Active probing provides predictive value for features extracted from both populations, but the effect is limited to increased stimulation amplitudes for population 2. (a) shows a simplified illustration of the model setting. For (b) variance, (c) skewness, (d) kurtosis and (e) lag-1 AC features, top rows show feature values extracted from the output of each population, as the slow inhibitory gain parameter (B1) is decreased. (f) shows the inter-population synchronization as measured by the mutual information feature. Bottom row shows Spearman’s rank correlation between each feature series and the shifted parameter (B1), as stimulus amplitude is increased.

https://doi.org/10.1371/journal.pcbi.1008377.g007

Similar feature profiles can be observed in Fig 8 with model setting II-K–sharp increases in variance, lag-1 AC, and mutual information and decrease in skewness and kurtosis, but only when probing stimuli are used.

[Figure omitted. See PDF.]

Fig 8. Feature series and probing efficiency for model setting II-K.

Probing provides predictive value for both stimulated populations. (a) shows a simplified illustration of the model setting. For (b) variance, (c) skewness, (d) kurtosis and (e) lag-1 AC features, top rows show feature values extracted from the output of each population, as the coupling gain parameter (K) is increased. (f) shows the inter-population synchronization as measured by the mutual information feature. Bottom row shows Spearman’s rank correlation between each feature series and the shifted parameter (K), as stimulus amplitude is increased.

https://doi.org/10.1371/journal.pcbi.1008377.g008

Testing for pro-ictal effects of active probing

To ensure that the probing stimuli have no pro-ictal effects, additional simulations were made to evaluate if the occurrence of ictal-like activity is precipitated by the probing stimuli. The same model settings were used, but parameters were shifted further to elicit sustained spike discharges–Fig 9 shows the onset times of this state for each population as a function of stimulus amplitude.

[Figure omitted. See PDF.]

Fig 9. Onset time of ictal-like activity as a function of stimulus amplitude.

Simulations with increased amplitudes did not affect the onset latencies of sustained spike discharges in any model setting.

https://doi.org/10.1371/journal.pcbi.1008377.g009

One-way ANOVA tests yielded overall significant differences in ictal-like activity onset times for model settings II-A (p = 0.028) and II-B (p = 0.000231). However, Tukey’s HSD of all stimulation amplitudes versus the passive observation case (amplitude = 0) did not result in any significant differences. That is, stimulation did not anticipate the occurrence of sustained discharges of spikes when compared to the passive observation case, implying that the probing approach did not have pro-ictal effects in any of the evaluated model settings and stimulation parameters.

Discussion

Using different model settings of coupled populations of a neural mass model, we have explored the ability of active probing to highlight early warning signs that reflect intrinsic parameter changes involved in triggering ictal states. In most model settings, the perturbational approach induced response changes that could be used as markers of changing parameters that lead to ictal activity–mainly the increase in variance and lag-1 autocorrelation and decreases in skewness and kurtosis of the simulated signals, as well as the increase in the inter-population synchrony. Importantly, a local effect of active probing was observed, as these changes are not observed in features extracted from the un-perturbed populations, except when unilateral stimuli are applied to the ictogenic population.

In model setting I, extracted features (namely variance and lag-1 AC) from the ictogenic population 1 are slightly sensitive to the decrease of slow inhibition (B1) and very sensitive to the increase in excitability (A1), even with passive observation. For features extracted from the probed population 2, the effects depended on which parameter was being altered. No detectable changes were found in setting I-B, but increases in the coupling factor (K) and excitability (A1) were correlated with features extracted from the perturbed population 2. In the case of I-A, the effects were subtle and restricted to high stimulation amplitudes. That is, the effects of active probing were local, being limited to the population being stimulated. This reaffirms the findings in Ref. [9], which highlighted the importance of suitably dense recordings for the active probing approach to be successful in extracting system information.

In model setting II, the predictive effects of probing the ictogenic population are observed regardless of which parameter shift was used to bring the model dynamics towards ictal activity–extracted features from stimuli responses are correlated with the increasing excitability (A1) or decreasing slow inhibition (B1), as well as the increase in coupling strength (K). Except for the increase in variance as A1 is increased, these changes are not visible with the passive observation of population activity. That is, the results indicate that the use of active probing to reveal approaching transitions from normal to ictal activity is very effective when stimuli are applied to the ictogenic population. Furthermore, the perturbational approach may benefit from the use of coordinated stimulation in different regions, since changes heralding the approach of ictal activity were visible not only in the ictogenic population activity, but also in the coupled neuronal populations with fixed parameters (as in the case of model setting II-A). This is especially relevant, considering that refractory epilepsy is increasingly being considered as a network disease [87]. The combination of stimulation of the ictal onset zone and remote structures has been shown previously to further increase DBS efficacy in epilepsy patients [88]. Network topology also plays a role in seizure generation and may be used for planning optimal surgery strategies–e.g., scale-free and rich-club networks display specific nodes that are critical for seizure generation and should be good candidates for removal [89]. Similar studies could show how different active probing configurations and parameters should be tailored for specific network topologies in order to leverage predictive power while assuring no pro-ictal effects occur.

In this work, external perturbation of a neural mass model highlighted common signatures of approaching critical transitions such as increased synchronization, variance, and lag-1 AC [38]. The increase in probing response variance and lag-1 autocorrelation is similar to the results found in Chang et al. [36], which examined the stability of the system’s response to stimulation of the Schaffer collaterals in isolated CA1 sections. The increased response amplitudes and longer recovery times from the perturbations were interpreted as markers of decreasing system resilience and critical slowing down. This is a common early warning signal of an impending critical transition, which are sudden shifts in dynamic states that have been reported in a myriad of complex dynamical systems [90] and are marked by common feature patterns [37]. It is still a matter of debate whether seizures would belong to this class of phenomena. Although it has been shown that human seizure termination can occur via critical transition [91], conflicting evidence has been found if this is true for the initiation of seizures. Comparing a model exhibiting self-organized criticality (SOC) with experimental data, deviations of a power-law were observed during human epileptic seizures and could indicate a shift of dynamics toward an ordered phase during ictal events [40]. Several other computational models involve critical transitions leading to seizure onsets [7,43,44,92] and experimental evidence in favor of this is found in Refs [36,45]. However, only a minority of patient seizures have been shown by Milanowski & Suffczynski [39] to present classical hallmarks or critical transitions before seizure onset. In another study, long-lasting iEEG recordings from 28 epilepsy patients did not reveal evidence for critical slowing down prior to epileptic seizures in humans [41]. On the other hand, a significant increase of lag1-AC was observed before seizure onset in 4 of 12 patients in Ref. [36], suggesting that critical slowing down should be evident in periods leading up to seizures in at least specific populations of epilepsy patients. This is also supported by a recent work which found that warning signals of critical slowing down (namely variance and autocorrelation) changed before seizure onset in 13 of 14 epilepsy patients over small temporal scales, but observed changes were also found in longer temporal scales (hours and days, as opposed to commonly observed scales of seconds to minutes) for 9 of 14 patients [42].

Critical slowing down can be evidenced by increased variance and longer recovery times from perturbations such as periodic stimuli or even by noise inputs [43]. This makes the use of probing stimuli an interesting tool to evaluate if seizure onsets are characterized by critical transitions in different computational and experimental seizure models, which would make its associated classical hallmarks promising markers for seizure forecasting, especially when combined with other features such as the rate of epileptiform spikes [42].

Active probing approaches may be viable but not needed in some cases where increased seizure susceptibility can be measured by intrinsic electrographic features, such as the increased signal variance of a population that has increased excitability parameter (A1), considering the model evaluated in this work. Intrinsic measures were used by Meisel et al. [93] to show how cortical excitability is reduced by antiepileptic drugs and increased across the wake period in patients. These measures were correlated with stimulation-evoked responses, supporting their viability as markers of cortical excitability and could play a role in seizure forecasting algorithms.

In model settings where stimuli response features were correlated with the shifted parameters, features evolved either gradually (as in I-K), or gradually until drastic changes appeared before ictal onset (II-K for high stimuli intensities). The latter changes are due to the appearance of sharp stimuli responses near the threshold to ictal activity and could be interpreted as interictal-like spikes (ILS) elicited by stimulation (see Fig 2). The results here are similar to the findings by Lopes et al. [86], which studied the dynamics of interictal spikes and early-warning signs of approaching critical transitions in a stochastic cellular automata neural network model (SCANNM). It was shown that the required external stimulation to evoke ILS became smaller as the model’s activity approached the ictal state, such that optimal stimuli intensity values could thus be tuned to elicit ILS in the minimally sufficient time for serving as clear early-warning signals of impending seizures. Electrical stimuli also appeared to induce ILS right before seizure onset in a PTZ seizure rat model in Medeiros et al. [45].

Probing stimuli were found to have no pro-ictal effect in the neural mass model used in this work. That is, sustained discharges of spikes were not precipitated by the use of low-frequency period stimulation. However, responses similar to ILS were evoked with increased stimuli intensities before the transition to ictal activity. While these responses can serve as clear early-warning signs, the role of interictal discharges in epilepsy and whether these events are anti or pro-ictal is still unclear [94–96]. Thus, despite the lack of pro-ictal effects of the active probing approach in this work and others [25,45,59], further experimental evidence is needed to ensure that the use of long-term probing stimuli would not induce pro-ictal effects while still providing seizure-forecasting value.

Conclusion

The lack of clinically validated and robust methods for seizure forecasting motivates the evaluation of active probing approaches, where the use of stimuli may be used to highlight electrographic features indicating impending ictal state transitions. By using a neural mass model with two coupled populations and different parameter changes leading up to seizures, this work showed that a low-frequency active probing approach provides predictive value. As parameters of a neuronal population were shifted from normal towards ictal activity, responses showed changes that were in general not visible with the passive observation of the populations’ activity, such as increased variance, lag-1 autocorrelation and synchrony between neuronal populations’ activities. A local effect of active probing was observed, if stimulation is applied to the ‘normal’ neuronal population. That is, activity dynamics heralding approaching transitions to ictal states are limited to the features being extracted from the stimulated population. The feature changes were consistent with the phenomenon of critical slowing down, where response variance and lag-1 autocorrelation increased as the model approached the ictal state. This is often related to the loss of a systems’ resilience and are common hallmarks of an impending critical transition. Furthermore, increasing stimulation amplitude did not cause significant changes in the onset times of sustained discharges of spikes, indicating the lack of pro-ictal effects of active probing in this neural mass model. Altogether, results indicate that low-frequency probing stimuli can provide predictive value and reveal underlying dynamics involved in transitions from normal to ictal states.

The optimistic projections that in the next few years part of the burden of epilepsy in some patients may be eased by seizure forecasting methods are motivated by current evidence of the feasibility of anticipating seizures and recent advances in the field [12]. New methods may leverage insights gained not only from active probing approaches discussed here, but also from increasingly available long-term iEEG databases [97], as well as increased adoption of new implantable devices that enable patient-tailored algorithms. Achieving this would be a crucial step that may enable to ultimately close the loop for targeted seizure-controlling interventions.

Supporting information

S1 Fig. Feature series and probing efficiency for unilateral probing of ictogenic population with increased excitability.

(a) shows a simplified illustration of the model setting. For (b) variance, (c) skewness, (d) kurtosis and (e) lag-1 AC features, top rows show feature values extracted from the output of each population, as the excitability gain parameter (A1) is increased. (f) shows the inter-population synchronization as measured by the mutual information feature. Bottom row shows Spearman’s rank correlation between each feature series and the shifted parameter (A1), as stimulus amplitude is increased. Black dots indicate significant differences from passive observation (stim. Amplitude = 0) according to Tukey’s HSD. For population 1, excitability increase (A1) is correlated with variance even with passive observation. Increased correlation of other features, including mutual information between population activities, with the ictogenic parameter A1 occurs only with stimulation. For features extracted from population 2, high stimulation amplitudes result in changes prior to ictal onset, even though this population was not stimulated–probably due to the ILS elicited in population 1.

https://doi.org/10.1371/journal.pcbi.1008377.s001

(TIF)

S2 Fig. Feature series and probing efficiency for unilateral probing of ictogenic population with decreased slow inhibition.

(a) shows a simplified illustration of the model setting. For (b) variance, (c) skewness, (d) kurtosis and (e) lag-1 AC features, top rows show feature values extracted from the output of each population, as the excitability gain parameter (A1) is increased. (f) shows the inter-population synchronization as measured by the mutual information feature. Bottom row shows Spearman’s rank correlation between each feature series and the shifted parameter (B1), as stimulus amplitude is increased. Black dots indicate significant differences from passive observation (stim. Amplitude = 0) according to Tukey’s HSD. For population 1, SDI decrease (B1) is slightly correlated with variance even with passive observation, but stimulation increases the correlation with this and other features, including mutual information between population activities. For features extracted from population 2, high stimulation amplitudes result in changes prior to ictal onset, even though this population was not stimulated–probably due to the ILS elicited in population 1.

https://doi.org/10.1371/journal.pcbi.1008377.s002

(TIF)

Acknowledgments

We would like to thank Daniel Medeiros, Flávio Mourão and Leonardo Guarnieri for helpful comments and discussions on the manuscript.

Citation: Carvalho VR, Moraes MFD, Cash SS, Mendes EMAM (2021) Active probing to highlight approaching transitions to ictal states in coupled neural mass models. PLoS Comput Biol 17(1): e1008377. https://doi.org/10.1371/journal.pcbi.1008377

References

1. Sackellares JC. Current Review in Clinical Science—Seizure Prediction. Epilepsy Curr. 2008;8: 55–59. pmid:18488065

2. Litt B, Echauz J. Prediction of epileptic seizures. Lancet Neurol. 2002;1: 22–30. Available: http://www.sciencedirect.com/science/article/pii/S1474442202000030 pmid:12849542

3. de Boer HM, Mula M, Sander JW. The global burden and stigma of epilepsy. Epilepsy Behav. 2008;12: 540–546. pmid:18280210

4. Cook MJ, O’Brien TJ, Berkovic SF, Murphy M, Morokoff A, Fabinyi G, et al. Prediction of seizure likelihood with a long-term, implanted seizure advisory system in patients with drug-resistant epilepsy: A first-in-man study. Lancet Neurol. 2013;12: 563–571. pmid:23642342

5. Nagaraj V, Lee ST, Krook-Magnuson E, Soltesz I, Benquet P, Irazoqui PP, et al. Future of seizure prediction and intervention: Closing the loop. J Clin Neurophysiol. 2015;32: 194–206. pmid:26035672

6. Mormann F, Andrzejak RG, Elger CE, Lehnertz K. Seizure prediction: the long and winding road. Brain. 2007;130: 314–33. pmid:17008335

7. Kalitzin SN, Velis DN, da Silva FHL. Stimulation-based anticipation and control of state transitions in the epileptic brain. Epilepsy Behav. Elsevier Inc.; 2010;17: 310–23. pmid:20163993

8. Schelter B, Timmer J, Schulze-Bonhage A. Seizure prediction in epilepsy: from basic mechanisms to clinical applications [Internet]. Neural Engineering. 2008. Available: http://link.springer.com/chapter/10.1007/0-306-48610-5_12

9. O’Sullivan-Greene E, Kuhlmann L, Nurse ES, Freestone DR, Grayden DB, Cook M, et al. Probing to Observe Neural Dynamics Investigated with Networked Kuramoto Oscillators. Int J Neural Syst. 2017;27: 1650038. pmid:27596927

10. Kalitzin S, Parra J, Velis DN, Lopes da Silva FH. Enhancement of phase clustering in the EEG/MEG gamma frequency band anticipates transitions to paroxysmal epileptiform activity in epileptic patients with known visual sensitivity. IEEE Trans Biomed Eng. 2002;49: 1279–86. pmid:12450358

11. Stacey WC, Litt B. Technology insight: neuroengineering and epilepsy-designing devices for seizure control. Nat Clin Pract Neurol. 2008;4: 190–201. pmid:18301414

12. Freestone DR, Karoly PJ, Cook MJ. A forward-looking review of seizure prediction. Curr Opin Neurol. 2017;30: 167–173. pmid:28118302

13. Lozano AM, Lipsman N. Probing and Regulating Dysfunctional Circuits Using Deep Brain Stimulation. Neuron. Elsevier Inc.; 2013;77: 406–424. pmid:23395370

14. Freestone DR, Karoly PJ, Peterson ADH, Kuhlmann L, Lai A, Goodarzy F, et al. Seizure Prediction: Science Fiction or Soon to Become Reality? Curr Neurol Neurosci Rep. 2015;15. pmid:26404726

15. Lopes da Silva FH, Blanes W, Kalitzin SN, Parra J, Suffczynski P, Velis DN. Dynamical diseases of brain systems: different routes to epileptic seizures. IEEE Trans Biomed Eng. 2003;50: 540–8. pmid:12769430

16. Penfield Wilder & Jasper H. Epilepsy and the functional anatomy of the human brain,. 1954.

17. O’Donnell BF, Vohs JL, Krishnan GP, Rass O, Hetrick WP, Morzorati SL. The auditory steady-state response (ASSR): A translational biomarker for schizophrenia. Suppl Clin Neurophysiol. 2013;62: 101–112. pmid:24053034

18. Oda Y, Onitsuka T, Tsuchimoto R, Hirano S, Oribe N, Ueno T, et al. Gamma Band Neural Synchronization Deficits for Auditory Steady State Responses in Bipolar Disorder Patients. Hashimoto K, editor. PLoS One. 2012;7: e39955. pmid:22792199

19. Matsubara T, Ogata K, Hironaga N, Uehara T, Mitsudo T, Shigeto H, et al. Monaural 40-Hz auditory steady-state magnetic responses can be useful for identifying epileptic focus in mesial temporal lobe epilepsy. Clin Neurophysiol. International Federation of Clinical Neurophysiology; 2019;130: 341–351. pmid:30669010

20. Casali AG, Gosseries O, Rosanova M, Boly M, Sarasso S, Casali KR, et al. A Theoretically Based Index of Consciousness Independent of Sensory Processing and Behavior. Sci Transl Med. 2013;5: 198ra105–198ra105. pmid:23946194

21. Kimiskidis VK, Tsimpiris A, Ryvlin P, Kalviainen R, Koutroumanidis M, Valentin A, et al. TMS combined with EEG in genetic generalized epilepsy: A phase II diagnostic accuracy study. Clin Neurophysiol. International Federation of Clinical Neurophysiology; 2017;128: 367–381. pmid:28007469

22. Spiegel A, Mentch J, Haskins AJ, Robertson CE. Slower Binocular Rivalry in the Autistic Brain. Curr Biol. Elsevier Ltd.; 2019; 1–6. pmid:31422885

23. Badawy R, MacDonell R, Jackson G, Berkovic S. The peri-ictal state: Cortical excitability changes within 24 h of a seizure. Brain. 2009;132: 1013–1021. pmid:19251759

24. Bauer PR, Helling RM, Perenboom MJL, Lopes da Silva FH, Tolner EA, Ferrari MD, et al. Phase clustering in transcranial magnetic stimulation-evoked EEG responses in genetic generalized epilepsy and migraine. Epilepsy Behav. Elsevier Inc.; 2019;93: 102–112. pmid:30875639

25. Freestone DR, Kuhlmann L, Grayden DB, Burkitt AN, Lai A, Nelson TS, et al. Electrical probing of cortical excitability in patients with epilepsy. Epilepsy Behav. Elsevier Inc.; 2011;22: S110–S118. pmid:22078511

26. Valentin A, Arunachalam R, Mesquita-Rodrigues A, Garcia Seoane JJ, Richardson MP, Mills KR, et al. Late EEG responses triggered by transcranial magnetic stimulation (TMS) in the evaluation of focal epilepsy. Epilepsia. 2008;49: 470–480. pmid:18028404

27. Kunieda T, Yamao Y, Kikuchi T, Matsumoto R. New approach for exploring cerebral functional connectivity: Review of cortico-cortical evoked potential. Neurol Med Chir (Tokyo). 2015;55: 374–382. pmid:25925755

28. Valentín A, Alarcón G, Honavar M, García Seoane JJ, Selway RP, Polkey CE, et al. Single pulse electrical stimulation for identification of structural abnormalities and prediction of seizure outcome after epilepsy surgery: a prospective study. Lancet Neurol. 2005;4: 718–726. pmid:16239178

29. Van’T Klooster MA, Zijlmans M, Leijten FSS, Ferrier CH, Van Putten MJAM, Huiskamp GJM. Time-frequency analysis of single pulse electrical stimulation to assist delineation of epileptogenic cortex. Brain. 2011;134: 2855–2866. pmid:21900209

30. Matsubara T, Ogata K, Hironaga N, Kikuchi Y, Uehara T, Chatani H, et al. Altered neural synchronization to pure tone stimulation in patients with mesial temporal lobe epilepsy: An MEG study. Epilepsy Behav. The Authors; 2018;88: 96–105. pmid:30243112

31. Pinto HPP, Carvalho VR, Medeiros D de C, Almeida AFS, Mendes EMAM, Moraes MFD. Auditory processing assessment suggests that Wistar audiogenic rat neural networks are prone to entrainment. Neuroscience. 2017;347: 48–56. pmid:28188855

32. Pinto HPP, de Oliveira Lucas EL, Carvalho VR, Mourão FAG, de Oliveira Guarnieri L, Mendes EMAM, et al. Seizure Susceptibility Corrupts Inferior Colliculus Acoustic Integration. Front Syst Neurosci. 2019;13: 105–121. pmid:31780904

33. Kalitzin S, Velis D, Suffczynski P, Parra J, Lopes Da Silva F. Electrical brain-stimulation paradigm for estimating the seizure onset site and the time to ictal transition in temporal lobe epilepsy. Clin Neurophysiol. 2005;116: 718–728. pmid:15721087

34. Freestone DR, Long SN, Frey S, Stypulkowski PH, Giftakis JE, Cook MJ. A method for actively tracking excitability of brain networks using a fully implantable monitoring system. Conf Proc. Annu Int Conf IEEE Eng Med Biol Soc IEEE Eng Med Biol Soc Annu Conf. 2013;2013: 6151–6154. https://doi.org/10.1109/EMBC.2013.6610957

35. de Castro Medeiros D, Raspante LBP, Mourão FAG, Carvalho VR, Mendes EMAM, Moraes MFD. Deep brain stimulation probing performance is enhanced by pairing stimulus with epileptic seizure. Epilepsy Behav. 2018; pmid:30352775

36. Chang WC, Kudlacek J, Hlinka J, Chvojka J, Hadrava M, Kumpost V, et al. Loss of neuronal network resilience precedes seizures and determines the ictogenic nature of interictal synaptic perturbations. Nat Neurosci. Springer US; 2018;21: 1742–1752. pmid:30482946

37. Scheffer M, Carpenter SR, Lenton TM, Bascompte J, Brock W, Dakos V, et al. Anticipating critical transitions. Science. 2012. pp. 344–348. pmid:23087241

38. Scheffer M, Bascompte J, Brock WA, Brovkin V, Carpenter SR, Dakos V, et al. Early-warning signals for critical transitions. Nature. Nature Publishing Group; 2009;461: 53–59. pmid:19727193

39. Milanowski P, Suffczynski P. Seizures Start without Common Signatures of Critical Transition. Int J Neural Syst. 2016;26: 1–15. pmid:27633895

40. Meisel C, Storch A, Hallmeyer-Elgner S, Bullmore E, Gross T. Failure of adaptive self-organized criticality during epileptic seizure attacks. PLoS Comput Biol. 2012;8. pmid:22241971

41. Wilkat T, Rings T, Lehnertz K. No evidence for critical slowing down prior to human epileptic seizures. Chaos. AIP Publishing LLC; 2019;29. pmid:31575122

42. Maturana MI, Meisel C, Dell K, Karoly PJ, D’Souza W, Grayden DB, et al. Critical slowing down as a biomarker for seizure susceptibility. Nat Commun. Springer US; 2020;11: 2172. pmid:32358560

43. Negahbani E, Steyn-Ross DA, Steyn-Ross ML, Wilson MT, Sleigh JW. Noise-Induced Precursors of State Transitions in the Stochastic Wilson–Cowan Model. J Math Neurosci. 2015;5: 1–27. pmid:26458900

44. Breakspear M, Roberts JA, Terry JR, Rodrigues S, Mahant N, Robinson PA. A unifying explanation of primary generalized seizures through nonlinear brain modeling and bifurcation analysis. Cereb Cortex. 2006;16: 1296–1313. pmid:16280462

45. Medeiros DDC, Oliveira LB, Mourão FAG, Bastos CP, Cairasco NG, Pereira GS, et al. Temporal Rearrangement of Pre-ictal PTZ Induced Spike Discharges by Low Frequency Electrical Stimulation to the Amygdaloid Complex. Brain Stimul. Elsevier Ltd; 2014;7: 170–178. pmid:24332185

46. Veraart AJ, Faassen EJ, Dakos V, Van Nes EH, Lürling M, Scheffer M. Recovery rates reflect distance to a tipping point in a living system. Nature. 2012;481: 357–359. pmid:22198671

47. Berzhanskaya J, Chernyy N, Gluckman BJ, Schiff SJ, Ascoli GA. Modulation of hippocampal rhythms by subthreshold electric fields and network topology. J Comput Neurosci. 2013;34: 369–389. pmid:23053863

48. Izhikevich EM. Neural Excitability, Spiking and Bursting. Int J Bifurc Chaos. 2000;10: 1171–1266.

49. Bastos AM, Usrey WM, Adams RA, Mangun GR, Fries P, Friston KJ. Canonical Microcircuits for Predictive Coding. Neuron. Elsevier Inc.; 2012;76: 695–711. pmid:23177956

50. Lytton WW. Computer modelling of epilepsy. Nat Rev Neurosci. 2008;9: 626–637. pmid:18594562

51. Kalitzin S, Koppert M, Petkov G, Velis D, da Silva FL. Computational model prospective on the observation of proictal states in epileptic neuronal systems. Epilepsy Behav. Elsevier Inc.; 2011;22: S102–S109. pmid:22078510

52. Taylor PN, Kaiser M, Dauwels J. Structural connectivity based whole brain modelling in epilepsy. J Neurosci Methods. Elsevier B.V.; 2014;236: 51–57. pmid:25149109

53. Kudela P, Franaszczuk PJ, Bergey GK. Model of the propagation of synchronous firing in a reduced neuron network. Neurocomputing. 1999;26–27: 411–418.

54. Tass PA. A model of desynchronizing deep brain stimulation with a demand-controlled coordinated reset of neural subpopulations. Biol Cybern. 2003;89: 81–88. pmid:12905037

55. Sinha N, Dauwels J, Kaiser M, Cash SS, Westover MB, Wang Y, et al. Predicting neurosurgical outcomes in focal epilepsy patients using computational modelling. Brain. 2017;140: 319–332. pmid:28011454

56. Proix T, Bartolomei F, Guye M, Jirsa VK. Individual brain structure and modelling predict seizure propagation. Brain. 2017;140: 641–654. pmid:28364550

57. Olmi S, Petkoski S, Guye M, Bartolomei F, Jirsa V. Controlling seizure propagation in large-scale brain networks. PLoS Comput Biol. 2019;15: 1–23. pmid:30802239

58. Wang Y, Trevelyan AJ, Valentin A, Alarcon G, Taylor PN, Kaiser M. Mechanisms underlying different onset patterns of focal seizures. PLoS Comput Biol. 2017;13: 1–22. pmid:28472032

59. Wendling F, Gerber U, Cosandier-Rimele D, Nica A, De Montigny J, Raineteau O, et al. Brain (Hyper)Excitability Revealed by Optimal Electrical Stimulation of GABAergic Interneurons. Brain Stimul. Elsevier Inc.; 2016;9: 919–932. pmid:27576186

60. Molaee-Ardekani B, Márquez-Ruiz J, Merlet I, Leal-Campanario R, Gruart A, Sánchez-Campusano R, et al. Effects of transcranial Direct Current Stimulation (tDCS) on cortical activity: A computational modeling study. Brain Stimul. 2013;6: 25–39. pmid:22420944

61. Nagaraj V, Lamperski A, Netoff TI. Seizure Control in a Computational Model Using a Reinforcement Learning Stimulation Paradigm. Int J Neural Syst. 2016;27: 1750012. pmid:28030999

62. Suffczynski P, Kalitzin S, Da Silva FL, Parra J, Velis D, Wendling F. Active paradigms of seizure anticipation: Computer model evidence for necessity of stimulation. Phys Rev E—Stat Nonlinear, Soft Matter Phys. 2008;78: 1–9. pmid:19113165

63. Carvalho VR, Moraes MFD, Mendes EMAM. On the Validity of Using Probing Stimuli for Seizure Prediction in the Epileptor Model. In: Cota V, Barone D, Dias D, Damázio L, editors. Communications in Computer and Information Science. Springer, Cham; 2019. pp. 269–281. https://doi.org/10.1007/978-3-030-36636-0_20

64. Mina F, Benquet P, Pasnicu A, Biraben A, Wendling F. Modulation of epileptic activity by deep brain stimulation: A model-based study of frequency-dependent effects. Front Comput Neurosci. 2013;7: 1–16. pmid:23355821

65. Wendling F, Bartolomei F, Bellanger JJ, Chauvel P. Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition. Eur J Neurosci. 2002;15: 1499–1508. pmid:12028360

66. Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol. 1952;117: 500–544. pmid:12991237

67. Wendling F, Bartolomei F, Mina F, Huneau C, Benquet P. Interictal spikes, fast ripples and seizures in partial epilepsies—combining multi-level computational models with experimental data. Eur J Neurosci. 2012;36: 2164–2177. pmid:22805062

68. Wilson HR, Cowan JD. Excitatory and Inhibitory Interactions in Localized Populations of Model Neurons. Biophys J. 1972;12: 1–24. pmid:4332108

69. Freeman WJ. Mass Action in the Nervous System [Internet]. New York: Academic Press; 1975. https://doi.org/10.1016/C2009-0-03145-6

70. Lopes da Silva FH, Hoeks A, Smits H, Zetterberg LH. Model of brain rhythmic activity. The alpha-rhythm of the thalamus. Kybernetik. 1974;15: 27–37. Available: http://www.ncbi.nlm.nih.gov/pubmed/4853232 pmid:4853232

71. Wright JJ, Liley DTJ. Simulation of electrocortical waves. Biol Cybern. 1995;72: 347–356. pmid:7748961

72. Medvedeva TM, Sysoeva M V., Lüttjohann A, van Luijtelaar G, Sysoev I V. Dynamical mesoscale model of absence seizures in genetic models. Charpier S, editor. PLoS One. 2020;15: e0239125. pmid:32991590

73. Taylor PN, Wang Y, Goodfellow M, Dauwels J, Moeller F, Stephani U, et al. A Computational Study of Stimulus Driven Epileptic Seizure Abatement. Bazhenov M, editor. PLoS One. 2014;9: e114316. pmid:25531883

74. Segal M. A correlation between hippocampal responses to interhemispheric stimulation, hippocampal slow rhythmic activity and behaviour. Electroencephalogr Clin Neurophysiol. 1978;45: 409–411. pmid:79479

75. Mirski MA, Rossell LA, Terry JB, Fisher RS. Anticonvulsant effect of anterior thalamic high frequency electrical stimulation in the rat. Epilepsy Res. 1997;28: 89–100. pmid:9267773

76. Cota VR, Drabowski BMB, de Oliveira JC, Moraes MFD. The epileptic amygdala: Toward the development of a neural prosthesis by temporally coded electrical stimulation. J Neurosci Res. 2016;94: 463–485. pmid:27091311

77. Wendling F, Hernandez A, Bellanger JJ, Chauvel P, Bartolomei F. Interictal to ictal transition in human temporal lobe epilepsy: Insights from a computational model of intracerebral EEG. J Clin Neurophysiol. 2005;22: 343–356. pmid:16357638

78. Harris CR, Millman KJ, van der Walt SJ, Gommers R, Virtanen P, Cournapeau D, et al. Array programming with NumPy. Nature. Springer US; 2020;585: 357–362. pmid:32939066

79. Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, et al. Scikit-learn: Machine Learning in Python. J Mach Learn Res. 2012;12: 2825–2830. https://dl.acm.org/citation.cfm?id=2078195

80. McKinney W. Data Structures for Statistical Computing in Python. 2010. pp. 56–61.

81. Hunter JD. Matplotlib: A 2D Graphics Environment. Comput Sci Eng. 2007;9: 90–95.

82. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al. SciPy 1.0—Fundamental Algorithms for Scientific Computing in Python. 2019.

83. Quian Quiroga R, Panzeri S. Extracting information from neuronal populations: information theory and decoding approaches. Nat Rev Neurosci. 2009;10: 173–185. pmid:19229240

84. Ross BC. Mutual Information between Discrete and Continuous Data Sets. Marinazzo D, editor. PLoS One. 2014;9: e87357. pmid:24586270

85. Zwillinger D, Kokoska S. CRC Standard Probability and Statistics Tables and Formulae. New York: Chapman & Hall; 2000.

86. Lopes MA, Lee KE, Goltsev A V. Neuronal network model of interictal and recurrent ictal activity. Phys Rev E. 2017;96: 1–8. pmid:29347379

87. Zweiphenning WJEM, Keijzer HM, Diessen E, van ‘t Klooster MA, Klink NEC, Leijten FSS, et al. Increased gamma and decreased fast ripple connections of epileptic tissue: A high-frequency directed network approach. Epilepsia. 2019; epi.16296. pmid:31329277

88. Vonck K, Sprengers M, Carrette E, Dauwe I, Miatton M, Meurs A, et al. A decade of experience with deep brain stimulation for patients with refractory medial temporal lobe epilepsy. Int J Neural Syst. 2013;23. pmid:23273130

89. Lopes MA, Richardson MP, Abela E, Rummel C, Schindler K, Goodfellow M, et al. An optimal strategy for epilepsy surgery: Disruption of the rich-club? PLoS Comput Biol. 2017;13: e1005637. pmid:28817568

90. Brett T, Ajelli M, Liu Q-H, Krauland MG, Grefenstette JJ, van Panhuis WG, et al. Detecting critical slowing down in high-dimensional epidemiological systems. Pitzer VE, editor. PLOS Comput Biol. 2020;16: e1007679. pmid:32150536

91. Kramer MA, Truccolo W, Eden UT, Lepage KQ, Hochberg LR, Eskandar EN, et al. Human seizures self-terminate across spatial scales via a critical transition. Proc Natl Acad Sci. 2012;109: 21116–21121. pmid:23213262

92. Jirsa VK, Stacey WC, Quilichini PP, Ivanov AI, Bernard C. On the nature of seizure dynamics. Brain. 2014;137: 2210–2230. pmid:24919973

93. Meisel C, Schulze-Bonhage A, Freestone D, Cook MJ, Achermann P, Plenz D. Intrinsic excitability measures track antiepileptic drug action and uncover increasing/decreasing excitability over the wake/sleep cycle. Proc Natl Acad Sci. 2015;112: 14694–14699. pmid:26554021

94. Karoly PJ, Freestone DR, Boston R, Grayden DB, Himes D, Leyde K, et al. Interictal spikes and epileptic seizures: Their relationship and underlying rhythmicity. Brain. 2016;139: 1066–1078. pmid:26912639

95. Avoli M, de Curtis M, Köhling R. Does interictal synchronization influence ictogenesis? Neuropharmacology. 2013;69: 37–44. pmid:22776544

96. Staley KJ, Dudek FE. Interictal Spikes and Epileptogenesis. Epilepsy Curr. 2006;6: 199–202. pmid:17260059

97. Reuben C, Karoly P, Freestone DR, Temko A, Barachant A, Li F, et al. Ensembling crowdsourced seizure prediction algorithms using long-term human intracranial EEG. Epilepsia. 2019; epi.16418. pmid:31883345