Parsing spatiotemporal dynamical stability in ECoG during seizure onset, propagation, and termination
Arian Ashourvan, Sérgio Pequito, Ankit N. Khambhati, Steven N. Baldassano, Kathryn A. Davis, Timothy Lucas, Jean M. Vettel, Brian Litt, George J. Pappas, Danielle S. Bassett
PParsing spatiotemporal dynamical stability in ECoGduring seizure onset, propagation, and termination
Arian Ashourvan , S ´ergio Pequito , Ankit N. Khambhati , Steven N. Baldassano ,Kathryn A. Davis , Timothy Lucas , Jean M. Vettel , Brian Litt , George J. Pappas ,and Danielle S. Bassett * Both authors contributed equally to this work. Department of Bioengineering, School of Engineering and Applied Science, University of Pennsylvania,Philadelphia, PA 19104 USA Department of Electrical and Systems Engineering, School of Engineering and Applied Science, University ofPennsylvania, Philadelphia, PA 19104 USA Penn Center for Neuroengineering and Therapeutics, University of Pennsylvania, Philadelphia, PA 19104 USA Department of Neurology, Hospital of the University of Pennsylvania, Philadelphia, PA 19104 USA Department of Neurosurgery, Hospital of the University of Pennsylvania, Philadelphia, PA 19104 USA U.S. Army Research Laboratory, Aberdeen Proving Ground, MD 21005 USA Department of Psychological & Brain Sciences, University of California, Santa Barbara, CA, 93106 USA † Correspondence to: Danielle S. Bassett, 210 S. 33rd Street, Philadelphia, PA 19104-6321,Email: [email protected]
ABSTRACT
Understanding brain dynamics in epilepsy is critical for establishing rigorous control objectives that enable new therapeuticmethods to mitigate seizure occurrence. In multichannel electrocorticography (ECoG) recordings acquired in 21 subjectsduring a total of 94 seizures, we apply dynamical systems stability analysis to assess the balance versus imbalance ofseizure dynamics across different timescales and brain regions. Specifically, we consider a sliding time window multivariateautoregressive linear approximation of the data captured by the ECoG channels, where eigendecomposition of the estimatedmatrix of coefficients describes the contribution of different regions to the spatiotemporal process (eigenvectors) associatedwith a particular timescale (eigenvalues). Interestingly, we observe a pattern of eigenvalue evolution and slowly changing (orapproximately time-invariant) eigenvectors across both seizures and subjects. The seizure-onset is marked by an increase inhigh frequency spatial information to which a few regions contribute for a long period. By contrast, the seizure termination ischaracterized by a sudden, small time period change in dynamics to which many regions contribute. As the seizure terminates,the relatively stable ictal dynamics rapidly transition into the post-ictal regime, marked by relatively fast-damping oscillations.Our methodology offers a subject-specific characterization of the spatiotemporal behavior of the seizure, providing new insightsinto the dynamic patterns and functional interactions between brain regions that occur over different timescales. More generally,our approach informs the development of engineering objectives that can be used to deploy new control strategies to preventseizure evolution or to hasten seizure termination.
Keywords: intracranial electrocorticography (ECOG); seizure-onset detection; multivariate time-series analysis; dynamicalstability analysis; eigenvalue-eigenvector structure
Introduction
Understanding the dynamic nature of the neurophysiological processes underlying seizure initiation is critical for developinginsight into epileptogenesis, seizure evolution, and seizure termination (Pitk¨anen et al., 2016). Fundamental knowledge aboutthese processes could augment burgeoning neurotechnologies, such as implantable devices (Morrell, 2011), by informingpatient-centric algorithms that modulate brain dynamics to abort seizures (Stacey and Litt, 2008a; Stanslaski et al., 2012; Afsharet al., 2013). Ideally, such an algorithm would account for both when and where to deliver a perturbative payload within thebrain’s physiological network. Indeed, spatiotemporal precision in stimulation protocols could improve quality of life for thenearly 1% of the world’s population affected by epilepsy (Mormann et al., 2007). Advances in stimulation protocols may befacilitated by theoretically-driven principles from engineering and control theory (Schiff, 2012), which can inform interventionsto mitigate the physiological effects of ictal activity (Schindler et al., 2007). Progress will critically depend on building a deeperunderstanding of the types of dynamical states that seizures represent in the context of broader neural systems. a r X i v : . [ q - b i o . N C ] J un ver the past decade or more, an important body of work has sought to connect seizure physiology to dynamical systemstheory (Schindler et al., 2007; Rummel et al., 2013). Several exploratory tools have been developed to identify features in EEGand ECoG signals that are associated with ictal processes and can serve as biomarkers for seizure detection and prediction(Baldassano et al., 2017; Brinkmann et al., 2016). These features are identified either by concepts in non-parametric statisticsand signal processing, or by multivariate analysis of parametric models (Schomer and Da Silva, 2012). The former approachcan provide classifications of seizure activity and seizure-onset types based on analysis of EEG or ECoG recordings (Peruccaet al., 2014). The latter approach relies on parametric models and has the added benefit of potentially capturing the substrateof the dynamic process itself – reflected in spatiotemporal properties of the EEG or ECoG activity (Schindler et al., 2007;Schiff et al., 2005; Nair et al., 2004; Kramer et al., 2010; Warren et al., 2010; Wilke et al., 2011; Burns et al., 2014; Khambhatiet al., 2015). The development of more accurate predictors and more physiologically interpretable models is an important andongoing area of research.These efforts are particularly challenging because epileptiform activity manifests across a range of temporal and spatialscales and displays marked inter-channel dependencies (Jirsa et al., 2014). Understanding the biological phenotype thereforerequires the ability to assess interactions between different channels across different timescales, with the potential of providingtestable causal mechanisms of seizure generation, progression and termination. To address this challenge, we exercise aparadigm that relies on the concept of dynamical systems stability to study the behavior of neural processes in the vicinity of aputative bifurcation between dynamically stable and unstable behaviors (Magnasco et al., 2009). The approach is built on thenotion that seizures can occur as a result of a dysfunction in regulatory mechanisms that might otherwise maintain stable braindynamics (Koepp et al., 2016). Our approach provides a compact representation and an intuitive visual mapping that capturesthe interactions between channels across different time-scales – a capability that can potentially provide a dynamics-basedinterpretation of epileptic neural signals in clinical neuroscience. Further, our approach offers a dynamical definition of the onsetand termination phases of a seizure, a critical step toward understanding the natural processes that govern seizure dynamics(Schiff et al., 2005).We apply this dynamical stability analysis approach to ECoG data acquired from twenty one subjects undergoing surgicaltreatment for medically refractory epilepsy that is likely of neocortical origin. Our results suggest that assessments basedon signal frequency alone may miss salient changes in the spatio-temporal evolution of the dynamical interaction processescharacteristic of epileptiform activity. Moreover, we demonstrate that spatio-temporal disruptions in dynamical stabilitymark periods of both seizure onset and propagation. Our data suggest that inter-ictal dynamics are an emergent macroscopicphenomenon characterized by coherent activity among different cortical regions, whereas the ictal regime is characterizedby the emergence of focal and persisting dynamics. Finally, we note that the proposed approach provides us with a criterionthat is useful for future efforts in stimulation control, supporting the development of new stimulation strategies to be used intherapeutic implantable devices (Osorio et al., 2005; Kim et al., 2009; Echauz et al., 2009; Morrell, 2006). Materials and Methods
Description of experimental data
In this study, we considered twenty one subjects undergoing surgical treatment for medically refractory epilepsy, likely ofneocortical origin. These subjects underwent the implantation of subdural electrodes to assess seizure location and evolutionafter noninvasive monitoring proved inconclusive. All data had previously been de-identified and was publicly available via theInternational Epilepsy Electrophysiology Portal (IEEG Portal) (Wagenaar et al., 2013). Patients exhibited a range of seizureetiology, location, type, and degree of severity, providing a rich dataset for the assessment of seizure dynamics. In Table 1, wereport important details about the subjects that were monitored at the Hospital of the University of Pennsylvania (Philadelphia,PA, USA) and at the Mayo Clinic (Rochester, MN, USA); subjects from these two cohorts were labeled as HUP and Study,respectively. Specifically, in Table 1, we report the gender (Sex) as well as the age at first reported onset and at phase IImonitoring (Age). Additionally, we report the localization of seizure-onset and the seizure etiology, which was clinicallydetermined through medical history, imaging, and long-term invasive monitoring. The different seizures observed (SeizureTypes) included simple-partial (SP), complex-partial (CP), and complex-partial with secondary generalization (CP+GTC).In this table, we also indicate the total number of seizures recorded in the epilepsy monitoring unit, as well as the clinicalimaging analysis (Imaging) that concludes if the seizure etiology is lesional (L) or non-lesional (NL). Finally, surgical outcome(Outcome) was based on either Engel score or ILAE score: seizure freedom to no improvement (I-V), no resection (NR), andno follow-up (NF).
Intracranial electroencephalogram recordings
ECoG signals were recorded and digitized at sampling rates of 512 Hz (Hospital of the University of Pennsylvania, Philadelphia,PA) and 500 Hz (Mayo Clinic, Rochester, MN). Surface electrode (Ad Tech Medical Instruments, Racine, WI) configurations,determined by a multidisciplinary team of neurologists and neurosurgeons, consisted of linear and two-dimensional arrays ubject (
Patient (IEEG Portal)
Sex Age (Years)(Onset/Surgery) seizure-onset Etiology Seizure Type Seizures (
Imaging Outcome
Table 1.
Patient Information . For each patient, we report gender (Sex), and age at first reported onset and at phase IImonitoring (Age). Additionally, we report the localization of seizure-onset and the seizure etiology, which wasclinically-determined through medical history, imaging, and long-term invasive monitoring. The different seizures observed(Seizure Types) include simple-partial (SP), complex-partial (CP), and complex-partial with secondary generalization(CP+GTC). We also indicate the total number of seizures recorded in the epilepsy monitoring unit, as well as the clinicalimaging analysis (Imaging) that concludes if the seizure etiology is lesional (L) or non-lesional (NL). Finally, surgical outcome(Outcome) was based on either Engel score or ILAE score: seizure freedom to no improvement (I-V), no resection (NR), andno follow-up (NF).(2.3 mm diameter with 10 mm inter-contact spacing) and sampled the neocortex for epileptic foci (depth electrodes were firstverified as being outside the seizure-onset zone and subsequently discarded from this analysis). Signals were recorded using areferential montage with the reference electrode, chosen by the clinical team, distant to the site of seizure-onset. Recordingspanned the duration of a patient’s stay in the epilepsy monitoring unit.A total of 94 seizure-onsets were determined and the seizure-onset time and localization were defined by the point of“earliest electrographic change” (EEC), which was annotated and marked by a team of practicing epileptologists (Litt et al.,2001). The EEC associated with a given seizure is found by moving backward in time from the unequivocal electrographiconset (UEO) until the earliest change from previous background activity can be observed visually on the ECoG traces. TheUEO is defined as the first unequivocal intracranial EEG sign of change from the background, which leads to a clear seizuredischarge without return to background activity (Spanedda et al., 1997). Using the EEC, we extracted multi-channel ECoGtime series from seizure onset to seizure termination. In addition, we also extracted ECoG time series 20 seconds directlypreceding (pre-ictal) and following (post-ictal) each seizure. Lastly, we extracted a total of 478 (22.7 ± . > . Dynamical stability approach
ECoG signals can be locally approximated at each time point by a linear system (Khalil, 2014). Let x ( k ) ∈ R n be the data from n channels, where the i th entry x i ( k ) corresponds to the data collected by channel i at time k . Then, we obtain x ( k ) = Ax ( k − ) + ε ( k ) , (1)where A is the n × n real matrix that results from fitting the data collected within the interval of time [ k − τ , k + τ ] using a leastsquares approach, and where ε ( k ) is the approximation error. As a consequence of the linear approximation, the results willdepend on the choice of the number n of channels and the size of the time window parametrized by τ . We will discuss theimpact of these parameters at greater length in the Results section and in the Supplementary Materials document.A significant advantage of using a linear approximation of sensor dynamics is that the dynamical properties of the underlyingprocess can be locally assessed through the eigendecomposition of A : that is, the n eigenvalue-eigenvector pairs. In otherwords, eigenvalue-eigenvector pairs capture linearly independent spatiotemporal dynamical processes. Specifically, for eachputative spatiotemporal process, the eigenvector weights the involvement of each sensor and the complex eigenvalue defines the requency of the oscillatory dynamics. More importantly, the stability of the dynamics can be captured by the absolute value ofthe eigenvalues, which forecast the exponential growth or decay along the associated eigenvector.Specifically, let A = V λ V (cid:124) be the eigendecomposition for a given time point k , where V = [ v , . . . , v n ] and λ = diag ( λ , . . . , λ n ) are the matrices of eigenvectors and eigenvalues (when the eigenvalues are all distinct complex numbers) respectively, where ( λ i , v i ) are the eigenvalues paired with the associated eigenvector, and where V (cid:124) is the transpose of V . Notice that someeigenvalues are complex numbers, which implies that no partial-order can be imposed on the eigendecomposition. Further, notethat after T time steps, one obtains from Eq. (1) that x ( k ) = A T x ( ) , which implies that x ( k ) = V λ T V (cid:124) x ( ) . Subsequently, byconsidering a linear combination of the original data z ( k ) = V ∗ x ( k ) , where z i ( k ) = v (cid:124) i x ( k ) is a weighted combination describedby the i th eigenvector associated with the i th eigenvalue.From these definitions, it follows that | z i ( k ) | = | λ i | T | z i ( ) | and three scenarios are possible: ( i ) | λ i | < | z i ( k ) | → t → ∞ ; ( ii ) | λ i | > | z i ( k ) | → ∞ as t → ∞ ; and | λ i | = | z i ( k ) | = | z i ( ) | for all times. In( i ), the process tends to vanish, and we therefore refer to these dynamics as asymptotically stable . In ( ii ), the process tendsto explode, and we therefore refer to these dynamics as unstable . Finally in ( iii ), the process oscillates between stability andinstability, and we therefore refer to these dynamics as stable . In practice, we consider that the stable regime is determinedby | λ i | ≈
1. Since these processes are associated with z ( k ) , it follows that they are associated with a specific eigenvector,and, consequently, we refer to a stable regime associated with a particular eigenvector. Thus, the interplay between thesethree different stages and different eigenvectors provides a dynamical stability characterization to identify seizure onset and tomonitor seizure evolution.To further characterize these dynamics, we note that the angle θ i associated with the polar coordinates of the i th complexeigenvalue provides a description of the frequency as follows: f i = θ i π δ t , where δ t corresponds to the sampling frequency; the timescale is given by ρ i = log ( | λ i | ) δ t which can be interpreted as the growth rate. Therefore, the dynamical process z ( k ) describes the spatiotemporal behavior ofthe dynamical system, where the timescale is encoded in the eigenvalue and the spatial scale is encoded in the eigenvector,indicating the relative contribution of a given channel or – by extension – the overall activity in a specific cortical region.In the main text, we use this approach to consider the time-evolution of eigenvalue-eigenvector pairs associated with higherand lower frequencies, and associated with fine and coarse timescales, to better understand seizure onset, propagation, andtermination. In the Supplementary Materials, we define and study synthetic examples of damping versus growing oscillationsand build the reader’s intuition for how dynamical stability analysis allows us to understand and quantitatively characterize theunderlying process. Statistical testing
Because the functional forms of the distributions of our statistics of interest were not known a priori , we applied non-parametricpermutation testing to assess statistical significance of our findings.As we will describe in greater detail in the Results section, we first sought to assess the extent of the time-invariance ofthe eigenvectors over the ictal period. We began by defining a summary statistic of interest. Specifically, we computed thesum of the element-wise squared differences between temporally contiguous eigenvectors; then, we calculated the averagestandard deviation in this value across all temporally contiguous pairs of eigenvectors across an ictal event. Intuitively, thissummary statistic quantifies the temporal variance of the given eigenvector across time: high values of this statistic indicatelarge temporal fluctuations, and small values of this statistic indicate small temporal fluctuations.Next, we assessed the frequency-dependence of this statistic at the group level by clustering eigenvectors into five non-overlapping frequency bands based on the eigenvectors’ average frequencies (across independent model fits over the timewindowed data) over each sample (see Supplementary Materials Fig. 11 and 12 and associated text for further details). Wecalculated subject-level summary statistics, averaged over all eigenvectors in each band, for (i) the ictal samples, averaged overall samples for each patient, and (ii) the inter-ictal samples, averaged over all samples for each patient. Then, we comparedthese summary statistics between ictal and inter-ictal samples at the group-level via bootstrap analysis. To create non-parametricnull distributions of the group-level average of the difference between the ictal and inter-ictal summary statistics computed atthe subject-level, we permuted the labels “ictal” and “inter-ictal” uniformly at random 50,000 times. Finally, we tested thehypothesis that the eigenvectors in the ictal samples at each frequency band on average displayed significantly less temporalfluctuations than inter-ictal samples across subjects. ollowing the examination of the time-dependence versus time-invariance of eigenvectors during the ictal period, and thefrequency dependence of our findings as a function of ictal versus inter-ictal samples, we finally turned to an assessment ofthe regional localization of the effects. To measure regional localization, we studied the absolute magnitude of the maximallycontributing elements of the eigenvectors. Then, we used a bootstrap analysis to test the hypothesis that the emergence of theslow-changing focal ictal dynamics was associated with an increase in the absolute magnitude of the maximally contributingelements of eigenvectors, particularly within a short time window (1, 3, 6, and 9 s) following the seizure onset. To determinethe specificity of the effects to the ictal period, we compared the subject-level summary statics in the ictal segments to thoseestimated from the pre-ictal segments (see Supplementary Materials Fig. 13 and associated text for further details).
Results
In ECoG traces obtained from twenty-one patients with medically refractory epilepsy, we used dynamical stability analysisto estimate Eq. (1) over a 1 s sliding-window for 100 ms shifts. Specifically, we estimated the matrices A using previouslydeveloped computational algorithms (Neumaier and Schneider, 2001), and then we calculated the eigenmodes or eigenvalue-eigenvector pairs of the estimated A matrices for each time point k . The eigenvectors are considered on the basis of theirspatiotemporal frequency: that is, the frequency determined by the corresponding eigenvalue. Because the eigenvector maycontain several non-zero entries, the frequency of the dynamic process along the direction determined by the eigenvector doesnot coincide with any specific ECoG channel. Instead, it corresponds to the spatiotemporal frequency along the directioncaptured by the eigenvector. For ease of visualization, we normalize the eigenvectors by dividing each entry by the maximumvalue obtained across the entries, and then we take the absolute value. Thus, all entries correspond to values between 0 and 1,and capture the contribution of different ECoG channels to the evolution of the underlying neurophysiological process at acertain spatiotemporal frequency. Dynamical stability analysis pinpoints seizure onset by distinct spatiotemporal patterns
Based on our preliminary observation of the sustained and focal nature of the ictal onset activity, we hypothesized thateigenmodes would evolve with slow spatio-temporal variance, characterized by eigenvectors that remained relatively consistentin their composition and magnitude over a period of time. Moreover, we anticipated distinct phases of a seizure could bedistinguished based on their frequency as given by the angle associated with the eigenvalues, and based on their stability asmeasured by the radius associated with the eigenvalues. To examine these hypotheses, we first consider the phenotypes in asingle subject, and then we examine the consistency of our results across the population.We begin by considering a randomly selected seizure from a single subject in our sample; this seizure is a complex-partialseizure with secondary generalization, as can be seen in the ECoG traces presented in Fig. 1A. As the seizure begins, weobserve the emergence of a sustianed regime whereby the eigenvectors associated with high frequency oscillatory dynamicsdisplay approximate time-invariance over several seconds (Fig. 1B-C: period between 20-27 s). Next, we observe that theeigenvector associated with the highest frequency (eigenvector 1) has high magnitude elements, in a small set of corticalregions: electrodes 4, 11,14, 74, and 79 annotated by red arrows in Fig. 1B at the period between 20-27 s. (For a thoroughcomparison between these marking and clinical markings, please see the supplementary materials section “Comparison toclinical annotations”.) This point characterized by regional specificity of the eigenvector loadings, as well as approximatetime-invariance of the eigenvectors marks seizure initiation. Following this point, we observe that the eigenvector associatedwith the highest frequency tracks the emergence of asymptotically stable behavior (i.e., fast-dampening oscillations), and againdisplays a period of approximate time-invariance in eigenvector loading patterns (Fig. 1B: period between 27-33 s). We notethat the loadings onto the eigenvector in this period are distinct from that in the previous seizure initiation period.Importantly, the intra-subject similarity of the eigenvectors’ temporal evolution at seizure-onset suggests that seizurestend to be led by a small set of cortical regions, consistent with the clinically-accepted notion that seizures originate fromspecific anatomical loci . These general trends were consistently observed across seizures within a single subject (Fig. 2,compare panels A-D drawn from four seizures experienced by the same subject). Notably, we observed that while the samequalitative phase-transition at seizure-onset occurs across seizures, the duration of the periods characterized by approximatetime-invariance of eigenvector loading patterns can differ across seizures.Following these qualitative observations, we next performed a quantitative assessment to statistically characterize eigenvectorevolution and to support our claim that – at seizure onset – the high frequency eigenvectors display approximate time-invarianceof loading patterns. Specifically, we calculated the standard deviation of the average pairwise differences between eigenvectorloading patterns in neighboring time points (see Methods). We consistently observed periods characterized by greaterspatiotemporal invariance than expected in inter-ictal data (Supplementary Materials Fig. 11). To show this, we extractedinter-ictal segments ( N = . ± .
35 per patient) from ECoG recordings that were taken more than 3 hours before or 3hours after a seizure (see Methods). Compared to inter-ictal segments, we observed that ictal onset is marked by a statisticallysignificant increase in spatio-temporal invariance across several eigenvectors with frequency content in the range 4 to 55 Hz .
10 20 30 40 50 60 70 80 90 100 110 12010 20 30 40 50 60 70 80 90 100 110 120Seconds10 20 30 40 50 60 70 80 90 100 110 12010 20 30 40 50 60 70 80 90 100 110 12010 20 30 40 50 60 70 80 90 100 110 120 E l e c t r o d e s B. E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r C.D.E. | θ | ( R a d ) | θ | ( R a d ) | θ | ( R a d ) | θ | ( R a d ) |λ||λ||λ||λ| E l e c t r o d e s E l e c t r o d e s E l e c t r o d e s E l e c t r o d e s Figure 1.
Evolution of eigenmodes over the ictal period. (A)
ECoG signals from 79 electrodes in subject 3 over a singleictal period, from 20 s before seizure onset to 20 s after seizure termination. Pre-ictal activity (yellow), earliest electrographicchange (orange), unequivocal seizure-onset (pink), ictal period (purple), and post-ictal period (gray) are delineated separately. (B-E)
The frequency (cyan trace) and stability (red trace) of eigenvalues associated with four representative eigenvectors – twoassociated with high frequencies and two associated with low-frequencies – whose evolution is displayed in heatmaps. Here weobserve that the onset of ictal activity is marked by the emergence of a period of approximate spatiotemporal-invariance in theloading patterns on eigenvectors associated with high frequency dynamics (see approximate spatiotemporal invariance from20-27 s and from 27-32 s in (B) ). We also observe an increase in the angle of the eigenvalues time-locked to the onset of theseperiods. Pointer arrows in panels (A) and (B) are offered as guides to the eye when comparing epileptiform activity in theECoG traces and in the eigenvector associated with the highest frequency dynamics. see Supplementary Results for more details). These comparisons between ictal and inter-ictal dynamics demonstrate thatthe temporal evolution of the eigenvectors reliably highlights the beginning of the seizure-onset period. These findings areconsistent with those from the single-subject analysis discussed previously.
Dynamical stability analysis tracks seizure progression and marks seizure offset
Our results thus far have demonstrated that the seizure-onset period is characterized by approximate time-invariance ofeigenvector loading patterns, and that these dynamics are unlike that observed in inter-ictal periods. It is next natural to examinethe dynamics of seizure propagation and termination. In both simple-partial and complex-partial seizures, we find that only afew electrodes have high loadings on high frequency eigenvectors, and that these loadings are maintained through the initialtime-invariant eigenvector period (Fig. 3). Interestingly, these regionally-specific, high eigenvector loadings are greater inmagnitude than expected in either inter-ictal or peri-ictal data (see Supplementary Materials Fig. 12 and Fig. 13, respectively).By contrast, the gradual emergence of a large number of electrodes across both hemispheres with high contribution to theeigenvectors marks seizure generalization (Fig. 3).Seizure offset is marked by sudden changes in the neurophysiological processes that manifest as a return to inter-ictaldynamics. Dynamical stability analysis can characterize these changes at seizure termination through transient changes in thestability metric | λ | (see red trace in northern insets of panels Fig. 1B-E, Fig. 2A-E, and Fig. 3B-E). We observed that in all ofthe CP and CP+GTC samples (after removing one CP+GTC sample with high artifact) the stability of several eigenmodesdrops immediately following clinically-marked seizure termination (Fig. 7). In fact, the average stability across all eigenmodesclearly captures this sudden decrease of the eigenmodes’ stability time-locked to seizure offset (Fig. 7A). Although the largestdecrease in the stability is observed in the CP+GTC samples, normalized average stability curves highlight the similaritybetween the profile of the stability drop across samples (Fig. 7B). We quantified the changes in the eigenmodes’ stability bycalculating the difference between in the average stability over a 3 s window before and after seizure offset. As seen in Fig. 7C,all the CP and CP+GTC samples show lower stability following seizure offset. Moreover we observed that this decrease instability is not unique to a single frequency band, but instead can observed across several frequency bands for the majority ofsamples. Nevertheless, the most robust results were observed for eigenmodes in the β -band (Fig. 7C). It is worth noting thatthis observation did not hold for the few SP seizures, where samples did not reveal a notable pattern of stability fluctuationsat the seizure offset. Following seizure termination, the angles and radii returned to inter-ictal values (less stability) quickly,and the eigenvector pattern becomes less time-invariant, and more unstructured (Fig. 1B). These observations reveal that thelong-lasting and stable epileptic activity ends abruptly at the seizure offset, which is quickly followed by a transition into aregime of relatively fast-damping dynamics. Inter-subject variability and seizure-type specificity
At the group-level, we observe (i) relative time-invariance of eigenvector loading patterns during the ictal period, (ii) this relativetime-invariance is unlike that observed in the inter-ictal phase, and is characteristic of eigenmodes representing dynamics inthe 4-55 Hz range, and (iii) seizure-onset is characterized by regional localization of eigenvector loading patterns, suggestingthe presence of anatomically localized epileptic foci . In addition to these group-level findings, we also observed interestingidiosyncrasies across subjects. Specifically, we observed that the relative time-invariance of loading patterns was better capturedin the eigenvectors associated with lower frequencies ( i.e. , 0-8 Hz) for patients with δ – θ onset epileptic activity as specifiedin the clinical annotations. Moreover, we observed that ictal samples varied in dynamical stability and in the patterns ofeigenvector evolution, suggesting pronounced variability in the spatiotemporal processes corresponding to the evolution andtermination of seizures. Finally, we also observed individual differences in the speed of the transition from seizure terminationto inter-ictal dynamics: some subjects showed a rapid transition (e.g., see subject 4 in Fig. 3) while others showed a moregradual transition (e.g., see subject 3 in Fig. 3). These findings point to individual variation in recovery mechanisms that drivebrain dynamics from unstable seizure states to normative and stable states.Importantly, we note that some of the inter-subject variability in the spatiotemporal evolution of the eigenvectors can beexplained by seizure type or characteristic frequency. Specifically, we observed that generalized seizures displayed a gradualincrease in the number and diversity of regions with high loadings onto the eigenvectors with the largest eigenvalues. In contrast,the contributions of different regions remained relatively unchanged over time in partial seizures ( e.g. , see subject 4 in Fig. 3).Further, we observed that the estimated angle and stability of eigenmodes during the ictal and inter-ictal periods differed inspatiotemporal frequency, as might be expected due to the different anatomical locations of ECoG sensors. For example, theangle of the highest frequency eigenmode (group average 60.4 ± ± > ± ± A. | θ | ( R a d ) |λ| E i g e n v e c t o r B. | θ | ( R a d ) |λ| E i g e n v e c t o r D. | θ | ( R a d ) |λ| E i g e n v e c t o r C. | θ | ( R a d ) |λ| E i g e n v e c t o r E i g e n v e c t o r E. Seconds |λ| | θ | ( R a d )
10 20 30 40 50 60 70 80 90 100 110 120 13010 20 30 40 50 60 70 80 90 100 110 120 130 14010 20 30 40 50 60 70 80 90 100 110 120 130 14010 20 30 40 50 60 70 80 90 100 110 120
Figure 2.
Evolution of the first eigenmode for subject 3 over several ictal periods, indicating cross-seizure consistency. (A)
The temporal evolution of eigenmode 1 for the seizure previously presented in Fig. 1. (B-E)
The temporal evolution ofeigenmode 1 across different seizures from the same subject. We consistently observe periods of approximate time-invarianceof eigenvector loading patterns immediately following the onset of the seizure. The seizure shown in panel (E) is a complexpartial seizure without secondary generalization; here we observe that the eigenmode properties during seizure propagationdiffer from that observed in seizures that secondarily generalize (panels (A-D) ). Pointer arrows in all panels are offered asguides to the eye for sensors with high loadings onto the highest frequency eigenvector. . E l e c t r o d e s B.C.D.E. | θ | ( R a d ) | θ | ( R a d ) | θ | ( R a d ) | θ | ( R a d ) |λ| E l e c t r o d e s E l e c t r o d e s
10 60 11010 60 110 E l e c t r o d e s E i g e n v e c t o r E l e c t r o d e s
10 60 110 E i g e n v e c t o r
10 30 50 10 50 9010 50 901020304050607080 10 50 9010203040506070801020304050607080 10 50 901020304050607080 10 50 901020304050607080 10 50 90 102030405060102030405060 1020304050102030405010203040501020304050102030405010 30 50102030405060 10 30 50102030405060 10 30 50Seconds102030405060 10 30 50 E i g e n v e c t o r E i g e n v e c t o r |λ||λ||λ|
10 50 9010 50 9010 50 9010 50 90 10.80.60.40.210.80.60.40.210.80.60.40.210.80.60.40.2 E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r Subject
Figure 3.
Evolution of sample eigenmodes of four subjects over ictal periods, indicating inter-subject variability andseizure-type specificity. (A)
ECoG signals from all electrodes in four subjects over a single ictal period, from 20 s beforeseizure onset to 20 s after seizure termination. (B-E)
The frequency (cyan trace) and stability (red trace) of eigenvaluesassociated with four representative eigenvectors – two associated with high frequencies and two associated withlow-frequencies – whose evolution is displayed in heatmaps. The seizures presented for subjects 3 and 16 are complex partialwith secondary generalization (compare to Fig. 2A-D), while the seizure for subject 18 is complex partial without secondarygeneralization (compare to Fig. 2E), and the seizure for subject 4 is simple partial. eliability and replicability across methodological choices
Finally, since the window length and the local approximation model can affect the estimated dynamics, we tested a varietyof window sizes (1, 2, 5, and 10 s) and higher-order linear models to assess the impact of these methodological choiceson our results. First, we demonstrate via synthetic examples that the linear model with a short window of 1s is able toaccurately capture broad timescales and overlapping dynamics, whereas longer windows lack temporal sensitivity (Fig. 5).These numerical studies are consistent with the proposed synthetic example constructed to assess the reliability of the proposedmethods (see Supplementary Materials for details). Specifically, longer time windows yield more reliable estimates of lowfrequency oscillatory modes (as measured by fewer time-points with unstable modes, | λ i | > i.e., Discussion
Here, we seek to understand the spatiotemporal evolution of seizure dynamics using dynamical stability features. Theproposed approach decomposes local dynamics into eigenmodes: eigenpairs in which the eigenvectors track the contributionof brain regions to the evolution of the modeled dynamical process and in which the eigenvalues monitor that process’stime-scale, frequency, and stability. We observe a relative time-invariance of eigenvector loading patterns during the ictal period.Interestingly, this relative time-invariance in the recruited spatiotemporal processes is unlike that observed in the inter-ictalphase, and is characteristic of eigenmodes representing dynamics in the 4-55 Hz range. We also find that seizure initiationis characterized by regional localization of eigenvector loading patterns, suggesting the presence of anatomically localizedepileptic foci . Finally, we observe that seizure termination is marked by a decrease in stability. Then, the angles and radii returnto inter-ictal values (less stability) quickly, and the eigenvector pattern becomes less time-invariant, and more unstructured. Ourapproach directly addresses the challenge that epileptiform activity manifests across a range of temporal and spatial scales, anddisplays marked inter-sensor dependencies. By assessing interactions between different channels across different time-scales,our approach provides putative mechanisms of seizure generation, progression and termination. More generally, our approachmay inform the development of new control strategies for seizure treatment.
Dynamic stability characterization of seizure-onset
A longstanding theory regarding the dynamical progression of seizures is that they evolve through different stages or dynamicalstates (Wulsin et al., 2013; Rummel et al., 2013; Burns et al., 2014; Khambhati et al., 2015). While previous studies havemotivated approaches to characterize dynamical and topological features of these states, a systems-level understanding ofhow these dynamical transitions occur has been elusive. Theoretical and experimental evidence points to a mechanismof seizure generation and progression in which the epileptic brain is more susceptible to transitions between stable andasymptotically stable dynamical states. However these studies undersample the spatiotemporal statistics of the network byrelying on information acquired only from a single channel (Jirsa et al., 2014; Worrell et al., 2002).In contrast, here we develop and apply a novel approach that enables us to characterize the onset and termination of seizuresin terms of the dynamical modes of the neural system by accounting for the epileptic network’s rich spatiotemporal structure.First, we find that the dynamics immediately following the EEC demonstrate an approximate spatio-temporal invariance ofeigenvectors associated with several frequency bands (ranging from 4-55 Hz), along with stereotypical shifts in the frequencyand stability of eigenmodes common across patients in peri-ictal events. These observations suggest that the emergence ofpersistent focal oscillatory sources eventually drives the underlying system into an ictal regime (Perucca et al., 2014; Curtis andAvoli, 2016; Khambhati et al., 2017). Second, we find that seizure termination is characterized by a transition from a stablestate in which a few regions have high eigenvector loadings, to a more asymptotically stable state in which many regions havecomparable eigenvector loadings, consistent with the dynamics observed during inter-ictal periods.Importantly, eigenmodes of the system distinguish between commonly diagnosed seizure types. In line with current clinicalclassification of seizure types (Fisher et al., 2017), we found that partial seizures – whose pathological dynamics remainfocal – present approximate time-invariance of several eigenvectors during the entire seizure with contributions from only a ew regions. In contrast, generalized seizures – whose pathological dynamics progress to bilateral tonic–clonic – present aheterogenous distribution of patterns as seizures spread to different regions and eventually across the hemispheres. Importantly,the observation that different sorts of seizures have different patterns of inter-sensor dynamics is not wholly new and has alsobeen reported in studies of time-resolved patterns of functional connectivity in ECoG (Khambhati et al., 2016; Schindler et al.,2006). Indeed, it is interesting to speculate that the previously reported seizure-type-specificity of the synchronizability of theECoG network estimated from functional connectivity may be related to differences in the underlying process, as revealedby our modeling approach. More generally, a fundamental understanding of the processes driving different seizure typescould inform the development of seizure-type-specific interventions in the form of stimulation via implantable devices perhapsthrough objective functions that seek to dynamically control parameters such as frequency and stability of the neural system.
Underlying drivers of seizure-onset and prevention
Our methodological approach enables an assessment of changes in dynamical stability between non-seizure and seizure periods.We find that the inter-ictal period is characterized by high temporal variance of the system’s eigenvectors, reminiscent offindings in prior work suggesting that the epileptic network traverses different dynamical states more frequently during theinter-ictal period than during the ictal period (Khambhati et al., 2015, 2017; Burns et al., 2014). We also find that the inter-ictalperiod is characterized by limited changes in frequency and stability. These insights from the dynamical stability-basedcharacterization suggest that a possible mechanism to mitigate seizure propagation would be to stimulate different regionstowards a homogenization of the contributions of these regions to the local dynamics. That is, rather than only targetingseizure-onset nodes, it may be important to stimulate nodes in the surrounding tissue to create heterogeneity in the dynamics,thereby increasing resilience to changes in the stability of the dynamical process. A similar suggestion has been made previously:that degree heterogeneity and lower synchronizability of the network relates to increased focalization of seizures, and that thelargest contributors to this relationship are found in the surrounding tissue, rather than in the seizure onset zone (Khambhatiet al., 2016; Weiss et al., 2013; Schevon et al., 2012). This study predicted that if these key electrodes in the surrounding tissuewere removed or virtually lesioned, one would expect a profound perturbation of the entire network’s synchronizability. Animportant area for future work is to understand the degree to which stimulation paradigms that alter synchronizability also alterdynamical stability, as studied here.Importantly, prior work has outlined the general aim of stimulating a region such that the surrounding tissue changesits properties, while ensuring that the stimulation protocol itself satisfies certain constraints (Alonso et al., 2015). Exampleconstraints include not targeting regions that are deemed unsafe or unethical to stimulate, and not using stimulation parametervalues (such as for frequency and magnitude) that are deemed ineffective, not technologically feasible, or unsafe for humantissue. In principle, our stability-based assessment could inform similar studies, by defining a time-varying stimulation strategythat alters the dynamical stability in such a way as to quiet or prevent seizures while minimizing the potential of pushing thesystem into a dynamical stability regime that is uncontrollable. Indeed, theoretical work based on similar principles to thosewe use here to model the ECoG data has begun to develop a framework to specify such a stimulation strategy (Pequito et al.,2017). In principle, should this theory prove effective, it could be used to inform ECoG sensor placement during installationof implantable devices, and sensor stimulation parameters by closed loop measurement and manipulation of system stability.The key advantage of this framework is that – if validated clinically – it could offer a thoroughly quantitative approach tostimulation protocols for implanted devices.More speculatively, the features that we obtain from eigenmode decomposition could be used as quantifiable performanceobjectives for optimizing stimulation parameters in closed-loop implantable devices (Stacey and Litt, 2008b; Johnson et al.,2013). Specifically, this approach would depart from the current practice of manual, intermittent, off-line tuning of deep brainstimulation parameters such as amplitude and frequency based on the count of seizures experienced by the patient. Instead,the features that characterize the seizure-onsets using the dynamical stability-based approach can be used to select channelsfor feature extraction in a closed-loop assessment and stimulation paradigm (Wu and Gotman, 1998). For example, sensorsincluded in the calculation of time-resolved seizure probability could be limited to those for which eigenvalue-eigenvectorfeatures provide reliable estimates of stability. This data reduction could simplify simplify and speed calculations, whilemaintaining optimal predictions about effective stimulation protocols.
Methodological Considerations
Several methodological considerations are pertinent to this work. First, the proposed analysis does not directly allow forthe assessment of underlying biological mechanisms. Second, although it is likely that the underlying dynamical processcaptured by ECoG data is non-linear, the dynamical stability -based approach allows for the analysis of global nonstationarityand non-linear dynamics by fitting a linear model to short temporal windows. Third, to assess reliability and reproducibility, weexamined the robustness of our results to different temporal windows and different models (see Supplement). Our reliablefindings therefore provide evidence for a more general phenomenon that cannot be explained by either noise or non-linearity(Netoff et al., 2004). Fourth, we have demonstrated that the dynamical stability analysis allows us to identify similar changes in he underlying process during the seizure onset period for different samples across patients. Despite our rather homogeneousdataset, which includes a large number of seizures with overlapping features, heterogeneity is still present in the evolutionof seizures across subjects and samples. It would be interesting in future to examine changes by brain region, assessing forexample whether a frontal onset pattern differs significantly from a temporal onset pattern, or whether a right-dominant onsetpattern differs significantly from a left dominant onset pattern. It would also be interesting to examine the duration of thereduced stability during the post-ictal period, and determine whether that duration differs by seizure type. Many factors suchas the location and number of electrodes, sampling rate, recording noise, and type of seizures partially explain the observedinter-subject variations. Nevertheless, the dynamical stability analysis provides an avenue for characterization of individualpatients’ idiosyncrasies and for establishing individualized features that can be used as biomarkers.
Future Directions
In summary, we provide a dynamical stability-based characterization of the onset of seizure activity that hinges on an eigenmodedecomposition of a computationally tractable model. Furthermore, the proposed approach uses linear stability analysis toreveal that seizures are a stereotyped behavior of brain networks that can be leveraged to aid the classification and identificationof seizures. This model is an important candidate for use in sensing and actuation design, as well as for new stimulationstrategies to be integrated into a future generation of implantable devices. The model provides fine-scale temporal assessmentof seizure dynamics, and can be extended to model the impact of control strategies on seizure evolution, with the goal of seizureprevention.
Acknowledgements
D.S.B., A.A., and A.N.K. would like to acknowledge support from the John D. and Catherine T. MacArthur Foundation,the Alfred P. Sloan Foundation, the Army Research Laboratory and the Army Research Office through contract numbersW911NF-10-2-0022 and W911NF-14-1-0679, the National Institute of Health (2-R01-DC-009209-11, 1R01HD086888-01,R01-MH107235, R01-MH107703, R01MH109520, 1R01NS099348, R21-M MH-106799, the Office of Naval Research, and theNational Science Foundation (BCS-1441502, CAREER PHY-1554488, BCS-1631550, and CNS-1626008). K.D. acknolwedgessupport from K23NS073801-01, and UH2-NS095495-01. S.P. and G.J.P. are supported in part by the TerraSwarm ResearchCenter, one of six centers supported by the STARnet phase of the Focus Center Research Program (FCRP) a SemiconductorResearch Corporation program sponsored by MARCO and DARPA. The content is solely the responsibility of the authors anddoes not necessarily represent the official views of any of the funding agencies. eferences
Pedram Afshar, Ankit N. Khambhati, Scott Stanslaski, David Carlson, Randy Jensen, Dave Linde, Siddharth Dani, MaciejLazarewicz, Peng Cong, Jon Giftakis, Paul Stypulkowski, and Tim Denison. A translational platform for prototypingclosed-loop neuromodulation systems.
Frontiers in Neural Circuits , 6:1—-15, 2013. ISSN 1662-5110. doi: 10.3389/fncir.2012.00117.Fabiola Alonso, Karin W˚ardell, and Simone Hemm-Ode. Influence on deep brain stimulation from lead design, operating modeand tissue impedance changes–a simulation study.
Brain Disorders and Therapy , 2015.Steven Baldassano, Benjamin Brinkmann, Hoameng Ung, Tyler Blevins, Erin Conrad, Kent Leyde, Mark Cook, AnkitKhambhati, Joost Wagenaar, Gregory Worrell, and Brian Litt. Crowdsourcing Seizure Detection: Algorithm Developmentand Validation on Human Implanted Device Recordings.
In Press: Brain , 2017.Benjamin H Brinkmann, Joost Wagenaar, Drew Abbot, Phillip Adkins, Simone C Bosshard, Min Chen, Quang M Tieng, JialuneHe, FJ Mu˜noz-Almaraz, Paloma Botella-Rocamora, et al. Crowdsourcing reproducible seizure forecasting in human andcanine epilepsy.
Brain , 139(6):1713–1722, 2016.Samuel P. Burns, Sabato Santaniello, Robert B. Yaffe, Christophe C. Jouny, and Nathan E. Crone. Network dynamics of thebrain and influence of the epileptic seizure onset zone.
Proceedings of the National Academy of Sciences of the United Statesof America , 111(49):E5321–5330, nov 2014. ISSN 1091-6490. doi: 10.1073/pnas.1401752111.Marco Curtis and Massimo Avoli. Gabaergic networks jump-start focal seizures.
Epilepsia , 57(5):679–687, 2016.Javier Echauz, Hiram Firpi, and George Georgoulas. Intelligent control strategies for neurostimulation. In
Applications ofintelligent control to engineering systems , pages 247–263. Springer, 2009.Robert S Fisher, J Helen Cross, Jacqueline A French, Norimichi Higurashi, Edouard Hirsch, Floor E Jansen, Lieven Lagae,Solomon L Mosh´e, Jukka Peltola, Eliane Roulet Perez, et al. Operational classification of seizure types by the internationalleague against epilepsy: Position paper of the ilae commission for classification and terminology.
Epilepsia , 58(4):522–530,2017.Viktor K Jirsa, William C Stacey, Pascale P Quilichini, Anton I Ivanov, and Christophe Bernard. On the nature of seizuredynamics.
Brain , 137(8):2210–2230, 2014.M. D. Johnson, H. H. Lim, T. I. Netoff, A. T. Connolly, N. Johnson, A. Roy, A. Holt, K. O. Lim, J. R. Carey, J. L. Vitek, andB. He. Neuromodulation for brain disorders: Challenges and opportunities.
IEEE Transactions on Biomedical Engineering ,60(3):610–624, March 2013. ISSN 0018-9294. doi: 10.1109/TBME.2013.2244890.H.K. Khalil.
Nonlinear Control . Always Learning. Pearson Education, Limited, 2014. ISBN 9781292060507.A N Khambhati, K A Davis, T H Lucas, B Litt, and D S Bassett. Virtual cortical resection reveals push-pull network controlpreceding seizure evolution.
Neuron , 91(5):1170–1182, 2016.Ankit N. Khambhati, Kathryn A. Davis, Brian S. Oommen, Stephanie H. Chen, Timothy H. Lucas, Brian Litt, and Danielle S.Bassett. Dynamic Network Drivers of Seizure Generation, Propagation and Termination in Human Neocortical Epilepsy.
PLOS Computational Biology , 11(12):e1004608, dec 2015. ISSN 1553-7358. doi: 10.1371/journal.pcbi.1004608.Ankit N. Khambhati, Danielle S. Bassett, Brian S. Oommen, H. Chen, Stephanie, Timothy H. Lucas, Kathryn A. Davis, andBrian Litt. Recurring functional interactions predict network architecture of interictal and ictal states in neocortical epilepsy. eNeuro , 4(1), 2017. ISSN 2373-2822. doi: 10.1523/ENEURO.0091-16.2017.JW Kim, JA Roberts, and PA Robinson. Dynamics of epileptic seizures: evolution, spreading, and suppression.
Journal oftheoretical biology , 257(4):527–532, 2009.Matthias J Koepp, Lorenzo Caciagli, Ronit M Pressler, Klaus Lehnertz, and S´andor Beniczky. Reflex seizures, traits, andepilepsies: from physiology to pathology.
The Lancet Neurology , 15(1):92–105, 2016.Mark A. Kramer, Uri T. Eden, Eric D. Kolaczyk, Rodrigo Zepeda, Emad N. Eskandar, and Sydney S. Cash. Coalescence andfragmentation of cortical networks during focal seizures.
The Journal of neuroscience : the official journal of the Society forNeuroscience , 30(30):10076–10085, jul 2010. ISSN 0270-6474. doi: 10.1523/JNEUROSCI.6309-09.2010. rian Litt, Rosana Esteller, Javier Echauz, Maryann D’Alessandro, Rachel Shor, Thomas Henry, Page Pennell, Charles Epstein,Roy Bakay, Marc Dichter, et al. Epileptic seizures may begin hours in advance of clinical onset: a report of five patients.
Neuron , 30(1):51–64, 2001.Marcelo O Magnasco, Oreste Piro, and Guillermo A Cecchi. Self-tuned critical anti-hebbian networks.
Physical review letters ,102(25):258102, 2009.Florian Mormann, Ralph G Andrzejak, Christian E Elger, and Klaus Lehnertz. Seizure prediction: the long and winding road.
Brain , 130(2):314–333, 2007.Martha Morrell. Brain stimulation for epilepsy: can scheduled or responsive neurostimulation stop seizures?
Current opinionin neurology , 19(2):164–168, 2006.Martha J. Morrell. Responsive cortical stimulation for the treatment of medically intractable partial epilepsy.
Neurology , 77(13):1295–1304, 2011. ISSN 00283878. doi: 10.1212/WNL.0b013e3182302056.Dileep R. Nair, Armin Mohamed, Richard Burgess, and Hans L¨uders. A critical review of the different conceptual hypothesesframing human focal epilepsy.
Epileptic disorders : international epilepsy journal with videotape , (2):77–83, 2004. ISSN1294-9361.Theoden I Netoff, Louis M Pecora, and Steven J Schiff. Analytical coupling detection in the presence of noise and nonlinearity.
Physical Review E , 69(1):017201, 2004.Arnold Neumaier and Tapio Schneider. Estimation of parameters and eigenmodes of multivariate autoregressive models.
ACMTransactions on Mathematical Software (TOMS) , 27(1):27–57, 2001.Ivan Osorio, Mark G Frei, Sridhar Sunderam, Jonathon Giftakis, Naresh C Bhavaraju, Scott F Schaffner, and Steven BWilkinson. Automated seizure abatement in humans using electrical stimulation.
Annals of neurology , 57(2):258–268, 2005.S Pequito, A Ashourvan, D S Bassett, B Litt, and G J Pappas. Spectral control of cortical activity.
Proceedings of the 2017American Control Conference , 2017, 2017.Piero Perucca, Franc¸ois Dubeau, and Jean Gotman. Intracranial electroencephalographic seizure-onset patterns: effect ofunderlying pathology.
Brain , 137(1):183–196, 2014.Asla Pitk¨anen, Wolfgang L¨oscher, Annamaria Vezzani, Albert J Becker, Michele Simonato, Katarzyna Lukasiuk, Olli Gr¨ohn,Jens P Bankstahl, Alon Friedman, Eleonora Aronica, et al. Advances in the development of biomarkers for epilepsy.
TheLancet Neurology , 15(8):843–856, 2016.Christian Rummel, Marc Goodfellow, Heidemarie Gast, Martinus Hauf, Fr´ed´erique Amor, Alexander Stibal, Luigi Mariani,Roland Wiest, and Kaspar Schindler. A systems-level approach to human epileptic seizures.
Neuroinformatics , 11(2):159–73,apr 2013. ISSN 1559-0089. doi: 10.1007/s12021-012-9161-2.Catherine A Schevon, Shennan A Weiss, Guy McKhann Jr, Robert R Goodman, Rafael Yuste, Ronald G Emerson, and Andrew JTrevelyan. Evidence of an inhibitory restraint of seizure activity in humans.
Nature communications , 3:1060, 2012.Steven J Schiff.
Neural control engineering: the emerging intersection between control theory and neuroscience . MIT Press,2012.Steven J Schiff, Tim Sauer, Rohit Kumar, and Steven L Weinstein. Neuronal spatiotemporal pattern discrimination: thedynamical evolution of seizures.
Neuroimage , 28(4):1043–1055, 2005.Kaspar Schindler, Howan Leung, Christian E Elger, and Klaus Lehnertz. Assessing seizure dynamics by analysing thecorrelation structure of multichannel intracranial eeg.
Brain , 130(1):65–77, 2006.Kaspar Schindler, Howan Leung, Christian E Elger, and Klaus Lehnertz. Assessing seizure dynamics by analysing thecorrelation structure of multichannel intracranial eeg.
Brain , 130(1):65–77, 2007.Donald L Schomer and Fernando Lopes Da Silva.
Niedermeyer’s electroencephalography: basic principles, clinical applications,and related fields . Lippincott Williams & Wilkins, 2012.F Spanedda, F Cendes, and J Gotman. Relations between eeg seizure morphology, interhemispheric spread, and mesial temporalatrophy in bitemporal epilepsy.
Epilepsia , 38(12):1300–1314, 1997. illiam C Stacey and Brian Litt. Technology insight: neuroengineering and epilepsy-designing devices for seizure control.
Nature clinical practice. Neurology , 4(4):190–201, 2008a. ISSN 1745-834X. doi: 10.1038/ncpneuro0750.William C Stacey and Brian Litt. Technology insight: neuroengineering and epilepsy?designing devices for seizure control.
Nature Clinical Practice Neurology , 4(4):190–201, 2008b.Scott Stanslaski, Pedram Afshar, Peng Cong, Jon Giftakis, Paul Stypulkowski, Dave Carlson, Dave Linde, Dave Ullestad,Al Thaddeus Avestruz, Timothy Denison, Peng Cong, Jon Giftakis, Paul Stypulkowski, Dave Carlson, Dave Linde, DaveUllestad, Al Thaddeus Avestruz, and Timothy Denison. Design and validation of a fully implantable, chronic, closed-loopneuromodulation device with concurrent sensing and stimulation.
IEEE Transactions on Neural Systems and RehabilitationEngineering , 20(4):410–421, jul 2012. ISSN 15344320. doi: 10.1109/TNSRE.2012.2183617.Joost B. Wagenaar, Benjamin H. Brinkmann, Zachary Ives, Gregory A. Worrell, and Brian Litt.
A multimodal platform forcloud-based collaborative research , pages 1386–1389. 2013. ISBN 9781467319690. doi: 10.1109/NER.2013.6696201.Christopher P. Warren, Sanqing Hu, Matt Stead, Benjamin H. Brinkmann, Mark R. Bower, and Gregory A. Worrell. Synchronyin normal and focal epileptic brain: the seizure onset zone is functionally disconnected.
Journal of neurophysiology , 104(6):3530–3539, dec 2010. ISSN 0022-3077. doi: 10.1152/jn.00368.2010.Shennan A Weiss, Garrett P Banks, Guy M McKhann Jr, Robert R Goodman, Ronald G Emerson, Andrew J Trevelyan, andCatherine A Schevon. Ictal high frequency oscillations distinguish two types of seizure territories in humans.
Brain , 136(12):3796–3808, 2013.Christopher Wilke, Gregory Worrell, and Bin He. Graph analysis of epileptogenic networks in human partial epilepsy.
Epilepsia ,52(1):84–93, jan 2011. ISSN 00139580. doi: 10.1111/j.1528-1167.2010.02785.x.Gregory A. Worrell, Stephen D. Cranstoun, Javier Echauz, and Brian Litt. Evidence for self-organized criticality in humanepileptic hippocampus.
Neuroreport , 13(16):2017–21, 2002. ISSN 0959-4965.L Wu and J Gotman. Segmentation and classification of eeg during epileptic seizures.
Electroencephalography and clinicalneurophysiology , 106(4):344–356, 1998.Drausin F Wulsin, Emily B. Fox, Brian Litt, and Emily B. Fox. Parsing Epileptic Events Using a Markov Switching ProcessModel for Correlated Time Series-Supplementary Materials.
ICML2013 , pages 1–15, 2013. upplementary Materials for“Parsing spatiotemporal dynamical stability in ECoG during seizure onset,propagation, and termination”
Supplementary Methods
In this supplementary methods section, we provide complementary information that aims to sharpen the intuition behind thedynamical stability-based approach exercised in the main text. First, we present an application of the dynamical stability-based approach to synthetic data. Second, we show that high-order linear models do not qualitatively change the proposedseizure-onset characterization.
Synthetic example
Here we provide a pedagogical example of the dynamical stability-based approach, and demonstrate its ability to accuratelyuncover a planted locus, frequency, and damping of oscillatory dynamics. The synthetic dataset that we study consists of fourchannels associated with four oscillatory sources, sampled at 512 Hz. Each source is modeled as an enveloped (exponential,Gaussian) sinusoid, and the four sources have a fixed lag between them to induce directionality of information flow acrosssensors. Specifically, the baseline activity of each channel is modeled as follows: y = sin ( ax + b ) e cx where a determines the frequency of the oscillation, b captures the phase lag between electrodes, and c determines the dampingof the oscillations. Negative values of c were used for early damping oscillations while positive values of c were used for lateemerging high frequency oscillations. In addition to this baseline activity, we also modeled transient activity that emerges anddisappears in brief intervals during the recordings. We model these transient dynamics as follows: y = sin ( ax + b ) √ σ π e − ( x − µ ) σ where µ is the mean and σ the standard deviation of a normal distribution. Finally, different simulated patterns of oscillatoryactivities over time were linearly superimposed to create the activity of individual channels. Additional details regarding thisprocess of sensor time series creation is depicted schematically in Fig. 4A.After constructing the sensor time series, we next applied the dynamic stability approach to determine its ability to accuratelyuncover the planted structure of the oscillatory dynamics. First, we show the eigenvector evolution associated with the highestfrequency dynamics in the signals (Fig. 4B), and then we show the eigenvector evolution associated with dynamics in a lowerfrequency (Fig. 4C). We observe that the locus of higher and lower frequency oscillations at every time point can be identifiedvia electrodes with relatively higher values in the eigenvectors 1 and 3, respectively. We also display the stability, | λ | associatedwith those same two eigenvectors (Fig. 4D-E). We notice that the changes in the angle match the changes of the underlyingfrequency of oscillations, while the intervals with radius ≈
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 181234
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 181234
32 Hz64 Hz1 Hz8 Hz 2 Hz
Frequency HZ -6-5-4-3-2-1 0 1 2 3 4 5 6 t au l og10 Windowsize SourceElectrods
A.B.C.D.E. E l e c t r o d e s E l e c t r o d e s F.G.
Time
ReIm ReIm
H. I. A n g l e | θ | R a d | λ | E i g e n v e c t o r E i g e n v e c t o r (Hz) E i g e n v a l u e E i g e n v a l u e Figure 4.
Dynamical stability-based characterization of synthetic time-series . (A) We constructed surrogate datacomposed of time-series that capture an underlying process: an enveloped (exponential, Gaussian) sinusoid. Briefly, thesimulated time-series capture various partially-overlapping spatiotemporal patterns of oscillations with varying frequencies.The colored bar at the top of the panel represents the duration over which the oscillation with that specific frequency is present.We used a 2-STD interval for oscillations with a Gaussian amplitude profile, and we used 1 − π /
2. The yellow shaded area on the right highlights the size on the moving window. (B-C)
The absolute values of the normalized eigenvectors 1 and 3 obtained using a window of 1 s duration moved by 100 msshifts. (D-E)
The evolution of the angle and radius associated with the four eigenvectors of the 4-sensor system. (F)
Thetemporal progression of the sliding window illustrated in changing hues of red; the yellow box represents the relative size of themoving window used to estimate the linear model. (G)
The temporal evolution of the frequency (Hz) and damping ( log ( τ ) )associated with the first eigenmode (‘dot’) and third eigenmode (‘star’). (H) The temporal evolution of the four eigenvalues inthe Argand complex plane. (I)
A zoomed-in version of panel (H) to show the changes in the angle and radius of the low angleeigenvalues in greater detail. Panels (G-I) use the same color-coding as in panel (F) . .25 0.5 1 2 4 8 16 32 64 128 Frequency (Hz) -6-5-4-3-2-1 0 1 2 3 4 5 6 t au l og10 Frequency (Hz) -6-5-4-3-2-1 0 1 2 3 4 5 6 t au l og10 Frequency (Hz) t au l og10 ReIm E. ReIm ReIm
F. G.H. I. J. W i n d o w s i z e ( S e c ) W i n d o w s i z e B. Time T a u ( L o g ) A.
2 4 6 8 10 12 14 16 181234
2 4 6 8 10 12 14 16 181234
2 4 6 8 10 12 141234
2 4 6 8 10 12 141234 E l e c t r o d s E l e c t r o d s C. D. E i g e n v e c t o r E i g e n v e c t o r Seconds SecondsSeconds
Figure 5.
Effect of window size on the dynamical stability-based characterization of synthetic time-series. (A)
Different sliding windows for the synthetic data shown in Fig. 4A are encoded by the yellow box, and the corresponding timeprogression is color-coded as in Fig. 4F. (B-D)
The absolute values of the normalized eigenvectors 1 and 3 obtained using awindow of 1, 5, and 10 s duration moved by 100 ms shifts (E-G)
The Argand diagram of the temporal evolution of the foureigenvalues for 1, 5, and 10 s windows, respectively. Insets: Zoomed-in version of the same plots to show the changes in theangle and radius of the low angle eigenvalues in greater detail. (H-J)
The temporal evolution of the frequency (Hz) anddamping ( log ( τ ) ) associated with the first eigenmodes (‘dots’) and third eigenmodes (‘stars’) associated with 1, 5, and 10 swindows. Panels (E-J) use the same color-coding as in panel (A) . ynamical stability-based approach using higher-order linear models The dynamical stability-based approach relies on a first-order approximation of a non-linear process, and therefore onlyprovides a local approximation of the underlying dynamics. Although computationally efficient and simple, it is importantto consider the drawbacks of this approach in comparison to approaches that employ higher-order models. Because the dataconsidered is in discrete-time and the linear model used can be understood as a coupled autoregressive model of first-order,we assess how a coupled autoregressive model of higher-order might impact the observed results. Towards this goal, we usedthe Schwarz’s Bayesian criterion (Neumaier and Schneider, 2001) to estimate the order of the coupled autoregressive model.For a 1 s time-window, the order of the coupled autoregressive model is generally lower than three, across time and subjects.Therefore, we contrast the results obtained with a coupled autoregressive of order two model (Fig. 6 in this Supplement) againstthose reported in Fig. 1 of the main manuscript.We observe notable differences between the dynamics of the two models, as expected from a theoretical point-of-view:higher-order models have extended memory that enables the monitoring of the spatiotemporal lag in the evolution of theunderlying dynamical process. First, the angles of the high frequency eigenmodes of the second-order autoregressive model aremuch higher than those of the first-order autoregressive model (compare Fig. 6A to Fig. 1A). Second, the eigenvectors of thefirst-order autoregressive model suggest that the second period of high frequency dynamics following the EEC onset is onlypresent between 27-33 s, while the eigenvectors of the second-order autoregressive model suggest that these dynamics extendfrom 27-80 s (compare Fig. 1A and Fig. 6A). These results support the presence of temporally overlapping high frequencydynamics that cannot be captured by the first-order autoregressive model. Despite these differences, both models presenteigenvectors associated with high frequencies that exhibit approximately time-invariant loading patterns at seizure onset andtermination. Moreover, they both show the contribution of only a few regions for the dynamic process associated with thespatiotemporal high frequencies at seizure onset, and at seizure termination they present a sudden transition to capture a largernumber of contributions with almost equal weight across different neocortical regions. These similarities between the twomodels corroborate the capability of the proposed dynamical stability-based approach to characterize seizure dynamics. A. E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r B.D. | θ | ( R a d ) | θ | ( R a d ) | θ | ( R a d ) E l e c t r o d e s E l e c t r o d e s E l e c t r o d e s E l e c t r o d e s | θ | ( R a d ) C. |λ||λ||λ||λ| Figure 6.
Evolution of eigenmodes estimated from a second order autoregressive model over the ictal period shownin Fig. 1A. (A-D)
The frequency (cyan trace) and stability (red trace) of eigenvalues associated with four representativeeigenvectors – two associated with high frequencies and two associated with low-frequencies – whose evolution is displayed inheatmaps. Naturally, these models are better able to capture spatiotemporal high frequency patterns. For instance, the period ofapproximately time-invariant eigenvector loading patterns shown in Fig. 1B between 27-33 s is here marked by high frequencytransient oscillatory activity present across a large number of electrodes, and extending over a long time period. upplementary Results
Seizure offset is marked by a sudden decrease in the stability of several egienmodes
Almost all ictal samples (with the exception of SP) reveal a similar decrease in the stability of several eigenmodes followingictal offset. The average stability of all eigenmodes at ictal offset – which we plot for all CP and CP+GTC samples (samplesfrom subjects 14,15, and 21 are not presented due to insufficient clinical annotations) in Fig. 7A (top) – highlights the suddenchange in the stability that is time-locked to the seizure offset. Although the magnitude of the decrease in the stability is muchlarger in most GTC samples (as seen in Fig. 7A&C), CP samples show comparable changes though at smaller magnitude. Thissimilarity is better captured in the normalized ( z -score) event-related average plots shown in Fig. 7B. Although the offset drop inthe stability is not limited to a single frequency band, calculating the difference between the average stability over a 3 s windowbefore and after ictal offset reveals that all CP and CP+GTC samples (after removing one CP+GTC sample with high recordingnoise) show a large offset stability drop in the β -band. These observations suggest that the estimated dynamics immediatelyfollowing ictal offset are more asymptotically stable – i.e., the post-ictal oscillations damp much faster than pre-offset (ictal)oscillations. Impact of the time-window length in the dynamical stability-based approach
Here we ask whether the results of the dynamic stability approach depend on the time window over which we estimate theARMA model. Our intuition is that if the time window is too large, we might not be able to accurately specify short time-scalechanges. We therefore consider 1, 2, 5 and 10 s time-windows, while maintaining 100 msec shifts. The calculations providesimilar intuitions to those we obtained from considering the synthetic data: we observe a loss of sensitivity to high frequencydynamics as the length of the time-window increases (Fig. 8). Despite these window-dependent effects, our main resultsstill hold: (i) the eigenvectors associated with high frequencies present approximate time-invariance of loading patterns for arelatively long period of time, with the contribution of a few regions at the onset of the seizure, and (ii) there is a sudden changeat the end of the seizure where a larger number of regions contribute in an approximately equal manner to the evolution of thespatiotemporal high frequency process. These consistent findings bolster confidence in the reported results. inter-ictal dynamical stability-based characterization
While our study focuses on ictal dynamics, it is also useful to gain an intuition for the dynamics observed during the inter-ictalperiod, as a natural comparator to the ictal state. We therefore apply the dynamical stability-based characterization to fourdifferent inter-ictal periods from subject 3. We observe that the angles of the eigenvalues display quite small fluctuations aroundthe mean (Fig. 9), in contrast to the large fluctuations observed over the ictal period (see Fig. 1 and Fig. 2). We do not observelong periods of approximately time-invariant eigenvector loading patterns with strong contribution from only a few regions.Instead, we observe high temporal variance in eigenvector loading patterns, and relatively homogeneous contributions from allelectrodes.
Comparison to clinical annotations
Seizure onset is marked by the emergence of slow-changing and approximately time-invariant dynamics, commonly with astereotypical spatial focus. Here, we provide examples that show that these observations accurately echo the clinical annotations,where seizure onset is associated with the emergence of high frequency activity lasting for a few seconds. We begin by studyingthe evolution of sample eigenmodes of a single patient over the ictal period presented in Fig. 1. Here, the clinical assessment ofthe time-series identifies the EEC at electrode 79, marked by a sharp wave followed by low voltage fast activity lasting ≈ ≈ ± ± Statistical testing demonstrating spatially constrained and slow changing dynamics
Here, we provide quantitative analyses supporting the relatively time-invariant nature of seizure onset dynamics in comparisonto inter-ictal dynamics. A time-invariant linear system, when estimated using a moving autoregressive model, will ideally yieldthe same A matrix and, consequently, the very same eigenvectors at every time-point. Nevertheless, in practice, physiologicaland measurement noise can lead to fluctuations of the estimated system.We first sought to assess the extent of the time-invariance of the eigenvectors over the ictal period at the group-level usingnon-parametric statistical testing. Specifically, we first defined summary statistic of interest as the sum of the element-wise - s c o r e | λ | z - s c o r e | λ | z - s c o r e | λ | z - s c o r e | λ | z - s c o r e | λ | z - s c o r e | λ || λ | - b a s e l i n e | λ | - b a s e l i n e | λ | - b a s e l i n e | λ | - b a s e l i n e | λ | - b a s e l i n e | λ | - b a s e l i n e C o u n t ( G T C ) C o u n t ( C P ) C o u n t ( G T C ) C o u n t ( C P ) C o u n t ( G T C ) C o u n t ( C P ) C o u n t ( C P ) C o u n t ( C P ) C o u n t ( C P ) A. B. C.
Seconds Seconds
Δ |λ| C o u n t ( G T C ) C o u n t ( G T C ) C o u n t ( G T C ) Figure 7.
Seizure offset is marked by a sudden decrease in the stability of several egienmodes. (A) (Top) The averagestability of all eigenmodes over a 40 s window centered around the ictal offset, and plotted for all CP (dashed line) andCP+GTC (solid line) samples, after removing a single CP+GTC sample with high recording noise and all samples fromsubjects 14,15, and 21 (insufficient clinical annotations). The y-axis shows the average stability values following thesubtraction of a 20 s pre-offset baseline. Note that large drops in the stability are present in all CP+GTC samples. The resultsfor the eigenmodes, after band-pass filtering based on the eigenmodes’ average frequency are presented from top to bottom toreflect δ , θ , α , β , and γ . (B) Average stability plots normalized calculating the z -score. The dashed black line represents theevent-related average of samples. (C) Histogram of the difference between the average stability over a 3 s window before andafter ictal offset. The blue histogram represents the CP+GTC results and the red histogram represents the CP results. Note thatall samples show higher average pre-offset average stability across all eigenmodes (Top). Although this effect is not unique toan individual frequency band, the most robust results are observed in the β -band. .B. θ E l e c t r o d e s E l e c t r o d e s E l e c t r o d e s E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r C.D.E. E i g e n v a l u e E i g e n v a l u e ( R a d ) | λ | Seconds
Figure 8.
Effect of window size on the dynamical stability-based characterization of ictal periods. (A-B)
For theexample seizure from subject 16 shown in Fig. 3A, we show the time evolution of the angles and stability of the eigenvaluesassociated with the different eigenmodes, for different window sizes (1, 2, 5, and 10 s; shown across different columns). Thered and cyan arrows highlight the onset and offset of the ictal period, respectively. (C-E)
The temporal evolution of threerepresentative eigenvectors – two associated with high frequencies and one associated with low-frequencies. . E l e c t r o d e s B.C.D.E. | θ | ( R a d ) | θ | ( R a d ) | θ | ( R a d ) | θ | ( R a d ) |λ| E l e c t r o d e s E l e c t r o d e s E l e c t r o d e s E i g e n v e c t o r E l e c t r o d e s |λ| |λ| |λ| E i g e n v e c t o r E i g e n v e c t o r E i g e n v e c t o r Sample
10 50 9010 50 9010 50 9010 50 9010 50 90 10 50 9010 50 9010 50 9010 50 9010 50 9010 50 9010 50 9010 50 9010 50 9010 50 9010 50 9010 50 9010 50 9010 50 9010 50 9010203040506070 10203040506070 1020304050607010203040506070102030405060701020304050607010203040506070 Seconds 102030405060701020304050607010203040506070102030405060701020304050607010203040506070102030405060701020304050607010203040506070 102030405060701020304050607010203040506070
Sample
Figure 9.
Evolution of sample eigenmodes for four sample inter-ictal periods of subject 3. (A)
ECoG signals from 79electrodes in subject 3 over four sample ictal periods (across the columns), from 20 s before seizure onset to 20 s after seizuretermination. (B-E)
The frequency (cyan trace) and stability (red trace) of eigenvalues associated with four representativeeigenvectors – two associated with high frequencies and two associated with low-frequencies – whose evolution is displayed inheatmaps. . B.C. D.
Eigenvector E l e c t r o d e s E l e c t r o d e s Eigenvector
Figure 10.
Evolution of sample eigenvectors of a subject over an ictal period. (A)
ECoG signals from 79 electrodes insubject 3 over a sample ictal period, from 20 s before seizure onset to 20 s after seizure termination, and applying a high-passfilter of 0 . (B) The same time-series after applying a high-pass filter at 60 Hz, highlighting the emergence of highfrequency activity at seizure onset. The clinical assessment identifies electrode 79 as the EEC, marked by a sharp wavefollowed by low voltage fast activity lasting approximately 3 s. The EEC is immediately followed by high frequency activity atelectrodes 11-14 and later followed by electrodes 3-4, and 53-54 (lasting approximately 7 s). These time-series are highlightedin red. (C-D)
The time evolution of one low and one high frequency eigenvector. Note that the same electrodes identifiedclinically as the seizure onset zone are also identifiable from the evolution of the eigenvectors (as indicated by red arrows). quared differences between temporally contiguous eigenvectors. Then, we calculated the average standard deviation in thisvalue across all temporally contiguous pairs of eigenvectors across an ictal event. Intuitively, this summary statistic quantifiesthe temporal variance of the given eigenvector across time: high values of this statistic indicate large temporal fluctuations,and small values of this statistic indicate small temporal fluctuations. To examine this temporal variance at the group levelacross frequency bands, we grouped the eigenmodes for each patient’s dataset based on the average estimated frequency ofthe eigenmodes into several frequency bands: delta (0-4 Hz), theta (4-8 Hz), alpha (8-12 Hz), beta (12-30 Hz), and gamma(30-55 Hz). We systematically removed the so-called relaxator eigenmodes, which are eigenmodes that are non-oscillatoryfor more than 50 % of the sample duration. Eigenmodes with average frequencies greater than 55 Hz were removed from theanalysis because only a handful of patient recordings displayed this high frequency dynamics. Next, we calculated subject-levelsummary statistics, averaged over all eigenvectors in each band, for (i) the ictal samples, averaged of all samples for eachpatient, and (ii) the inter-ictal samples, averaged of all samples for each patient. Then, we compared these summary statisticsbetween ictal and inter-ictal samples at the group-level via bootstrap analysis. To create non-parametric null distributions ofthe group-level average of the difference between the ictal and inter-ictal summary statistics computed at the subject-level,we permuted the labels “ictal” and “inter-ictal” uniformly at random 50000 times. Finally, we tested the hypothesis that theeigenvectors in the ictal samples at each frequency band on average displayed significantly less temporal fluctuations thaninter-ictal samples across subjects. We observed significantly less eigenvector fluctuations across several frequency bands overthe ictal samples in comparison to the inter-ictal samples ( p < .
05 by permutation testing with N =
50 000; Fig. 11). Weobserve that the difference between the ictal and inter-ictal samples was clear even for the 1 s window following seizure onsetin high frequency γ . However, this effect is best captured with larger window sizes for eigenmodes with slower dynamics.Following the examination of the time-dependence versus time-invariance of eigenvectors during the ictal period, andthe frequency dependence of our findings as a function of ictal versus inter-ictal samples, we turned to an assessment of theregional localization of the effects. To measure regional localization, we studied the absolute magnitude of the maximallycontributing elements of the eigenvectors. Then, we used a bootstrap analysis to test the hypothesis that the emergence of theslow-changing focal ictal dynamics was associated with an increase in the absolute magnitude of the maximally contributingelements of eigenvectors, particularly within a short time window (1, 3, 6, and 9 s) following the seizure onset. We observedsystematic increases in the average maximum absolute value of the eigenvector elements following seizure onset across severalfrequency bands (Fig. 12). Interestingly, the high frequency γ did not show any increase while α and β bands did whenconsidering periods longer than 6 s following the seizure onset. Finally, to determine the specificity of the effects to the ictalperiod, we compared the subject-level summary statics in the ictal segments to those estimated from the inter-ictal and pre-ictalsegments. The comparison between the ictal and preictal samples revealed analogous results as seen in (Fig. 13). Together, ourcomplementary analyses reveal the relatively slow-changing nature of seizure onset dynamics across subjects. p = 0.06734 -0.01 -0.005 0 0.005 0.010500100015002000 p = 0.00098 -0.02 -0.01 0 0.01 0.0205001000150020002500 p = 0.00026 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.00128 -0.15 -0.1 -0.05 0 0.05 0.10500100015002000 p = 2e-05 -0.02 -0.01 0 0.01 0.020500100015002000 p = 0.37648 -0.01 -0.005 0 0.005 0.010500100015002000 p = 0.067 -0.02 -0.01 0 0.01 0.020500100015002000 p = 0.32238 -0.01 -0.005 0 0.005 0.01050010001500 p = 0.04014 -0.02 -0.01 0 0.01 0.020500100015002000 p = 0.22138 -0.02 -0.01 0 0.01 0.0205001000150020002500 p = 0.03582 -0.04 -0.02 0 0.02 0.040500100015002000 p = 0.05642 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.0149 -0.1 -0.05 0 0.05 0.10500100015002000 p = 2e-05 -0.15 -0.1 -0.05 0 0.05 0.10500100015002000 p = 2e-05 - H z - H z - H z - H z Group average of the difference of the standard deviation of eigenvectors’ fluctuations between each subjects’ ictal and interictal samples. c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t A. B. C. - H z c o u n t c o u n t c o u n t -0.02 -0.01 0 0.01 0.020500100015002000 p = 0.01462 -0.02 -0.01 0 0.01 0.0205001000150020002500 p = 0.00018 -0.02 -0.01 0 0.01 0.02050010001500 p = 2e-05 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 2e-05 -0.1 -0.05 0 0.05 0.10500100015002000 p = 2e-05 c o u n t c o u n t c o u n t c o u n t c o u n t D. * ** **** *** ** Figure 11.
Approximate spatio-temporal invariance of eigenvector loading patterns at seizure onset.
Here we plot thenull distributions for the group average of the summary statistic of interest: taking the sum of the element-wise squareddifferences between temporally contiguous eigenvectors, we then calculated the average standard deviation in this value acrossall temporally contiguous pairs of eigenvectors across an ictal event. (A-D)
Columns indicate different time window sizes (1, 3,6 and 9 s); rows indicate groupings of eigenmodes into several frequency bands based on their average frequency over samples.The null distribution (shown in blue) was generated by shuffling the ictal and inter-ictal subject-level means and re-calculatingthe group-level average of subject-level difference between the subjects’ ictal and inter-ictal sample means. The dashed greenlines show the 95% confidence interval ( N =
50 000 permutations) and the red vertical bar shows the empirical group averageof the difference between the subjects’ ictal and inter-ictal sample means. Plots with significant results are marked by a red staron the empirical group average. p = 0.22034 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.37518 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.25468 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.0596 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.30316 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.19456 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.05986 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.01348 -0.04 -0.02 0 0.02 0.04 0.06050010001500 p = 0.31826 -0.04 -0.02 0 0.02 0.04050010001500 p = 0.03476 -0.06 -0.04 -0.02 0 0.02 0.04050010001500 p = 0.00094 -0.05 0 0.05050010001500 p = 2e-05 -0.1 -0.05 0 0.05 0.1050010001500 p = 0.2304 -0.1 -0.05 0 0.05 0.1050010001500 p = 0.04008 -0.1 -0.05 0 0.05 0.10500100015002000 p = 0.00154 -0.1 -0.05 0 0.05 0.10500100015002000 p = 8e-05 -0.4 -0.2 0 0.2 0.40500100015002000 p = 0.19 -0.4 -0.2 0 0.2 0.40500100015002000 p = 0.25356 -0.4 -0.2 0 0.2 0.4050010001500 p = 0.2578 -0.4 -0.2 0 0.2 0.4050010001500 p = 0.24616 - H z - H z - H z - H z Average difference between individual patient’s ictal and interictal samples’ maximum (window average) eigenvector’s values per timepoint. c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t A. B. C. - H z c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t D. * *** * Figure 12.
Increased eigenvector loading following seizure onset compared to the interical period.
Here we plot thenull distributions for the group average obtained as a result of computing the difference between the average maximum entry ofthe eigenvectors following the seizure onset and that of the inter-ictal samples at the subject-level. (A-D)
Columns indicatedifferent time window sizes (1, 3, 6 and 9 s); rows indicate groupings of eigenmodes into several frequency bands based ontheir average frequency over samples. The null distribution was generated by shuffling the ictal and inter-ictal subject-levelmeans and re-calculating the group-level average of subject-level differences between the sample means. The dashed greenlines show the 95% confidence interval ( N =
50 000 permutations) and the red vertical bar shows the empirical group averageof the difference between the subjects’ ictal and inter-ictal sample means. Plots with significant results are marked by a red staron the empirical group average. p = 0.22258 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.3744 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.25344 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.06112 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.30094 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.19082 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.05886 -0.04 -0.02 0 0.02 0.0405001000150020002500 p = 0.01382 -0.05 0 0.05050010001500 p = 0.31686 -0.04 -0.02 0 0.02 0.04050010001500 p = 0.03458 -0.04 -0.02 0 0.02 0.04 0.06050010001500 p = 0.00114 -0.05 0 0.05050010001500 p = 6e-05 -0.1 -0.05 0 0.05 0.1050010001500 p = 0.23002 -0.1 -0.05 0 0.05 0.1050010001500 p = 0.0401 -0.1 -0.05 0 0.05 0.10500100015002000 p = 0.00172 -0.1 -0.05 0 0.05 0.10500100015002000 p = 0.00012 -0.4 -0.2 0 0.2 0.40500100015002000 p = 0.18828 -0.4 -0.2 0 0.2 0.40500100015002000 p = 0.2515 -0.4 -0.2 0 0.2 0.40500100015002000 p = 0.25562 -0.4 -0.2 0 0.2 0.4050010001500 p = 0.24236 - H z - H z - H z - H z Average difference between individual patient’s ictal and perictal samples’ maximum (window average) eigenvector’s values per timepoint. c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t A. B. C. - H z c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t c o u n t D. ** * ** Figure 13.
Increased eigenvector loading following seizure onset compared to the preictal period.
Here we plot the nulldistributions for the group average obtained as a result of computing the difference between the average maximum entry of theeigenvectors following the seizure onset and that of the preictal samples at the subject-level. (A-D)
Columns indicate differenttime window sizes (1, 3, 6 and 9 s); rows indicate groupings of eigenmodes into several frequency bands based on their averagefrequency over samples. The null distribution was generated by shuffling the ictal and preictal subject-level means andre-calculating the group-level average of subject-level differences between the sample means. The dashed green lines show the95% confidence interval ( N =
50 000 permutations) and the red vertical bar shows the empirical group average of thedifference between the subjects’ ictal and preictal sample means. Plots with significant results are marked by a red star on theempirical group average. Note the similarity between these results and those presented in Fig. 12.50 000 permutations) and the red vertical bar shows the empirical group average of thedifference between the subjects’ ictal and preictal sample means. Plots with significant results are marked by a red star on theempirical group average. Note the similarity between these results and those presented in Fig. 12.