Analysis of EEG data using complex geometric structurization
AAnalysis of EEG data using complex geometricstructurization
E. A. Kwessi ∗ L. J. Edwards † Abstract
Electroencephalogram (EEG) is a common tool used to understand brain ac-tivities. The data are typically obtained by placing electrodes at the surface ofthe scalp and recording the oscillations of currents passing through the electrodes.These oscillations can sometimes lead to various interpretations, depending on thesubject’s health condition, the experiment carried out, the sensitivity of the toolsused, human manipulations etc. The data obtained over time can be considereda time series. There is evidence in the literature that epilepsy EEG data may bechaotic. Either way, the embedding theory in dynamical systems suggests thattime series from a complex system could be used to reconstruct its phase spaceunder proper conditions. In this paper, we propose an analysis of epilepsy elec-troencephalogram time series data based on a novel approach dubbed complexgeometric structurization. Complex geometric structurization stems from the con-struction of strange attractors using embedding theory from dynamical systems.The complex geometric structures are themselves obtained using a geometry tool,namely the α -shapes from shape analysis. Initial analyses show a proof of conceptin that these complex structures capture the expected changes brain in lobes underconsideration. Further, a deeper analysis suggests that these complex structurescan be used as biomarkers for seizure changes. Keywords:
EEG, time series, complex structure, morphometry, alpha-shape
In neuroscience, to understand brain functions and brain related diseases, it is verycommon to place electrodes on a subject’s head and record the electrical activities thatresult. These activities, when represented, are often oscillations that change in shape,frequency, and range. They can be analyzed to see if meaningful information can beextracted that gives a clue about the brain state of the subject at a certain time, orin a certain region of the brain. The accuracy of the recording highly depends on theinstrument used, the region of interest (ROI), the experimenter’s experience, the time of ∗ Corresponding author, Trinity University, Department of Mathematics, Visiting the University ofAlabama at Birmingham, Department of Biostatistics, [email protected], [email protected] † Department of Biostatistics, University of Alabama at Birmingham, [email protected] a r X i v : . [ ee ss . SP ] F e b he day, the subject’s discipline during the recording and many other factors. This is tosay that uncertainty in the accuracy of the information recorded is always present. Evenwhen the information is being recorded on the same subject, on the same day, on thesame ROI, but at different intervals, changes are bound to occur. Electroencephalogramor EEG, a term coined by Berger (1929), are electrical activities recorded on humans oranimals that display prominent oscillatory behavior subject to important changes duringvarious behavioral states. These changes show a high degree of nonlinearity in the signalsthat may be important. Indeed, in the field of biomedical signal processing (analysis ofheart rate variability, electrocardiogram, hand tremor, EEG), the presence or absenceof nonlinearity often conveys information about the health condition of a subject. Inparticular, EEG signals are often examined using nonlinearity analysis techniques orby comparing signals that are recorded during different physiological brain states (e.g.epileptic seizure). The differences observed during these analyses can either be due togenuine differences in dynamical properties of the brain or due to differences in recordingparameters. The EEG are often analyzed as times series and there are many methods foranalysis of times series in the literature. The methods can be grouped into two categories:univariate measures and multivariate measures. A good review on the topic can be foundin Carney et al. (2011).Among the methods that have been touted as more efficient at providing an insightinto the real dynamics of EEG is the famous Embedding Theorem of Takens (1981). ThisEmbedding Theorem has been instrumental to understanding how to reconstruct the truedynamics of systems based on the times series measured on these systems. Essentially,this reconstruction theory, in layman’s terms, suggests that a times series measured oversufficiently long period of times contains enough information to reconstruct the phasespace in which the associated system normally evolves. This allows also to show that thereare other intricate subunits that influence the changes observed in the measured variablesthat are represented in the time series. This theorem was used for example by Grassbergerand Procaccia (1983) to propose a measure called Correlation Dimension which was inturn instrumental to Lehnertz and Elger (1998) for the prediction of epilepsy seizures.Despite all the methods proposed in the literature, there is no agreement on which methodconstitutes the best tool at extracting the most meaningful information that could beuseful for the physician for the prediction of seizure-like diseases. Moreover, with theincreasing use of the concepts of chaos and complexity in health sciences, it is becomingmore and more difficult to distinguish their adequate application. For example, therehave been evidence of chaos in EEG data, see Destexhe (1992); Destexhe and Babloyantz(1986); Destexhe et al. (1998), and since chaotic systems are inherently complicated, theymay look complex. Likewise, complex systems may also look chaotic. Distinguishingthese two notions is important in applications, especially in health sciences. Henceforth,we will adopt a terminology along the lines of Rickles et al. (2007) for understandingcomplexity and chaos.In this paper, we propose a new method for analyzing EEG times series data, whichwe call Complex Geometric Structurization (CGS). The complex nature of the methodstems from the fact that we have multiple subunits interacting together resulting in arich collective behavior feeding back into the behavior of individual parts. The methodis inspired by the Embedding Theorem of Takens (1981) for the construction of a ge-ometric structure whose volume can be evaluated from shape analysis technique. The2olume of this geometric structure behaves as a key statistic akin to a biomarker for thephenomenon or ROI of interest. Using data driven approaches to study brain pathologiesis a very active field of study nowadays due to improvement in life expectancy across theworld with its cohorts of problems such as brain disorders as illustrated in Zheng et al.(2019), David et al. (2020) and the references therein. Moreover, the push to use dataand methodological driven approaches to brain pathologies is also evidenced by the nu-merous grants offered by the National Institute of Health and private foundations suchas the Michael J. Fox Foundation for Parkinson Disease, the Bill and Melinda GatesFoundations, just to name some.The remainder of the paper is organized as follows: In Section 2, we make a briefreview on how to use the Embedding Theorem to construct strange attractors; in Section3, we review important notions of statistical morphometry; in Section 4, we introduce thecomplex geometric structurization method; in Section 5, we discuss some applications ofthe CGS method on real data; in Section 6, we discuss the pros and cons of the proposedmethod in different contexts. As its name suggests, a dynamical system is a system whose variables evolve over time.Its phase space is a geometric representation of the trajectory of its variables over time.The values taken by the system’s variables at an instant describe the system’s states.To understand how to reconstruct the phase space of a dynamical system based onobservations (times series) of one of its variables, we need to revisit the EmbeddingTheory of Takens (1981), which is essentially a high dimension transformation of thetime series. Consider the n -dimensional space R n . We recall that a manifold M in thespace R n is a topological space that locally looks like a Euclidian space near each point.Topology here means that bending is ignored. For example, the surface of the globe is atopological manifold in the space R . Now, consider a dynamical system with a system’sstate x ( t ) lying on a manifold M of R n . Let ρ be a sampling interval and let the timeseries s ( t ) = g ( x ( t )) be given as a one-dimensional observation of the system dynamicsthrough an observation function g . Takens (1981) embedding theory states that foralmost every smooth function g , the delay coordinate map defined as F : R n → R m withF( x ( t )) = [ s ( t ) , s ( t − ρ ) , · · · , s ( t − ( m − ρ )] T is an embedding, that is, it is a one-to-oneimmersion of the state space attractor with dimension d when m > d . In other words,the result states that F( x ( t )) is a representative of x ( t ), even if the true state space M has not been observed, see Figure 1. The quantity m is referred to in the literature asthe embedding dimension and ρ as the time delay (or lag). The embedding theory ispredicated on the observation that a time series observed over a long period of time mayshow an internal structure. In fact, considering a time series s ( t ) = g ( x ( t )), we onlyobserve an incomplete picture of x ( t ) since s ( t ) is a scalar. However, if we observe it fora long period of time, a more precise description will emerge, which will help understandits dynamics. In practice, most of the focus is given into how to find appropriate valuesfor the time delay ρ and the embedding dimension m , see Appendix 8.2. We observethat this type of reconstruction technique has been applied successfully in the case ofepilepsy EEG data before, see for instance Destexhe and Babloyantz (1986), Destexhe(1992), and Destexhe et al. (1998). 3igure 1: An illustration of the embedding mechanism. The Takens’ embedding theory can be used to reconstruct theLorenz, a famous attractor often mentioned in dynamical systems. The Lorenz (1963)systems of differential equations is given as ˙ x = s ( y − x )˙ y = rx − y − xz ˙ z = xy − bz . It is known for example that the Lorenz system is chaotic for s = 10 , r = 28, and b = 8 / x = x ( t ) , y = x ( t − ρ ), and z = x ( t − ρ ). The time step used is ∆ t = 0 .
005 foran interval time of [0, 75] for the Lorenz system. We also note that the time lag or delayis ρ = 31∆ t , and was found using either the autocorrelation (ACF) or average mutualinformation (AMI), see for instance Section 8.2 below.4igure 2: Lorenz attractor (red) with its reconstructed counterparts (blue) and plottedin the same coordinates system. We can observe the topological equivalence between theoriginal phase space and its reconstructed counterpart. Real Data:
There have been many applications of Takens’ embedding theory with realdata. For example, in Fisher et al. (2009),Carney et al. (2011), an attractor is constructedfrom a publicly available (at ) epilepsy data set called here EDATA for simplicity. The dataconsist of five sets A, B, C, D, and E. Each containing 100 single-channel EEG segmentsof 23.6 seconds, each of which was selected after visual inspection for artifacts (such asacoustic and electrical shielding, separation of earth ground for laboratory, interconnec-tivity of devices on the same phase and ground centrally and locally) and has passed aweak stationarity criterion. Sets A and B were obtained from surface EEG recordingsof five healthy subjects with eyes open and closed, respectively. Data were obtained inseizure-free intervals from five patients in the epileptogenic zone for set D and from thehippocampal formation of the opposite hemisphere of the brain for set C. Set E containsseizure activity, selected from all recording sites exhibiting ictal activity. Sets A and Bhave been recorded extracranially, whereas sets C, D, and E have been recorded intracra-nially. All EEG signals were recorded with the same 128-channel amplifier system, usingan average common reference [omitting electrodes containing pathological activity (C,D,and E) or strong eye movement artifacts (A and B)]. After 12 bit analog-to-digital con-version, the data were written continuously onto the disk of a data acquisition computersystem at a sampling rate of 173.61 Hz. Band-pass filter settings were 0.53–40 Hz (12dB/oct.), see Andrzejak et al. (2001). Figure 3 below is an illustration of the datasetEDATA. Each row represents one time series from sets A, B, C, D, and E respectively.Clearly, time series in the seizure set E has a much pronounced amplitude synonymouswith more brain activities. In Figure 4 below, we represent the reconstructed phasespaces based on the time selected times series from set A – E . We note that the axesare x = x ( t ) , y = x ( t − ρ ) , z = x ( t − ρ ) where ρ = 1∆ t , with ∆ t = fs = 5 .
76 ms.5 − A − B − − C − D − E Figure 3: This figure represents one time series selected at random from each set A–E.We observe that the amplitude is much more pronounced in the set E, which representsthe seizure prone patients. 6
A) (B)(C) (D)(E)
Figure 4: Reconstructed phase space diagram, restricted to the space xyz .7 Statistical Morphometry
Given a set of points in a two or three-dimensional space, statistical morphology (ormorphometrics) amounts to finding an appropriate geometric characterization of thevariability of the shape and size of the set of points. This characterization includes butis not limited to volumes and surface area, surface roughness, deviation from convex-ity, porosity, and permeability. We note that in general, the set of points may not beconvex in the strictest sense, so there is a need for a method that relaxes the convexityrestriction. The notions of α -convexity and α -shape represent alternatives that relax thestrict convexity assumption. These new α -shaped objects are even more flexible than α -convex objects in that α now controls the spatial scale of the estimator. note that α isa unitless quantity. In fact, as α decreases, the α -shape shrinks and more space appearsamong the sample points, whereas as α increases, the α -shape object converges to theconvex hull of the sample. In other words, in α -shaped objects, α controls the amount ofporosity between the sample points. These α -shaped objects have been used in variousfields for the characterization of biological systems, see for instance Lafarge et al. (2014)and the references therein, or also Gardiner et al. (2018). We can mention the pivotalwork of Edelsbrunner and M¨ucke (1994) where the main algorithm for the constructionof α -shaped objects can be found. Example
In the example below, we illustrate the α -shape construction respectively in two andthree dimension, based on a random sample of data taken from the original object S .We consider 2500 points in 3D, obtained from the object S which is the objectdelimited by the curve and z = x + y where x and y are random points selected in theinterval [ − , α -shape object (red) for α = 0 (a), α = 0 . α = 0 . α = 2 (d). (a) (b) (c) (d) Figure 5: α -shape object (red) of S ( z = x + y ) for α = 0 (a), α = 0 . α = 0 . α = 2 (d). In view of the apparent shape that can be observed from the reconstructed phase spaceabove, the question of how to compare these complex structures arises. In other words,this question is related to the question of how to compare strange attractors. Among the8ethods proposed, we can mention the work of Grassberger and Procaccia (1983) andits many variants. The method we propose is borne out of the observation that thereseems to be a complex geometric structure whose shape changes from healthy patientsto seizure-prone ones. So we will rely on the notion of α -shape to construct the complexstructure related to each situation. We will use the following algorithm to construct our“complex structure” (CGS). The motivation comes from the fact that, given a time series,if it has been observed long enough, it carries the signature of the original phase spacediagram in which the true system it originated from evolves. If we have many such timeseries from the same system, we should be able to capture enough information about thetrue phase space diagram. The large number of these time series should be enough toeventually eliminate noise or undesired artefacts from the reconstruction. We then expecteach reconstructed phase diagram to be topologically equivalent to any other, thereforeforming a structure compact enough to be the revolving center of all other reconstructedphase space diagram. In general, the the dimension of the space in which the true systemevolves is greater than three, making it impossible to visualize with the naked-eye. Ourempirical approach is to consider only a cross section of the reconstructed phase space indimension three in this case. the reason for this choice is two-fold: one, we can actuallyvisualize a part of the true phase space diagram. Second, we can use existing method toevaluate the volume of the structure in lower dimensions. This is to say that the volumeis preserved in reconstruction. In doing this, we are trying to find a measurable identifieror markup for this group of time series that will vary from groups to groups and fromindividual to individual. Given N time series, we use the Takens reconstruction technique to obtain the embeddingdimension m n (using ACF) and time delay (or sampling interval) ρ n for n = 1 , · · · , N (using the method of false nearest neighbors). Let m = min { m , · · · , m N } and ρ =min { ρ n , n = 1 , , · · · , N } . If m ≥
3, then for each time series, we obtain the complexstructure CGS α ( n ) . Let α = min { α ( n ) : Vol D (CGS α ( n ) ) is maximized } . We use the α -shape technique to obtain the volume of the CGS α in 3D. This step is crucial sincewe choose to represent the complex structure only using the first three delay-coordinates x = x ( t ) , x = x ( t − ρ ) , x = x ( t − ρ ) even if the actual space has dimension m >
3. Themotivation for this selection is that the volume of the CGS based on any combinationof three coordinates would not significantly be different from any other volume obtainedfrom any other combination of three coordinates. The reasoning behind the choice of α is that the volume of the complex structure is bounded by the volume of its convex hull,as α increases, so we select the smallest alpha that maximizes the volume, see section4.4.3 below for an illustration. This leads to the following algorithm
1. Use all the times series collected on patients in different groups.2. For each time series, reconstruct a 3D “strange attractor” from the first threedelay-coordinates. 9. Use the α -shape technique to construct a “Complex Geometric Structure” (CGS)related to all strange attractors.4. Define a measure related to each CGS that that can be statistically analyzed. Thiscould be the volume, the surface area, the center, etc.5. Repeat this procedure for all replicates (if there are any) and thus obtain new data.6. (a) For comparison:Use a statistical test to see if there is a significant difference between measuresof different groups. This could be a parametric or a nonparametric test.(b) For Prediction:Use the data obtained for training in machine learning, see Kwessi and Ed-wards (2020) (preferred to standard generalized linear model (GLM) modelsbecause of the absence of a specific model) and test it on potential new data.
1. The CGS terminology stems from the fact the structures obtained after reconstruc-tion do not have of classical geometric shapes. Their shapes are rather “complex”in nature.2. We propose to use all time series available for a particular region of the brain, ata specific instant; for example before seizure, during seizure, or after seizure. It iscommon to use one time series, see for instance Fisher et al. (2009). We note thatrepeated measures on the same subject do not yield the same values and thus inthe reconstruction, there could be individual times series whose excursions in thephase space are wider than others yielding a bigger geometric structure in size andvolume. In this case, individual times series could be used and the subsequent largevolumes could be trimmed out if necessary. This is particularly important if one isinterested for example in obtaining a richer data set to analyze.3. The embedding dimension m often obtained is greater than 3, so in the absenceof a visualization mechanism for data in spaces of dimensions higher than 3, wepropose to focus on only the first three delay coordinates.4. Using the R-package Alphashape3d , measures such as volumes or surface areas canbe obtained for the α -shape object under consideration.5. As we saw in the example above, the choice of α is critical in the α -shape con-struction. We propose to select the minimum value of α for which the volume ismaximized.6. Finally, we note that in practice, the experimenter can used data from chaoticor complex systems. One can check for chaos in data by calculating Lyapunovexponents. If the data are sensitive to initial conditions and chaotic, the attractorwould be referred to as a “strange attractor”, see Celso et al. (1987). Even if on theother hand the data are sensitive to initial conditions but non chaotic, we would stillkeep the same “strange attractor” denomination because in this case (the Lyapunov10xponents are non-positive), the attractor obtained is still strange, see for instanceCelso et al. (1984); Paladin and Vulpiani (1987). It also noteworthy to observe thatnon chaotic systems that are sensitive to initial conditions are sometimes referredto in the literature as complex systems, see for instance the comparative reviewbetween chaotic and complex systems by Rickles et al. (2007). In this section, we examine some technical considerations to keep in mind when imple-menting this method. α As we mentioned above, the choice of α is critical in this process. In what follows, weshow that the optimal value of α is smallest that maximizes the volume of the CGS.Moreover, if one is interested in comparing CGS among groups, it would be adequate toselect a common value of α , for example as the largest α value, among the values thatmaximize the CGS for each group. For instance, using the EDATA set, we obtained m = 10 and ρ = 1∆ t , with ∆ t = fs = 5 .
76 ms. In Figure 6 below, we observe that fora value of α ≈ α ≈ α optim = 580 as the optimal valueof α if we want to compare the two CGSs. a V o l u m e 0123456 0 . . . . . . . . CS for ACS for E
Figure 6: Evolution of the volume CGS as function of α . We note that the embedding dimension is closely related to the time lag. We are notgoing to discuss which should be estimated first. However, we note that there are twomethods for estimating the time lag and they do not always yield the same embeddingdimension. We propose to select the embedding dimension as the smallest value of thetwo embedding dimensions when they are different.11 .4.3 Choice of the delay-coordinates
In Section 4.4.2, we mentioned that we will select the first three delay-coordinates toconstruct the complex structure. However, it is worthwhile to consider the question ofwhether the volume would change if a different combination of delay-coordinates is used.Our choice is based on the conjecture that it does not matter which combination of delay-coordinates is selected. To emphasize that point, we will discuss the case of subsets Aand E. In the case of subsets A and E, the embedding dimension found is m = 10. Sothere are (cid:0) (cid:1) = 120 possible different combinations of three delay-coordinates, selectedamong 10. So, for each combination, we will construct the corresponding CGS andassess whether the volume changes significantly across all the different CGS’s. Figure 7shows the volume of the CGS by combination of three delay-coordinates for subsets Aand E. The boxplots suggest that the distribution of volumes for each set is reasonablyconcentrated around its median (thick red and blue lines) with a small range and nooutliers. This is to say that taking the first three delay coordinates seems reasonabledespite minor variations otherwise. (A) (E) . . . . . . V o l u m e o f t he C S . . . . . V o l u m e o f t he C S Figure 7: Evolution of the volume the CGS for the 120 different combinations of 3 delay-coordinates for subsets A and E.The main takeaway from the plots above is that the variation of the volume as a functionof the combination of delay-coordinates appear not to be significant. Even for seizuredata like subset E, the variation in volume appears to be mild, which is our impetus forconjecturing that a similar observation could be made for other datasets. Obviously wecannot predict what will happen for all datasets, and ultimately that endeavor wouldrequire a mathematical proof that may go beyond the scope of the present manuscript.
In Figure 4, we have shown a representation of the reconstructed phase space diagramfor one time series in the space x = x ( t ) , y = x ( t − ρ ) , z = x ( t − ρ ) and for each set12–E. In Figures 8–12, we construct the phase space diagram and the complex structuresfor the EDATA set, for all the 100 time series in each subset. We observe that for eachset, we obtain a compact structure whose volume we can now evaluate, see ( A1) – (E1) .We selected α = 580 for each case, because of all sets A-E, 580 is the value of alpha thatmaximizes the volume of the complex structure of E, the largest of all of them (see theselection criterion given above), see figure 6 above. The CGS structure for each set isthen obtained, see (A2) – (E2) . We observe that the CGS for set E appears to be biggerthan all other CGSs, which is a sign of larger excursions in the phase space and thereforesynonymous with a much intense brain activity during seizure. (A1) (A2) Figure 8: Complex geometric structure for set A in EDATA.13
B1) (B2)
Figure 9: Complex geometric structure for set B in EDATA. (C1) (C2)
Figure 10: Complex geometric structure for sets C in EDATA.14
D1) (D2)
Figure 11: Complex geometric structure for set D in EDATA. (E1) (E2)
Figure 12: Complex geometric structure for set E in EDATA.Pooling the data may seem to inflate the volume of the CGS. So instead, we couldconstruct the complex structure related to individual channel and evaluate their volumes.This will give a richer data to analyze, where outliers can be removed. Figure 13 belowshows the boxplot of the volume of the complex structure for each channel. The volumefor subset E is very large compared to A-E and obscures the distribution for subsets A-D.15s a result, in the 2nd panel in Figure 13 we remove subset E so that the distributionsof subsets A-D can be more clearly seen.
A B C D E . . . . . V o l u m e o f t he C S ABCDE A B C D . . . . . . Figure 13: Boxplots for the CGS for subsets A–E in EDATA, obtained from individualtimes series.This Figure is further evidence of what already observable in the raw data in Figure 3and Figure 4. Sets A–D have comparable amplitude at the time series level confirmed bythe fact that the volume of CGS are in similar range. The amplitude at the time serieslevel of Set E is much higher than that of sets A–D, which is also confirmed by a largervolume at CGS level. This is further evidence that during seizure, these volumes increasesubstantially when compared to seizure-free intervals. More importantly, we can extractmeaningful statistics from this data as suggested in Table 1 below:Min Q Median Q Max Mean SDA 12.30 63.63 115.87 154.17 317.88 109.84 66.08B 37.89 63.63 173.50 291.31 479.06 364.25 266.63C 7.73 68.38 161.30 262.70 351.60 262.10 333.70D 9.45 62.78 174.10 4569.00 17800.00 692.30 206.54E 516.3 3277.30 9279.40 36270.50 134300.00 26422.50 33682.48Table 1: Table of volumes of CGS for each set A–E. The values are of order 10 − . In this example, we discuss data collected for the Brain Core Initiative at UAB. Thisdata was collected on 20 individuals in two regions of their brain, namely the Auditory16nd Visual cortex, see Figure 14 below. The patients were subject to three “tasks”:Auditory, Visual, and Rest. Measurements of brain activity was obtained as EEG timesseries, after removing unnecessary artefacts. For this data set, the optimum value of α is α optim = 300. We will consider the following subsets: Auditory cortex-Auditory task,Auditory cortex-Visual task, Auditory cortex-Rest, Visual cortex-Visual task, Visualcortex-Visual task, Visual cortex-Rest.Figure 14: Functional diagram of the Brain lobes. Image credit Chen (2011). Macro-level analysis
Here all the time series are used. Figure 15 below represents the CGS for each of thesubsets above. 17 uditory–Auditory Task Auditory–Visual task Auditory–RestVol=7.8181 Vol=7.8911 Vol=11.449Visual-Auditory Task Visual–Visual task Visual–RestVol=8.2875 Vol=7.3373 Vol=14.745
Figure 15: Complex geometric structures for the above subsets.Comparing the first and second row, the volumes seem to differ by cortex. Comparing thefirst, second, and third column, the difference in volumes of tasks are respectively 0.449,0.5538, and 3.296. These numbers show similarity between auditory and visual tasksbut they both differ from rest. The increase of volume during rest may not be counter-intuitive. In fact, there is vast literature linking resting brain activities to underlyinghigh cognitive processes such as moral reasoning, self-consciousness, remembering pastexperiences, or planning for the future, see for instance Bruckner et al. (2008).
Micro-level analysis
However, if we use individual time series, since they represent individual patients, thenwe can obtain their CGS by cortex and by task. This would enable having a richerdataset and a more in-depth comparison. In Figure 16 and 17 below, we plot the densityand boxplot of the volumes of the CGS by cortex and/or by task. Plots (A) representthe plots of CGS by cortex, plots (B) the plots of cortexes by task, plots (C) the plotsby task, and plots (D) the plots tasks by cortex. At first glance, it seems as if a Beta orLognormal distribution would provide a good fit to these densities.18 . . . . . Volume D en s i t y Auditory cortexVisual cortex(A) 0 2 4 6 8 10 . . . . . . Volume D en s i t y Auditory cortex−Auditory taskAuditory cortex−Visual taskAuditory cortex−RestVisual cortex−visual TaskVisual cortex−Auditory taskVisual cortex−Rest(B)0 2 4 6 8 10 . . . . . Volume D en s i t y Auditory taskVisual taskRest(C) 0 2 4 6 8 10 . . . . . Volume D en s i t y Auditory Task−Auditory cortexAuditory Task−Visual cortexVisual task−Auditory cortexVisual Task−Visual cortexRest−Auditory cortexRest−Visual cortex(D)
Figure 16: Density plots for the volume of CGS obtained from individual times series.19
Auditory cortexVisual cortex (A)
13 1112 919 56 19
Auditory cortex−Auditory taskAuditory cortex−Visual taskAuditory cortex−RestVisual cortex−Auditory taskVisual cortex−Visual taskVisual cortex−Rest (B)
Auditory taskVisual taskRest (C)
Auditory task−Auditory cortexAuditory task−Visual cortexVisual task−Auditory cortexVisual task−Visual cortexRest−Auditory cortexRest−Visual cortex
Auditory task−Auditory cortexAuditory task−Visual cortexVisual task−Auditory cortexVisual task−Visual cortexRest−Auditory cortexRest−Visual cortex (D)
Figure 17: Boxplots for the volume of CGS obtained from individual times series.We note there are 20 patients and three tasks per cortex for a total of 120 observations.When grouping by cortex, we labelled the observations 1 through 60 per cortex. Whengrouping by task, we labelled the observations 1 through 40, 1 through 20 when groupingby task within cortex. The numbers shown in the boxplots in Figure 17 represent thelabels for the corresponding outlying observations. Boxplots (A) does not suggest adifference between the auditory and visual cortexes, since the notches do overlap for themost part. It shows that there are two common outliers (37 and 57) in both cortexes.This means that there are two patients whose brain activities are more pronounced thanothers in both cortexes, causing wide excursions in the phase space and therefore largevolumes of their CGS. Boxplots (B) and (C) suggest a similarity between the tasks sinceall the notches do overlap. Boxplots (D) suggest that the difference observed in Boxplot(A) is due to the difference between the auditory tasks in both cortexes. Boxplots (B)basically confirm the observations in Boxplots (A) that all the tasks are similar. A more20n-depth analysis is needed to make more meaningful conclusions.In the tables below, we report the intrinsic discrepancy, that is, the minimum be-tween d KL ( F, G ) and d KL ( G, F ), where d KL ( F, G ) is defined as the Kullback-Leibler (KL)directed divergences between the probability distributions F and G (see Appendix Sec-tion 8.4). So the smaller the intrinsic discrepancy, the bigger the difference between F and G . Auditory cortex Visual cortex P-valueAuditory cortex 0 0.160 0.021Table 2: Statistical analysis of the intrinsic discrepancy between the cortexes.Auditory task Visual task Rest P-valueAuditory task 0.000 0.104 0.029– (0.98) (0.98) 0.451Visual task 0.104 0.000 0.080(0.98) – (0.98)Table 3: Statistical analysis of the intrinsic discrepancy between the tasks.To make a meaningful use of Table 2 and 3, a reference point is needed. However, thedensities above allows us to make the hypothesis that the volumes are not normallydistributed. A nonparametric test such as a Wilcoxon-Mann-Whitney would be neces-sary to compare for example the volumes of the two cortexes, controlling for the tasksperformed. Such a test yields a p-value of 0.021, which is small enough to suggest a statis-tically significant difference between the CGS volumes for Auditory and Visual cortex. Italso suggests that the intrinsic discrepancy 0.160 obtained above is a sign of a statisticalsignificant difference between the two cortexes. Likewise, the nonparametric Kruskall-Wallis test is used to compare the tasks and it yields a p-value of 0.451, which is large,suggesting that the volumes of the CGS are not significantly different. The p-values forthe pairwise Wilcoxon rank sum test (in parenthesis) with Bonferroni correction are alsolarge. This could be further interpreted as the brain activities of individuals during theauditory task are not significantly different from their brain activities during the visualtask. Now, whether these differences or lack thereof are corroborated at the biologicallevel remains to be proved. Comparison with a Cross-Correlation Function
In this section, we will compare our method with the cross-correlation function (CCF),see Appendix, Section 7.5. In the heatmaps below, represented are distances ρ ( X, Y )where X and Y are the EEG data (times series) collected on the 20 different patients indifferent regions of the brain, when subjected to the three tasks above. The heat gradient21epresents the values of ρ ( X, Y ). Figure 18 (A) below represents the comparison betweenauditory and visual cortexes. Its symmetric nature is indicative of the similarity betweenthe two cortexes. This is in agreement with Figure 17 (A). From Figure 18 (B) below,we observe that the heatmaps are very similar. In particular, there are four areas thatare redder than the majority and indicative of the strongest correlation. A classificationby a k-means algorithm confirms the existence of these two clusters, see Figure 19 (B)in the Appendix. This is similar to saying that there is no significant difference betweenthe tasks. The same conclusion can be made with boxplots in Figure 17 (C), where theoutliers in the volume of CGS match precisely the regions with high correlation. Lookingat the first three heatmaps in 18 (C) and (D) , there is not a huge difference betweenthe heatmaps, confirming that there is no significant difference between the tasks withinthe auditory cortex, an observation also made from Figure 17 (B). The same observationcan be made about the last three maps in the visual cortex (A) (B)
Auditory cortex vs Visual cortex
Distance
A task vs V task
A task vs Rest
V task vs Rest
Distance (C) (D)
AA vs AR
AA vs AV
AV vs AR
Distance
VV vs VR
VV vs VA
VA vs VR
Distance
Figure 18: Plot of ρ ( X, Y ) between the auditory and visual cortexes.
We offer some items for discussion as a result of our findings:1. The method we are proposing has the potential to discriminate brain activities indifferent parts of the brain. Further, changes in volume could be an indication of changes22f activity level in the part of the brain under study. This predictive potential could beused to predict or monitor potential brain disorder.2. This method has the potential to be extended to experiments in which EEG dataare suitable such as epilepsy, Alzheimers, attention deficit disorder, learning disabilities,anxiety disorders, fetal alcohol syndrome, autism, chronic pain, insomnia, dyslexia, etc,see for instance Satheesh Kumar and Bhuvaneswari (2012) and the references therein.3. Preliminary analysis of Local Field Potential (LFP) data suggests that this methodcan be extended beyond EEG. In fact LFP differs from EEG in that the electrodes areinserted in the brain tissue rather than at the surface of the scalp. However, we did notobtain the authorization to release the results on the LFP data at this time.4. There have been studies suggesting that group averaging of neuroimaging data arenot clinically relevant. One possible explanation for this issue is the lack of focus onindividuals. The method we propose can remedy that by providing a measure that isindividual-focused.5. We have not addressed the weak stationarity of the EEG data used here and howmuch stationarity or lack thereof plays a role in the method we are proposing. The dataset EDATA was checked for weak stationarity but the data collected for the Brain CoreInitiative was too small for stationarity. This overall is an issue in multiple studies wheresamples tend to be small.6. One of the drawbacks of the “distance” above is that its definition may change andtherefore this could affect the interpretation of results. In fact, given two times series X and Y , the distance could also be defined as ρ ( X, Y ) = 1 − max ≤ τ ≤ m l {| CCF ( X, Y, τ ) |} orby ρ ( X, Y ) = mean { CCF ( X, Y, τ ) } with CCF ( X, Y, τ ) = E [( X t + τ − X ) · ( Y t − Y )] √ E [( X t − X ) ] E [( Y t − Y ) ] , where X t + τ is the time-shifted version of X t , and τ is the time lag separating the two timesseries X and Y , X, Y are the respective means of the times series X and Y . The readercan refer to Arbabshirani et al. (2012) for more information. In this paper, we have proposed a method called complex geometric structurization tohelp analyze EEG signals. The method is based on the embedding theory of dynamicalsystems and shape analysis. The method works well on epilepsy data. The method hasthe advantage to discriminate individual EEG and also discriminates groups of EEG sig-nals. More importantly, the method performs better when compared to Cross-CorrelationFunction. It is important to note that the results and the method, though empirical, offeran important proof of concept relating time series and the volume of their reconstructedcomplex in three dimension that need to be explored further by validating them withmore solid mathematical concepts. If successful, this method could be an importantaddition to the literature of EEG signal analysis and can be used to explore other brainpathologies such as sleep disorder, schizophrenia, or Alzheimers. Its discrimination po-tential also can be used to improve our understanding of functional brain connectivity ingeneral. 23
Appendix8.1 α -convex and α -shapeDefinition 1. A set A ⊆ R m is said to be α -convex, for α > if A = C α ( A ) = ∪ B B α ( x ) , (8.1) where B = { B α ( x ) : A ∩ B α ( x ) (cid:54) = 0 } , B α ( x ) = { y : (cid:107) y − x (cid:107) ≥ α } , and (cid:107)·(cid:107) is a norm in R m . The quantity C α ( A ) in equation (8.1) is called the α -convex hull of A . Now suppose wehave a sample S n = { X , · · · , X n } obtained from an object S in R m , then an estimatorof S is C α ( S n ) if S is assumed α -convex. Estimators of α -convex objects are constructedfrom arc of circles in two-dimension (2D) and spherical caps in three-dimension (3D). Definition 2.
Let A ⊆ R m be α -convex. An α -shaped estimator of A is an estimator ofits α -convex hull C α ( A ) obtained by approximating arc of circles with polygonal curves in2D and spherical caps with polyhedral surfaces in 3D. There are two popular methods for estimating the timedelay ρ : the autocorrelation function (ACF) and the average mutual information (AMI).Indeed, consider N measurements of a time series s ( t ). Then the sample ACF is definedas ρ ( t ) = n (cid:88) n =1 ( s ( n + t ) − s ) ( s ( n ) − s ) n (cid:88) n =1 ( s ( n ) − s ) , with s = N − N (cid:88) n =1 s ( n ) . (8.2)The time delay is chosen as the ρ = min t> { ρ ( t ) < } .Now define the AMI as I ( t ) = 1 N N (cid:88) n =1 P r ( s ( n ) , s ( n + t )) log (cid:18) P r ( s ( n ) , s ( n + t )) P r ( s ( n )) P r ( s ( n + t )) (cid:19) , (8.3)where P r ( s ( n )) and P r ( s ( n ) , s ( n + t )) are respectively the probability of observing s ( n )and the probability of observing s ( n ) and s ( n + t ). The time delay is estimated as thefirst local minima of I ( t ). The most popular method for estimating the embedding dimension m is the so-calledFalse Nearest Neighbors technique, see for example Kennel et al. (1992).24 .4 Kullback-Leibler divergence function Given two distributions F and G of a continuous random variable, with respective densityfunctions f and g , the Kullback-Leibler divergence function is defined as d KL ( F, G ) = (cid:90) ∞−∞ f ( x ) log (cid:18) f ( x ) g ( x ) (cid:19) dx . (8.4)This quantity measures the degree to which F diverges from G . We will use it to assessthe difference between the densities of the different cortex per task and thus the differencebetween their brain activities, see Section 5.2 The cross-correlation function (CCF) is defined as : given a time lag τ , and two time-series X = { X t } and Y = { Y t } , the CCF is defined as CCF ( X, Y, τ ) = E [( X t − τ − X )( Y t − Y )] (cid:113) E [( X t − τ − X ) ] E [( Y t − τ − Y ) ] , where X and Y are the respective mean of the time series X and Y .The CCF is then calculated over range of temporal lags and a “distance” ρ ( · , · ) isdefined as the maximum absolute CCF over the interval [1 , m l ] ( for a maximum lag valueof m l to be selected) ρ ( X, Y ) = max ≤ τ ≤ m l {| CCF ( X, Y, τ ) |} . Note that difference defines the similarity between the two time series.25
A) (B)
Auditory cortex vs Visual cortex c l u s t e r c l u s t e r Distance
A task vs V task c l u s t e r c l u s t e r A task vs Rest
V task vs Rest
Distance (C) (D)
AA vs AR c l u s t e r c l u s t e r AA vs AV
AV vs AR
Distance
VV vs VR c l u s t e r c l u s t e r VV vs VA
VA vs VR
Distance
Figure 19: Plot of ρ ( X, Y ) between the Auditory and Visual cortexes by cluster. Theheat gradient represents the values of ρ ( X, Y ). Bibliography
R. G. Andrzejak, K. Lehnertz, F. Mormann, C. Rieke, P. David, and C. E Elger. Indica-tions of nonlinear deterministic and finite-dimensional structures in time seriesof brainelectrical activity: Dependence on recording region and brain state.
Physical ReviewE , 2001. doi: 10.1103/PhysRevE.64.061907.M. R. Arbabshirani, M. Havlicek, K. A. Kiehl, G. D. Pearlson, and V. D. Calhoun.Functional network connectivity during rest and task conditions: a comparative study.
Human brain mapping , 34(11):2959–2971, 2012.H. Berger. ¨Uber das elektroenkephalogramm des menschen.
Arch Psychiatr Nervenkr ,87(1):527–570, 1929.R. L. Bruckner, J. R. Andrews-Hanna, and D. L. Schacter. The brain’s default network:Anatomy, function, and relevance to disease.
Ann. N.Y. Acad. Sci , 1124:1–38, 2008.P. R. Carney, S. Myers, and J. D. Geyer. Seizure prediction: methods.
Epilepsy &behavior , 22:S94–S101, 2011. 26. Celso, E. Ott., S. Pelikan, and J. A. Yorke. Strange attractors that are not chaotic.
Physica D: Nonlinear Phenomena. Elsevier BV , 13(1–2):261–268, 1984.G. Celso, E. Ott., and J. A. Yorke. Chaos, strange attractors, and fractal basin boundariesin nonlinear dynamics.
Science , 238(4827):632–638), 1987.P. Chen. Principles of biological science. (Last accessed February 19, 2021), 2011. URL http://bio1152.nicerweb.com/ .S. A. David, J.A.T Machado, C.M.C. I´nacio, and C.A. Valentin. A combined measureto differentiate eeg signals using fractal dimension and mfdfa-hurst.
Communicationsin Nonlinear Science and Numerical Simulation , 84:1, 2020.A Destexhe.
Nonlinear Dynamics of the RhythmicalActivity of the Brain: . PhD thesis,Universit´e Libre de Bruxelles, Brussels, Belgium, 1992.A. Destexhe and A. Babloyantz. Low-dimensional chaos in an instance of epilepsy.
Proc.Natl. Acad. Sci , 83:3513–3517, 1986.A. Destexhe, J. A. Sepulchre, and A. Babloyantz. A comparative study of the experi-mental quantification of deterministic chaos.
Physics Letters , A 132:101–106, 1998.H. Edelsbrunner and E.P. M¨ucke. Three-dimensional alpha shapes.
ACM Transactionson Graphics , 13(1):43–72, 1994.N. K. Fisher, S. S. Talathi, A. Cadotte, and P. R. Carney. Epilepsy detecting andmonitoring. In S. Tong and N. V. Thakor, editors,
Quantitative EEG analysis methodsand clinical applications , chapter 6, pages 141–165. Artech House, 2009.J. D. Gardiner, J. Behnsen, and C. A. Brassey. Alpha shapes: determining 3d shapecomplexity across morphologically diverse structures.
BMC Evolutionary Biology , 18(184), 2018. doi: 10.1186/s12862-018-1305-z.P. Grassberger and I. Procaccia. Measuring the strangeness of strange attractors.
PhysicaD. Nonlinear Phenomena , 9(1-2):189–208, 1983.M. Kennel, R. Brown, and H. Abarbanel. Determining embedding dimension for phase-space reconstruction using a geometrical construction.
Physical Review A , 45(6):3403–3411, 1992.E. A. Kwessi and L. J. Edwards. Artificial neural networks with a signed-rank objectivefunction and applications.
Communication in Statistics-Simulations and Computa-tions , 2020. doi: 10.1080/03610918.2020.1714659.T. Lafarge, B. Pateiro-Lopez, A. Possolo, and J. P. Dunkers. R implementation of apolyhedral approximation to a 3d set of points using the α -shape. J. Stat. Software ,54(4):1–19, 2014.K. Lehnertz and C. E. Elger. Can epileptic seizures be predicted? evidence from nonlineartime series analysis of brain electrical activity.
Physical Review Letters , 80:5019–5022,1998. 27. N. Lorenz. Deterministic nonperiodic flow.
Journal of the Atmospheric Sciences , 20(2):130–141, 1963.G. Paladin and A. Vulpiani. Anomalous scaling laws in multifractals objects.
PhysicsReports , 156(4):147–225, 1987.D. Rickles, P. Hawe, and A. Shiell. A simple guide to chaos and complexity.
Journal ofepidemiology and community health , 61(11):933–937, 2007.J. Satheesh Kumar and P. Bhuvaneswari. Analysis of electroencephalography (eeg) sig-nals and its categorization-a study.
Procedia Engineering , 38:2525–2536, 2012.F. Takens. Detecting strange attractors in turbulence dynamical systems and turbulence(lecture notes in mathematics), vol. 898, 1981.J. Zheng, H. Fushing, and L. Ge. A data-driven approach to predict and classify epilepticseizures from brain-wide calcium imaging video data.