Conex-Connect: Learning Patterns in Extremal Brain Connectivity From Multi-Channel EEG Data
CConex-Connect: Learning Patterns inExtremal Brain Connectivity From
Multi-Channel EEG Data
Matheus B Guerrero, Rapha¨el Huser, and Hernando Ombao
King Abdullah University of Science and Technology (KAUST)Statistics Program, CEMSE Divison
January 26, 2021
Abstract
Epilepsy is a chronic neurological disorder affecting more than 50 million peopleglobally. An epileptic seizure acts like a temporary shock to the neuronal system,disrupting normal electrical activity in the brain. Epilepsy is frequently diagnosedwith electroencephalograms (EEGs). Current methods study the time-varying spec-tra and coherence but do not directly model changes in extreme behavior. Thus, wepropose a new approach to characterize brain connectivity based on the joint tailbehavior of the EEGs. Our proposed method, the conditional extremal dependencefor brain connectivity (Conex-Connect), is a pioneering approach that links the as-sociation between extreme values of higher oscillations at a reference channel withthe other brain network channels. Using the Conex-Connect method, we discoverchanges in the extremal dependence driven by the activity at the foci of the epilepticseizure. Our model-based approach reveals that, pre-seizure, the dependence is no-tably stable for all channels when conditioning on extreme values of the focal seizurearea. Post-seizure, by contrast, the dependence between channels is weaker, and de-pendence patterns are more “chaotic”. Moreover, in terms of spectral decomposition,we find that high values of the high-frequency Gamma-band are the most relevantfeatures to explain the conditional extremal dependence of brain connectivity.
Keywords:
Epilepsy, Extreme-value theory, Conditional extremes, Penalized likelihood,Non-stationary time series. 1 a r X i v : . [ q - b i o . N C ] J a n Introduction
Electroencephalograms (EEGs) are multidimensional spatio-temporal signals that measurebrain electrical activity from electrodes placed on the scalp. EEGs capture changes in brainsignals following a shock to the neuronal system, such as an external stimulus, stroke,or epilepsy. These shocks show a profound impact on the neuronal system, includingchanges in frequency content, wave amplitudes and connectivity structure of the networkof signals. Often these shocks have an impact on the distributions (in particular, at thetails). Here, our goal is to develop a new statistical approach for investigating changes in theextremal dependence structure, i.e., the dependence structure prevailing in the tail of thedistribution, between EEG channels during an epileptic seizure. While classical methodsmostly rely on the behavior in the bulk (or around the center) of the distribution, ourextreme-value method is natural and theoretically justified for modeling extremely largesignal amplitudes that are observed during the onset of an epileptic seizure. Therefore,the proposed method, which provides new insights into “extremal” brain connectivity, mayhave a potential impact on public health, given that epilepsy affects nearly 50 millionpeople worldwide (WHO, 2019). The anticipated impact will be deeper understanding onthe etiology of seizure and refined diagnosis brought about by the ability to differentiatebetween seizure subtypes and hence develop more targeted treatments. Figure 1 displaysEEG traces of the left and right sides of the brain during an epileptic seizure from a patientpreviously diagnosed with left temporal lobe epilepsy. The seizure onset is believed to beat channel T3 (in red).The most frequently used measure of dependence is cross-correlation. Consider a bivari-ate random vector (
X, Y ) T with joint density f ( x, y ) with support on the domain D withmarginal means E( X ) = E( Y ) = 0 and variances var( X ) = var( Y ) = 1. Then, the cross-2 igure 1: Recorded EEG signals (10–20 electrode system) from a female patient diagnosed with left temporal lobe epilepsyduring an epileptic seizure. The seizure onset happens roughly at t = 350 seconds on the reference channel T3 (red). Thedataset consists of 50 ,
000 time points at a sampling rate of 100 Hz. The central temporal (Fz), central (Cz) and centralparietal (Pz) channels are not displayed in the figure, albeit they are included in the analysis. X and Y is ρ XY = E( XY ) = (cid:82) D xyf ( x, y ) d x d y , where the integral iscalculated across the entire domain D . Indeed, classical time-domain and frequency-domainmeasures of dependence, in multivariate time series, are all derived by some averaging acrossthe entire domain D . The main limitation of these measures is that they cannot capturelocal, and sometimes conflicting, variations in the dependence structure that may happenacross the domain D . That is, the dependence structure at the “center” of the distributionmight be different from that at the “tails.” This is especially important when analyzingEEGs that are recorded during an epileptic seizure when there is an abnormal disruption inelectrical activity in the brain that results in increased amplitude and changes in spectraldecomposition at the localized focal region. In this paper, we develop a new statisticalprocedure for studying how an extreme signal amplitude caused by an epileptic seizure atthe focal region (left temporal lobe T3) may trigger a change in the signal amplitude ofanother channel of the brain network. In other words, the proposed method will be utilizedto study the conditional extremal brain connectivity of an epileptic person.Our proposed method, conditional extremal dependence for brain connectivity ( Conex-Connect , in short), jointly utilizes extreme-value theory and spectral analysis to modelthe conditional extremal dependence of brain connectivity. Conex-Connect adapts theconditional extremes model (Heffernan and Tawn, 2004) for multivariate extremes to time-varying brain signals and further extends it to the spectral setting where we study theimpact of extreme oscillations in one channel on the behavior of other channels; see Davisonand Huser (2015) for a general overview of the extreme-value theory.Consider a random vector (
X, Y ) T with marginal cumulative distribution functions(CDFs) F X and F Y , respectively. The Heffernan and Tawn (H&T) model is an asymp-totically motivated model for the conditional distribution of X given that Y is large, i.e.,4 | Y > F − Y ( p ), as p ↑
1. In finite samples, we condition on Y being larger than its p -quantile, for some (fixed) high non-exceedance probability p ∈ (0 , X and Y : asymp-totic independence, which arises when χ = lim p → Pr (cid:0) X > F − X ( p ) | Y > F − Y ( p ) (cid:1) = 0, andasymptotic dependence, which arises when χ >
0. The distinction between these tworegimes is key for extrapolating to higher levels of the conditioning variable Y . However,at the onset of a seizure, EEGs display sudden bursts or increases in amplitude, which maybe attributed to some (but not all) frequency bands. One limitation of the H&T model isthat it is unable to distinguish such features and hence needs to be adapted accordingly.Our proposed Conex-Connect method overcomes these limitations: it examines extremaldependence behavior in oscillatory components of brain signals for specific frequency bands.Moreover, the H&T model is designed for stationary data with respect to covariates.Thus, it needs to be adapted to seizure EEG signals which are highly non-stationary – bothin their marginal and dependence behavior. Here, we further modify the H&T model toflexibly capture the time dynamics of extremal brain connectivity based on an improvedversion of the penalized piecewise constant (PPC) approach of Ross et al. (2018).The remainder of the paper is organized as follows. Section 2 presents our proposedextreme-value approach to analyze EEG data. Section 3 illustrates our novel methodologyby an application to a real dataset from a patient with left temporal lobe epilepsy. Section4 concludes and discusses perspectives on future research.5 Conditional extremal dependence for connectivity(Conex-Connect)
Let ˙ Y dt be the EEG amplitude, in absolute value, of channel d ∈ { , . . . , D } at time t ∈ { , . . . , T } . Alternatively, it may also represent the amplitude of filtered signal compo-nent at a given frequency band; see Section 2.4 for more details on spectral decomposition.Figure 1 displays the behavior of these multi-channel time series. In a nutshell, our pro-posed method proceeds as follows. The first step is to select a reference channel (i.e., the conditioning variable ). In our case, it is natural to select the left temporal channel T3 (inred; Figure 1) as the reference channel because it is believed to be the focal point of seizuredischarges for this particular patient. We label it as ˙ Y = { ˙ Y t } Tt =1 . The remaining channelsare called the associated channels ( associated variables ) and are labeled as { ˙ Y d } Dd =2 , where˙ Y d = { ˙ Y dt } Tt =1 . Next, since the dependence structure in the brain network is thought toevolve across time, we segment the multi-channel EEG dataset into B (approximately)time-homogeneous blocks. Observations within the same time block are assumed to havecommon extremal characteristics. Finally, we model the changes in both the marginal anddependence properties.More precisely, for each channel d ∈ { , . . . , D } , we fit a PPC marginal model by as-suming that high threshold exceedances follow a generalized Pareto (GP) distribution witha time-varying scale parameter ν db > b ∈ { , . . . , B } . Using the empirical CDF belowthe threshold, we then transform the data from their original scale to a common marginalscale chosen to be the Laplace distribution, based on the probability integral transform.The original data are denoted by ˙ Y d , while the transformed Laplace-scale data are denoted6y Y d , d ∈ { , . . . , D } . The GP marginals are transformed to the same standard scale sothat they become comparable. Moreover, the Laplace distribution is chosen because theH&T model expressed on this scale yields a wide range of extremal dependence structures,from asymptotic independence with negative association to asymptotic dependence withpositive association (Keef et al., 2013).To estimate the conditional extremal dependence of Y d given that Y exceeds a highthreshold, we finally fit the H&T model with a time-varying dependence parameter α db ∈ [ − , α db , arediscussed in Section 2.3. The penalized term in the PPC approach stems from roughness-penalization parameters included in the likelihood functions to control the extent of vari-ation of the time-varying parameter estimates for the marginal and dependence models.We assess the uncertainty of the estimates through a bootstrap technique that accounts forthe temporal dependence present in the data. The optimal penalty parameter is obtainedusing a walk-forward cross-validation procedure designed to keep the temporal dependencein the data. The next subsections provide details on the two stages of the Conex-Connectmethod: the marginal modeling of EEG channels and their connectivity characterization. For each channel ˙ Y d = { ˙ Y dt } Tt =1 , d ∈ { , . . . , D } , and for some high non-exceedance proba-bility τ d ∈ T d ⊂ (0 , ψ db ( τ d ) as the empirical τ d -quantile for the b -th time block( b = 1 , . . . , B ). Then, motivated by extreme-value theory, we model the upper tail of ˙ Y dt F GP ( ˙ y dt ; ξ d , ν db , ψ db ( τ d )) = Pr (cid:16) ˙ Y dt ≤ ˙ y dt | ˙ Y dt > ψ db ( τ d ) (cid:17) = 1 − (cid:20) ξ d ν db { ˙ y dt − ψ db ( τ d ) } (cid:21) − ξd , for ˙ y dt ∈ (cid:0) ψ db ( τ d ) , ˙ y + db (cid:3) , where ˙ y + db = ψ db ( τ d ) − ν db /ξ d if ξ d <
0, and ˙ y + db = ∞ , other-wise. Here, ξ d ∈ R is the shape parameter, assumed constant across time blocks, while { ν db } Bb =1 ∈ (0 , ∞ ) B are the block-wise scale parameters. In the above expression, it is im-plicitly understood that time t is contained within the b -th time block, but in practice weneed to identify the correct time block for each time point, and assign parameter valuesaccordingly. For each channel d , the shape parameter ξ d controls the heaviness of the tailof the GP distribution (when compared to the tail of an exponential distribution), and wehere keep it time-constant for reasons of parsimony and because this parameter is usuallydifficult to estimate. Depending on the value of ξ d , there are three types of upper tail:bounded ( ξ d < ξ d = 0), and heavy ( ξ d > ξ d = − ξ d = 0), and the (shifted)Pareto ( ξ d > τ d , indicating the flexi-bility of the model in capturing the channel-specific extremal characteristics. In addition,despite τ d being invariant over time blocks, it does not imply threshold invariance (overtime) since the τ d -quantile is specific to each time block.Let θ d := (cid:0) { ν db } Bb =1 , ξ d (cid:1) T ∈ Θ := (0 , ∞ ) B × R be the parameter vector of the marginalGP model for channel d . Under the working assumption of independence, the likelihoodfunction is L τ d ( θ d ) = B (cid:89) b =1 (cid:89) t ∈ T b ˙ y dt >ψ db ( τ d ) ν db (cid:20) ξ d ν db { ˙ y dt − ψ db ( τ d ) } (cid:21) − ξd − , T b is the set of time points within block b , and { ˙ y dt } Tt =1 are the observed data. Tocontrol the variance of temporal fluctuations in the estimated GP scale parameters, ν db , weadd a roughness-penalization parameter λ d >
0. The negative penalized log-likelihood is (cid:96) τ d ,λ d ( θ d ) = − log L τ d ( θ d ) + λ d B B (cid:88) b =1 ν db − (cid:32) B B (cid:88) b =1 ν db (cid:33) . (1)Marginal parameter estimates are obtained by minimizing (1), while λ d is selected by cross-validation; see Section 2.5. Note that larger λ d values imply an increased smoothness ofestimates of the GP scale parameter ν db across time blocks. Also, given λ d , each channel hasa different penalized log-likelihood, i.e., each marginal fit has a different set of parameters,consisting of ( B + 1) parameters from the GP fit.Now, using the probability integral transform, we transform data to the standardLaplace scale. For observations below the threshold ψ db ( τ d ), we use the empirical CDF,denoted here by F E ( · ). First, the raw data, { ˙ y dt } Tt =1 , is transformed to the uniform scale, { u dt } Tt =1 , as follows: u dt = τ d F E ( ˙ y dt ), if ˙ y dt ≤ ψ db ( τ d ); and u dt = τ d + (1 − τ d ) F GP ( ˙ y dt ),if ˙ y dt > ψ db ( τ d ). Then, we use the inverse of the standard Laplace CDF, F L ( · ), i.e., y dt = F − L ( u dt ) = sign(0 . − u dt ) log (2 min { − u dt , u dt } ), to obtain common standardLaplace margins, { y dt } Tt =1 . After fitting the marginal models for all D channels, we obtain a standard Laplace-scalesample { y t , y t , . . . , y Dt } Tt =1 , with y t representing the transformed time series of the ref-erence channel; here, T3. Here, the lowercase { y t , . . . , y Dt } Tt =1 denote (transformed) real-ized values, while the uppercase { Y t , . . . , Y Dt } Tt =1 denote the corresponding random vari-ables. The next goal is to study the conditional dependence of the associated variables,9 Y t , . . . , Y Dt } , given that the conditioning variable, { Y t } , takes on a large value (in theupper tail of the distribution). Define ˜ τ ∈ ˜ T ⊂ (0 ,
1) to be a non-exceedance probabil-ity such that φ (˜ τ ) is the ˜ τ -quantile of the standard Laplace distribution for the referencechannel Y . Let ˜ θ = (cid:16) { α db } D,Bd =2 ,b =1 , { β d } Dd =2 , { µ d } Dd =2 , { σ d } Dd =2 (cid:17) T ∈ ˜ Θ := [ − , ( D − B × ( −∞ , ( D − × R ( D − × (0 , ∞ ) ( D − be the parameter vector of the H&T model for allchannels. Thus, according to the H&T model, for all t ∈ T b (within time block b ) such that y t > φ (˜ τ ), conditional on Y t = y t , Y dt may be expressed as Y dt = α db y t + y β d t W dt , d = 2 , . . . , D, b = 1 , . . . , B, (2)where, for model estimation purposes, the components of the random variable W dt areassumed to be mutually independent and normally distributed with mean µ d and standarddeviation σ d , both time-constant.The parameters { α db } D,Bd =2 ,b =1 are the “first-order” dependence parameters, and theycapture the extent of block-wise linear extremal dependence between channels Y d and large Y , which is allowed to change over time (blocks). When α db ∈ (0 , Y d and large Y within time block b . This association becomesstronger as α db increases, with (positive) asymptotic dependence corresponding to α db = 1and β d = 0. For α db ∈ [ − , α db → −
1. As α db →
0, the linear dependence weakens. In particular, when Y d and large Y are independent, then α db = β d = 0. On the other hand, the time-constant parameters { β d } Dd =2 may be considered as capturing the “second-order” dependence characteristicssince they specify the spread of the data around the linear relationship given by { α db } forincreasing values of Y . As β d → −∞ , the distribution of the data (around the linearrelationship dictated by the α ’s) becomes tighter for higher values of the conditioningvariable. If β d >
0, we have the opposite behavior; the distribution of the data has a wider10pread around the linear relationship (given by the α ’s) for higher values of the conditioningvariable. Notice that among all dependence parameters, α db is the only one that is allowedto vary across time-blocks b = 1 , . . . , B . The other parameters ( β d , µ d , σ d ) are intentionallykept constant over time in order to reduce the overall uncertainty, and avoid interferenceswith the estimation of α db , which drives the main dependence feature.From (2), we may rewrite the model as Y dt | Y t = y t ∼ N (cid:16) α db y t + µ d y β d t , ( σ d y β d t ) (cid:17) , (3)for all t ∈ T b (within time block b ) such that y t > φ (˜ τ ), d = 2, . . . , D , and b = 1, . . . , B .The parameters for each channel are estimated by minimizing the negative log-likelihoodbased on (3), i.e.,˜ (cid:96) ˜ τ,d ( ˜ θ ) = B (cid:88) b =1 (cid:88) t ∈ T b y t >φ (˜ τ )
12 log (2 π ) + log ( σ d y β d t ) + 12 (cid:32) y dt − ( α db y t + µ d y β d t ) σ d y β d t (cid:33) . (4)Here, there are ( B + 3) parameters per channel d . Keef et al. (2013) proposed additionalconstraints on ˜ Θ , leading to a smaller set of feasible parameters, which reduce the varianceof the estimators and overcome complications on the modeling of negatively associatedvariables and parameter identifiability. Also, these additional constraints avoid drawingconditional inferences inconsistent with the marginal distributions. Here, we use a prag-matic approach and restrict { β d } Dd =2 to the interval [0 , . Y dt giventhat Y t equals its own 0 . . α db across time blocks, which stabilizestheir fluctuations and reduces uncertainty. The joint negative penalized log-likelihood is˜ (cid:96) ˜ τ, ˜ λ ( ˜ θ ) = D (cid:88) d =2 ˜ (cid:96) ˜ τ,d ( ˜ θ ) + ˜ λ D (cid:88) d =2 B B (cid:88) b =1 α db − (cid:32) B B (cid:88) b =1 α db (cid:33) , where ˜ λ > λ d , performing D − .4 Extremal Spectral Analysis One major novelty of our proposed Conex-Connect approach is that we apply it not onlyto the original time series (i.e., from a time-domain perspective) but also to the time series“dissected” into several frequency bands (i.e., from a frequency-domain perspective). Thisreveals hidden features of the extremal dependence structure between channels that onlyaffect certain frequency bands but may not be visible by looking at the entire time series.Most of the current work on spectral analysis of electrophysiological data focus on thespectral estimation and its association with behavioral measures and brain states. None ofthese examine the downstream effect of extremal dependence in a reference channel on theentire network (Krafty et al., 2011; Krafty and Collinge, 2013; Fiecas and Ombao, 2016;Krafty et al., 2017; Ombao et al., 2018; Scheffler et al., 2020). Similarly, none of the work onconditional extreme-value theory (e.g., Heffernan and Tawn (2004), or Keef et al. (2013))examine extremal dependence from a spectral perspective. Here, fill this gap and study theconditional extremal dependence in a brain network from a frequency-domain perspective.EEGs are zero-mean signals that can be expressed as a mixture of frequencies oscillatingat the following standard frequency bands: Delta-band (Ω : 1–4 Hz), Theta-band (Ω : 4–8Hz), Alpha-band (Ω : 8–12 Hz), Beta-band (Ω : 12–30 Hz), and Gamma-band (Ω : 30–50Hz); assuming a sampling rate of 100 Hz (Ombao et al., 2016; Nunez and Srinivasan, 2007).The spectral decomposition of channel T3 is displayed in Figure 2 for both the pre- andpost-seizure onset phases.Following Ombao and Van Bellegem (2008) and Gao et al. (2020), two EEG signals, Y ( t ) and Y ( t ), may be decomposed into the 5 rhythms as Y ( t ) ≈ (cid:80) k =1 Y , Ω k ( t ) and Y ( t ) ≈ (cid:80) k =1 Y , Ω k ( t ), where Y d, Ω k ( t ) denotes the Ω k -waveform in channel Y d ( t ). In practicethese rhythms are obtained by applying a linear filter such as the Butterworth filter (see13 igure 2: Spectral decomposition of channel T3 signal during 1 second of both pre- and post-seizure onset phases from theEEG recording of a patient diagnosed with left temporal lobe epilepsy. The frequency bands are: Delta-band (Ω : 1–4 Hz),Theta-band (Ω : 4–8 Hz), Alpha-band (Ω : 8–12 Hz), Beta-band (Ω : 12–30 Hz), and Gamma-band (Ω : 30–50 Hz). Ω k ( t ) = max (cid:96) | cor ( Y , Ω k ( t ) , Y , Ω k ( t + (cid:96) )) | . Another major contribution of this paper is a novel measure of dependence based onthe extremal behavior of the oscillations of brain channels. We propose a measure thatexamines the impact of unusually large frequency band amplitudes in the T3 channel(the reference channel projecting from the seizure foci on the temporal lobe) on the otherchannels. In order to study the effect of the extremal behavior in the Ω k -waveform, westudy the conditional distribution of | Y d, Ω i ( t ) | given | Y , Ω j ( t ) | > τ , i, j = 1 , . . . ,
5. Thenew measure produces new interesting results that give us deeper insights into the highlynon-linear interactions and dependence between channels during an epileptic seizure event.
This section gives details on both the block bootstrap and walk-forward cross-validationprocedures designed to assess the uncertainty of estimated parameters and to select theoptimal roughness penalty parameters respectively, while handling the auto-correlation ofthe time series.In both the marginal and the H&T models, we need to select the optimal value ofthe roughness penalty parameter, λ >
0, in an objective manner. In both cases, for eachchannel and within each time block, we rely on a 10-fold walk-forward cross-validationprocedure (Hyndman and Athanasopoulos, 2019) since the data are auto-correlated acrosstime. For this procedure, we divide each time block b into k = 10 adjacent non-overlappingfolds. Then, we consider the first r folds as the training set and the last k − r folds as the15est set. For a grid of λ values and for each split configuration, we then fit the model tothe training set. We use the parameter estimates to obtain the value of the log-likelihoodfunction on the test set, and we finally add up these log-likelihood values of each splitconfiguration to get an overall score. The value of λ that maximizes the overall score istaken as the optimal value for the penalization parameter. Figure 3 shows a schematicillustration of this cross-validation procedure. Figure 3:
Schematic for a 10-fold walk-forward cross-validation procedure for the auto-correlated data within time block 1of the pre-seizure phase of the EEG application. The blue dots represent the training sets. The red dots are the test sets.
To assess the uncertainty of parameter estimates, we employ a block bootstrap pro-cedure which preserves local temporal dependence in the data. As described in Lahiri(2003), we resample entire blocks of consecutive observations (rather than single observa-tions) keeping the correspondence between the conditioning variable Y and the associatedvariables Y , . . . , Y D , and then refit the marginal and H&T dependence models. From the16ootstrap samples, we compute confidence intervals for the parameters that evolve withtime and histograms for those that are constant over time. Notice that the bootstrap blocksdiffer from the B time blocks used in the PPC model formulation.Specifically, we split each of these B time blocks into C non-overlapping (sub-)blocks,of size s , which are resampled with replacement to generate M bootstrap samples. Theblock size s is determined with the aid of the sample autocorrelation function to make sureit is large enough to keep the overall temporal dependence structure.Simultaneously, for each bootstrap sample, the non-exceedance probability for bothmarginal and H&T models is also randomly sampled from a uniform distribution givena range of reasonable values to account for thresholding uncertainty: τ d ∈ T d ⊂ (0 , d ∈ { , . . . , D } , and ˜ τ ∈ ˜ T = ⊂ (0 , We examine the changes in connectivity (or dependence between channels) using the pro-posed Conex-Connect method. In this paper, the scalp EEG recording is from a femalepatient diagnosed with left temporal lobe epilepsy. The EEG data is recorded from a 10–20 system, and D = 19 channels are available for analysis. Figure 1 shows the electrodeplacement and EEG traces for some channels. In this specific case, the physician knew inadvance that the seizure onset would be on the left temporal lobe region, which would be17aptured by the recording at the T3 channel ( d = 1). Thus, T3 is set to be the referencechannel (i.e., the conditioning variable in the H&T model). The remaining 18 channelsare treated as the associated variables ( d = 2 , . . . , ,
000 time points perchannel.The first stage of the Conex-Connect method is to model the marginal extremes foreach channel d independently. Each channel d is split into pre- and post-seizure onsetphases at time t = 350 seconds and thus T = 15 ,
000 points for each of the post-seizureonset and pre-seizure phase. Each phase is segmented into B = 12 time blocks of 1 , M = 500 block bootstrap samples aregenerated within each time block b = 1 , . . . , B = 12, with bootstrap block size s = 25data points. For all channels d = 1 , . . . , D , a common non-exceedance probability interval, T d = (0 . , . T = (0 . , . Figure 4 displays a dashboard with the results of the Conex-Connect method for channel F7given high values of the reference channel T3, both for pre- and post-seizure onset phases.Panel I) provides the individual behavior of the H&T parameter estimates: { ˆ α ,b } Bb =1 , ˆ β ,ˆ µ , and ˆ σ . Here, the subscript d = 2 refers to channel F7 and b ∈ { , . . . , } denotes thetime block index. The 95% confidence bands for ˆ α ,b and the histograms for ˆ β , ˆ µ , andˆ σ are drawn out from the block bootstrap procedure; see Section 2.5. In the pre-seizurephase, the estimates of the first-order dependence parameter, ˆ α ,b , b = 1 , . . . ,
12, are veryhigh and stable, close to 1, indicating a strong positive linear relationship between F7 andlarge values of T3. This makes sense as both channels F7 and T3 are located on the leftside of the brain. Note that the confidence bands are narrow, suggesting a high level ofcertainty concerning this dependence structure. The histogram of ˆ β for the pre-seizurephase is centered around 0 . α , being much smaller than ˆ α , , indicating weaker extremal dependence. Also, the behaviorof ˆ α ,b , b = 13 , . . . ,
24, are much more erratic with wider confidence bands, suggesting ahigher level of uncertainty after the seizure onset. This is interesting because it reveals themore chaotic nature of the seizure process. For some blocks, the estimates are close to 0 . β , there is a shift to the left in the histogramwith a high concentration around 0, pointing to a possible change towards independence.By contrasting the findings above to panel I) of Figure 5, we can see that a differentstory prevails on the right side of the brain. For the right frontal channel F8 ( d = 3),given high values of the reference channel T3, it appears that seizure does not affect theextremal dependence structure as much, since there is practically no discernible change inthe estimates ˆ α ,b before and after the seizure. Also, in the pre-seizure phase, the linearrelationship between F8 and large values of T3 is weaker when compared to F7. This ispartly because F7 is closer to the seizure focus location, which is around the referencechannel T3. In addition, the histograms for ˆ β , ˆ µ , and ˆ σ , both for pre- and post-seizureonset phases, are relatively stable, indicating the absence of changes in the second-orderextremal dependence characteristics.Back to Figure 4, panel II) displays the time evolution of the conditional distribution ofthe conditional 0.975-quantile in (3) given that T3 reaches its own 0.975-quantile, estimatedsemi-parametrically as explained in Section 2.3. In this plot, we can analyze the parametersof the H&T model jointly. Also, with the different violin plots, we investigate how the effectof large values of T3 at time t impacts channel F7 at time t + (cid:96) , for time lags (cid:96) = 0, 0 . .
25, and 0 .
50 seconds. The Conex-Connect method shows that the extremal dependenceis stronger at lag (cid:96) = 0 and weakens for other lags, indicating a stronger contemporaneousextremal dependence than lagged extremal dependence. Moreover, when contrasting lag (cid:96) = 0 to the other lags, the discrepancy between the violins, both in terms of mediansand shapes, is more pronounced after the seizure onset. This suggests that seizure has animpact on the conditional extremal dependence of the brain network. The same panel of20 igure 4:
A dashboard with results of the Conex-Connect method for pre- and post-seizure onset phases. Channel F7(blue) given high values of T3 (red) are highlighted in the EEG scalp cartoon.
Panel I)
In the first line, the evolution ofthe estimated first-order dependence parameter α db (solid blue line) through time with its bootstrap mean (solid black line)and 95% confidence bands (dashed lines). In the second line, the histograms of the bootstrap estimates for scale exponentparameter β d , and for residual mean µ d and scale σ d . Panel II)
Bootstrap violin plots for the 0 . . Panel III)
Effect of the different frequency bands (Ω : 1–4, Ω : 4–8, Ω :8–12, Ω : 12–30, and Ω : 30–50 Hz) on the first-order dependence parameter. Darker blue pixels indicate higher dependence. igure 5: A dashboard with results of the Conex-Connect method for pre- and post-seizure onset phases. Channel F8(blue) given high values of T3 (red) are highlighted in the EEG scalp cartoon.
Panel I)
In the first line, the evolution ofthe estimated first-order dependence parameter α db (solid blue line) through time with its bootstrap mean (solid black line)and 95% confidence bands (dashed lines). In the second line, the histograms of the bootstrap estimates for scale exponentparameter β d , and for residual mean µ d and scale σ d . Panel II)
Bootstrap violin plots for the 0 . . Panel III)
Effect of the different frequency bands (Ω : 1–4, Ω : 4–8, Ω :8–12, Ω : 12–30, and Ω : 30–50 Hz) on the first-order dependence parameter. Darker blue pixels indicate higher dependence. (cid:96) = 0 is indistinguishable from those of higherlags, and at a lower level overall than for channel F7. This denotes less synchrony (inthe extremal dependence) between the channel corresponding to the seizure foci (referencechannel) and the channel on the contra-lateral side of the foci.The Conex-Connect method produces interesting results regarding how the extreme val-ues of the different oscillations of the reference channel T3 impact brain connectivity. Figure4, panel III), shows the estimated first-order dependence coefficient ˆ α ,b ( b = 1 , . . . , B ) be-tween the (absolute) Ω k -waveform in channel F7 given high (absolute) amplitudes of thesame waveform in channel T3. Values of ˆ α ,b closer to 1 are darker. In the pre-seizure phase,the extreme values in the high-frequency Gamma-band exhibit the lowest level of extremaldependence. This seems to be consistent across the entire pre-seizure phase. However, thedependence pattern changes in the post-seizure phase. First, the extremes of Gamma-bandfrom T3 shows the highest level of extremal dependence with the Gamma-band from F7.This is quite interesting because sudden outbursts of high-frequency oscillations typicallycharacterize seizure onset as shown by Medvedev et al. (2011). Moreover, since the post-seizure onset is typically non-stationary, we see that this dependence structure also evolvesover time blocks. In the right side of the brain, Figure 5, panel III), shows that beforethe seizure occurs, the higher values of the mid-frequency Beta-band lead the changes inthe extremal dependence structure. In the post-seizure phase, we notice that, immediatelyafter the seizure onset, the high-frequency Gamma-band becomes more prominent, similarto the left side of the brain (channel F7).Figure 6, presents a dashboard with the results of classical methods based on cross-correlation and cross-coherence for comparison. The left column shows the results for F7,23hile the right column shows the results for F8. In all panels, we display both pre- andpost-seizure onset phases. Panel I) displays the evolution of the cross-correlation overtime blocks. Panel II) shows how the cross-correlation between T3 at time t evolves whencomputed for future values of channel F7 and F8 at time t + (cid:96) , for time lags (cid:96) = 0, 0 . .
25, and 0 .
50 seconds. Finally, panel III) exhibits the impact of high values of the differentfrequency bands of T3 in its cross-coherence with F7 and F8.
Figure 6:
A dashboard with results of the classical analysis for pre- and post-seizure onset phases. Left column displaysresults for channel F7 (brain left size) while the right column show the results for F8 (brain right side).
Panel I)
Evolutionof the cross-correlation estimates (black solid line), along with its confidence bands (dashed lines), between channels F7 andT3 (reference channel) and between channels F8 and T3.
Panel II)
Violin plots to access the changes in the distributionof the cross-correlation both over time and for different lag values of the associated channels.
Panel III)
Effect of thedifferent frequency bands on the cross-coherence. Darker blue pixels indicates higher values of cross-coherence. For all panels,uncertainty is obtained trough 500 bootstrap samples. (cid:96) = 0) againstother higher lags. These similarities suggest either (a.) that the conditional extremaldependence dominates the global dependence (as measured by correlation) or (b.) that thephenomenon that we see in the joint tail is similar to other less extreme quantiles of thedistribution. Regarding the spectral decomposition, our method shows that for both sidesof the brain, immediately after the seizure, the high-frequency Gamma-band becomes themost relevant frequency in explaining the conditional extremal dependence. This findingdoes not agree with the results in terms of the classical coherence; see Figure 6, panelIII). We further investigate the conditional extremal dependence in terms of frequencyoscillations by decomposing the associated channels in their canonical frequency bands.We refit the model for all pairs of Ω k -waveforms. Figure 7 displays the impact of extremevalues of the different frequency bands of the reference channel T3 on the different frequencybands of the associated channels F7 and F8. Note that the heatmap in panel III) of Figures4 and 5 corresponds to the diagonals of Figure 7. Here, beyond the previous findings interms of seizure phases and sides of the brain, we notice that the medium-frequency bands,Beta and Alpha, of the reference channel T3, impact the first-order dependence parameterestimates, mainly on the Gamma-band of the associated channels. This may indicate25hat the Gamma-band can be used for feature engineering to improve the performance ofmachine learning algorithms for epilepsy detection. In this paper, we present the novel Conex-Connect method, the first extreme-value model-based approach to learn patterns in the extremal dependence during periods of high volatil-ity in brain signals. The method extends the Heffernan and Tawn (2004) model to capturetime-varying extremal dependence features, both from time-domain and frequency-domainperspectives, while adequately assessing estimation uncertainty using a block bootstrapprocedure. We here give a full characterization of the conditional extremal dependence ofbrain connectivity, shedding light on how the brain network responds to an epileptic seizureevent. We also study the extremal brain connectivity based on the spectral decompositionof the different channels. To the best of our knowledge, our method is a pioneer in linkingthe association between extreme values of frequency oscillations in a reference channel withoscillations in other channels of the brain network.Essentially, our methodology relies on “dissecting” the original brain signal into variouscomponents at different frequency bands, and to study the extremal dependence structure ofeach frequency component separately based on distinct conditional extreme-value models.While our statistical analysis is motivated by a very concrete applied problem and ourempirical results reveal salient features of brain connectivity during an epileptic seizure,it would also be interesting in future research to further investigate the theoretical linkbetween the (conditional) extremal dependence structure of each frequency component tothat of the overall signal. Specifically, we expect that the frequency component with the26 igure 7:
Effect of the different frequency bands of the reference channel T3 on the different frequency bands of the associatedchannels, F7 and F8, pre- and post-seizure onset phases. The solid blue line represents the evolution of the estimated first-order dependence parameter α through time with its bootstrap 95% confidence bands (shaded area). The vertical red lineis the seizure onset. The frequency bands are: Delta-band (Ω : 1–4 Hz), Theta-band (Ω : 4–8 Hz), Alpha-band (Ω : 8–12Hz), Beta-band (Ω : 12–30 Hz), and Gamma-band (Ω : 30–50 Hz). References
Cohen, M. X. (2014).
Analyzing neural time series data . USA: The MIT Press.28rowder, M. (2001). On repeated measures analysis with misspecified covariance structure.
J. R. Statist Soc. B 63 , 55–62.Davison, A. and R. Huser (2015). Statistics of extremes.
Annual Review of Statistics andIts Application 2 , 203–235.Fiecas, M. and H. Ombao (2016). Modeling the evolution of dynamic brain processes duringan associative learning experiment.
J. Am. Stat. Assoc. 111 , 1440–1453.Gao, X., W. Shen, B. Shahbaba, N. Fortin, and H. Ombao (2020). Evolutionary state-spacemodels with applications to time-frequency analysis of local field potentials.
StatisticaSinica 30 , 1561–1582.Heffernan, J. E. and J. A. Tawn (2004). A conditional approach for multivariate extremevalues.
J. R. Statist Soc. B 66 , 497–546.Hyndman, R. J. and G. Athanasopoulos (2019).
Forecasting: principles and practice (3rded.). Australia: OTexts. https://OTexts.com/fpp3 . Accessed on 08/20/2020.Keef, C., J. Papastathopoulos, and J. A. Tawn (2013). Estimation of the conditionaldistribution of a multivariate variable given that one of its components is large: additionalconstraints for the heffernan and tawn model.
J. Multivar. Anal. 115 , 396–404.Krafty, R. T. and W. O. Collinge (2013). Penalized multivariate whittle likelihood forpower spectrum estimation.
Biometrika 100 , 447–458.Krafty, R. T., M. H. Hall, and W. Guo (2011). Functional mixed effects spectral analysis.
Biometrika 98 , 583–598.Krafty, R. T., O. Rosen, D. S. Stoffer, D. J. Buysse, and M. H. Hall (2017). Conditionalspectral analysis of replicated multiple time series with application to nocturnal physi-ology.
J. Am. Stat. Assoc. 112 , 1405–1411.Lahiri, S. N. (2003).
Resampling methods for dependent data . USA: Springer-Verlag.29edvedev, A. V., A. M. Murro, and K. J. Meador (2011). Abnormal interictal gammaactivity may manifest a seizure onset zone in temporal lobe epilepsy.
InternationalJournal of Neural Systems 21 , 103–114.Nunez, P. L. and R. Srinivasan (2007). Electroencephalogram.
Scholarpedia 2 , 1348.Ombao, H., M. Fiecas, C. M. Ting, and Y. F. Low (2018). Statistical models for brainsignals with properties that evolve across trials.
NeuroImage 180 , 609–618.Ombao, H., M. Lindquist, T. W., and J. Aston (2016).
Handbook of neuroimaging dataanalysis . USA: Chapman and Hall/CRC.Ombao, H. and S. Van Bellegem (2008). Coherence analysis: a linear filtering point ofview.
IEEE Transactions on Signal Processing 56 (6), 2259–2266.Ross, E., S. Sam, D. Randell, G. Feld, and P. Jonathan (2018). Estimating surge in extremeNorth Sea storms.
Ocean Engineering 154 , 430–444.Scheffler, A. W., A. Dickinson, C. DiStefano, S. Jeste, and D. S¸ent¨urk (2020). Covariate-adjusted hybrid principal components analysis. In M. J. Lesot, S. Vieira, M. Z. Refor-mat, J. P. Carvalho, A. Wilbik, B. Bouchon-Meunier, and R. R. Yager (Eds.),
Informa-tion processing and management of uncertainty in knowledge-based systems . Switzerland:Springer.Wadsworth, J. L. and J. Tawn (2019, December). Higher-dimensional spatial extremes viasingle-site conditioning. arXiv e-prints , arXiv:1912.06560.WHO (2019). Epilepsy: a public health imperative - summary.