Cortical chimera states predict epileptic seizures

Chimera states exhibit a hybrid structure of coexisting synchronous (coherent) and asynchronous (incoherent) behavior. This phenomenon was first described by Kuramoto and Battogtokh in 2002¹ and named “chimera states” by Abrams and Strogatz in 2004.² Panaggio and Abrams³ published a systematic review of the field and collected several references to “bump states,” as they were called before the invention of chimera states. In fact, in 1967, Winfree⁴ discussed synchronization and mixed states. Experimental chimeras were reported by Tinsley et al.,⁵ Hagerstrom et al.,⁶ Nkomo et al.,⁷ Abrams and Strogratz,² and Dudkowski et al.,⁸ among others. Outside laboratories, many real world chimeras have been reported. These include unihemispheric sleep, ventricular fibrilation, power grids, social systems, neural networks, and the brain (see Refs. 3 and 9 for references). Tognoli and Kelso¹⁰ report the coexistence of metastable and stable states in electroencephalography (EEG) data for cognitive tasks. Chimera states have also been reported by Showalter and his group and a special issue in Chaos was dedicated to him in 2019.¹¹ 

The connection between human epilepsy and chimeras has been of recent interest, and networks of model neurons have been studied.^12–17 The present study investigated this question and hypothesized that the asynchrony-synchrony transition that leads to epileptic seizures can be explained by the induction and dynamics of chimera states. We analyzed electrocorticography (ECoG) data from patients with intractable focal epilepsy, and found that in one patient, each seizure was consistently preceded by a chimera state. Therefore, in this patient, chimera states were 100% predictive of seizures before their onsets. Additionally, our analyses revealed that the brain stayed in these chimera states for up to 2 hours prior to a seizure. This is, to the best of our knowledge, the first direct evidence of chimera states in the human epileptic brain.

Delay differential analysis (DDA) is a nonlinear computational method that combines differential embeddings with linear and nonlinear nonuniform functional delay embeddings^18–20 to relate the current derivatives of a system to the current and past values of the system variables.^21,22 Inspired by Max Planck’s “natural units”,²³ the DDA model maps experimental data onto a set of natural embedding coordinates.

The general nonlinear DDA model is

for $τn,mn,i∈N0$⁠, where $N$ is the number of delays (usually 2), $I$ is the number of terms (typically around 3), and $xτn=x(t−τn)$⁠, relating the signal derivative $x˙$ to the signal nonuniformly shifted in time. We then used the coefficients $ai$ and the least square error $ρ$ as features. To restrict the complexity of the DDA model, most of the terms in Eq. (1) were set to zero. We, therefore, considered here DDA models with two delays $τn$⁠, three terms, and one degree $∑imn,i≤4$ of nonlinearity.

DDA is a computational framework that (1) uses unprocessed data so as not to disturb the nonlinear properties of the data, (2) uses sparse models that match the macroscopic architecture of the underlying dynamical system, (3) disregards amplitude information to focus on the dynamical aspects of the data, and (4) can be extended for the detection of dynamical causality in the data to understand information flow in the system. Previous applications to brain data have demonstrated that DDA efficiently captured essential properties of brain dynamics and classified brain states with a high level of performance.^24–26 

The patient demographics and characteristics are described in Ref. 25.

A DDA model with minimum error was chosen using a genetic algorithm (GA) to maximally discriminate between preictal and postictal recordings. Around one million one second data segments were used from one hour periods centered around the seizure onset times.²⁵ 

The selected DDA model is

with $xi=x(t−τi)$ and $τ=(7,10)$$δt$⁠, where $δt=1fs$ with the sampling rate $fs=500$ Hz. This model (2) bifurcates at the seizure onset, as shown in Ref. 25: It has a low error after seizure onset, but not before. The delays were chosen in the same way: $τ2=10$ is a good delay (low error) before seizures and $τ1=7$ is a good delay after seizures. The combination $τ=(7,10)δt$ then proved to be a good choice and we, therefore, used this model for seizure characterization. As shown in Refs. 25 and 27, the parameter $a1$ by itself is sufficient for seizure analysis.

All of the analyses in this article were performed on 100 h of recording from one patient (electrode locations are shown in Fig. 1) with a sampling rate of 500 Hz and 77 implanted electrodes (one bad electrode was excluded from the analysis). The sliding data windows were half a second long and the window shift was a quarter second. The data were not preprocessed, except for normalizing each data window to zero mean and unit variance. Figure 2 shows the DDA feature $a1$⁠, and Fig. 3 shows the phases of the DDA feature $a1$ as a function of time for all 77 implanted electrodes. The channel numbers were rearranged such that the first 32 channels are the onset channels as determined by the neurologist and the remaining 45 channels are the nononset channels. Each of the seven seizures has a preceding time segment of chimera states. These states last for around 2 hours (2.1951, 2.0418, 2.2509, 0.7108, 0.5645, 1.2334, and 2.6969 hours) before the seizure onsets.

Figure 4 shows the phases of $a1$ at two time points: 10 min before and 10 min after the fourth seizure. Before the seizure onset, the onset channels (the first 32 channels) were desynchronized, while the nononset channels were synchronized. In contrast, all channels became synchronized after the seizure offset.

An expanded version of seizure 4 is shown in Fig. 5 (upper panel). Microevents of seizurelike behavior can be observed before the seizure onset. Interestingly, after seizure onset the mixed synchrony/asynchrony behavior switches during the seizure. At seizure termination most of the channels are immediately synchronized, but a few need more time to synchronize with the rest. A movie for this seizure from 5 min before to 5 min after the seizure onset is shown in the supplementary material (S4.mp4). This behavior is consistent across all seven seizures as can be seen in Fig. 6 where the phases of all seven seizures are plotted and centered around the seizure onsets (one hour of data) as marked by the neurologist.

We calculated the order parameter as introduced by Kuramoto,¹ 

from $a1(t)$ for all 77 implanted electrodes, where $k$ is the rearranged channel number. In Fig. 5 (lower panel), we show this order parameter for a three minute window around seizure 4. Figure 7(a) shows the order parameter for the entire recording period, and Fig. 7(b) shows the zoomed time frame around the seizures. The seizures have the lowest order parameters followed by the time segments before the seizures during which chimera states are observed, consistent with Refs. 28 and 29.

In addition to the order parameter, which was computed from phases, we introduced a new DDA order parameter that can be calculated directly from DDA feature $a1$ (see Fig. 8): all pairwise L2 differences in $a1$ between channels were computed in a sliding window 5 s in duration (20 $a1$ values) and a shift of 250 ms. A smaller pairwise difference indicates higher degree of similarity and synchronization between channels. The channels in the lower 30% of the L2 distances were displayed as a function of time (see Fig. 8). In a chimera state, the synchronized channels have values that are between 0.5 and 1 [dark to light brown colors in Figs. 8(c) and 8(d)] and in the desynchronized state, the channels have values that are between 0 and 0.5 [light to dark gray colors in Figs. 8(c) and 8(d)]. In a nonchimera state, where all channels are either synchronized or desynchronized the value for all channels should be one or zero [light colors in Figs. 8(c) and 8(d)]. The DDA order parameter is then the sum of these values for each time window divided by $N$⁠, where $N$ is the number of channels.

We found chimera states in human ECoG data for one patient that reliably preceded each of this patient’s seven seizures. These chimeras can be seen in the DDA feature $a1$ and in the phases of $a1$⁠. We have so far analyzed 15 patients and only found chimeras in one patient. This might be because we used the same delays for all patients. We restricted our analysis to this one delay pair because of the amount of data to be analyzed. We, therefore, developed a DDA order parameter that detects chimera states in a completely automatic manner. With such an automatic procedure, we hope to find more chimera states in other patients by looking at more delay pairs.

Up to now, chimera states were reported in neuronal and epilepsy models. The findings in this paper are the first evidence of chimeras in the human brain.

See the supplementary material for a movie of the phases of seizure 4 from 5 min before to 5 min after the seizure onset.

This work was supported by the National Institute of Health (NIH)/NIBIB (Grant No. R01EB026899-01), NIH/NINDS (Grant No. R01NS104368), Office of Naval Research (ONR) (Grant No. N0014-15-1-2328), NIH/NINDS (Grant Nos. NS062092 and NS088568), and mission funding from the U.S. Army Research Laboratory. The authors would like to thank Robert Kim and Christopher Gonzalez for valuable discussions.