Lokalisatie van epileptische aanvallen met tensor decompositie technieken

Lokalisatie van epileptische aanvallen met tensor decompositie technieken Daan Camps

Thesis voorgedragen tot het behalen van de graad van Master of Science in de ingenieurswetenschappen: wiskundige ingenieurstechnieken Promotoren: Prof. dr. ir. Lieven De Lathauwer Prof. dr. ir. Sabine Van Huffel Assessoren: Prof. dr. ir. Marc Van Barel Prof. dr. ir. Johan Suykens Begeleiders: ir. Borbála Hunyadi ir. Laurent Sorber

Academiejaar 2012 – 2013

© Copyright KU Leuven Zonder voorafgaande schriftelijke toestemming van zowel de promotoren als de auteur is overnemen, kopiëren, gebruiken of realiseren van deze uitgave of gedeelten ervan verboden. Voor aanvragen tot of informatie i.v.m. het overnemen en/of gebruik en/of realisatie van gedeelten uit deze publicatie, wend u tot het Departement Computerwetenschappen, Celestijnenlaan 200A bus 2402, B-3001 Heverlee, +32-16327700 of via e-mail [email protected]. Voorafgaande schriftelijke toestemming van de promotoren is eveneens vereist voor het aanwenden van de in deze masterproef beschreven (originele) methoden, producten, schakelingen en programma’s voor industrieel of commercieel nut en voor de inzending van deze publicatie ter deelname aan wetenschappelijke prijzen of wedstrijden.

Voorwoord Het begin van deze tekst is een geschikt moment om enkele personen te bedanken die me met raad en daad bijgestaan hebben bij het onderzoek voor deze thesis. Aanvangen doe ik vanzelfsprekend met mijn beide promotoren, prof. Lieven De Lathauwer en prof. Sabine Van Huffel. Hun uitgebreide expertise in respectievelijk multilineaire algebra en biomedische dataverwerking heeft regelmatig tot gegronde opmerkingen en feedback op dit werk geleid en me nieuwe, verbeterde inzichten gegeven. De in dit werk gebruikte methode is bovendien door prof. Lieven De Lathauwer geïntroduceerd, zonder hem zou dit onderzoek niet mogelijk geweest zijn. Ten tweede wil ik ook mijn begeleiders, Borbála Hunyadi en Laurent Sorber bedanken om me verder te helpen met mijn vragen. Zonder het werk van Laurent Sorber, die de gebruikte Tensorlab toolbox ontwikkeld heeft, zou deze thesis niet tot stand gekomen zijn. De assessoren, prof. Marc Van Barel en prof. Johan Suykens, wil ik alvast bedanken voor het kritisch nalezen van de tekst. Afsluitend wil ik ook mijn familie en vriendin bedanken. Niet alleen hebben ze me quasi onafgebroken gesteund, ook heb ik hen verveeld met het nalezen van de tekst en hebben ze me op, onderhand talloze, typo’s en taalfouten gewezen. Ware het niet voor hun, dan zou u zich wellicht aan heel wat meer taalfouten geërgerd hebben. Daan Camps

i

Inhoudsopgave Voorwoord

i

Samenvatting

iv

Lĳst van figuren en tabellen

v

Lijst van afkortingen en symbolen 1 Literatuurstudie 1.1 Een inleiding tot tensor decompositie technieken 1.1.1 Matrixalgebra . . . . . . . . . . . . . . . . 1.1.2 Tensoralgebra . . . . . . . . . . . . . . . . 1.1.2.1 MLSVD . . . . . . . . . . . . . . 1.1.2.2 CPD . . . . . . . . . . . . . . . . 1.1.2.3 BTD . . . . . . . . . . . . . . . . 1.2 Epilepsie: eigenschappen, monitoring en analyse

viii . . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

2 Simulatiestudie 2.1 Het genereren van artificiële epileptische insulten . . . . . . . . . . 2.2 Constructie van de tensor met een continue wavelet transformatie 2.3 De ruiscomponent in simulaties . . . . . . . . . . . . . . . . . . . . . 2.4 Dipoolsimulaties met statisch gesitueerde bronnen . . . . . . . . . . 2.4.1 Stationaire frequenties en CPD . . . . . . . . . . . . . . . . . 2.4.1.1 Invloed van het type ruis . . . . . . . . . . . . . . . 2.4.1.2 Vergelijking met SOBI, een andere BSS methode 2.4.1.3 Lokalisatie van de bron . . . . . . . . . . . . . . . . 2.4.1.4 Meerdere statisch gesitueerde bronnen met stationaire frequenties . . . . . . . . . . . . . . . . . 2.4.2 Veranderlijke frequenties en BTD-(L,L,1) . . . . . . . . . . . 2.4.2.1 Het bepalen van de rang (Lr ) in de BTD-(L,L,1) ontbinding . . . . . . . . . . . . . . . . . . . . . . . . 2.4.2.2 Lokalisatie van de bron (reprise) . . . . . . . . . . 2.4.2.3 Meerdere statisch gesitueerde bronnen met veranderlijke frequenties . . . . . . . . . . . . . . . 2.5 Dipoolsimulaties met dynamisch gesitueerde bronnen . . . . . . . . 2.6 Tensoren geconstrueerd uit Hankel matrices en BTD-(L, L, 1) . . . 3 Validatie op klinische data ii

. . . . . . . . . . . . . . .

. . . . . . .

1 1 4 5 6 7 9 10

. . . . . . . .

15 15 18 18 21 21 25 28 29

. . 32 . . 33 . . 38 . . 41 . . 42 . . 43 . . 47 53

Inhoudsopgave 3.1 3.2 3.3

MA 1 - Patiënt 22 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 MA 2 - Patiënt 11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 MA 3 - Patiënt 26 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

4 Conclusie

63

A Gebruikte Software

69

B Artikel: Decomposition of a tensor constructed from ictal scalp EEG in rank-(Lr , Lr , 1) terms determines the evolution in location and frequency of the seizure

71

C Poster

79

Bibliografie

81

iii

Samenvatting Tensor decompositie technieken, met name een ontbinding in rang-1 termen of CPD, zijn reeds eerder succesvol toegepast om epileptische activiteit uit ruizige EEG data te scheiden en een lokalisatie van het insult uit te voeren. De tensor wordt hier typisch geconstrueerd door de CWT te berekenen van ieder kanaal. Tekortkomingen van deze methode komen naar voren als de epileptische activiteit van frequentie of lokalisatie verandert. Indien de frequentie veranderlijk is, zijn de temporele en frequentie component van de epileptische tensor verschillend van rang-1. Indien de lokalisatie verandert, zijn de ruimtelijke en temporele component verschillend van rang-1. In dit onderzoek is een ontbinding in termen van lage multilineaire rang of BTD-(L, L, 1) gebruikt om epileptische activiteit uit een ruizige EEG tensor, geconstrueerd met CWT, te scheiden. Afzonderlijk zijn een verandering in frequentie en een veranderlijke lokalisatie succesvol gemodelleerd met een BTD(L, L, 1) ontbinding en kunnen beide tot op lage SNR bepaald worden. De methode is uitvoerig getest in een uitgebreide simulatiestudie en vergeleken met SOBI als referentie BSS techniek. Nadien zijn er drie voorbeelden uit een klinische dataset besproken die een verdere validatie van de methode geven. Voor deze voorbeelden is de lokalisatie bepaald BTD-(L, L, 1) steeds in overeenstemming met de beschikbare SISCOM data. De resultaten van de BTD-(L, L, 1) ontbinding duiden in één van de drie voorbeelden op een frequentieverandering gedurende het insult. In een ander voorbeeld is er een kleine verandering in de lokalisatie van het insult opgemerkt met BTD-(L, L, 1). Naast de tensorisatie met CWT is er een tweede methode onderzocht die de tensor construeert door uit elk kanaal een Hankel matrix te vormen. Dit geeft aanleiding tot matrices van lage rang en bijgevolg een (L, L, 1) structuur in de tensor met zuivere epileptische activiteit. Deze methode laat in combinatie met BTD-(L, L, 1) opnieuw een betrouwbare scheiding toe van de epileptische activiteit uit de ruizige tensor tot op lage SNR. Bij de klinische validatie zijn de resultaten van deze methode in goede overeenstemming met de andere methode en de SISCOM data.

iv

Lĳst van figuren en tabellen Lĳst van figuren 1.1 1.2 1.3 1.4 1.5 1.6

Kolommen, rijen en tubes van een derde orde tensor . . . . . . . . mode-1, mode-2 en mode-3 uitvouwing van een derde orde tensor CPD van een tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . BTD-(Lr , Lr , 1) van een tensor . . . . . . . . . . . . . . . . . . . . Plaatsing van de 21 elektroden op het hoofd . . . . . . . . . . . . Tensoriseren van EEG data naar een derde orde tensor via CWT

2.1

Dipoollocatie in het hoofd, samen met elektrodeposities en -benamingen, en de lokatie van de neus en oren . . . . . . . . . . . . . . . . . . . . . . . . Twee voorbeelden van een dipoolsimulatie met respectievelijk stationaire en veranderlijke frequentie . . . . . . . . . . . . . . . . . . . . . . . . . . . . Artefacten geëxtraheerd met CCA . . . . . . . . . . . . . . . . . . . . . . . MLSVD van de gentensoriseerde artefacten . . . . . . . . . . . . . . . . . . Evidentie voor rang-1 structuur bij stationaire frequentie . . . . . . . . . . Voorbeeld van een CPD analyse . . . . . . . . . . . . . . . . . . . . . . . . . Gereconstrueerde ictale activiteit . . . . . . . . . . . . . . . . . . . . . . . . Correlatie gerelateerd aan de vorm en amplitudeverdeling in functie van SNR (R = 2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Correlatie gerelateerd aan de vorm en amplitudeverdeling in functie van SNR (hogere R) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Relatieve fout van de CPD . . . . . . . . . . . . . . . . . . . . . . . . . . . . Vergelijking tussen CPD en SOBI . . . . . . . . . . . . . . . . . . . . . . . . Lokalisatiefout in functie van SNR bij stationaire frequentie . . . . . . . . Evolutie van de geschatte dipoollocatie bij stationaire frequentie . . . . . Simulatie met twee dipolen met stationaire frequentie . . . . . . . . . . . . Correlatie betreffende tijd-frequentie matrices en lokalisatiefout in functie van SNR bij twee dipolen met stationaire frequentie . . . . . . . . Evolutie van de twee geschatte dipoollocaties . . . . . . . . . . . . . . . . . Evidentie voor (L, L, 1) structuur bij veranderlijke frequenties . . . . . . . Tijd-frequentie matrices bij veranderlijke frequenties met CPD en BTD-(L, L, 1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10 2.11 2.12 2.13 2.14 2.15 2.16 2.17 2.18

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. 2 . 3 . 8 . 10 . 11 . 13

17 17 19 20 22 24 25 26 27 27 28 30 31 33 34 34 35 36 v

Lĳst van figuren en tabellen 2.19 Correlatie tussen tijd-frequentie matrices en amplitudeverdeling in functie van de rang van de beschrijvende term . . . . . . . . . . . . . . . . 2.20 Artefacten in tijd-frequentie matrix en signaal bij BTD-(L, L, 1) . . . . . 2.21 Details van originele, BTD-(L, L, 1) en CPD tijd-frequentie matrices, alsook van de gereconstrueerde signalen . . . . . . . . . . . . . . . . . . . . 2.22 Vergelijking tussen geclusterde BTD-(L, L, 1) en CPD . . . . . . . . . . . 2.23 Lokalisatiefout in functie van SNR bij veranderlijke frequenties . . . . . . 2.24 Correlatie tussen tijd-frequentie matrices in functie van SNR in het geval van meerdere bronnen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.25 Tijdreeks en afstand tot het pad van een bewegende bron . . . . . . . . . 2.26 Evolutie van de bronlokalisatie voor een bewegende bron . . . . . . . . . . 2.27 MLSVD en statistieken van een signaal met stationaire frequentie en veranderlijke locatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.28 Resultaten van de veranderlijke lokalisatie met BTD-(L, L, 1) . . . . . . . 2.29 Structuur van tensoren geconstrueerd met Hankel matrices . . . . . . . . . 2.30 Resultaten met tensoren geconstrueerd uit Hankel matrices . . . . . . . . 2.31 Resultaten met gecomprimeerde tensoren geconstrueerd uit Hankel matrices 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10

Eerste 20s van de EEG opname van patiënt 22 met MA 1 na de start van het epileptisch insult . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . MLSVD van patiënt met MA 1 . . . . . . . . . . . . . . . . . . . . . . . . Topografische kaart en tijdfrequentie matrix van de epileptische activiteit bij de patiënt met MA 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . Evolutie in lokalisatie van de epileptische aanval en vergelijking van de methodes bij de patiënt met MA 1 . . . . . . . . . . . . . . . . . . . . . . Uitvergroting van de verschillende lokalisaties en evolutie van de veranderingen per coördinaat voor het geval MA 1 . . . . . . . . . . . . . Eerste 20s van de EEG opname van patiënt 11 met MA 2 na de start van het epileptisch insult . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Topografische kaart en tijdfrequentie matrix van de epileptische activiteit bij de patiënt met MA 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . Vergelijking tussen de lokalisatie op basis van SISCOM data en op basis van BTD-(L, L, 1) voor patiënt met MA 2 . . . . . . . . . . . . . . . . . . Vergelijking tussen de lokalisatie op basis van SISCOM data en op basis van BTD-(L, L, 1) voor patiënt met MA 3 . . . . . . . . . . . . . . . . . . Evolutie van de veranderingen per coördinaat voor het geval MA 3 . . .

37 38 39 40 41 43 44 45 46 47 50 51 52

. 55 . 56 . 56 . 57 . 58 . 59 . 60 . 60 . 62 . 62

A.1 Vergelijking van verschillende methodes om de CWT en zijn inverse uit te voeren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 C.1 Poster voor de masterproefbeurs . . . . . . . . . . . . . . . . . . . . . . . . . 80 vi

Lĳst van figuren en tabellen

Lĳst van tabellen 2.1 2.2 2.3

Eigenschappen van het hoofdmodel . . . . . . . . . . . . . . . . . . . . . . . 16 Minimale SNR waarden voor bruikbare lokalisatie bij stationaire frequentie 30 Minimale SNR waarden voor bruikbare lokalisatie bij veranderlijke frequentie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

vii

Lijst van afkortingen en symbolen Afkortingen BSS SVD PCA MLSVD CPD BTD ALS NLS ICA CCA SOBI CWT EEG EMG SISCOM SPECT MRI SNR

viii

Blind Source Separation Singular Value Decomposition Principal Component Analysis Multilinear Singular Value Decomposition Canonical Polyadic Decomposition Block Term Decomposition Alternating Least Squares Nonlinear Least Squares Idenpendent Component Analysis Canonical Correlation Analysis Second Order Blind Identification Continuous wavelet transform Electoroencephalogram Electromyogram Subtraction ictal SPECT co-registered to MRI Single Photon Emission Computerized Tomography Magnetic Resonance Imaging Signal to Noise Ratio

Symbolen a, b, ... a, b, ... A, B, ... A, B, ... aijk , (A)ijk A(n) ∥A∥F K ⊗ ⊙ ×n ○

scalairen vectoren matrices tensoren element van A op positie ijk mode-n matrix uitvouwing van A Frobenius norm van A een algemeen veld (R of C) Kronecker product Khatri-Rao product mode-n product uitproduct

ix

Hoofdstuk 1

Literatuurstudie De literatuurstudie voor deze thesis is tweeledig. Enerzijds wordt het meer wiskundige aspect betreffende tensor decomposities toegelicht, maar anderzijds zijn ook de reeds behaalde resultaten op het gebied van het lokaliseren van epileptische aanvallen in de hersenen bekeken. We geven een overzicht hiervan in de volgende twee secties.

1.1

Een inleiding tot tensor decompositie technieken

Tensoren. Het zijn de wiskundige structuren waarvan we enkele eigenschappen in dit werk willen uitbuiten. We vangen deze sectie dus aan met enkele definities en begrippen. Als eerste de belangrijkste: een tensor is een meerdimensionale array. Een vector a ∈ KI is een eerste speciaal geval van een tensor, het is namelijk een eerste orde tensor. Hier wordt K gebruikt voor het neerschrijven van een algemener veld, het kan R of C zijn. Een tweede speciaal geval is dit van een tweede orde tensor, of een matrix A ∈ KI×J . Meer algemeen is A ∈ KI1 ×I2 ×⋯×IN een tensor van orde N . Tensoren van orde 3 of hoger noemen we tensoren van hogere orde. Voor het onderzoek in deze thesis is gebruik gemaakt van tensoren van orde 3. Voor we een inleiding geven tot enkele tensor decompositie technieken, geven we eerst nog een overzicht van de gebruikte notatie en enkele definities betreffende tensoren. Deze zijn allen geïnspireerd op de literatuur [16, 17, 20, 21]. Scalairen of getallen worden genoteerd met kleine letters (a), vectoren worden, zoals eerder, aangeduid met vette kleine letters (a), matrices met vette hoofdletters (A) en tensoren van orde drie of hoger met kalligrafische letters (A). Deze notatie wordt ook doorgetrokken naar lagere orde delen van een hogere orde structuur. Het element op de i-de rij en j-de kolom van een matrix A, ook wel (A)ij , wordt genoteerd als aij vermits het een scalair is. Analoog hieraan hebben we (a)i = ai en (A)ijk = aijk , met A een derde orde tensor. Omdat we meestal met tensoren van orde drie zullen werken, veronderstellen we dit vanaf nu en vermelden we de orde enkel nog expliciet indien deze van drie verschilt. De ide kolom vector van een matrix A wordt genoteerd als ai , dus A = [a1 a2 ⋯]. We geven nu enkele nuttige definities, 1

1. Literatuurstudie beginnende met een veralgemening van de begrippen kolommen en rijen van een matrix: Definitie 1.1. Een mode-n vector van een N de orde tensor A, is een vector die bekomen wordt uit A door de nde index te variëren terwijl de andere indices constant gehouden worden. Voor tensoren van orde drie hebben de drie verschillende types van mode-n vectoren ook nog een specifieke benaming zoals bij matrices: kolommen, rijen en tubes. Deze worden geïllustreerd in Figuur 1.1.

Figuur 1.1: Illustratie van (A) kolommen, (B) rijen en (C) tubes van een derde orde tensor. Figuur overgenomen uit [20].

Vervolgens definiëren we twee types van matrixproducten: Definitie 1.2. Het Kronecker product van twee matrices wordt genoteerd met ⊗ en gedefinieerd als, ⎡a11 B a12 B ⋯⎤ ⎢ ⎥ ⎢ ⎥ A ⊗ B = ⎢a21 B a22 B ⋯⎥ . ⎢ ⎥ ⎢ ⋮ ⎥ ⋮ ⎣ ⎦

Definitie 1.3. Het Khatri-Rao of kolomsgewijze Kronecker product van twee matrices wordt genoteerd met ⊙ en gedefinieerd als, A ⊙ B = [a1 ⊗ b1 a2 ⊗ b2 ⋯]

We definiëren ook het Tucker mode-n matrix product van een tensor [47]: Definitie 1.4. Het mode-n product van een tensor A ∈ KI1 ×I2 ×⋯×IN met een matrix U ∈ KJ×In wordt genoteerd met A ×n U en heeft als grootte I1 × ⋯ × In−1 × J × In+1 × ⋯ × IN . De elementen van het mode-n product zijn dan gedefinieerd als In

(A ×n U )i1 ⋯in−1 jin+1 ⋯iN = ∑ ai1 i2 ⋯in ⋯iN ujin . in =1

2

1.1. Een inleiding tot tensor decompositie technieken Definitie 1.5. Het uitproduct A ○ B van een tensor A ∈ KI1 ×⋯×IM met een tensor B ∈ KJ1 ×⋯×JN is de tensor gedefinieerd door (A ○ B)i1 ⋯iM j1 ⋯jN = ai1 ⋯iM bj1 ⋯jM ,

voor alle verschillende waarden van de indices.

De volgende definitie geeft een methode om een tensor door een matrix voor te stellen, dit wordt ook het uitvouwen van een tensor genoemd en komt erop neer dat de mode-n vectoren van de tensor de kolommen in de resulterende matrix worden: Definitie 1.6. We beschouwen een tensor A ∈ KI1 ×I2 ×I3 . De mode-1 matrix uitvouwing van deze tensor wordt genoteerd met A(1) ∈ KI1 ×I2 ⋅I3 , de mode-2 uitvouwing met A(2) ∈ KI2 ×I1 ⋅I3 en de mode-3 uitvouwing met A(3) ∈ KI3 ×I1 ⋅I2 . Ze zijn gedefinieerd door: (A(1) )i1 ,(i2 −1)I3 +i3 = (A)i1 i2 i3 ,

(A(2) )i2 ,(i3 −1)I1 +i1 = (A)i1 i2 i3 ,

(A(3) )i3 ,(i1 −1)I2 +i2 = (A)i1 i2 i3 ,

voor alle verschillende waarden van de indices.

Conceptueel is deze manier om een tensor om te zetten in een matrix eenvoudig, de notatie maakt het op het eerste zicht echter een stuk ingewikkelder. Duidelijker wordt het wanneer de operatie grafisch wordt afgebeeld zoals in Figuur 1.2.

Figuur 1.2: Grafische voorstelling van respectievelijk de mode-1, mode-2 en mode-3 uitvouwing van een tensor A. Figuur overgenomen uit [18].

De Frobenius norm van een tensor is vaak handig wanneer men met tensoren werkt. Deze is gedefinieerd als: 3

1. Literatuurstudie Definitie 1.7. De Frobenius norm van een tensor A ∈ KI1 ×I2 ×⋯×IN wordt genoteerd als ∥A∥F en is gedefinieerd als 2 IN ⎛ I1 I2 2⎞ ∥A∥F = ∑ ∑ ⋯ ∑ ∣ai1 i2 ⋯iN ∣ . ⎝i1 =1 i2 =1 iN =1 ⎠ 1

1.1.1

Matrixalgebra

Om er nu toe te komen waarom we voorgaande kader rond tensoren opgetrokken hebben, kijken we even terug naar het doel dat we voor ogen hebben en bekijken we de mogelijke oplossingen die een aanpak met matrices kan bieden. Het uiteindelijke doel van het onderzoek is om op een betrouwbare manier data van een elektro-encefalogram (EEG) van een patiënt te ontbinden om zo de verschillende bronnen aanwezig in de data te scheiden. Alle details betreffende het uitvoeren van EEG opnames met patiënten worden in het tweede deel van de literatuurstudie besproken, alsook enkele typische eigenschappen van epileptische aanvallen. Op dit moment volstaat het te weten dat een EEG opname met verscheidene elektroden wordt gemaakt zodat de resulterende datasets matrices zijn met dimensies: kanalen × tijd. Ook is het reeds interessant te vermelden dat een epileptische aanval vaak gekarakteriseerd wordt als een periodische golf in het EEG met een karakteristieke frequentie waarvan de bron vaak sterk gelokaliseerd is in de hersenen. De periodische golf kan van frequentie veranderen gedurende de aanval en de lokalisatie van de bron kan wijzigen. Het scheiden van de verschillende bronnen die aanwezig zijn in een gemengd signaal onder zo weinig mogelijk voorwaarden die opgelegd dienen te worden, wordt in de literatuur met de term Blind Source Separation (BSS) [9] aangeduid. Deze scheiding in bronnen is in feite een factorisatie van de data en hiervoor zijn verschillende technieken ontwikkeld, zowel voor het ontbinden van matrices als, meer algemeen, voor het ontbinden van tensoren. Als we de datamatrix van de EEG opname noteren met A, dan wensen we deze te ontbinden in rang-1 termen, als volgt1 : A = F ×2 G,

met F de mixing matrix en G de matrix die de bronsignalen als kolommen bevat. Het probleem met deze ontbinding in rang-1 termen is echter dat deze niet uniek is, we hebben namelijk:

1

A = (F M ) ⋅ (M −1 GT ) ˜ ×2 G. ˜ =F

Merk op dat we hier het mode-n product uit Definitie 1.4 gebruikt hebben om de vermenigvuldiging van een matrix met een getransponeerde matrix te noteren

4

1.1. Een inleiding tot tensor decompositie technieken ˜ en G ˜ alternatieve mixing- en bron-matrices. Een manier om de ontbinHierin zijn F ding van een matrix uniek uit te voeren, is met een singuliere waarden ontbinding (SVD). Deze wordt gedefinieerd als volgt: Definitie 1.8. Beschouw een matrix A ∈ KI×J , de SVD van deze matrix is een factorisatie van de vorm A = U ⋅ Σ ⋅ V T = Σ ×1 U ×2 V ,

met U ∈ KI×I en V ∈ KJ×J orthogonale matrices met als kolommen respectievelijk de linkse singuliere vectoren en de rechtse singuliere vectoren. Σ ∈ KI×J een nietnegatieve diagonaal matrix met als diagonaal elementen σi de singuliere waarden van de matrix A. Deze ontbinding is uniek dankzij de orthogonaliteitsvoorwaarden die opgelegd worden, en we kunnen ze interpreteren als een ontbinding in rang-1 termen als we ze herschrijven als R

A = ∑ σr ur v Tr , r=1

waarin R de rang is van de matrix A en ui en v i de i-de kolom is van respectievelijk de matrix U en V . We geven nu twee equivalente definities van de rang van een matrix A: Definitie 1.9. De rang van een matrix A is R als: 1. Het aantal lineair onafhankelijke rijen van A, dat ook gelijk is aan het aantal lineair onafhankelijke kolommen van A, gelijk is aan R. 2. Het minimaal aantal rang-1 termen waarvan de som gelijk is aan A gelijk is aan R. De SVD heeft enkele interessante eigenschappen, zo is de SVD beperkt tot ρ termen, met ρ ≤ R de beste rang-ρ benadering van A. Dit resultaat staat bekend als de stelling van Eckart-Young [23]. De beperkingen om een unieke factorisatie te bekomen zijn vrij streng en voor sommige problemen niet zinvol. Toch bestaan er BSS methoden die hierop gebaseerd zijn, zoals de principale component analyse (PCA) [23]. Zoals we zullen zien in wat volgt, kan dit probleem verholpen worden door over te schakelen op tensoren.

1.1.2

Tensoralgebra

Alvorens we enkele tensor decompositie technieken introduceren, zijn er twee problemen die we dienen te bespreken. Ten eerste, heeft het in het kader van een EEG opname zin om tensor decomposities te gebruiken? De data die we verkrijgen uit een EEG is namelijk tweedimensionaal. Het blijkt echter dat het heel nuttig is 5

1. Literatuurstudie met het oog op het lokaliseren van epileptische aanvallen in de hersenen om tensor decomposities te beschouwen [1, 22, 21, 20]. Hiervoor moet de tweedimensionale EEG data dan wel eerst getensoriseerd worden. Een mogelijke methode om een tensor te construeren is met behulp van een continue wavelet transformatie (CWT) die de kanalen × tijd data omzet in kanalen × tijd × frequentie data, zoals in [1, 22]. Een andere mogelijke methode is het Hankelizeren van de data, hierbij wordt van elk kanaal een Hankel matrix geconstrueerd die vervolgens in een derde orde tensor geplaatst wordt [16]. Een tweede probleem is het veralgemenen van het begrip rang naar tensoren. Dit is minder eenduidig te definiëren als voor matrices vermits in het algemeen niet geldt dat elke mode-n rang gelijk is, vandaar een eerste definitie van de rang [16] Definitie 1.10. Een derde orde tensor A ∈ KI1 ×I2 ×I3 heeft multilineaire rang (R1 , R2 , R3 ) als en slechts als rang(A(1) ) = R1 , rang(A(2) ) = R2 en rang(A(3) ) = R3 . De waarden R1 , R2 en R3 zijn respectievelijk de mode-1, mode-2 en mode-3 rang. Een tweede veralgemening van het begrip rang naar tensoren is gegeven door: Definitie 1.11. De rang van een tensor A is het minimaal aantal rang-1 tensoren die A vormen als een lineaire combinatie. Hierbij is een rang-1 tensor een tensor die te schrijven is als het uitproduct van niet-nul vectoren. We gaan nu over tot het bespreken van enkele tensor decompositie technieken, en zullen ons toespitsen op drie verschillende ontbindingen: de multilineaire SVD (MLSVD), de Canonical Polyadic Decompositie (CPD) en de ontbinding in bloktermen (BTD). 1.1.2.1

MLSVD

Definitie 1.12. De MLSVD [18] van een derde orde tensor A ∈ KI1 ×I2 ×I3 van rang-(R1 ,R2 ,R3 ) is een ontbinding van A van de vorm

met

A = S ×1 U (1) ×2 U (2) ×3 U (3)

• de orthogonale of unitaire matrices U (1) ∈ KI1 ×I1 , U (2) ∈ KI2 ×I2 en U (3) ∈ KI3 ×I3 • de tensor S ∈ KI1 ×I2 ×I3 , die

– geheel-orthogonaal of geheel-unitair is, dit wil zeggen: ⎛ I2 I3 ⎞2 (1) ∗ = σi1 als j1 = j2 = i1 ∑ ∑ sj1 i2 i3 sj2 i2 i3 ⎝i2 =1 i3 =1 ⎠ 1

6

= 0 als j1 ≠ j2

1.1. Een inleiding tot tensor decompositie technieken

⎛ I1 I3 ⎞2 (2) ∗ = σi2 als j1 = j2 = i2 ∑ ∑ si1 j1 i3 si1 j2 i3 ⎝i1 =1 i3 =1 ⎠ 1

= 0 als j1 ≠ j2

⎞2 ⎛ I1 I2 (3) ∗ = σi3 als j1 = j2 = i3 ∑ ∑ si1 i2 j1 si1 i2 j2 ⎠ ⎝i1 =1 i2 =1 1

= 0 als j1 ≠ j2

– geordend is: (1)

σ1

(2)

σ1

(3)

σ1 (1)

(2)

≥ σ2

(1)

≥ σ2

(2)

≥ σ2

(3)

(3)

≥ ⋯ ≥ σR1 > σR1 +1 = ⋯ = σI1 = 0, (1)

(1)

(1)

≥ ⋯ ≥ σR2 > σR2 +1 = ⋯ = σI2 = 0, (2)

(2)

(2)

≥ ⋯ ≥ σR3 > σR3 +1 = ⋯ = σI3 = 0. (3)

(3)

(3)

De waarden σi1 , σi2 en σi3 zijn respectievelijk de mode-1, mode-2 en mode-3 singuliere waarden. De kolommen van U (1) , U (2) en U (3) zijn respectievelijk de mode-1, mode-2 en mode-3 singuliere vectoren. 1.1.2.2

CPD

De oorsprong van de Canonical Polyadic Decompositie (CPD) begint in 1927 wanneer Hitchcock de polyadische vorm van een tensor introduceert [28]. De polyadische vorm van een tensor is de voorstelling van een tensor als de som van een eindig aantal rang-1 tensoren. Het duurde echter nog tot de jaren 1970 tot het concept aanhang vond in het psychometrisch onderzoek onder de naam Candecomp [7] en Parafac [26]. Nadien vond het concept meerdere toepassingen, onder andere in de chemometrie [6]. De CPD van een tensor wordt als volgt gedefinieerd: Definitie 1.13. Een CPD van een rang-R tensor A ∈ KI1 ×I2 ×I3 is een ontbinding van A in een lineaire combinatie van R rang-1 termen: R

A = ∑ dr ○ er ○ f r . r=1

Indien de rang-R tensor niet ontbonden wordt in R rang-1 termen, maar in minder dan wordt een deel van de variatie niet verklaard en in een residu tensor bevat: Q

A = ∑ dr ○ er ○ f r + E. r=1

7

1. Literatuurstudie

Figuur 1.3: CPD van een tensor A in R rang-1 componenten

De factor matrices van de ontbinding zijn de combinatie van de vectoren: D = [d1 d2 ⋯ dR ], analoog voor E en F . Met deze notatie kunnen we de CPD van een tensor ook schrijven voor de mode-n uitvouwing (Definitie 1.6) van deze tensor als: A(1) = D (E ⊙ F )T .

Het voordeel van een CPD ontbinding in vergelijking met matrix factorisatie is dat ze uniek is onder milde voorwaarden. In [41] wordt een kort overzicht gegeven van de resultaten betreffende de uniciteit van de CPD ontbinding. De eerste gepubliceerde uniciteitsresultaten zijn deze van Harshman in 1970 [26, 27]. Andere resultaten volgden hier snel op [40], uiteindelijk culmineerde dit in het bekende resultaat van Kruskal [31] dat stelt dat de CPD ontbinding essentieel uniek is als: rangk (D) + rangk (E) + rangk (F ) ≥ 2R + 2.

Hierin is rangk de Kruskal rang van de betreffende matrix. Deze is gedefinieerd als de maximale r waarvoor geldt dat elke verzameling van r kolommen van de matrix lineair onafhankelijk is. De essentiële uniciteit houdt in dat een tensor van rang R maar als een enkele combinatie van rang-1 tensoren geschreven kan worden op uitzondering van een permutatie en schaling operatie. Deze operaties noemt men ook wel elementaire onbepaaldheden. In de literatuur zijn er meerdere uniciteitsvoorwaarden over CPD terug te vinden, die elk andere veronderstellingen maken om het resultaat van Kruskal te verbeteren. Zoals bijvoorbeeld uniciteitsresultaten in [41, 46, 33, 13]. Voor de berekening van de CPD van een tensor zijn verschillende algoritmes ontwikkeld. De klassieke aanpak om de CPD van een tensor te berekenen is een Alternating Least Squares (ALS) algoritme. Deze methode werd reeds voorgesteld in [7] en bestaat eruit dat elke factor matrix in de CPD afwisselend wordt herberekend via een Least Square (LS) update terwijl de anderen constant gehouden worden. Andere methodes om de CPD van een tensor te berekenen zijn voorgesteld in de literatuur. Zo is er een Nonlinear LS (NLS) procedure [44] die op een robuuste en efficiënte manier de ontbinding berekent. De CPD kan ook berekend worden onder verschillende randvoorwaarden zoals bijvoorbeeld niet-negativiteit van de 8

1.1. Een inleiding tot tensor decompositie technieken factoren. Verschillende algoritmes voor de berekening van de CPD ontbinding zijn geïmplementeerd in de Tensorlab toolbox [45] die gebruikt wordt voor de data-analyse in deze thesis, zowel voor de CPD als voor de tensor decompositie die we in wat volgt introduceren. 1.1.2.3

BTD

De CPD van een tensor is soms moeilijk te berekenen doordat het probleem slecht geconditioneerd is. Dit heeft tot gevolg dat er degeneratie ontstaat wat betekent dat de termen in de ontbinding naar oneindig gaan maar dat in de som dit effect gedeeltelijk opgeheven wordt zodat de fit toch verbetert. Dit probleem kan ontstaan als de onderliggende componenten in de tensor niet rang-1 zijn. Een oplossing voor dit probleem is bijvoorbeeld een randvoorwaarde van niet-negativiteit opleggen bij de ontbinding. Beter is om de afwijking van rang-1 in de componenten te brengen. Dit leidt ons tot een zogenaamde Block Term Decomposition (BTD) van de tensor. We bekijken meer bepaald de ontbinding van een tensor in rang-(Lr , Lr , 1) termen, geïntroduceerd in [14, 15, 19]. Zoals de naam doet vermoeden, is dit een ontbinding waarbij één van de componenten rang-1 is, terwijl de andere twee componenten van lage rang zijn, hetgeen uiteindelijk een ontbinding oplevert in componenten van lage multilineaire rang. Wanneer er slechts één component is in de ontbinding, dan reduceert deze zich tot de MLSVD. Als voor elke component Lr gelijk is aan 1, dan hebben we te maken met een CPD. Een BTD in rang-(Lr , Lr , 1) termen is gedefinieerd als [16]: Definitie 1.14. Een ontbinding van een tensor A ∈ KI1 ×I2 ×I3 in een som van rang(Lr , Lr , 1) termen (1 ≤ r ≤ R) is een ontbinding van de vorm A = ∑ (D r ⋅ E Tr ) ○ f r , R

r=1

waarbij de matrix Gr = D r ⋅ E Tr ∈ KI1 ×I2 van rang ≤ Lr is en de vector f r verschillend is van nul. Het aantal componenten R wordt minimaal verondersteld. De decompositie is schematisch weergegeven in Figuur 1.4. In de literatuur zijn enkele voorwaarden afgeleid waaronder de ontbinding essentieel uniek is [16]. Indien de ontbinding uitgevoerd wordt in minder dan R componenten dan kan een deel van de data ook nog in een residu tensor blijven zoals eerder bij de CPD. Algoritmes om de BTD-(Lr , Lr , 1) te berekenen op basis van ALS en NLS procedures zijn beschikbaar in de Tensorlab toolbox [45]. Een nog algemenere ontbinding in Block Termen van lage multilineaire rang is een ontbinding in rang-(L, M, N ) termen [15]. Deze ontbinding is echter niet beschikbaar in Tensorlab. In [30] wordt een meer uitvoerig overzicht gegeven van een groot deel van de hier voorgestelde theorie, ook worden een aantal andere tensorontbindingen toegelicht. Een ander standaardwerk in het domein van tensoren is [42].

9

1. Literatuurstudie

Figuur 1.4: Ontbinding van een tensor A in rang-(Lr , Lr , 1) termen. Figuur overgenomen uit [16].

Zowel CPD als BTD worden gebruikt als BSS technieken voor het scheiden van de aanwezige componenten in EEG data. Het doel is om op een betrouwbare manier de component, te wijten aan epileptische activiteit in de hersenen, te extraheren uit de data zodat lokalisatie en karakterisering van de epileptische aanval vergemakkelijkt wordt. Het voordeel van de ontwikkeling van een betrouwbare lokalisatiemethode op basis van EEG data is ten eerste dat EEG metingen relatief gemakkelijk uit te voeren zijn. Ook hebben ze een nauwkeurige tijdsresolutie en kunnen de metingen over langdurige periodes uitgevoerd worden.

1.2

Epilepsie: eigenschappen, monitoring en analyse

Het brein is opgebouwd uit een enorme verzameling van zenuwcellen of neuronen. Deze neuronen staan via een ingewikkeld netwerk met elkaar in verbinding en communiceren met elkaar door chemische signalen uit te wisselen aan hun synapsen. Deze chemische uitwisseling zorgt ervoor dat het ontvangende neuron anders gepolariseerd wordt. Als het effect groot genoeg is, genereert het neuron een actiepotentiaal dat via zijn axon, wat een soort geleider is die de verbinding legt met andere neuronen, verspreid wordt naar andere regionen in de hersenen. EEG metingen hebben net als doel om deze elektrische signalen te registreren in functie van de tijd en van plaats. Ze kunnen echter enkel een verandering gewaarworden wanneer een groot aantal neuronen gelijktijdig gepolariseerd worden. De spanning die gegenereerd wordt door een groep neuronen, moet wel nog doordringen door het hersenvocht, de schedel en de hoofdhuid, maar kan op de hoofdhuid nog steeds gedetecteerd worden. Een typische EEG meetopstelling wordt in Figuur 1.5 voorgesteld. Hier worden 21 elektrodes op welbepaalde plaatsen op de hoofdhuid geplaatst. Dit is de elektrodepositionering gebruikt voor de metingen die de beschikbare dataset uitmaken. Deze opstelling is gebaseerd op het bekende 10-20 systeem [36] met T 1 en T 2 als twee extra elektrodes. Merk op dat de naamgeving van de elektrodes steeds gebaseerd is 10

1.2. Epilepsie: eigenschappen, monitoring en analyse op de onderliggende hersenkwab en dat de linker-hersenhelft een oneven nummering heeft, terwijl de rechter-hersenhelft met even nummers aangeduid wordt.

Figuur 1.5: Plaatsing van de 21 elektroden op het hoofd volgens de opstelling gebruikt bij de beschikbare dataset.

Een EEG tracht de onderliggende golfvorm op te meten die de hersenactiviteit karakteriseert. De interessante frequentieband gaat ruwweg van 0 − 30 Hz en is opgedeeld in vier delen: • β : 13 ≤ f < 30 Hz. Dit type golven is gerelateerd aan actieve gedachten en treden op wanneer een persoon geconcentreerd met een taak bezig is.

• α : 8 ≤ f < 13 Hz. Deze golven indiceren dat een persoon op ontspannen wijze alert is. • θ : 4 ≤ f < 8 Hz. Deze golven komen voor wanneer je in een lichte slaap bent. • δ : 0.5 ≤ f < 4 Hz. Deze golven duiden op een diepe slaap.

Het belangrijkste type hersenactiviteit van toepassing is deze gerelateerd aan epileptische aanvallen. Epilepsie is een neurologische aandoening die gekarakteriseerd wordt door abnormale synchrone ontladingen van neuronen in de hersenen [4]. Dit is geen heel duidelijke definitie van het fenomeen wat te wijten is aan de grote diversiteit waarin het zich manifesteert. Een EEG opname van een epilepsie patiënt kan ruwweg opgedeeld worden in twee delen. Ten eerste is er de ictale periode, dit is de periode waarin er een epileptische aanval plaatsvindt. De periode tussen aanvallen wordt de interictale periode genoemd. De ictale periode is herkenbaar in het EEG als een periode waarin er continu ritmische golven of herhaaldelijk scherpe pieken waar te nemen zijn. De ritmische golven zijn eerder typisch voor volwassen patiënten, terwijl de pieken meer voorkomen bij pasgeborenen met epilepsie. De ritmische activiteit kan sterk veranderen 11

1. Literatuurstudie gedurende de epileptische aanval, zowel in frequentie, amplitude en ruimtelijke situering. Typisch is dat ritmische activiteit toeneemt in amplitude en daalt in frequentie. Als EEG opnames gedurende een epileptische aanval geregistreerd worden, dan zal niet enkel de ictale activiteit vastgelegd worden maar tegelijkertijd worden ook een deel artefacten geregistreerd. Dit is waar BSS technieken van toepassing zijn om de ictale activiteit zo goed mogelijk te scheiden van de artefacten. Vooraleer een overzicht te geven van de BSS methodes die reeds toegepast zijn om het probleem op te lossen, wordt er een overzicht gegeven van enkele types artefacten die gedurende epileptische aanvallen op kunnen treden. Vermits een epileptische aanval vaak gepaard gaat met het onbeheersbaar samentrekken van verschillende spieren in het gelaat en elders, zijn musculaire artefacten vaak voorkomend in EEGs. Dit type artefacten staat ook bekend als ElectroMyoGram (EMG) artefacten en het heeft een brede spectrale signatuur daar het van 0 tot 200 Hz kan gaan [24]. Het spectrum van de artefacten overlapt dus met de frequentieband waarin we geïnteresseerd zijn. Een gewone filter zal dus niet in staat zijn om de componenten te scheiden. In [12] wordt Canonical Correlation Analysis (CCA) [5] succesvol toegepast om het BSS probleem op te lossen en de component gelinkt aan spieractiviteit te verwijderen. Een tweede type artefacten is gerelateerd aan bewegingen van de ogen. Zowel het knipperen van de oogleden als het bewegen van de oogbol heeft een effect op de EEG meting en dan vooral aan de Fp1 en Fp2 elektroden. Dit type artefacten manifesteert zich typisch als een kortstondige grote uitwijking bij de frontale elektroden. Ook hier zijn verschillende technieken ontwikkeld [29, 35, 48] om de ongewenste component uit de data te verwijderen opdat de interpretatie vergemakkelijkt wordt. Andere types artefacten spelen een minder significante rol en kunnen verder opgedeeld worden in twee categorieën, namelijk artefacten van fysiologische oorsprong (zoals de twee voorgaande) en artefacten van niet-fysiologische oorsprong (zoals bijvoorbeeld interferentie van het spanningsnet). In [1, 22] wordt de aanpak om het BSS probleem op te lossen omgekeerd. De ongewenste artefacten worden niet verwijderd uit de data, maar de component gerelateerd aan de epileptische activiteit wordt rechtstreeks uit de dataset gehaald met behulp van CPD. Zoals reeds eerder gesteld, is EEG data tweedimensionaal, terwijl tensorontbindingen zoals CPD en BTD van toepassing zijn op hogere orde tensoren. In [22] lossen ze dit op door de EEG data te tensoriseren via een CWT. Dit proces is schematisch weergegeven in Figuur 1.6. Er zijn ook andere methoden om van een EEG matrix een tensor te construeren, zo is het mogelijk om van ieder kanaal apart een Hankel matrix te construeren en deze in een tensor te plaatsen. Naar andere tensorisaties is echter nog niet veel onderzoek gedaan in de literatuur, dit omdat de CWT een duidelijke interpretatie van de componenten toelaat en een rang-1 structuur oplevert voor sinusoïdale signalen. Na de ontbinding bestaat iedere term uit een ruimtelijke-, een tijds- en een frequentiecomponent. In [16] is echter aangetoond dat de BTD-(Lr , Lr , 1) ontbinding essentieel uniek is indien de tensor uit Hankel matrices is opgebouwd en de onderliggende bronnen lineaire combinaties zijn van exponentiëlen of algemener exponentiële polynomen. Dit zet ertoe aan om 12

1.2. Epilepsie: eigenschappen, monitoring en analyse dit type tensorisatie nader te bekijken in het kader van dit onderzoek.

Figuur 1.6: Tensoriseren van een EEG matrix naar een derde orde tensor via CWT.

In [22] wordt er, zowel door een simulatiestudie als door onderzoek op een klinische dataset, aangetoond dat een CPD ontbinding van een EEG tensor de correcte lokalisatie van de epileptische aanval kan bepalen. Voor analyse van de klinische data werd er steeds gewerkt met een kort tijdssegment van vlak na de start van de ictale periode. Het is dus belangrijk om op te merken dat de CPD analyse geen aanvallen detecteert in tijd, maar dat het een handige manier is om de aanval beter te identificeren. Er kwamen wel enkele tekortkomingen van de CPD analyse naar voren bij [21]. Zo blijkt uit de simulatiestudie dat bij een simulatie met een sinusoïdaal signaal alle componenten perfect het onderliggende bronsignaal beschrijven. Wanneer de simulatie uitgevoerd wordt met een signaal dat van frequentie verandert in de tijd, dan blijkt dat de lokalisatie correct blijft, maar de andere componenten minder goed gemodelleerd worden met CPD. Dit is te verklaren doordat in dit geval de CWT van het bronsignaal geen rang-1 matrix meer is, waardoor het gebruik van één component in de CPD decompositie als ictale term niet meer gerechtvaardigd is. In dit geval is er een verbetering mogelijk door een decompositie in termen van lage multilineaire rang te beschouwen. Dat dit een realistisch scenario is, blijkt uit het feit dat epileptische aanvallen vaak van frequentie veranderen in de tijd. Een ander geval waar BTD een voordeel kan bieden in vergelijking met CPD, is als de aanval van lokalisatie verandert. Een verspreiding van ictale activiteit over de hersenen in de loop van de aanval is niet ongewoon en dit is een van de redenen waarom er in [22] slechts met beperkte tijdssegmenten werd gewerkt. Dit zijn enkele van de onderzoeksvragen die in deze thesis onderzocht zullen worden.

13

Hoofdstuk 2

Simulatiestudie De verschillende simulaties die uitgevoerd zijn om de methodes, uitgelijnd in het vorige hoofdstuk, te testen, vormen het onderwerp van dit hoofdstuk. Vooreerst wordt in Sectie 2.1 de methode toegelicht voor het genereren van artificiële EEG data die een epileptische aanval of insult nabootsen. De constructie van de tensor met CWT wordt in meer detail besproken in Sectie 2.2. De ruis die gebruikt is om de simulaties meer realistisch te maken, vormt het onderwerp van Sectie 2.3. De methode wordt in Secties 2.4 en 2.5 uitgebreid getest in verschillende situaties. Het verschil tussen stationaire en veranderlijke frequenties wordt bekeken, alsook de invloed van het type ruis, de aanwezigheid van meerdere bronnen en het resultaat voor bewegende bronnen. Ook wordt de methode vergeleken met een andere BSS techniek. In Sectie 2.6 wordt ten slotte een andere methode gebruikt voor het construeren van de tensor uit EEG data, deze gebruikt Hankelmatrices. In het volgende hoofdstuk worden de methodes met CWT en Hankel tensoren, die reeds getest zijn op artificiële EEG data, uitgevoerd op enkele voorbeelden uit een beschikbare klinische EEG dataset.

2.1

Het genereren van artificiële epileptische insulten

Het genereren van artificiële EEG data met ictale activiteit gedurende een epileptisch insult is uitgevoerd met de FieldTrip toolbox [38] voor Matlab. Een beschrijving van de gebruikte functies en software is in Bijlage A gegeven. Conceptueel werkt het als volgt: een dipool wordt in een meerlagig sferisch hoofdmodel geplaatst en het voorwaartse probleem wordt opgelost om de intensiteiten aan de elektrodelocaties te bepalen. De dipool wordt gedefinieerd aan de hand van de locatie in het hoofd, het dipoolmoment en de gebruikte golfvorm. De invloeden van de locatie van de epileptische bron in het hoofd, alsook van de vorm van de tijdreeks (sinusoïdaal, frequentieveranderend, ...) worden besproken in wat volgt. Het gebruikte hoofdmodel is een volumegeleider bestaande uit drie concentrische sferen. Dit zijn de drie lagen van het model. Hierbij komt de binnenste en grootste laag met de hersenen, de tweede met de schedel en de buitenste met de hoofdhuid overeen. De straal van deze buitenste laag valt samen met de straal, relek , waarop de elektroden geplaatst 15

2. Simulatiestudie zijn. In Tabel 2.1 is een overzicht gegeven van de verdeling van het hoofd in drie lagen [22] en de geleidbaarheid van elke laag [37]. De precieze straal van het hoofd, relek , is onbelangrijk voor de simulatie, maar is bruikbaar om afstanden met elkaar te vergelijken. In navolging van [22] wordt relek = 9.2 cm als standaardwaarde gekozen.

Regio Hersenen Schedel Hoofdhuid

Straal (relek ) 0.88 0.92 1.00

Geleidbaarheid (× 3.3 10−4 /Ω/mm) 1 1 16

1

Tabel 2.1: De eigenschappen van het gebruikte hoofdmodel.

De enige noodzakelijke keuze die ons rest alvorens het model opgesteld kan worden, is de locatie van de elektroden op de hoofdhuid. Hiervoor zijn 21 elektrodes gebruikt, geplaatst volgens het 10-20 systeem [36] met bijkomende elektroden T1 en T2 geplaatst ter hoogte van de temporale kwab. Deze opstelling is reeds vermeld in Sectie 1.2 en het is de gebruikte opstelling voor de beschikbare klinische dataset. Verder wordt dezelfde bemonsteringsfrequentie, fs = 250 Hz, als voor de klinische data verondersteld. Ter illustratie zijn in Figuur 2.2 drie voorbeelden van een simulatie van ictale activiteit weergegeven. Voor elke simulatie is de (x, y, z)-locatie van de dipool gelijk aan (0.3, −0.3, 0.3), deze positie is in Figuur 2.1 met rood aangeduid. Het assenstelsel is zo gedefinieerd dat de oorsprong zich recht onder de Cz elektrode bevindt, meer bepaald centraal in het hoofd. De neus is in de +x richting gesitueerd, het linkeroor in de +y richting en het rechteroor in de −y richting. Het dipoolmoment is voor de simulatiestudie steeds in de +x richting gehouden, zoals in [22]. De segmenten in Figuur 2.2 zijn respectievelijk gegenereerd met een 6 Hz sinus en een frequentieveranderend signaal (ook chirp genoemd) waarbij de frequentie gedurende het segment afvalt van 12 tot 6 Hz. Het eerste is reeds uitgebreid bestudeerd in een eerdere studie met CPD [22], bij signalen die van frequentie veranderen neemt de correctheid van de modellering met CPD af. De hypothese die het uitgangspunt vormt voor dit hoofdstuk is dat een modellering met BTD meer gerechtvaardigd is, zoals reeds geargumenteerd in Hoofdstuk 1. Merk op dat ritmische ictale activiteit gedurende epileptische insulten zelden een perfect stationaire frequentie heeft. Typisch neemt de frequentie af tijdens het insult. Een halvering in frequentie over een tijdspanne van 2s is echter wel weinig waarschijnlijk in klinische EEG opnames.

16

2.1. Het genereren van artificiële epileptische insulten

1

F7 Left Ear T3 P3 Pz

0.5

Fz

F3 Cz

C3

Fp1 NoseF4 C4Fp2

z

T2

F8

T5

P4

0

O1 T6

O2

−0.5 1

Right Ear T4 1 T1 0.5 0

0.5 0

−0.5 −0.5

x

−1 −1

y

Figuur 2.1: De dipoollocatie in het hoofd (rood) gebruikt voor de simulaties in Figuur 2.2. De elektrodeposities en -benamingen zijn aangegeven, alsook de locatie van de neus, en het linker- en rechteroor.

b: Chirp (12-6 Hz)

a: Sinus (6 Hz) Fp2 F8 T4 T6 O2 F4 C4 P4 Fz Cz Pz Fp1 F7 T3 T5 O1 F3 C3 P3 T2 T1

7 uV

0

Fp2 F8 T4 T6 O2 F4 C4 P4 Fz Cz Pz Fp1 F7 T3 T5 O1 F3 C3 P3 T2 T1

0.5

1

Time (sec)

1.5

2

7 uV

0

0.5

1

1.5

2

Time (sec)

Figuur 2.2: De resultaten van een dipoolsimulatie in een sferisch hoofdmodel. De dipool is steeds op positie (0.3, −0.3, 0.3) geplaatst zoals aangegeven in Figuur 2.1. Twee verschillende tijdreeksen zijn gebruikt.

17

2. Simulatiestudie

2.2

Constructie van de tensor met een continue wavelet transformatie

Zoals reeds eerder aangehaald in Sectie 1.2 is een eerste methode om de derde orde tensor X te construeren uit de datamatrix X het berekenen van een continue wavelet transformatie (CWT) van elk kanaal in de datamatrix. Alvorens deze transformatie en de tensorontbinding uitgevoerd worden, wordt de EEG matrix genormaliseerd. Dit gebeurt door het gemiddelde van elk kanaal af te trekken en het te delen door zijn standaardafwijking. Zo heeft elk kanaal als gemiddelde 0, en een eenheidsstandaardafwijking. Na de tensorontbinding wordt de ruimtelijke component opnieuw vermenigvuldigd met de vector van standaardafwijkingen, zodat de amplitude-informatie bewaard blijft. Een CWT is een redundante voorstelling van het oorspronkelijke signaal in het tijd-frequentie domein. Om de inverse transformatie te kunnen definiëren wordt gebruik gemaakt van de integraalvoorstelling van Morlet [10]. De gekozen schalen of frequenties in de transformatie zijn steeds lineair verdeeld tussen fmin en fmax . Voor fmax is als standaardwaarde 30 Hz gekozen omwille van twee redenen. Vooreerst is de klinische EEG data digitaal gefilterd met een band-pass filter met een passband van 0.3 − 35 Hz [22]. Ten tweede is 30 Hz reeds een zeer hoge frequentie om realistisch nog ictale activiteit op terug te vinden, daar het zich typisch op lagere frequenties voordoet. Verder is 1 Hz als standaard gekozen voor zowel fmin als de frequentieresolutie. Er worden dus 30 schalen tussen 1 en 30 Hz gebruikt. Als wavelet is er voor de mexican hat wavelet geopteerd, die zijn naam ontleent aan de treffende gelijkenis met een typisch Mexicaans hoofddeksel, en gedefinieerd is als de tweede afgeleide van de Gaussische functie. Deze keuze is gemaakt omdat het ten eerste een reële wavelet is, wat de interpretatie van de verschillende componenten in de tensorontbinding duidelijk bepaalt. En ten tweede omdat het oscillerende signalen duidelijk benadrukt, het is ook reeds in eerder onderzoek naar epilepsie gebruikt [1, 32]. In Bijlage A is een overzicht gegeven van de gebruikte methode om de CWT en inverse CWT te bepalen.

2.3

De ruiscomponent in simulaties

In klinische EEG data opgemeten gedurende epileptische insulten is er nooit enkel ictale activiteit aanwezig. Zoals besproken in Sectie 1.2 van Hoofdstuk 1 zijn er vaak verschillende types artefacten aanwezig gedurende de ictale periode. De twee belangrijkste types artefacten van fysiologische oorsprong benoemen we als musculaire artefacten en oculaire artefacten. Om een zo realistisch mogelijke simulatiestudie uit te voeren moet er dus nog ruis toegevoegd worden aan de gegenereerde ictale activiteit. Hiervoor zijn twee segmenten uit de beschikbare klinische dataset geëxtraheerd. In navolging van [11, 12] wordt 18

2.3. De ruiscomponent in simulaties Canonical Correlation Analysis (CCA) [5] als BSS techniek gebruikt om de gewenste component uit de data te scheiden. CCA lost het BSS probleem op door bronnen te construeren die onderling minimaal gecorreleerd zijn, maar een maximale autocorrelatie bezitten. Spieractiviteit kan op deze manier van de aanwezige hersenactiviteit gescheiden worden vermits de bronnen gerelateerd met spieractiviteit typisch een veel lagere autocorrelatie hebben dan bronnen gerelateerd aan hersenactiviteit [12]. In Figuur 2.3a wordt het gebruikte segment musculaire activiteit weergegeven. Het is geconstrueerd uit de 9 bronnen met laagste autocorrelatie uit de EEG meting van patiënt 27 van de klinische dataset (zie Hoofdstuk 3) die sterk gecontamineerd is met musculaire artefacten na het begin van het insult. In Figuur 2.3b wordt het segment met oculaire artefacten weergegeven. Ook dit is geëxtraheerd met behulp van CCA uit een EEG meting, ditmaal van patiënt 39 uit de klinische dataset. Het is geconstrueerd uit één enkele bron uit de CCA analyse die gerelateerd is aan het knipperen van de oogleden.

a: Musculaire artefacten Fp2 F8 T4 T6 O2 F4 C4 P4 Fz Cz Pz Fp1 F7 T3 T5 O1 F3 C3 P3 T2 T1

b: Oculaire artefacten

179 uV

0

Fp2 F8 T4 T6 O2 F4 C4 P4 Fz Cz Pz Fp1 F7 T3 T5 O1 F3 C3 P3 T2 T1 2

4

6

Time (s)

8

10

345 uV

0

2

4

6

8

10

Time (s)

Figuur 2.3: De gebruikte artefact segmenten in de simulatiestudie, beide geëxtraheerd met CCA.

De matrix met musculaire artefacten is van rang 9, terwijl de oculaire artefacten een matrix met een rang 1 structuur vormen. Dit is rechtstreeks het gevolg van het aantal bronnen waaruit ze opgebouwd zijn. Om een idee te krijgen wat de structuur van de resulterende tensor is, wordt in Figuur 2.4 de MLSVD van de getensoriseerde versie van de matrices weergegeven. De tensor is hierbij gevormd zoals beschreven in de vorige sectie. In de ruimtelijke component zijn er opnieuw respectievelijk 9 en 1 significante singuliere waarde(n). Bij de oculaire artefacten (Figuur 2.4b) observeren we dat deze voldoet aan een (L,L,1) structuur. De MLSVD van de musculaire artefacten vertoont minder structuur. Gedurende de simulatiestudie wordt het effect van het toevoegen van de volgende types ruis beschouwd: 19

2. Simulatiestudie

a: Musculaire artefacten

b: Oculaire artefacten

1200

3000

mode−1 − Spatial mode−2 − Frequency mode−3 − Temporal

800

600

400

200

0


2500

Singular value

Singular value

1000

2000

1500

1000

500

2

4

6

8

10

12

14

16

18

0

20

2

4

6

8

#

10

12

14

16

18

20

#

Figuur 2.4: MLSVD van de getensoriseerde artefacten.

1. zuivere musculaire artefacten (MA) (Figuur 2.3a); 2. zuivere oculaire artefacten (OA) (Figuur 2.3b); 3. Gaussisch verdeelde random ruis (G); 4. Uniform verdeelde random ruis met gemiddelde 0 (U). Ook worden er combinaties van deze types ruis gebruikt. De geselecteerde ruismatrix N wordt steeds eerst genormaliseerd tot dezelfde Frobenius norm als de gegenereerde simulatiematrix S, N =N⋅

∣∣S∣∣F , ∣∣N ∣∣F

alvorens de geobserveerde matrix X wordt samengesteld als: X(λ) = S + λ ⋅ N .

Het variëren van de parameter λ verandert de aanwezige hoeveelheid ruis in X. Een belangrijke grootheid om de aanwezige hoeveelheid ruis in het signaal te kwantificeren is de signaal-ruisverhouding (SNR), berekend als: SNR(λ) =

RMS(S) . RMS(λ ⋅ N )

Hierin stelt RMS het kwadratisch gemiddelde voor van de respectievelijke matrix, hetgeen voor een signaalmatrix M met K kanalen die elk bestaan uit L datapunten, gedefinieerd is als: ¿ Á 1 K L−1 2 À RMS(M ) = Á ∑ ∑ (M (k, l)) . K ⋅ L k=1 l=0 20

2.4. Dipoolsimulaties met statisch gesitueerde bronnen

2.4

Dipoolsimulaties met statisch gesitueerde bronnen

Het uitvoeren van een diepgaande studie waarbij dipoolsimulaties met één of meerdere statisch gesitueerde bronnen gebruikt zijn, vormt het onderwerp van dit deel. In Sectie 2.4.1 geven we een analyse voor het geval dat de golfvorm van het genererende signaal een stationaire frequentie heeft, met andere woorden als het een sinusoïdaal signaal is. Eerst wordt er geargumenteerd dat in dit geval een beschrijving met een rang-1 term een goede keuze is voor de ictale activiteit, nadien worden de resultaten van de methode in functie van de hoeveelheid en het type van de ruis besproken (Sectie 2.4.1.1). In Sectie 2.4.1.2 wordt er een vergelijking gemaakt met een andere BSS methode, om in Sectie 2.4.1.3 de resultaten van een nauwkeurige lokalisatie van de epileptische bron te bespreken. Tot nog toe is er een enkele dipool met stationaire frequentie beschouwd. In Sectie 2.4.1.4 wordt er besproken wat het gedrag van de methode is indien er twee of meerdere statische bronnen met epileptische activiteit aanwezig zijn. In Sectie 2.4.2 is de analyse gegeven indien de golfvorm van het genererende signaal niet één stationaire frequentie maar meerdere frequenties bevat. Het standaardvoorbeeld hiervan is het reeds besproken chirp-signaal. Eerst wordt er aangetoond dat in dit geval een beschrijving als een term van lage multilineaire rang gerechtvaardigd is, er wordt bekeken welke rang L optimaal is, en er worden ook enkele kanttekeningen geplaatst bij het gedrag van de methode bij lage SNR en grotere L. In Sectie 2.4.2.1 wordt een methode geanalyseerd voor een systematische bepaling van de rang L van de ictale term in de ontbinding, terwijl we in Sectie 2.4.2.2 de lokalisatie opnieuw nader bekijken. In Sectie 2.4.2.3 wordt ten slotte de invloed van de aanwezigheid van meerdere bronnen op de methode besproken.

2.4.1

Stationaire frequenties en CPD

Het eerste geval dat we behandelen is dat van ictale activiteit op basis van een sinusoïdaal signaal. Hier wordt een segment voor gebruikt dat gebaseerd is op een sinus met een frequentie van 6 Hz, zoals in Figuur 2.2a, maar dan bestaande uit 2000 samples of een tijdsduur van 8s. Zoals reeds aangegeven in Hoofdstuk 1, is uit eerder onderzoek [22, 20] gebleken dat een CPD een goede methode is voor dit type data: zowel het signaal als de potentiaalverdeling van de lokalisatie wordt hiermee goed gemodelleerd. In Figuur 2.5 wordt aangetoond waarom dit zo is. Figuur 2.5a toont de CWT van de sinus waarmee de dipoolsimulatie is uitgevoerd, deze CWT is een rang-1 matrix en kan dus geschreven worden als het uitproduct van twee vectoren. Figuur 2.2 suggereert bovendien dat in een simulatie met een enkele dipool, elk kanaal een veelvoud is van de andere kanalen. Dit wordt bevestigd doordat ook deze matrix rang-1 is. De tensor gevormd door een puntbron met zuivere ictale activiteit op een stationaire frequentie is bijgevolg een rang-1 tensor, zoals zichtbaar is in het bovenste gedeelte van Figuur 2.5b. De MLSVD toont slechts één significante 21

2. Simulatiestudie singuliere waarde in elke mode. Wanneer er ruis aan de simulatie toegevoegd is, gaat de rang-1 structuur van de tensor duidelijk verloren (Figuur 2.5b, onder). In dit geval is er MA ruis toegevoegd aan de data en we observeren eenzelfde structuur als de superpositie van Figuur 2.4a en het bovenste gedeelte van Figuur 2.5b. We besluiten dat een CPD in dit geval een goede methode is aangezien we aangetoond hebben dat de component die we willen extraheren rang-1 is.

b: MLSVD van de zuivere, getensoriseerde ictale activiteit (boven), en van de getensoriseerde signaal-ruis combinatie (SNR: 0.32) (onder)

a: CWT van de genererende sinus

Singular value

Original data (S)

10


2000 1000 0

2

4

6

8

10

12

14

16

# Original + Noise (X)

15

Singular value

Frequency samples

5

3000

20

25

30 200

400

600

800

1000

1200

Time samples

1400

1600

1800

2000

1500 1000 500 0

2

4

6

8

10

12

14

16

18

20

#

Figuur 2.5: Evidentie voor de rang-1 structuur van de CWT van een sinusoïdaal signaal en de hieruit geconstrueerde tensor.

Nu we weten dat in dit geval een beschrijving als rang-1 termen gerechtvaardigd is voor het modelleren van de ictale component, bekijken we een voorbeeld in Figuur 2.6. In Figuur 2.6a is de originele, ruisloze data (S) gegeven. Hieraan is OA ruis toegevoegd om tot een datamatrix X te komen met een SNR van 0.32, deze is weergegeven in Figuur 2.6e. Dit is reeds een zeer lage SNR waarde, de oorspronkelijke ictale activiteit is met het oog al moeilijker te herkennen in X en het vormt, zelfs voor een geoefend persoon, een grote uitdaging om op basis hiervan een lokalisatie van het insult uit te voeren. Hier komen BSS technieken van pas om de ictale activiteit van de artefacten te scheiden. Voor onze methode tensoriseren we X zoals beschreven in Sectie 2.2 met een CWT over 30 frequenties van 1 − 30 Hz tot een tensor X . Nadien voeren we een CPD ontbinding uit. Hiervoor dient het aantal noodzakelijke termen R in de ontbinding (Definitie 1.13) nog bepaald te worden. Een eerste optie is R = 2, dit omdat er twee bronnen in de data X aanwezig zijn. Dit is echter niet de beste keuze vermits de OA tensor niet rang-1 is, maar wel rang-(L, L, 1). Gezien in dit voorbeeld de ruis sterk domineert ten opzichte van de ictale activiteit en bij R = 2 de tensor beschreven wordt als de best mogelijke benadering met twee rang-1 termen, is er geen term met zuivere 22

2.4. Dipoolsimulaties met statisch gesitueerde bronnen ictale activiteit in deze ontbinding. De term waar ictale activiteit in waar te nemen is, is gecontamineerd door de artefacten. We dienen R dus te verhogen, als tweede optie bekijken we R = 5. In dit geval zijn vier termen in de ontbinding gerelateerd aan het knipperen van de ogen, en één aan de ictale activiteit. De vier OA termen hebben elk eenzelfde ruimtelijke component, de maximale Euclidische afstand tussen de vier ruimtelijke vectoren van lengte 21 is van grootteorde 10−6 . Bijgevolg vindt de CPD ontbinding essentieel een rang-(L, L, 1) term zonder deze randvoorwaarde op te leggen. De term gerelateerd aan de ictale activiteit is nu niet langer gecontamineerd door de artefacten. In het geval van een klinische dataset met vermoedelijke ictale activiteit, zal R typisch hoger gekozen worden om zoveel mogelijk bronnen in de data te identificeren. Een MLSVD is in dit geval handig voor het bepalen van het aantal benodigde termen R in een CPD of voor de structuur van een BTD-(L, L, 1) ontbinding. Een CPD ontbinding met R = 5, beschrijft het grootste deel van de variatie in de data. Om dit te kwantificeren berekenen we de relatieve fout, erel , tussen X en de CPD als: erel =

∥X − XCPD,R ∥F . ∥X ∥F

In dit voorbeeld is erel slechts 0.04. In Figuur 2.6 worden de drie componenten van twee termen uit deze CPD getoond. Figuren 2.6b, 2.6c en 2.6d tonen respectievelijk componenten f , e en d uit Definitie 1.13 voor de ictale term uit de CPD. De ruimtelijke component d is weergegeven als een topografische kaart in functie van de elektrodeposities, zodat er ook een ruwe schatting van de lokalisatie kan gebeuren. De ruimtelijke component bevat zowel amplitude informatie (absolute waarde), als fase informatie (het teken). Figuren 2.6f, 2.6g en 2.6h tonen dezelfde componenten maar dan voor de dominante OA term uit de vier resterende OA termen. We zien dat het signaal in Figuur 2.6b sterk op een sinus gelijkt. Uit de frequentie component kan de frequentie van de ictale activiteit bepaald worden als 6 ± 1 Hz. De tijd-frequentie matrix van term r wordt bepaald als er ○ f r en is een rang-1 matrix. Om de overeenkomst van de geëxtraheerde ictale term met het oorspronkelijke signaal te kwantificeren, gebruiken we de Pearson correlatie tussen de twee. In dit geval is deze 1.0, wat een perfecte overeenkomst betekent. Figuur 2.7 toont het resultaat nadat de geselecteerde ictale term uit de CPD terug tot een matrix is getransformeerd door een inverse CWT. Het is duidelijk dat er een grote overeenkomst is met de originele ictale activiteit en dat de interpretatie eenvoudiger is in vergelijking met Figuur 2.6e. In wat volgt bekijken we eerst de invloed van de ruis op het resultaat. Vervolgens vergelijken we de methode met een andere BSS-methode en ten slotte bekijken we de resultaten van een precieze lokalisatie van de bron van de ictale activiteit in functie van de SNR.

23

2. Simulatiestudie

7 uV

1

2

3

Time (s)

4

5

Original simulation

a: Originele data S Fp2 F8 T4 T6 O2 F4 C4 P4 Fz Cz Pz Fp1 F7 T3 T5 O1 F3 C3 P3 T2 T1 0

1

2

Time (s)

4

5

Simulation with noise added

3

6

6

7

7

8

8

e: Data met OA ruis toegevoegd X (SNR: 0.32) 25 uV

0


0

1

2

3

Time (s)

4

5

6

7

8

3.5 3

2

2.5

1.5 1 0.5 0

5

10

20

Frequency (Hz)

15

25

30

1

2

3

4

Time (s)

5

6

7

8

2.5

2

1.5

1

0.5

5

10

15

20

Frequency (Hz)

25

30

T3

T1

F7

T5

T5

F7

FP1

F3

C3

P3

O1

FP1

F3

C3

P3

O1

Fz

Cz

Pz

Pz

Cz

Fz

Fp2

F4

C4

P4

02

Fp2

F4

C4

P4

02

F8

T6

T6

F8

T2

T4

T4

T2

h: Topografische kaart van de ruimtelijke component van een OA term

T3

T1

d: Topografische kaart van b: Temporele component van de ic- c: Frequentie component van de ic- ruimtelijke component van tale term tale term de ictale term 4 3 2 1 0 −1 −2 −3 −4

20 15 10 5 0 −5 −10 −15 0

f: Temporele component van een g: Frequentie component van een OA term OA term

Intensity Intensity

Figuur 2.6: Voorbeeld van de analyse van ruizige ictale activiteit met stationaire frequentie aan de hand van een CPD ontbinding.

24

Signal Signal



7 uV

0

1

2

3

4

5

6

7

8

Time (s)

Figuur 2.7: Gereconstrueerde ictale EEG matrix

2.4.1.1

Invloed van het type ruis

Of de aard van de ruismatrix N al dan niet van belang is voor de resultaten van de CPD van de tensor X is de onderzoeksvraag van deze sectie. Om dit te onderzoeken is een eerste test uitgevoerd in navolging van [20]. Daar is bevonden dat het optimaal aantal termen in de CPD ontbinding meestal twee is. Deze schatting is uitgevoerd op basis van de Core Consistency Diagnostic (Corcondia) [42]. Als startpunt beschouwen we dus een CPD ontbinding met R = 2. De term gerelateerd aan de epileptische activiteit is dan geselecteerd en vergeleken met de eigenschappen van de oorspronkelijke simulatie. Het resultaat hiervan is in Figuur 2.8 weergegeven in functie van de SNR en het type ruis. Zo toont Figuur 2.8a de correlatie tussen de CWT van de genererende sinus en de geselecteerde rang-1, tijd-frequentie matrix geconstrueerd uit de temporele en frequentie component. Figuur 2.8b toont de correlatie tussen de amplitudeverdeling. Deze amplitudeverdeling over de verschillende kanalen is bepaald voor de uit de epileptische rang-1 tensor gereconstrueerde datamatrix en vergeleken met die van de oorspronkelijke simulatie. Vermits er slechts 21 kanalen zijn, is de monstergrootte waarvoor de correlatie bepaald wordt veel beperkter dan bij de tijd-frequentie matrix. Dit betekent dat de correlatie hier groter moet zijn om significant te zijn. Uit Figuur 2.8 besluiten we dat een CPD ontbinding met R = 2 tot betrekkelijk lage SNR waarden de epileptische activiteit uit de data kan halen. Eén van de twee rang-1 termen in de CPD komt voor SNR > 0.6 steeds tot op een correlatie van minstens 0.95 overeen met de oorspronkelijke epileptische activiteit en dit zowel qua vorm en amplitudeverdeling. Voor lagere SNR wordt de ruisbijdrage belangrijker. Wat tot gevolg heeft dat in een beste beschrijving van de tensor als de som van twee rang-1 termen (CPD), de term gerelateerd aan de epileptische activiteit verstoord wordt door de ruis. Dit komt omdat de ruisbijdrage niet perfect te beschrijven is als een enkele rang-1 term. Voor MA en OA ruis is dit reeds aangetoond in 25

2. Simulatiestudie

1

1

0.95

0.95

0.9

0.9

0.85

0.85

Correlation

Correlation

a: Correlatie tussen de tijd-frequentie ma- b: Correlatie gerelateerd aan de amplitudetrices verdeling over de verschillende kanalen

G U MA OA MA+OA

0.8 0.75 0.7 0.65

0.7 0.65

0.6

0.6

0.55

0.55

0.5 0.1

G U MA OA MA+OA

0.8 0.75

0.2

0.3

0.4

0.5

0.6

SNR

0.7

0.8

0.9

1

0.5 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

SNR

Figuur 2.8: Correlaties gerelateerd aan de vorm en amplitudeverdeling van het geëxtraheerde signaal in functie van de SNR (R = 2).

Figuur 2.4. Voor G en U ruis is dit ook het geval, deze ruis matrices zijn steeds van volle rang, wat zich in de MLSVD van de overeenkomstige tensor vertaalt als een tiental significante singuliere waarden in de frequentie component, 21 in de ruimtelijke component en nog een stuk meer in de temporele component. Als we de relatieve fout tussen de tensor en deze CPD ontbinding bekijken, merken we al snel op dat deze veel hoger is dan in het voorbeeld van Figuur 2.6. In het geval van G ruis is erel = 0.47 bij een SNR van 1 en dit stijgt zelfs tot erel = 0.94 bij een SNR van 0.1.

Dit alles zet ertoe aan om een ontbinding te beschouwen met een hoger aantal rang-1 termen. Uit het voorbeeld weten we dat R = 5 een goede keuze is voor OA ruis. Voor de MA ruis, alsook voor de G en U ruis is er geen duidelijke structuur zoals (L, L, 1) waar te nemen. Hierdoor is het niet éénduidig te bepalen wat een goede keuze is voor het aantal termen in de CPD. Om een schatting te maken, beschouwen we de relatieve fout erel van een CPD ontbinding van een zuivere ruistensor. Voor OA ruis is deze 0.04 indien R = 4, MA ruis heeft meer significante singuliere waarden in elke component, dit zorgt ervoor dat de fout hier 0.14 is voor R = 30. Bij G ruis is er nog minder structuur aanwezig, hierdoor is er een nog hogere fout van 0.51 bij R = 30. Het resultaat voor U ruis is vergelijkbaar. Gebaseerd op dit resultaat is ervoor gekozen om de analyse met een hoger aantal rang-1 termen uit te voeren, namelijk R = 5 voor OA ruis, R = 12 voor MA ruis en R = 22 voor G en U ruis. Het resultaat hiervan is weergegeven in Figuur 2.9. Het is duidelijk dat het verhogen van het aantal termen in de ontbinding de overeenkomst tussen de geëxtraheerde epileptische activiteit en de oorspronkelijke simulatie sterk vergroot. De verbetering is het meest opvallende voor OA ruis. Bij een beschrijving met 2 termen zijn de resultaten daar het minst. Bij hogere R is er in dit geval een nagenoeg perfecte beschrijving mogelijk van de epileptische activiteit en 26


a: Correlatie tussen de tijd-frequentie ma- b: Correlatie gerelateerd aan de amplitudetrices verdeling over de verschillende kanalen 1

1

0.98

0.98

0.96

0.96 0.94

0.92

Correlation

Correlation

0.94

G U MA OA

0.9 0.88

0.92

0.88

0.86

0.86

0.84

0.84

0.82

0.82

0.8 0.1

G U MA OA

0.9

0.8 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.1

0.2

0.3

0.4

SNR

0.5

0.6

0.7

0.8

0.9

1

SNR

Figuur 2.9: Correlaties gerelateerd aan de vorm en amplitudeverdeling van het geëxtraheerde signaal in functie van de SNR voor hogere R.

dit voor nagenoeg alle SNR waarden. Dit is zo omdat er een perfecte overeenkomst is tussen de structuur van de tensor en een ontbinding in vijf rang-1 termen. Voor MA, G en U ruis is deze perfecte overeenkomst er niet, maar zorgt het verhogen van R wel voor een aanzienlijk betere overeenkomst. In Figuur 2.10 is erel weergegeven in functie van SNR. Dit bevestigt de vorige verklaring: voor OA ruis is de fout onafhankelijk van de SNR, wat aantoont dat de ruiscomponent goed beschreven is. Voor de andere ruistypes stijgt de fout voor dalende SNR, dit geeft aan dat de ruiscomponent goed, maar niet volledig gemodelleerd wordt.

0.8

G U MA OA

0.7

Relative error

0.6 0.5 0.4 0.3 0.2 0.1 0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

SNR

Figuur 2.10: De relatieve fout erel van de CPD in functie van SNR voor de verschillende situaties weergegeven in Figuur 2.9.

27

2. Simulatiestudie 2.4.1.2

Vergelijking met SOBI, een andere BSS methode

In deze sectie vergelijken we de resultaten bekomen met de CPD methode met wat mogelijk is met een andere BSS methode. Vermits CCA reeds gebruikt is voor MA en OA ruis te extraheren, valt nu de keuze op een algoritme voor Second Order Blind Identification (SOBI) [2, 3, 8]. SOBI maakt, net zoals CCA, gebruik van de temporele karakteristieken van de signalen. Er worden, net als bij CCA, bronnen verondersteld die onderling niet gecorreleerd zijn maar een maximale autocorrelatie hebben. Vertrekkende vanuit deze veronderstelling is SOBI opgebouwd als een ICA algoritme gebaseerd op tweede orde statistieken in plaats van hogere orde statistieken. De veronderstellingen gemaakt betreffende de onderliggende bronnen zijn bij SOBI dus sterk gelijkaardig aan die voor CCA, de BSS methode gebruikt voor het extraheren van het MA segment. Niet geheel onverwacht dus blijkt uit Figuur 2.11a dat SOBI tot op zeer lage SNR waarden de ictale activiteit uit de MA ruis kan scheiden. Op lagere SNR is het steeds zo dat SOBI 10 bronnen in de data terugvindt, waarbij de eerste 9 overeenkomen met de 9 CCA bronnen van het MA segment en de laatste de ictale activiteit is. Betreffende de amplitudeverdeling (Figuur 2.11b) is er bij lage SNR toch een mindere overeenkomst met de oorspronkelijke data. Het is echter duidelijk dat in dit geval SOBI betere resultaten geeft dan de CPD methode. De CPD methode is in deze analyse steeds uitgevoerd met R = 10 voor de verschillende types ruis. Deze keuze is gemaakt om een zo goed mogelijke vergelijking te kunnen maken, vermits de epileptische activiteit steeds uit de eerste 10 SOBI bronnen geselecteerd is.

b: Correlatie met de oorspronkelijke intensiteitsverdeling

a: Correlatie met de genererende sinus 1

1

0.95 0.9

0.95

Correlation

Correlation

0.85 0.8 0.75

CPD − MA SOBI − MA CPD − G SOBI − G CPD − G+MA SOBI − G+MA

0.7 0.65 0.6 0.55 0.5 0.1

0.2

0.3

0.4

0.5

0.6

SNR

0.7

0.8

0.9

0.9

CPD − MA SOBI − MA CPD − G SOBI − G CPD − G+MA SOBI − G+MA

0.85

0.8

1

0.75 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

SNR

Figuur 2.11: Vergelijking van de CPD methode met de resultaten bekomen met SOBI.

In het geval van MA ruis heeft SOBI echter wel een concurrentieel voordeel omwille van de constructie van de ruis. Als er random, Gaussische ruis aan de data toegevoegd wordt, zijn de resultaten met de CPD methode beter dan deze van SOBI. De correlatie op de amplitudeverdeling bij SOBI is min of meer in overeenkomst 28

2.4. Dipoolsimulaties met statisch gesitueerde bronnen met het oorspronkelijk signaal en vergelijkbaar met CPD. De overeenkomst met de genererende sinus is voor alle SNR waarden echter een stuk minder met SOBI. Het met SOBI geëxtraheerd signaal bevat steeds nog aanzienlijk veel ruis terwijl de CPD methode bij lage SNR steeds een vervormd maar ruisvrij signaal oplevert omwille van de rang-1 structuur (Figuur 2.7). De meest realistische vergelijking van de twee methodes wordt gegeven door het derde geval in Figuur 2.11. Dit is een combinatie van random Gaussische ruis en de gelokaliseerde MA ruis (G+MA). In dit geval is er gelokaliseerde ruis aanwezig, maar heeft SOBI niet het voordeel dat er slechts 10 bronnen in de data aanwezig zijn. In dit geval is het resultaat met de CPD methode opnieuw beter. De lokalisatie is vergelijkbaar, maar de vorm van het signaal wordt met CPD nauwkeuriger bepaald om dezelfde reden als bij de Gaussische ruis. We besluiten dat, vermits de CPD methode geen veronderstellingen maakt over de bronnen, behalve dat ze een rang-1 structuur hebben, deze in het algemeen beter in staat is om ictale activiteit met stationaire frequentie uit de data te scheiden. Het is echter ook niet gegarandeerd dat bij klinische insulten de ictale activiteit een stationaire frequentie heeft en met een rang-1 tensor gemodelleerd kan worden. Daar wordt verder in dit hoofdstuk dieper op ingegaan.

2.4.1.3

Lokalisatie van de bron

De ruimtelijke component van de ictale term in de CPD ontbinding laat een lokalisatie van de bron toe (Figuur 2.6d). Een andere mogelijkheid is om de ictale term terug tot een matrix te transformeren (Figuur 2.7). Een specialist kan zo een schatting van de locatie van het epileptisch insult maken op basis van dit resultaat. Het is ook mogelijk om een dipool te fitten aan deze matrix om zo de dipoollocatie te schatten. In dit werk is hiervoor de FieldTrip toolbox [38] gebruikt. Vertrekkende vanuit hetzelfde hoofdmodel en dezelfde elektrodeposities, als eerder gebruikt voor het opstellen van de simulatie, wordt een optimalisatieprobleem opgelost. Zo wordt vanuit een initiële dipoollocatie een optimale dipoollocatie gevonden om de data te verklaren. De methode is gevalideerd door de oorspronkelijke gesimuleerde ictale activiteit te fitten. Het blijkt dat dan steeds de werkelijke dipoollocatie teruggevonden wordt, onafhankelijk van de gekozen initiële dipoollocatie. In Bijlage A wordt de gebruikte methode besproken. Om de resultaten van de dipoollokalisatie te vergelijken is de CPD uitgevoerd bij verschillende SNR waarden en voor verschillende types ruis. Deze analyse gebruikt hetzelfde aantal termen R in de ontbinding als eerder in Figuur 2.9. De ictale term in de CPD ontbinding is steeds geïdentificeerd en teruggetransformeerd. Het fitten van een dipool gebeurt op basis van dit resultaat. De Euclidische afstand tussen de werkelijke locatie, (xdip , ydip , zdip ) = (0.3, −0.3, 0.3), en de locatie gefit op het resultaat, (xfit , yfit , zfit ), is in Figuur 2.12 weergegeven voor verschillende types ruis en bij verschillende SNR waarden. Zoals vermeld in het onderschrift van de figuur is 29

2. Simulatiestudie de afstand of fout gegeven in eenheden van de straal waarop de elektroden op de schedel geplaatst zijn. Een lokalisatiefout van 0.1 relek komt overeen met een afstand van 9.2 mm in SI eenheden. Localization error 1

Localization error (relek)

0.9 0.8 0.7 0.6

MA G UNI OA

0.5 0.4 0.3 0.2 0.1 0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

SNR

Figuur 2.12: De fout in de lokalisatie uitgevoerd met de geëxtraheerde ictale activiteit bij verschillende SNR en voor verschillende types ruis in vergelijking met de werkelijke locatie van de dipool (0.3, −0.3, 0.3). De eenheid van de y-as is relek = 9.2 cm.

Het blijkt uit deze Figuur dat de lokalisatie voor ieder type ruis nauwkeurig uitgevoerd kan worden tot op zeer lage SNR waarden. Er zijn enkele verschillen te bemerken maar deze zijn eerder beperkt. Als we als limiet voor een bruikbare lokalisatie een fout opleggen van maximaal 20 mm ten opzicht van de werkelijke dipoollocatie, dan verkrijgen we de resultaten uit Tabel 2.2 als ondergrenzen op de SNR. Gemiddeld genomen is de lokalisatie nauwkeurig tot op een SNR waarde van 0.3, wat reeds een zeer grote hoeveelheid ruis is. Deze analyse is ook uitgevoerd met een CPD ontbinding die de tensor benadert als een som van twee rang-1 termen. In dat geval groeit de fout sneller, vooral bij MA en OA ruis: de grens is dan gelijk aan 0.7. Bij U en G ruis is het verschil kleiner daar ligt de grens op respectievelijk een SNR van 0.35 en 0.4. We besluiten dat het afstemmen van het aantal termen op de rang van de tensor tot een meer nauwkeurige lokalisatie leidt.

Type ruis MA OA G U

SNRmin 0.40 0.25 0.16 0.28

Tabel 2.2: Minimale SNR waarden voor een bruikbare lokalisatie voor verschillende types ruis.

Het resultaat uit Tabel 2.2 is positief om verschillende redenen. Ten eerste, is het precies lokaliseren van een epileptisch insult in de hersenen aan de hand van 30

2.4. Dipoolsimulaties met statisch gesitueerde bronnen een EEG meting moeilijk vermits de ruimtelijke sampling bij een EEG opname met 21 elektroden eerder beperkt is. Bijgevolg is het lokalisatie probleem ill-posed, een kleine verandering aan de data kan reeds een groot effect hebben op het resultaat. Dat de fouten beperkt blijven tot lage SNR waarden toont dus aan dat CPD een goede methode is om de ictale component uit de data te extraheren. Ten tweede is het positief vermits er een goede lokalisatie mogelijk blijkt tot een SNR van 0.3. Bij een correct opgenomen klinische ictale EEG meting ligt de SNR waarde typisch veel hoger [20]. Vervolgens gaan we in Figuur 2.13 nauwkeuriger in op de precieze fouten in de lokalisatie van de dipool. Deze toont de verschillende locaties voor de verschillende ruistypes samen met de 21 elektrodelocaties, de aanduiding voor de neus en twee oren en de exacte dipoollocatie. De gebruikte SNR waarden zijn dezelfde als deze in de vorige figuren. De hoogste SNR waarde komt steeds overeen met de best gelokaliseerde schatting en de laagste SNR waarde met de slechtst gelokaliseerde schatting. Uit deze Figuur is het duidelijk dat de meeste locaties zich in de omgeving van de exacte locatie bevinden, maar dat de richting waarin de fouten evolueren verschillend is voor ieder ruistype.

Evolution in source localization

F3

C3

1

Fz

Cz

F4

0.8 F7

0.6

LeftP3 Ear T3

0.4

z

0.2

Fp1

Nose C4

Pz

Fp2 P4

T5

F8

T1

Right Ear T4

0 O1

−0.2

T6

O2

T2

−0.4 −0.6

1

−0.8

0.5

−1 1

0 0.5

−0.5

0 −0.5 −1

y

−1

x

Figuur 2.13: De evolutie van de geschatte locatie van de dipool in het hoofd voor verschillende types ruis en bij variërende SNR. De betekenis van de kleuren is gelijkaardig aan Figuur 2.12: MA (rood), G (zwart), U (groen) en OA (blauw). De werkelijke locatie van de dipool is aangegeven in cyaan.

31

2. Simulatiestudie 2.4.1.4

Meerdere statisch gesitueerde bronnen met stationaire frequenties

Indien er niet één maar meerdere bronnen van epileptische activiteit in de simulatie aanwezig zijn, is de doelstelling dat beide door de methode gedetecteerd en gelokaliseerd kunnen worden. In dit deel beperken we ons tot twee aanwezige bronnen van ictale activiteit, maar de resultaten kunnen eenvoudig veralgemeend worden naar meerdere bronnen. Indien beide bronnen dezelfde stationaire frequentie hebben, blijkt dat het scheiden van de twee bronnen met deze methode niet mogelijk is. Er kan dan slechts één bron geëxtraheerd worden met CPD, die min of meer halfweg tussen de locatie van de twee bronnen gesitueerd is. Het verschil in locatie alleen is dus niet voldoende om de bronnen van elkaar te scheiden. Indien de twee bronnen een verschillende frequentie hebben, kunnen ze wel van elkaar gescheiden worden met een CPD ontbinding. Nadien kunnen ze elk apart gelokaliseerd worden zoals al reeds eerder gebeurd is. Om dit te illustreren, bespreken we een voorbeeld. In Figuur 2.14a is de resulterende tijdreeks van een dipoolsimulatie met twee dipolen weergegeven. De eerste dipool is gesitueerd in (−0.3, 0.3, 0.5) met een sinus van 6 Hz als tijdreeks, de tweede dipool is gesitueerd in (0.3, −0.3, 0.3) en heeft een cosinus van 10 Hz als tijdreeks. In dit geval is er dus een frequentie- en een faseverschil. Het signaal in Figuur 2.14a is niet geschikt om de lokalisatie van de twee bronnen op uit te voeren, vermits het een superpositie is van de twee bijdragen. Uit de vorige resultaten is het geweten dat deze twee bijdragen elk een rang-1 structuur hebben, omwille van deze reden identificeren we de twee bronnen als twee verschillende termen in een CPD ontbinding. Figuur 2.14b, boven toont de MLSVD van de getensoriseerde data. We zien dat CPD niet geschikt is om de originele data te modelleren in één term. Een beschrijving met twee rang-1 termen is wel mogelijk, alsook een beschrijving met een enkele BTD-(2, 2, 2) term. De eerste optie is echter te verkiezen, vermits we de twee aparte bronnen willen terugvinden. Figuur 2.14b, onder toont aan dat het toevoegen van ruis de eenvoudige structuur opnieuw verstoort. Hier is er MA ruis met een SNR van 0.32 aan de data toegevoegd. De CPD is uitgevoerd met één term meer dan bij een enkele bron van ictale activiteit: er zijn 13 termen gebruikt bij MA ruis en 23 bij G ruis. In Figuur 2.15a wordt de correlatie tussen de tijd-frequentie matrix en de CWT van het genererende signaal in functie van de SNR waarde getoond. We zien dat de correlatie voor MA ruis typisch groter is als voor G ruis en bemerken ook een verschil op tussen de twee bronnen: de bron met de laagste frequentie heeft steeds een hogere correlatie dan de bron met de hogere frequentie. Dit effect is het grootst bij Gaussische ruis, maar is ook waarneembaar bij MA ruis. De precieze reden hiervoor is onzeker, er is geen grote frequentieafhankelijkheid opgemerkt bij enkelvoudige bronnen, noch is er een afhankelijkheid van de locatie waargenomen. Vermoedelijk is de oorzaak dus te vinden in de aanwezigheid van twee bronnen.

32


a: Tijdreeks

b: MLSVD Original simulation

0

Singular value

13 uV

2

3

4

5

6

7

8

mode−1 matrix − Spatial mode−2 matrix − Frequency mode−3 matrix − Temporal

2000 1000

2

4

6

8

10

12

14

16


1500 1000 500 0

1

Original data (S)

3000

0

Singular value


5

10

Time (s)

15

20

25

30

#

Figuur 2.14: Tijdreeks en MLSVD van een simulatie met twee dipolen. Deze hebben (−0.3, 0.3, 0.5) en (0.3, −0.3, 0.3) als locatie en respectievelijk een sinus van 6 Hz en een cosinus van 10 Hz als tijdreeks.

Figuur 2.15b toont de lokalisatiefout in functie van SNR voor de verschillende gevallen. Deze lokalisatie is steeds uitgevoerd op de gereconstrueerde ictale data, opnieuw door een dipool aan de data te fitten. Bij lage SNR is het signaal met hoge frequentie steeds minder precies gelokaliseerd dan het signaal met lage frequentie. Bij hogere SNR is het resultaat vergelijkbaar voor MA en G ruis en voor de twee verschillende bronnen. Alhoewel de lokalisatiefout nog steeds klein blijft tot lage SNR waarden, zijn de fouten wel groter dan bij een enkelvoudige bron. De grens van 20 mm wordt hier voor de twee bronnen en onafhankelijk van het ruistype overschreden bij een SNR van 0.5. Bijgevolg is vooral het resultaat bij G ruis significant minder goed dan voor een enkele epileptische bron. Figuur 2.16 toont de evolutie van de lokalisatie in functie van SNR. De hoogste SNR waarde komt steeds overeen met de best gesitueerde lokalisatie. Het besluit is dat de methode gebaseerd op een CPD ontbinding ook in staat is om meerdere bronnen van ictale activiteit met stationaire frequentie te lokaliseren, al is het mogelijk dat dit minder precies is.

2.4.2

Veranderlijke frequenties en BTD-(L,L,1)

In deze sectie worden de resultaten besproken die voortgekomen zijn uit de analyse van een dipoolsimulatie met een frequentie-veranderend signaal. Eerst bekijken we opnieuw een enkele bron en houden de locatie ervan op dezelfde positie als in de vorige sectie. Het gebruikte chirp signaal is sterk gelijkaardig aan dat weergegeven in Figuur 2.2b, enkel beschouwen we hier een signaal van 8s en verloopt de frequentiedaling van 12 Hz tot 6 Hz over de volledige duur van het signaal. Het is in klinische insulten niet ongewoon dat deze van frequentie veranderen, zoals reeds eerder aangegeven. 33

2. Simulatiestudie

a: Correlatie

b: Lokalisatiefout Localization error

1

0.6


0.98

Correlation

0.96

SIN − MA COS − MA SIN − G COS − G

0.94

0.92

0.9

0.88

0.86

0.5

SIN − MA COS − MA SIN − G COS − G

0.4

0.3

0.2

0.1

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0.3

0.4

0.5

0.6

SNR

0.7

0.8

0.9

1

SNR

Figuur 2.15: (a): Correlatie tussen de tijd-frequentie matrix en de CWT van het genererende signaal in functie van de SNR. De 10 Hz cosinus is aangeduid met ’COS’, de 6 Hz sinus met ’SIN’. In het geval van MA ruis zijn er 13 termen in de CPD gebruikt, bij G ruis zijn het er 23. (b): De lokalisatiefout van beide dipolen in functie van SNR en type toegevoegde ruis.

Evolution in source localization

F3

C3

1

Fz

Cz

F4

0.8 F7

0.6

LeftP3 Ear T3

0.4

z

0.2

Fp1

Nose C4

Pz

Fp2 P4

T5

F8

T1

Right Ear T4

0 O1

−0.2

T6

O2

T2

−0.4 −0.6

1

−0.8

0.5

−1 1

0 0.5

−0.5

0 −0.5 −1

y

−1

x

Figuur 2.16: De evolutie van de geschatte locaties van de dipool in het hoofd voor verschillende types ruis en bij variërende SNR. De betekenis van de kleuren is: MA (rood) en G (zwart). De werkelijke locaties van de dipolen zijn aangegeven in cyaan.

34


In Figuur 2.17 wordt opnieuw de CWT van het signaal, alsook de MLSVD van de derde orde tensor met en zonder toegevoegde ruis weergegeven. Uit deze Figuren kunnen twee nuttige conclusies getrokken worden. Ten eerste, de CWT van het signaal (Figuur 2.17a) geeft een duidelijk beeld van de frequentieverandering in functie van de tijd. Dit is een nuttige aanpak voor de karakterisering van het insult. Ten tweede, de MLSVD van de ruisvrije tensor S (Figuur 2.17b, boven) toont dat hier de rang-1 structuur behouden blijft in de ruimtelijke component. Dit terwijl de temporele en frequentie component, die samen de CWT opmaken, van lage rang zijn. Zo is aangetoond dat de bijdrage van de ictale activiteit beschreven kan worden als een rang-(L, L, 1) term en dus is een BTD in rang-(Lr ,Lr ,1) termen, zoals gedefinieerd in Definitie 1.14, in dit geval gerechtvaardigd. De onderste helft van de Figuur toont dat bij het toevoegen van ruis, de rang-(L,L,1) structuur verloren gaat.

b: MLSVD van de zuivere, getensoriseerde ictale activiteit (boven), en van de getensoriseerde signaal-ruis combinatie (SNR: 0.32) (onder)

a: CWT van de genererende chirp

Singular value

Original data (S)

10

mode−1 matrix − Spatial mode−2 matrix − Frequency mode−3 matrix − Temporal

2000 1000 0

5

10

15

20

25

30

25

30


15

Singular value

Frequency samples

5

3000

20

25

30 200

400

600

800

1000

1200

Time samples

1400

1600

1800

2000

1500 1000 500 0

5

10

15

20

#

Figuur 2.17: Evidentie voor de (L, L, 1) structuur van de CWT van een chirp signaal en de hieruit geconstrueerde tensor.

Op basis van voorgaande Figuur 2.17, is het aannemelijk dat indien een BTD(L, L, 1) de ictale term uit de tensor kan extraheren, er een verbeterde overeenkomst is wat betreft de vorm van het signaal in vergelijking met een beschrijving als een rang-1 term. In de ruimtelijke component verwachten we een goede overeenkomst tussen CPD en BTD-(L, L, 1). Om dit te controleren voeren we als eerste test een CPD en een BTD-(L, L, 1) uit op deze data. Hiervoor is er MA ruis aan het chirp-signaal toegevoegd om tot een SNR van 1 te komen. De CPD ontbindt de resulterende tensor in twee rang-1 termen, voor de BTD-(L, L, 1) kiezen we, gebaseerd op Figuur 2.17b, voor één rang-(3, 3, 1) term en één rang-1 term. In de term van lage multilineaire rang is de ruimtelijke 35

2. Simulatiestudie component de rang-1 term, zodat de tijdfrequentie-matrix gemodelleerd wordt door een matrix van lage rang. Het resultaat hiervan is in Figuur 2.18b weergegeven. Er valt meteen een groot verschil op tussen de twee gevallen. In de beschrijving met een rang-1 term is de frequentieverandering van 12 Hz tot 6 Hz over het verloop van het signaal, zoals te zien is in Figuur 2.17a, niet meer waar te nemen. In plaats daarvan is gedurende het volledige signaal een band waar te nemen die zich uitstrekt over het volledige frequentiebereik. Een beschrijving in termen van lage multilineaire rang is in staat om de oorspronkelijke tijdfrequentie matrix beter te benaderen. We nemen in Figuur 2.18b de frequentiedaling gedurende het verloop van het signaal duidelijk waar. Dit geeft een indicatie dat een BTD-(L, L, 1) ontbinding een verbeterde karakterisering van het insult toelaat.

b: Tijdfrequentie-matrix van het chirp siga: Tijdfrequentie-matrix van het chirp signaal, geëxtraheerd met BTD-(L, L, 1) (rangnaal, geëxtraheerd met CPD (3, 3, 1) term) 5

Frequency samples

Frequency samples

5

10

15

20

25

10

15

20

25

30

30 200

400

600

800

1000

1200

Time samples

1400

1600

1800

2000

200

400

600

800

1000

1200

1400

1600

1800

2000

Time samples

Figuur 2.18: Structuur van de tijd-frequentie matrix van het chirp signaal, bepaald met respectievelijk CPD en BTD-(L, L, 1). De data heeft een oorspronkelijke SNR van 1.

Om deze hypothese verder te onderzoeken, is er een eerste test uitgevoerd. Deze vergelijkt de correlatie tussen de tijdfrequentie-matrices, alsook de amplitudeverdelingen, in functie van de rang L van de BTD-(L, L, 1). We kiezen voor L = 1 (CPD), en 2 tot en met 5, en vergelijken ook de invloed van de SNR waarde, alsook de invloed van het type ruis. De resultaten van deze test zijn weergegeven in Figuur 2.19. Enkele conclusies en bemerkingen kunnen uit deze Figuur gemaakt worden. Ten eerste, het is duidelijk dat, onafhankelijk van de SNR waarde en het type ruis, er een significant grotere correlatie is tussen de geëxtraheerde en de originele tijdsfrequentiematrix bij L = 2 ten opzichte van een rang-1 term. Ten tweede blijkt het dat verder verhogen van L geen grotere correlatie met de tijdfrequentie-matrix oplevert. Integendeel, bij MA ruis daalt de correlatie voor grotere L, dit effect is bovendien groter bij kleinere SNR. In het geval van Gaussische ruis is dit effect kleiner, maar aanwezig, zeker bij de laagste SNR waarde. Als we de correlatie met de amplitudeverdeling 36


a: MA ruis

b: G ruis Amplitude

Time−frequency

Amplitude

1

1

0.99

0.995

0.99

0.995

0.98

0.99

0.98

0.99

0.97

0.985

0.96

0.98

0.95

0.975

0.94 1

2

3

Rank

4

5

0.97

SNR : 1 SNR : 0.56 SNR : 0.32 1

2

3

Rank

4

5

Correlation

1

Correlation

1

Correlation

Correlation

Time−frequency

0.97

0.985

0.96

0.98

0.95

0.975

0.94 1

2

3

Rank

4

5

0.97

SNR : 1 SNR : 0.56 SNR : 0.32 1

2

3

4

5

Rank

Figuur 2.19: De correlatie in functie van de rang (L) tussen de originele en geëxtraheerde tijdfrequentie-matrix (links) en tussen de originele en geëxtraheerde amplitudeverdeling (rechts). Dit bij verschillende SNR waarden en voor respectievelijk MA ruis (a) en G ruis (b).

bekijken, zien we dat bij de twee hoogste SNR waarden, dit min of meer onafhankelijk is van de rang. Bij een SNR van 0.32 daalt de correlatie steeds voor stijgende L. Deze daling is groter voor MA ruis als voor Gaussische ruis. Hieruit besluiten we dat met een BTD-(L, L, 1) ontbinding de ictale term het best uit de data geëxtraheerd kan worden indien L = 2 gekozen wordt. Tijdens deze test is het ook opgevallen dat de ictale tijdfrequentie-matrices bij een BTD-(L, L, 1) ontbinding vaak extra artefacten bevatten. Dit gebeurt onafhankelijk van het type ruis, maar is meer aanwezig bij de twee laagste SNR waarden en de hogere rang matrices (L > 2). Een voorbeeld hiervan is te zien in Figuur 2.20a. Deze toont de tijdfrequentie-matrix van een ontbinding met L = 5 bij data met MA ruis en een SNR van 0.56. In deze matrix zijn duidelijk enkele artefacten aanwezig die in de oorspronkelijke CWT van het signaal niet aanwezig zijn. Waarschijnlijk wordt dit effect veroorzaakt doordat de grotere vrijheidsgraad in een BTD-(L, L, 1) ontbinding ook gemakkelijker toelaat dat een deel van de ruiscomponent mee gemodelleerd wordt in de ictale term van de ontbinding. Dit verklaart ook mee waarom het prominenter aanwezig is voor grotere L, en respectievelijk lagere SNR. Het effect vertaalt zich ook naar de ICWT van de tijdfrequentie-matrices. Dit is zichtbaar in Figuur 2.20b, die een segment van de ICWT toont. Er zijn duidelijk afwijkingen ten opzichte van de originele golfvorm van het chirp-signaal zichtbaar. Deze afwijkingen zijn minder aanwezig bij een CPD ontbinding, zoals aangetoond wordt in wat volgt. In Figuur 2.21 wordt een segment van de originele tijd-frequentie matrix vergeleken met het overeenkomstig segment uit de BTD-(L, L, 1) en de CPD ontbinding. Op deze schaal observeren we dat er duidelijk vervormingen waar te nemen zijn in de 37

2. Simulatiestudie a: Tijdfrequentie-matrix uit een BTD(L, L, 1) ontbinding met L = 5 b: Segment van de ICWT van (a) 1

0.8 5

0.4 0.2

Signal

Frequency samples

0.6 10

15

0 −0.2

20 −0.4 −0.6

25

−0.8 30

−1 200

400

600

800

1000

1200

Time samples

1400

1600

1800

2000

5.8

5.9

6

6.1

6.2

6.3

6.4

6.5

6.6

6.7

Time (s)

Figuur 2.20: Tijdfrequentie-matrix (a) en een segment van het signaal dat ermee gereconstrueerd is (b). De oorspronkelijke data bestaat uit het chirp-signaal met toegevoegde MA ruis zodat er een SNR van 0.56 bekomen is.

BTD-(L, L, 1). CPD is robuuster tegen dit soort vervormingen op kleine schaal, maar heeft als nadeel dat er geen karakterisering van het frequentieverloop van het insult mogelijk is zoals wel mogelijk is bij BTD-(L, L, 1). In Figuur 2.21d is het effect van deze vervorming zichtbaar op het overeenkomstige segment van het gereconstrueerde signaal.

2.4.2.1

Het bepalen van de rang (Lr ) in de BTD-(L,L,1) ontbinding

Een volgend probleem is hoe we in een BTD-(L, L, 1) ontbinding de rang Lr bepalen, die van term tot term kan veranderen. Bij een CPD is de enige noodzakelijke keuze het aantal rang-1 termen R in de ontbinding. De methode gebruikt uit Tensorlab [45] voor de BTD-(L, L, 1) verwacht een specificatie van de te gebruiken rang van iedere term in de ontbinding. In deze simulatiestudie is het mogelijk om hiervoor een goede schatting te maken aan de hand van de originele data. In het geval van de analyse van klinische EEG data, hebben we deze echter niet beschikbaar en is het net de bedoeling deze uit de data te halen. Hiervoor bekijken we ook een andere methode om een BTD-(L, L, 1) ontbinding uit te voeren, deze berekent de BTD-(L, L, 1) op basis van een CPD en clustert de termen samen [43]. Het is gebaseerd op [15, 19]. De methode ondersteunt eveneens ALS en NLS algoritmen, het verschil is echter dat hier de BTD-(L, L, 1) ontbinding gespecificeerd wordt op basis van het totale aantal rang-1 termen R in de ontbinding en het aantal termen T waarin deze samen geclusterd dienen te worden. Deze clustering gebeurt met behulp van een k-means++ algoritme. We vergelijken in Figuur 2.22 de methode van een geclusterde BTD-(L, L, 1) ontbinding met een CPD ontbinding. Het eerste wat uit Figuur 2.22 af te leiden is, is dat ze in overeenstemming is 38


a: Segment van de originele tijdfrequentie- b: Segment van de BTD-(L, L, 1) met L = 4 matrix tijdfrequentie-matrix 5

Frequency samples

Frequency samples

5

10

15

20

25

10

15

20

25

30

30 1420

1440

1460

1480

1500

1520

1540

1560

1580

1600

1420

1440

1460

Time samples

1480

1500

1520

1540

1560

1580

1600

Time samples

c: Segment van de CPD tijdfrequentie- d: Segment van de reconstructie uit zowel matrix CPD als BTD-(L, L, 1) 0.8 0.6 5

0.2

Signal

Frequency samples

0.4 10

15

0 −0.2

20 −0.4 25

BTD − L=4 CPD Original

−0.6 −0.8

30 1400 1420 1440 1460 1480 1500 1520 1540 1560 1580 1600

Time samples

5.6

5.7

5.8

5.9

6

6.1

6.2

6.3

6.4

Time (s)

Figuur 2.21: a: Segment van de originele tijdfrequentie-matrix, b: de ictale tijd-frequentie matrix bepaald met BTD-(L, L, 1), en c: de ictale tijd-frequentie matrix bepaald met CPD. In d is het overeenkomstige segment van de ICWT van b en c, alsook het originele signaal zichtbaar. Aan de data is Gaussische ruis toegevoegd tot een SNR van 0.32.

39

2. Simulatiestudie a: Correlatie tussen originele en geëxtraheerde tijd-frequentie matrices in functie van b: Correlatie betreffende de amplitudeverdeSNR ling in functie van SNR 1

22

2

2 2

2 2 2 2

2 2

1

2 2 2 2

0.95

0.95

Correlation

Correlation

1 0.9

0.85

1

0.9

0.85

0.8

0.8

MA − BTD−(L,L,1)c


MA − CPD (2 terms) 0.75

0.1

0.2

0.3

0.4

0.5

0.6

SNR

0.7

0.8

0.9

1

MA − CPD (2 terms) 0.75 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

SNR

Figuur 2.22: Vergelijking van het gedrag van een geclusterde BTD-(L, L, 1) ontbinding, aangeduid met het subscript c, met dat van een CPD ontbinding. Er is steeds MA ruis aan de data toegevoegd om tot de gewenste SNR waarde te komen. De CPD ontbinding is steeds uitgevoerd in twee termen, waaruit de ictale term geïdentificeerd is. De geclusterde BTD-(L, L, 1) methode berekent steeds 10 rang-1 termen in de ontbinding en clustert deze samen tot 3 termen. Enkel voor de twee laagste SNR waarden van 0.1 en 0.16 was uit dit resultaat de ictale term niet éénduidig te identificeren, hiervoor is R = 5 en T = 4 gekozen. De rang van de uit de BTD-(L, L, 1) geïdentificeerde ictale term is in Figuur 2.22a aangeduid, voor Figuur 2.22b is deze hetzelfde.

met onze bevinding uit Figuur 2.19. Het blijkt ook hier namelijk dat een rang L = 2 meestal optimaal is als ictale term. Dit onafhankelijk van het vorige resultaat en enkel op basis van een clustering methode. Enkel bij de twee laagste SNR waarden is er een rang-1 term gebruikt voor de ictale term, maar hier zijn de parameters voor de geclusterde BTD-(L, L, 1) dan ook aangepast. Uit Figuur 2.22a besluiten we ook dat de correlatie voor de tijd-frequentie matrix steeds hoger is bij een BTD-(L, L, 1) ontbinding, en dat, wanneer L = 1 wordt, deze vergelijkbaar is met het CPD resultaat. We bemerken dat er ook hier vervormingen op kleine schaal waar te nemen zijn, zeker voor lagere SNR waarden. In de volgende sectie onderzoeken we wat hiervan het effect op de lokalisatie is. Uit de analyse van de amplitudeverdelingen in Figuur 2.22b blijkt dat er voor een SNR groter dan 0.4 geen verschil is tussen BTD-(L, L, 1) en CPD. Merk op dat de CPD hier slechts uitgevoerd is in twee termen, terwijl we in Sectie 2.4.1.1 aangetoond hebben dat er een verbeterd resultaat is wanneer R groter gekozen wordt. Daar is het echter zo dat de ictale term rang-1 is, vermits dat in dit geval niet meer is, bemerken we bij het verhogen van R hier hetzelfde effect op als eerder bij de oculaire artefacten. Meerdere termen in de CPD zijn dan te wijten aan de ictale activiteit, deze hebben allen een gelijkaardige ruimtelijke component. Al loopt de maximale afstand tussen de ruimtelijke componenten hier op 40

2.4. Dipoolsimulaties met statisch gesitueerde bronnen tot 0.05, wat groter is dan eerder. De BTD-(L, L, 1) ontbinding gebruikt meer dan twee termen omdat hier de ictale activiteit wel in een enkele term beschreven kan worden. We besluiten dat de randvoorwaarde van een (L, L, 1) structuur in de ontbinding hier voordelig is omdat een enkele ruimtelijke component dan gegarandeerd is. 2.4.2.2

Lokalisatie van de bron (reprise)

Gelijkaardig aan Sectie 2.4.1.3 fitten we een dipool aan de gereconstrueerde ictale term die we uit de data gehaald hebben. We gebruiken hier echter wel de chirp simulatie en onderzoeken de verschillen tussen CPD en BTD-(L, L, 1). Hiervoor voegen we twee types ruis (MA en G) aan de data toe en bekijken de lokalisatie voor SNR waarden tussen 0.25 en 1. De resultaten zijn weergegeven in Figuur 2.23. We merken een verschil op in de manier waarop de analyse is uitgevoerd voor MA en G ruis. In het geval van MA ruis zijn er twee termen in de CPD ontbinding gekozen en is R = 10 en T = 3 gekozen voor de geclusterde BTD-(L, L, 1) methode. Dit levert voor de BTD-(L, L, 1) ontbinding steeds L = 2 op voor de ictale term. In het geval van G ruis is de analyse anders uitgevoerd. Er is slechts 1 term in de CPD ontbinding uitgerekend. Dit omdat wanneer er twee termen bepaald worden, de ictale activiteit opnieuw opgesplitst wordt over de termen. Indien we bij BTD-(L, L, 1) voor de clustering methode opteren en R, T kiezen zoals eerst, wordt de rang van de ictale term zodanig hoog dat het resultaat er minder door wordt. Verlagen van deze parameters heeft weinig effect, daarom is hier de methode gebruikt waar de rang van iedere term in de BTD-(L, L, 1) ontbinding gespecificeerd wordt. Er is steeds 1 term bepaald met L = 2. Deze verschillen zijn te wijten aan de eerder besproken variaties in de aard van de ruis. Localization error 0.5

MA − CPD MA − BTD G − CPD G − BTD


0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

SNR

Figuur 2.23: De fout in de lokalisatie uitgevoerd met de geëxtraheerde ictale activiteit bij verschillende SNR, voor twee types ruis en twee tensorontbinding methodes in vergelijking met de werkelijke locatie van de dipool (0.3, −0.3, 0.3). De eenheid van de y-as is relek = 9.2 cm.

Uit Figuur 2.23 blijkt dat de hypothese dat de lokalisatie in dit geval onafhankelijk is van het feit of er een CPD dan wel een BTD-(L, L, 1) ontbinding gebruikt wordt 41

2. Simulatiestudie te rechtvaardigen is. Deze hypothese is gemaakt omdat de ruimtelijke component van de ictale term rang-1 is (Figuur 2.17b). Ondanks dat de fout in de lokalisatie vergelijkbaar is, is ze voor G ruis en SNR > 0.4 toch systematisch kleiner indien een BTD-(L, L, 1) ontbinding gebruikt is. Het verschil in de fout is hier ongeveer 1 tot 2 mm. Voor MA ruis is de lokalisatiefout (meestal) kleiner voor SNR > 0.3 in het geval van een BTD-(L, L, 1) ontbinding. Hier is de trend echter wel minder systematisch. We besluiten dat het aanpassen van de methode aan de data met een BTD-(L, L, 1) ontbinding een verbeterde lokalisatie tot gevolg heeft en dit ondanks dat de ruimtelijke structuur van de tensor rang-1 is. Voor de laagste SNR waarden is een CPD ontbinding preciezer, vermits deze robuuster is tegen het modelleren van ruis in de ictale term dan een BTD-(L, L, 1). Verder valt op wanneer we Figuur 2.23 vergelijken met Figuur 2.12 de lokalisatiefouten bij MA ruis kleiner zijn bij vergelijkbare SNR en dit zowel voor CPD als BTD-(L, L, 1). Het enige verschil is het gebruikte signaal. Een verklaring is niet duidelijk, maar het blijkt dat wanneer het insult verschillend is van een stationaire sinus, het tot op lagere SNR waarden nauwkeurig gelokaliseerd kan worden. Deze test is herhaald om dit resultaat te valideren, wat een vergelijkbaar resultaat geeft. We kunnen dus analoog aan Tabel 2.2 grenzen opstellen waarvoor de fout kleiner is dan 20 mm. Dit is weergegeven in Tabel 2.3.

Type ruis MA G

SNRmin (CPD) 0.27 0.27

SNRmin (BTD-(L,L,1)) 0.32 0.32

Tabel 2.3: Minimale SNR waarden voor een bruikbare lokalisatie voor verschillende types ruis bij een signaal met veranderlijke frequenties.

2.4.2.3

Meerdere statisch gesitueerde bronnen met veranderlijke frequenties

Indien we opnieuw meerdere bronnen beschouwen maar ditmaal beide met veranderlijke frequentie is het, gelijkaardig aan Sectie 2.4.1.4 opnieuw het doel om beide bronnen van elkaar en de ruis te scheiden. Bij dit probleem blijkt het in het algemeen moeilijker om de ictale termen met BTD-(L, L, 1) uit de data te scheiden. Om dit te illustreren is een simulatie uitgevoerd met twee bronnen op verschillende locaties en verschillende tijdreeksen: de eerste bron is een chirp signaal met een frequentieverloop van 12 tot 8 Hz gedurende 8s, de tweede heeft een frequentieverloop van 6 tot 4 Hz in dezelfde tijdspanne. In Figuur 2.24 zijn de correlaties tussen de tijd-frequentie matrices weergegeven in functie van SNR en dit voor OA en G ruis. Het geval van OA ruis heeft een duidelijke (L, L, 1)-structuur, zowel de twee ictale bijdragen als de ruis en bijgevolg geeft een BTD-(L, L, 1) een goed resultaat voor 42

2.5. Dipoolsimulaties met dynamisch gesitueerde bronnen

1

0.95

Correlation

0.9

0.85

0.8

128 − OA 64 − OA 128 − G 64 − G

0.75

0.7

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

SNR

Figuur 2.24: Correlatie tussen de tijd-frequentie matrices in functie van SNR voor OA en G ruis. De bron met hoge frequenties is aangeduid met “128”, de bron met lage frequenties met “64”. In het geval van oculaire artefacten is een BTD-(L, L, 1) uitgevoerd in drie bloktermen met L = 2, 2 en 4. De termen van rang-2 komen overeen met de twee ictale bronnen, de term van rang-4 met de OA ruis. In het geval van G ruis is er een CPD ontbinding uitgevoerd waarbij termen te wijten aan eenzelfde bron handmatig samengevoegd zijn.

alle SNR waarden. De lokalisatie op basis van dit resultaat gebeurt ook nauwkeurig. Indien er G ruis aanwezig is in de data, is het probleem moeilijker. De ruis vertoont geen structuur en bij een modellering met BTD-(L, L, 1) is het dus niet mogelijk om een structuur uit te buiten. Het blijkt zelfs niet mogelijk om een goede keuze voor L voor de verschillende termen te vinden. Hierdoor is een CPD gebruikt als oplossing, in de ontbinding zijn de term(en) gerelateerd aan de twee bronnen geïdentificeerd. Voor de “64” bron zijn er steeds twee termen in de CPD teruggevonden, deze hebben een sterk gelijkaardige ruimtelijke component. Bijgevolg zijn de temporele en frequentie component van deze twee termen samengevoegd tot een tijd-frequentie matrix van rang-2. Dit verhoogt steeds de correlatie met de oorspronkelijke tijd-frequentiematrix: voor de rang-2 matrices is deze steeds groter dan 0.8 (zie Figuur 2.24), terwijl deze voor de twee afzonderlijke rang-1 termen meestal rond 0.7 ligt en steeds lager is dan de BTD-(L, L, 1) term. Voor de “128” bron is er steeds slechts één term in de CPD teruggevonden onafhankelijk van R, deze is dus steeds rang-1.

2.5

Dipoolsimulaties met dynamisch gesitueerde bronnen

In deze Sectie bootsen we het effect na van een aanval die van locatie verandert gedurende het insult. Dit wordt gedaan door korte segmenten van een FieldTrip simulatie, ieder met een andere locatie, samen te voegen tot een langere tijdreeks. De veranderende locatie is zo gekozen dat de dipool langs een rechte beweegt. Deze aanpak is gevolgd omdat de toolbox geen veranderende locaties ondersteunt. We bekijken twee verschillende voorbeelden met een bewegende dipool. Ten eerste een chirp signaal van 8s met een frequentieverloop van 12 tot 8 Hz, waarbij de locatie 43

2. Simulatiestudie iedere 2s verandert langs een rechte. De startpositie is (−0.2, 0.1, 0.5), de eindpositie (0.2, 0.2, 05), bijgevolg legt de dipool een afstand af van ongeveer 0.4 relek . Het resultaat van de samengestelde simulatie is weergegeven in Figuur 2.25a. De veranderingen in amplitude te wijten aan de veranderende locatie zijn subtiel en best waar te nemen in de kanalen met weinig activiteit, zijnde T4, C4, Cz, T3 en C3. We komen hier in het tweede voorbeeld uitgebreider op terug. De tijd en frequentie component van de getensoriseerde data zijn opnieuw van lage rang omwille van de veranderlijke frequenties, de ruimtelijke component is in dit geval echter ook van lage rang. Vermits het doel is om steeds een enkele locatie te schatten, voegen we G ruis aan de dataset toe en segmenteren we de 8s in 7 segmenten van 2s die elk 1s met elkaar overlappen. Op ieder van deze 7 segmenten is een BTD-(L, L, 1) uitgevoerd en er is steeds een term van L = 2 met de epileptische activiteit geassocieerd. De lokalisatie is uitgevoerd zoals eerder en in Figuur 2.25b is de loodrechte afstand van het resultaat voor elk segment tot het pad van de bewegende bron weergegeven voor drie verschillende SNR waarden. In Figuur 2.26 is het resultaat weergegeven in een 3D plot. De vier oorspronkelijke locaties van de bron liggen op een rechte en zijn in cyaan aangeduid, de locaties van het resultaat zijn voor de verschillende SNR waarden aangeduid in dezelfde kleur als Figuur 2.25b.

b: Loodrechte afstand tot het pad van de bron

a: Tijdreeks van de bewegende bron Original simulation

0.2

7 uV

0

SNR: 1.7 SNR: 1 SNR: 0.56

0.18 0.16 0.14

d (relek)


0.12 0.1 0.08 0.06 0.04

1

2

3

4

Time (s)

5

6

7

8

0.02

1

2

3

4

5

6

7

Segment

Figuur 2.25: Gesimuleerde tijdreeks van een bron van epileptische activiteit die van frequentie en locatie verandert, en de loodrechte afstand van het resultaat van de lokalisatie uit een gepartitioneerde BTD-(L, L, 1) tot het werkelijke pad van de bron.

Uit deze Figuren besluiten we dat door de ruizige data te segmenteren en op ieder segment een BTD-(L, L, 1) uit te voeren, de evolutie van de locatie van de bron gedurende het insult geschat kan worden. Deze schatting verloopt nauwkeuriger bij hogere SNR waarden: bij een SNR van 0.56 is de trend moeilijker waar te nemen, terwijl dit bij 1 en 1.78 wel mogelijk is. In het tweede voorbeeld bekijken we het geval van een dipool met stationaire frequentie en veranderlijke locatie. Hiertoe wordt er opnieuw een samengestelde 44

2.5. Dipoolsimulaties met dynamisch gesitueerde bronnen

Cz

1 0.8

C3

Pz

0.6 P4

Fz

C4 P3

F3

F4

0.4

z

0.2 0

O2 T6 O1

Left Ear T3

RightT5 Ear T4

F8

Nose F7 Fp2 Fp1

−0.2 T1

−0.4

T2

−0.6 −0.8 −1 −1

−0.5

0

0.5

−1 1

1

0

y x

Figuur 2.26: De evolutie van de locatie van de bron in de oorspronkelijke dataset, aangegeven in cyaan, en de resultaten bij verschillende SNR: 1.78 (rood), 1 (blauw) en 0.56 (groen).

simulatie geconstrueerd: een sinusoïdaal signaal verandert iedere periode, of om de 0.1s, lichtjes van locatie om zo over een tijdspanne van 8s, of in 80 stappen langs een rechte te bewegen van (−0.4, −0.4, 0.2) tot (0.4, 0.4, 0.5). De dipool legt bijgevolg een afstand af van 1.17 relek . Er is nu eenzelfde strategie mogelijk als in het eerste voorbeeld door de data te segmenteren, uit ieder segment de ictale activiteit te extraheren (in dit geval met CPD) en elk segment te lokaliseren. Uit Figuur 2.27a die de MLSVD van de bewegende dipool toont, blijkt dat er ook een andere aanpak mogelijk is. De data kan, in goede benadering, beschreven worden met een BTD-(L, L, 1) term maar dan wel met een frequentiecomponent van rang-1, de rang van de spatiale en temporele component is 2. Deze methode benut de kracht van de BTD-(L, L, 1) ontbinding ten volle en op een andere manier dan voordien. Door de spatiale component van lage rang te kiezen, is het mogelijk de veranderende locatie te modelleren in een enkele term. Alvorens de behaalde resultaten te bespreken, komen we even terug op het probleem dat reeds in Figuur 2.25a opgevallen is. Het is namelijk zo dat het effect van de veranderlijke locatie van de bron de grootste invloed heeft op de kanalen met de laagste activiteit, analoog aan het eerste voorbeeld. Om de veranderingen te maximaliseren is er een betrekkelijk grote verplaatsing gebruikt. Figuur 2.27b geeft enkele statistieken weer in verband met de grootte van de relatieve verandering in activiteit gedurende de verplaatsing (zwart) en de amplitude (wit) in de verschillende kanalen. De relatieve verandering in activiteit is berekend als het verschil tussen de maximale en minimale amplitude in een kanaal gedeeld door de maximale amplitude. De verhoudingen in amplitude tussen de kanalen is bepaald door de maximale amplitude van elk kanaal te delen door het maximum van deze maxima. 45

2. Simulatiestudie

b: Statistieken van het bewegende sinusoïdaal signaal

a: Ruisvrije MLSVD Original data (S) mode−1 matrix − Spatial mode−2 matrix − Frequency mode−3 matrix − Temporal

2000

Singular value

Relative activity change / Relative amplitude (%)

2500

1500

1000

500

0

1

2

3

4

5

6

7

8

100 90 80 70 60 50 40 30 20 10 0


#

Figuur 2.27: Links: MLSVD van de ruisvrije getensoriseerde data van het signaal met stationaire frequentie en veranderlijke locatie. Deze suggereert een BTD-(L, L, 1) met de frequentiecomponent als rang-1. Rechts: statistieken van het signaal: zwart duidt de relatieve verandering in activiteit in een kanaal aan, wit duidt de relatieve maximale amplitude van het kanaal aan in vergelijking met de maximale amplitude over alle kanalen, de pijl duidt de richting van de verandering in het kanaal aan.

We zien dat de vijf eerder vernoemde kanalen hier ook de grootste relatieve verandering in activiteit ondergaan, alsook dat deze elk slechts een beperkte amplitude hebben. Dit zorgt voor een dubbel effect. Ten eerste, is het absoluut effect van de positieverandering klein, waardoor het moeilijker te bepalen is. Ten tweede, wanneer er ruis aan de data toegevoegd wordt, gebeurt dit uniform, maar het effect op de kanalen met kleinste amplitude is relatief groter. Als de gehele datamatrix bijvoorbeeld een SNR heeft van 10.0, dan is de gemiddelde SNR over de 5 voornoemde kanalen nog slechts 0.5. Het blijkt dat het extraheren mogelijk is tot op zekere SNR, indien er teveel ruis aanwezig is gaat het subtiele effect van de veranderende lokalisatie teniet. De ictale activiteit wordt nog geëxtraheerd, maar de lokalisatie gaat steeds meer naar een gemiddelde positie en de trend in de data wordt minder gevolgd. In Figuur 2.28 zijn de resultaten van de methode gepresenteerd. Er zijn vier verschillende SNR waarden gebruikt die steeds als gemiddelde over alle kanalen is gegeven (TOT) alsook voor de vijf kanalen met laagste activiteit (LA). Er is G ruis aan de data toegevoegd en na de BTD-(L, L, 1) is de lokalisatie op 50, evenredig verspreide posities in het signaal, uitgevoerd. Deze veranderende lokalisatie is uit het resultaat van BTD-(L, L, 1) bepaald door de matrix van lage rang te construeren uit de ruimtelijke en temporele component en op verschillende tijdstippen een lokaal gemiddelde van de ruimtelijke component uit deze matrix te bepalen. Figuur 2.28a toont de loodrechte afstand tot het pad en toont dat deze klein blijft voor alle beschouwde SNR waarden. Om 46

2.6. Tensoren geconstrueerd uit Hankel matrices en BTD-(L, L, 1) te illustreren dat de posities niet steeds verspreid blijven langs het pad, wordt in Figuur 2.28b de afstand getoond tussen de werkelijke startpositie en de geschatte startpositie, hetzelfde voor de eindpositie. De observatie is dat de geschatte locatie voor de beschouwde SNR waarden dicht bij het pad blijven, maar steeds meer gecentreerd geraken en de trend minder goed volgen. Bij lagere SNR (< 1) wordt een deel van de ruis met vergelijkbare frequentie als het signaal mee gemodelleerd en evolueert het resultaat naar een gemiddelde positie.

a: Loodrechte afstand tot het pad van de bron 0.07

b: Afstand tussen werkelijk begin en geschat begin, idem voor einde 0.45

SNR: TOT: 32, LA: 1.7 SNR: TOT: 18, LA: 0.95 SNR: TOT: 10, LA: 0.53 SNR: TOT: 3, LA: 0.17

0.06

0.4 0.35

d (relek)

d (relek)

0.05

0.04

0.03

0.3


0.25 0.2

0.02 0.15 0.01 0.1 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

start

end

Relative position

Figuur 2.28: Links: de loodrechte afstand tot het pad op 50 verschillende tijdstippen. Rechts: de afstand tussen de startpositie en de eerste geschatte positie, alsook de afstand tussen de eindpositie en de laatste geschatte eindpositie.

We besluiten dat het voorbeeld aangetoond heeft dat BTD-(L, L, 1) in staat is om een veranderlijke lokalisatie uit de data te bepalen, al mag de SNR niet te laag zijn. Dit is voornamelijk omwille van het subtiele effect van de locatie verandering op de simulatie, bovendien is het effect ook nog geconcentreerd op de kanalen met lage activiteit. Indien de effecten in klinische data groter zijn dan voor de simulaties, is het ook mogelijk om de veranderlijke lokalisatie beter te bepalen. Het voordeel van de methode waarbij de spatiale component van lage rang is, is ontegensprekelijk dat met een enkele term uit de tensorontbinding de lokalisatie op elk tijdstip bepaald kan worden. Een methode die de data segmenteert zal meerdere ontbindingen dienen uit te voeren.

2.6

Tensoren geconstrueerd uit Hankel matrices en BTD-(L, L, 1)

Als laatste onderdeel van de simulatiestudie voor deze thesis bekijken we een alternatieve manier om een derde orde tensor te vormen uit EEG data en geven een overzicht van de resultaten hiermee behaald. De reeds gebruikte methode op basis van een 47

2. Simulatiestudie CWT is al meermaals gebruikt in eerder onderzoek, onder andere in [1, 22, 34]. Ze biedt enkele voordelen. Ten eerste is er een duidelijke interpretatie mogelijk van de toegevoegde component als het frequentiespectrum van de data. Ten tweede is een CWT van een sinusoïdaal signaal een rang-1 matrix, en verandert dit in een matrix van lage rang voor een signaal met veranderende frequentie. Ten derde is het mogelijk om aan de hand van een tensorontbinding de epileptische activiteit uit de tensor te scheiden, zoals meermaals aangetoond in de vorige secties en in eerder onderzoek. De methode die in deze sectie wordt gebruikt, vervangt elk kanaal, of rijvector [a1 a2 ⋯ aN ], door een (J × K) Hankel matrix met J + K − 1 = N . Deze matrices worden in een derde orde tensor X ∈ KI×J×K opgeslagen als: (X )ijk = (X)i,j+k−1 ,

met 1 ≤ i ≤ I, 1 ≤ j ≤ J, 1 ≤ k ≤ K. Hierin is I het aantal kanalen. Concreet wordt de vector dus tot de volgende matrix omgevormd: ⎡ a1 ⎢ ⎢ a2 ⎢ ⎢ ⎢ a3 ⎢ ⎢ ⋮ ⎢ ⎢a ⎣ J

a2 a3 a3 ⋯ ⋯ aK ⋰ ⋰ aJ+1 ⋯

⋯ aK ⎤⎥ aK aK+1 ⎥⎥ ⎥ aK+1 aK+2 ⎥ . ⎥ ⋰ ⋰ ⎥⎥ aN −1 aN ⎥⎦

Indien N oneven is, kan er een vierkante matrix bekomen worden met J = K = N2+1 . Na de tensorontbinding kan het signaal gereconstrueerd worden door uit te middelen langs de anti-diagonaal van de verkregen Hankel matrix. In [16] wordt aangetoond dat de Hankel matrices van lage rang zijn, dus ook dat de verkregen tensor een (L, L, 1) structuur heeft en dat de BTD-(L, L, 1) ontbinding uniek bepaald is indien de bronsignalen combinaties van exponentiëlen zijn of exponentiële polynomen. Via de gelijkheid van Euler is het dan duidelijk dat de uniciteit geldig is voor sinusoïdale signalen. Een opmerking hierbij is dat indien N groot wordt en dus J = K ≈ N2 ook groot wordt, de resulterende tensor X van grootteorde 21 × N2 × N2 is. Dit is potentieel veel groter dan de tensoren gebruikt met de CWT methode en is computationeel intensiever. Om dit probleem op te lossen, volgen we twee methodes. Ten eerste, beschouwen we kortere segmenten dan voor de CWT methode, indien N tweemaal kleiner is, zal de resulterende tensor viermaal minder elementen bevatten. Een tweede methode comprimeert X , dit gebeurt op eenzelfde manier als uitgelijnd in [16]. Hiervoor wordt de multilineaire rang van X afgeschat aan de hand van ˜ 2 en R ˜ 3 een bovengrens zijn op respectievelijk de mode-2 en de MLSVD. Laat R ˜ ˜ mode-3 rang van X . Beschouw dan U 2 ∈ KJ×R2 en U 3 ∈ KK×R3 die gevormd worden ˜ 2 mode-2 singuliere vectoren en de eerste R ˜ 3 mode-3 door respectievelijk de eerste R ˜ 2 ×R ˜3 I×R singuliere vectoren van X . Dan kunnen we de tensor T ∈ K van gereduceerde 48

2.6. Tensoren geconstrueerd uit Hankel matrices en BTD-(L, L, 1) grootte beschouwen, T = X ×2 U T2 ×3 U T3 .

Door BTD-(L, L, 1) uit te voeren op T wordt de rekentijd verkleind, maar merk wel op dat het uitrekenen van de MLSVD van een grote tensor ook computationeel intensief zal zijn. Het beperken van de grootte van het segment is nog steeds zinvol. Uit de BTD-(L, L, 1) van T kan de benadering voor de Hankel matrix van volle grootte Gr ∈ KJ×K berekend worden als Gr = U 2 ⋅ H r ⋅ U T3 ,

met H r ∈ KR2 ×R3 de matrix van lage rang uit de BTD-(L, L, 1) van T . Uit Gr kan dan het signaal afgeschat worden door de eerder vernoemde uitmiddeling langs de anti-diagonaal. Een laatste opmerking is dat voor een tensorisatie met Hankel matrices opnieuw elk kanaal door zijn standaardafwijking gedeeld wordt, na de tensorontbinding wordt hier terug voor gecorrigeerd. Het gemiddelde wordt hier niet van de data afgetrokken, vooreerst omdat het vaak dicht bij 0 ligt en ook omdat dit een extra pool kan veroorzaken. ˜

˜

Om deze methode te testen, beschouwen we een signaal met stationaire en een signaal met veranderlijke frequenties. Het eerste is een sinus met een frequentie van 6 Hz, het tweede een chirp met een frequentieverloop van 12 tot 6 Hz. Beide signalen hebben een lengte van 4s of 1000 samples, de laatste sample wordt steeds weggelaten in de tensorisatie om tot een tensor van grootte 21 × 500 × 500 te komen. Figuur 2.29 toont de MLSVD voor de twee ruisvrije tensoren. We merken op dat er in geen geval meer sprake is van een rang-1 term bij stationaire frequentie. De mode-2 en 3 component zijn door de tensorconstructie van dezelfde rang. Bijgevolg heeft de zuivere ictale activiteit wel steeds een (L, L, 1) structuur. Voor de stationaire frequentie is L = 2 en dit neemt toe tot 18 indien er een signaal met veranderlijke frequenties gebruikt wordt. In een eerste test passen we een BTD-(L, L, 1) ontbinding toe zonder de tensoren eerst te comprimeren. In de BTD-(L, L, 1) ontbinding van het sinus signaal wordt een term met L = 2 gebruikt voor de beschrijving van de ictale activiteit, voor het chirp signaal is L = 18 gebruikt. Aan de data is steeds een zekere hoeveelheid G+MA ruis toegevoegd, vermits we eerder opgemerkt hebben dat dit het meest realistische geval is. De resultaten in functie van SNR zijn weergegeven in Figuur 2.30. Figuur 2.30a toont aan dat de Hankel methode de temporele component zeer goed uit de data kan scheiden tot op zeer lage SNR. De resultaten zijn vergelijkbaar of beter dan voor de CWT methode (Figuren 2.9a en 2.22a). We nemen wel waar dat de correlatie lager ligt bij het chirp signaal dan bij de sinus. Dit is voor de CWT methode echter ook zo. Figuur 2.30b toont aan dat ook de correlatie tussen de originele en geëxtraheerde ruimtelijke component voor alle SNR zeer hoog blijft. Bij lage SNR is het resultaat voor het sinus signaal opnieuw iets beter dan voor het chirp signaal. Als we dit resultaat vergelijken met Figuren 2.9b en 2.22b valt op dat de correlatie hier significant hoger ligt, vooral voor lage SNR. Het blijkt dus 49

2. Simulatiestudie a: MLSVD van getensoriseerd sinusoïdaal signaal

b: MLSVD van getensoriseerd chirp signaal

Original data (S) 2500

mode−1 − Spatial mode−2 mode−3

2000

Singular value

Singular value

2000

1500

1000

500

0

Original data (S)

2500


1500

1000

500

1

2

3

4

5

#

6

7

8

0

2

4

6

8

10

12

14

16

18

20

#

Figuur 2.29: De MLSVD van tensoren geconstrueerd met de methode gebaseerd op Hankel matrices. Links: het geval van een sinusoïdaal signaal, rechts: een chirp signaal.

dat de Hankel tensorisatie methode samen met BTD-(L, L, 1) er beter in slaagt om ictale activiteit uit ruizige data te scheiden dan de CWT methode in combinatie met CPD/BTD-(L, L, 1). Om verder te illustreren dat de correcte ruimtelijke component uit de ruizige data bepaald is, toont Figuur 2.30c de Euclidische afstand tussen de oorspronkelijke en geëxtraheerde ruimtelijke component. Het is duidelijk dat deze klein blijft voor alle SNR. De verwachting is dus dat de lokalisatie op basis van een dipoolfit van deze resultaten een kleine fout geven. Figuur 2.30d toont de lokalisatiefout in functie van SNR en als we dit vergelijken met Figuren 2.12 en 2.23 valt op dat deze voor alle SNR waarden toch aanzienlijk hoger ligt dan voor de CWT methode. De reden hiervoor is niet duidelijk, vermits de andere resultaten wijzen op een betere bepaling van de epileptische activiteit met deze methode en gezien voor klinische data de lokalisatie enkel op basis van de ruimtelijke component gebeurt zonder dipoolfit, is het ook van ondergeschikt belang. In een tweede test wordt het effect van het comprimeren van de tensor, zoals eerder beschreven, op het resultaat van de ontbinding onderzocht. De resultaten van deze test zijn weergegeven in Figuur 2.31. Voor het sinus signaal wordt de mode-2 en 3 component gecomprimeerd tot verschillende waarden tussen 400 en 2, vertrekkende van een oorspronkelijke grootte van 500. Voor het chirp signaal zijn verschillende waarden tussen 400 en 18 onderzocht. Deze ondergrens is gekozen vermits het chirp signaal in de ontbinding beschreven wordt met een term van L = 18. Figuur 2.31a toont de temporele correlatie in functie van de gecomprimeerde grootte. Voor de sinus is de correlatie onafhankelijk van de gebruikte compressie, en ook grotendeels onafhankelijk van de SNR (zoals in Figuur 2.30a). Voor het chirp signaal is er een afhankelijkheid waarneembaar van de compressiegrootte, indien de tensor te klein is, wordt de correlatie significant lager. Toch is deze nog steeds betrekkelijk groot. Figuur 2.31b toont de ruimtelijke correlatie in functie van de gecomprimeerde grootte. Deze is onafhankelijk van de gecomprimeerde grootte voor zowel de sinus als het 50

2.6. Tensoren geconstrueerd uit Hankel matrices en BTD-(L, L, 1)

a: Temporele correlatie

b: Spatiale correlatie

1

1

0.995 0.995

Spatial correlation

Temporal correlation

0.99 0.985 0.98 0.975 0.97 0.965 0.96

Chirp Sine 0.985

0.98

0.975

Chirp Sine

0.955 0.95

0.99

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.97

0.3

0.4

0.5

SNR

0.35

0.8


0.9

Euclidean distance

0.3 0.25 0.2 0.15 0.1

0.4

0.5

0.6

SNR

0.9

1

0.7

0.8

0.9

1

0.7 0.6 0.5 0.4 0.3 0.2 0.1

Chirp Sine 0.3

0.8

d: Lokalisatiefout

0.4

0

0.7

SNR

c: Euclidische afstand

0.05

0.6

0.7

0.8

0.9

1

0

Chirp Sine 0.3

0.4

0.5

0.6

SNR

Figuur 2.30: Resultaten in functie van SNR: temporele correlatie tussen het oorspronkelijk en geëxtraheerd signaal, ruimtelijke correlatie tussen de originele en geëxtraheerde ruimtelijke component, Euclidische afstand tussen de originele en geëxtraheerde ruimtelijke component en de fout op een lokalisatie uitgevoerd met een dipoolfit.

chirp signaal. De afhankelijkheid in functie van SNR uit Figuur 2.30b valt ook hier opnieuw op. We besluiten dat het verkleinen van de grootte van de tensor met de beschreven methode, nog toelaat om de ictale component uit de data te scheiden met behulp van een BTD-(L, L, 1) ontbinding. Dit kan tot op zeer kleine grootte en bij verschillende SNR. Indien het onderliggend signaal een sinus is die met een term van lagere rang beschreven kan worden, zijn de resultaten iets beter, vooral betreffend de temporele correlatie.

51

2. Simulatiestudie

a: Temporele correlatie

b: Spatiale correlatie

1.01

1

1

0.998 0.996

0.98 0.97 0.96

Chirp − SNR: 1 Chirp − SNR: 0.55 Chirp − SNR: 0.32 Sine − SNR: 1 Sine − SNR: 0.55 Sine − SNR: 0.32

0.95 0.94 0.93 0.92

50

100

150

200

250

Compressed size

300

350

400

Spatial correlation


0.99

0.994 0.992


0.99 0.988 0.986 0.984

50

100

150

200

250

300

350

400

Compressed size

Figuur 2.31: Resultaten in functie van gecomprimeerde tensorgrootte: temporele correlatie tussen het oorspronkelijk en geëxtraheerd signaal en ruimtelijke correlatie tussen de originele en geëxtraheerde ruimtelijke component. Steeds voor zowel het sinus, als het chirp signaal bij drie verschillende SNR waarden.

52

Hoofdstuk 3

Validatie op klinische data Voor het onderzoek op klinische data is een dataset met 40 EEG opnames van 40 verschillende patiënten beschikbaar. Hieruit zijn er 34 geselecteerd waarvoor de starttijd van het epileptisch insult nauwkeurig bepaald is. Deze 34 opnames vormen de gebruikte klinische dataset en zijn opgedeeld in vier categorieën afhankelijk van de aanwezige hoeveelheid musculaire en andere artefacten (MA). Deze categorieën zijn MA 1, 2, 3 en 4. Hierbij staat 1 voor geen of amper aanwezigheid van artefacten, 2 voor enige aanwezigheid van artefacten zonder dat deze evenwel een grote invloed hebben op de interpretatie door een geoefend oog, 3 voor een redelijke hoeveelheid artefacten die een invloed hebben op de interpretatie van het EEG, en ten slotte 4 voor een sterke aanwezigheid van artefacten die de interpretatie van het EEG betrekkelijk moeilijk maakt. Onder deze 34 resterende EEG opnames is er één met MA 1 gedurende het insult, acht met MA 2, zeventien met MA 3 en acht met MA 4. In totaal zijn 24 van de patiënten vrouwelijk, 10 patiënten zijn mannelijk. Voor 16 patiënten is er zogenaamde subtraction ictal SPECT co-registered to MRI (SISCOM) data beschikbaar. Hierbij staat SPECT voor Single Photon Emission Computerized Tomography. Dit is een medische beeldvormingstechniek waarbij de epileptische zone in de hersenen gedurende de ictale periode nauwkeurig opgemeten kan worden. In tegenstelling tot EEG heeft deze techniek een goede ruimtelijke resolutie, maar een slechtere tijdsresolutie. Deze SISCOM data is, indien beschikbaar, steeds als referentie gebruikt om het resultaat van de tensorontbinding van de EEG data mee te vergelijken. De lokalisatie van de epileptische aanval op basis van een tensorontbinding is in dit hoofdstuk op een andere manier uitgevoerd dan voordien. De reden hiervoor is dat een lokalisatie op basis van een dipoolfit inconsistente resultaten geeft voor de klinische data. Een vermoedelijke verklaring is dat de simulaties in het vorige hoofdstuk een specifiek en vereenvoudigd hoofdmodel gebruiken, zowel voor het opstellen van de simulatie als het fitten van een dipool aan het resultaat van de tensorontbinding. De klinische data voldoet dus niet aan dit model. De gebruikte methode is een eenvoudige weging met de elektrodeposities. Ze is in overeenstemming 53

3. Validatie op klinische data met de topografische kaart van de ruimtelijke component en met de SISCOM data. Indien P de matrix is die per rij een elektrodepositie bevat, en a de ruimtelijke vector gerelateerd aan de ictale activiteit bepaald met een tensorontbinding en reeds terug gecorrigeerd voor de geëxtraheerde standaardafwijking. Dan wordt de geschatte ˆ gegeven door: locatie x a ˆ= x ⋅P. ∥a∥1

Deze methode geeft dan weer inconsistente resultaten voor simulaties, wat de reden is dat er twee verschillende methodes gebruikt zijn. In wat volgt geven we de resultaten van drie uitgewerkte voorbeelden uit de klinische dataset in volgorde van toenemende hoeveelheid ruis en artefacten. Indien beschikbaar is SISCOM data gebruikt ter referentie.

3.1

MA 1 - Patiënt 22

De eerste patiënt heeft een MA van 1, en zoals in Figuur 3.1 te zien is, is er duidelijke ictale activiteit waarneembaar in verschillende kanalen gedurende het tijdsverloop van de opname. De amplitude van het insult neemt toe met de tijd en is het hoogste aan de Fp1, F7 en T3 elektroden. De ictale activiteit is ook duidelijk waarneembaar in de meeste andere kanalen. Voor deze patiënt is er geen SISCOM data beschikbaar om het resultaat van een tensorontbinding mee te vergelijken, maar vermits de amplitudeverdeling duidelijk waarneembaar is, is het ook zo mogelijk om te controleren of het een zinvol resultaat betreft. Strikt genomen heeft een geoefend persoon hier geen BSS technieken nodig om een schatting van de locatie uit te voeren omwille van het gebrek aan ruis. De meest opvallende artefacten zijn de uitwijkingen in Fp1 en Fp2 kanalen te wijten aan het knipperen van de ogen, deze activiteit wordt in mindere mate ook in naburige kanalen opgevangen. Kanaal T5 lijkt wel iets meer ruis te bevatten, alsook een artefact net na 14s. De aanwezige artefacten belemmeren echter de interpretatie van het EEG niet. Hierdoor vormt dit voorbeeld ook een goede eerste validatie test van de methodes uit het vorige hoofdstuk, en laten ze ook een meer nauwkeurige en kwantitatieve analyse van het insult toe in vergelijking met een observatie van het EEG. Concreet vergelijken we de resultaten van vier methodes. De eerste drie construeren de tensor met een CWT over 30 frequenties en de volledige 20s of 5000 samples, de laatste gebruikt de methode met Hankel matrices en beschouwt twee opeenvolgende segmenten van 999 samples. Voor de CWT methode bekijken we ten eerste een ontbinding in rang-1 termen, of CPD, van de tensor. Uit deze CPD ontbinding wordt de ictale term geëxtraheerd en op basis daarvan de positie van het insult geschat. De geschatte positie zal de gemiddelde positie van het insult over de volledige tijdreeks zijn. De tweede methode gebruikt een BTD-(L, L, 1) ontbinding van het volledige segment, waarbij de ruimtelijke component rang-1 is. Dit geeft opnieuw een gemiddelde positie over het insult en zal dus een eventuele 54

3.1. MA 1 - Patiënt 22


204 uV

0

2

4

6

8

10

12

14

16

18

20

Time (s) Figuur 3.1: Eerste 20s van de EEG opname van patiënt 22 met MA 1 na de start van het epileptisch insult.

positieverandering van de bron niet opmerken. Wel is het in dit geval mogelijk om een eventuele frequentieverandering van het insult te detecteren. De derde methode gebruikt ook een BTD-(L, L, 1) ontbinding maar neemt als rang-1 component de frequentie component. Dit laat toe om een eventuele verandering van lokalisatie op te merken, zoals eerder ook uitgevoerd in het vorige hoofdstuk. De vierde methode gebruikt dus een tensorisatie met Hankel matrices en een BTD-(L, L, 1) ontbinding met de ruimtelijke component als rang-1 component. Zoals reeds in vorig hoofdstuk opgemerkt geeft dit aanleiding tot termen van lage rang. Figuur 3.2 toont de MLSVD van zowel de tensor geconstrueerd met CWT als voor deze met Hankel matrices. Bij de CWT methode is de frequentie component van laagste rang, dit suggereert dat een BTD-(L, L, 1) met de frequentiecomponent als rang-1 een goed idee is. De MLSVD van de Hankel tensor toont zoals eerder dat de mode-2 en mode-3 componenten dezelfde structuur bezitten. De vier verschillende methodes zijn elk in staat om de epileptische activiteit uit de data te scheiden. Figuur 3.3a toont de topografische kaart van de ruimtelijke component van de epileptische activiteit bepaald met een BTD-(L, L, 1) ontbinding met de ruimtelijke component rang-1 en L = 2 voor de ictale term. Figuur 3.3b toont de bijhorende tijdfrequentie matrix van deze term in de ontbinding. Deze matrix is dus van rang 2. De topografische kaart laat toe om een eerder ruwe schatting van de lokalisatie van het insult te krijgen. Als we dit vergelijken met de amplitudeverdeling in Figuur 3.1, dan lijkt het resultaat overeen te komen. De tijd-frequentie matrix 55

3. Validatie op klinische data

a: CWT tensor

b: Hankel tensor

3500

1800


3000


1600 1400

Singular value

Singular value

2500

2000

1500

1200 1000 800 600

1000 400 500

0

200

2

4

6

8

10

12

14

16

18

0

20

1

2

3

4

5

6

#

7

8

9

10

11

12

#

Figuur 3.2: MLSVD van respectievelijk het met CWT getensoriseerde segment van 5000 samples en de MLSVD van een segment van 999 samples getensoriseerd met Hankel matrices.

toont aan dat de epileptische aanval in intensiteit toeneemt gedurende het insult alsook dat de frequentie tot op de gebruikte resolutie constant blijft gedurende het insult en gelijk is aan 5 ± 1 Hz. Ook toont dit aan dat deze term in de ontbinding ontdaan is van de aanwezige artefacten, want de aanwezige oculaire artefacten in Figuur 3.1 zijn niet aanwezig in Figuur 3.3b. a: Topografische kaart van de ruimtelijke component van de epi- b: Tijd-frequentie matrix van de epileptische activiteit leptische activiteit

F3

Fz

F4

F8

T1

T2

T3

Cz

C3

T5

P3

O1

Pz

C4

P4

02

T4

Frequency samples

F7

5

Fp2

FP1

10

15

20

T6 25

30 500

1000

1500

2000 2500 3000 Time samples

3500

4000

4500

5000

Figuur 3.3: De topografische kaart en tijdfrequentie matrix van de term gerelateerd aan de epileptische activiteit in de BTD-(L, L, 1) ontbinding van de CWT tensor met de ruimtelijke component rang-1 en L = 2.

De lokalisatie van het epileptisch insult is het hoofddoel van de methode en daarom vergelijken we de geëxtraheerde locatie voor de vier verschillende methodes in Figuur 3.4. Wat opvalt is dat op deze schaal de resultaten van de vier methodes 56

3.1. MA 1 - Patiënt 22 moeilijk te onderscheiden vallen en dus dicht bij elkaar gelegen zijn Dit is een bewijs van de correctheid van elk van de vier methodes.

Cz Fz

1 C3

F3

0.8

C4 F4 Pz

0.6 P3

z

0.4 0.2

Left Ear

Nose Fp1P4

F7

Fp2

T3

F8 Right Ear

T5

T4

0 T1 O1

T6

O2

−0.2

1 T2

−0.4 1

0.5

0 0.5

−0.5

0 −0.5 −1

y

−1

x

Figuur 3.4: Lokalisatie op basis van een CPD term van de CWT tensor (rood), op basis van een BTD-(L, L, 1) ontbinding met de ruimtelijke component rang-1 (zwart). De verscheidene lokalisaties aan de hand van een BTD-(L, L, 1) ontbinding met de frequentie component rang-1 en het gemiddelde van deze posities zijn aangegeven in blauw, de twee lokalisaties op basis van een BTD-(L, L, 1) ontbinding van twee tensoren die geconstrueerd zijn uit Hankel matrices van twee EEG segmenten en het gemiddelde zijn weergegeven in groen.

Om toch een betere vergelijking mogelijk te maken is in Figuur 3.5a een vergroting van dit resultaat weergegeven. Hier zijn de verschillende lokalisaties wel waarneembaar. Vooral het resultaat van de BTD-(L, L, 1) ontbinding met de frequentiecomponent als rang-1 component is interessant. Uit deze ontbinding is uit een enkele term de lokalisatie van het insult om de seconde geschat. Deze posities zijn als blauwe punten weergegeven en hun gemiddelde als een grotere blauwe bol. We observeren dat de positie van de bron gedurende het insult slechts weinig veranderlijk is, de maximale afstand tussen alle deze blauwe punten is slechts 0.08 relek ofwel ongeveer 0.8cm als we opnieuw een straal relek van 9.2cm veronderstellen. Vermits dit zeer klein is, suggereert het dat het insult niet van locatie verandert. De groene punten zijn geschat op basis van de ontbinding van twee Hankel tensoren uit twee opeenvolgende segmenten van 999 samples waarvan het eerste 5s na de start van het insult begint. De tensor is via een MLSVD gecomprimeerd tot grootte 100 in mode-2 en mode-3 en de epileptische activiteit is beschreven door een term met L = 6. We observeren dat de lokalisatie van de twee segmenten verschillend is, maar ook deze punten liggen slechts 0.04 relek van elkaar af. Als de afstanden tussen de vier 57


b: Verandering in de verschillende coördinaten

a: Vergroting van Figuur 3.4

0.6

0.5 0

0.4

X coordinate Y coordinate Z coordinate

Coordinate

−0.01 −0.02

z

−0.03 −0.04 −0.05

0.2

0.1

−0.06 0.6

−0.07 −0.08 0.56

0.3

0.55 0.54

0

0.5 0.52

y

0.5

0.48

0.45

−0.1

x

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Relative position

Figuur 3.5: Links: een uitvergroting van Figuur 3.4, dezelfde kleuren zijn behouden. Voor BTD-(L, L, 1) met de frequentie component rang-1 zijn de 20 schattingen als kleine blauwe punten weergegeven, het gemiddelde als een grote blauwe bol. Voor het resultaat van de Hankel tensor zijn de lokalisaties op de twee segmenten met kleine groene punten weergegeven, het gemiddelde als een grote groene bol. Rechts: de evolutie in de verschillende coördinaten gedurende het verloop van het segment bepaald met een BTD-(L, L, 1) met de frequentie component rang-1.

(gemiddelde) lokalisaties vergeleken worden, dan liggen deze maximaal 0.04 relek van elkaar af. Figuur 3.5b toont het verloop voor de drie coördinaten van de blauwe punten in functie van de positie in het tijdsverloop van het insult. Ook hier valt het op dat de verschillen beperkt zijn en is geen trend waarneembaar. We besluiten dat dit voorbeeld heeft aangetoond dat de lokalisaties op basis van vier verschillende methoden om de tensor te ontbinden in overeenstemming met elkaar zijn. Ook is er aangetoond dat het insult noch van frequentie, noch van locatie verandert.

3.2

MA 2 - Patiënt 11

In het tweede voorbeeld van deze klinische validatie studie is de EEG opname iets meer door ruis verstoord en heeft een MA van 2. De ruizige opname is weergegeven in Figuur 3.6. Op deze schaal is het al moeilijker om de epileptische aanval waar te nemen in de data, een uitvergroting van deze data laat echter nog wel een interpretatie door een geoefend persoon toe. Voor deze patiënt is er ook SISCOM data beschikbaar dus het resultaat van de tensorontbindingen kan hiermee vergeleken worden. Als eerste bekijken we opnieuw het resultaat van een BTD-(L, L, 1) ontbinding 58

3.2. MA 2 - Patiënt 11


323 uV

0

2

4

6

8

10

12

14

16

18

20

Time (s) Figuur 3.6: Eerste 20s van de EEG opname van patiënt 11 met MA 2 na de start van het epileptisch insult.

van de met CWT getensoriseerde data waarbij de ruimtelijke component rang-1 is. De ontbinding is uitgevoerd in 4 termen, ieder met L = 2. De epileptische activiteit in de data is geassocieerd met de tweede term in deze ontbinding. Figuur 3.7a toont de topografische kaart van de ruimtelijke component, Figuur 3.7b de tijdfrequentie matrix van rang 2. Dit resultaat laat een eerste ruwe schatting van de locatie van de aanval toe en laat vooral een zichtbare frequentiestijging in de tijdfrequentie matrix zien gedurende het verloop van het insult. Inspectie van de geëxtraheerde tijdreeks bevestigt dat de frequentie stijgt van ongeveer 5 Hz tussen seconde 4 en 8 tot ongeveer 8 Hz in de laatste 4s van het signaal. Naast deze eerste methode, die in Figuur 3.8b in zwart is weergegeven, gebruiken we nog twee andere methodes om de lokalisatie uit te voeren. Ten eerste, een BTD(L, L, 1) ontbinding van de CWT tensor met de frequentie component rang-1, de ictale term heeft L = 2. Dit laat opnieuw toe om de lokalisatie in functie van de tijd te schatten. De schatting is opnieuw iedere seconde uitgevoerd en zijn weergegeven met blauwe punten in Figuur 3.8b. Het gemiddelde van deze posities is met een grotere blauwe bol weergegeven. Ten tweede bekijken we een BTD-(L, L, 1) ontbinding van de Hankel tensor. Vermits hier meer ruis in de data aanwezig is, is het niet meer zo dat de epileptische activiteit de belangrijkste component is. Hierdoor is een enkel segment van 999 samples gebruikt dat 13s na het begin van de aanval begint. De epileptische term is de derde term in de ontbinding met L = 2 en de geschatte positie is in Figuur 3.8b weergegeven in groen. 59


a: Topografische kaart van de ruimtelijke component van de epileptische activiteit

Fz

F4

F8

T1

T2

T3

Cz

C3

T5

P3

O1

Pz

C4

P4

T4

Frequency samples

F3

5

Fp2

FP1 F7

b: Tijd-frequentie matrix van de epileptische activiteit

10

15

20

T6

02

25

30 500

1000

1500


3500

4000

4500

5000

Figuur 3.7: De topografische kaart en tijdfrequentie matrix van de term gerelateerd aan de epileptische activiteit in de BTD-(L, L, 1) ontbinding van de CWT tensor met de ruimtelijke component rang-1 en L = 2.

a: SISCOM lokalisatie van het insult

b: Lokalisatie op basis van BTD-(L, L, 1) 1

1

0.5

0.5

0

LE

−0.5 −1

RE

−0.5

0

0.5

RE

0

1

−0.5 1

0.5

0

N

−0.5

−1

1

N

0.5

CWT BTD s=1 CWT BTD f=1 (sep)

RE 0

LE

CWT BTD f=1 (mean)

−0.5

1

0.5

0

−0.5

Hankel BTD

−1 −1

Figuur 3.8: Links: lokalisatie van het epileptisch insult op basis van SISCOM data. Rechts: lokalisatie van het insult op basis van BTD-(L, L, 1) met drie verschillende methoden. De posities van de neus, het linker- en rechteroor zijn in de figuren respectievelijk aangeduid met N, LE en RE. De figuren tonen beide respectievelijk het achteraanzicht, het rechterzijaanzicht en het bovenaanzicht van de locatie.

60

3.3. MA 3 - Patiënt 26 Figuur 3.8a toont de drie aanzichten van de lokalisatie met SISCOM. Wanneer deze vergeleken worden met de drie aanzichten van het resultaat van de drie BTD(L, L, 1) methoden is het duidelijk dat er een goede overeenkomst is. Dit geeft een onafhankelijke validatie van de gebruikte methode. In Figuur 3.8b merken we wel iets grotere verschillen op tussen de verschillende methoden onderling dan in het vorige voorbeeld. Het maximale verschil in afstand tussen de drie gemiddelde geschatte posities is hier 0.39 relek , maar vermits de lokalisatie met SISCOM ook geen eenduidige positie weergeeft, zijn ze allen in overeenstemming. De methode waarbij een veranderende lokalisatie geschat kan worden op basis van een enkele BTD-(L, L, 1) term geeft hier een groter verschil in maximale verplaatsing van de locatie gedurende het insult. Deze is namelijk 0.14 relek , maar is nog steeds te klein om door een veranderlijke lokalisering veroorzaakt te zijn.

3.3

MA 3 - Patiënt 26

Het derde en laatste voorbeeld uit deze klinische validatie heeft een iets grotere hoeveelheid ruis en is onderverdeeld in de categorie MA 3. De tijdreeks van de EEG opname is niet weergegeven, maar is kwalitatief vergelijkbaar met de vorige twee voorbeelden, enkel met een iets grotere hoeveelheid ruis. De opname duurt na de detectie van het begin van het insult iets minder dan 20s. Hierdoor is er voor de CWT methode een segment gebruikt van 15s in tegenstelling tot de vorige twee voorbeelden. De lokalisatie op basis van een ontbinding van een tensor geconstrueerd uit Hankel matrices gebruikt opnieuw een segment van 999 samples dat 7s na het begin van het insult aanvangt. De lokalisatie is uitgevoerd op dezelfde drie manieren als voor het vorige voorbeeld. De CWT BTD-(L, L, 1) methode met de ruimtelijke component rang-1 ontbindt de tensor in 6 termen met L = 2, waarbij de eerste gerelateerd is aan de epileptische activiteit. In dit resultaat is geen wezenlijke frequentieverandering waargenomen, de frequentie van het insult is gedurende de volledige tijdsduur 4±1 Hz. De CWT BTD-(L, L, 1) methode met de frequentie component rang-1 is uitgevoerd met 4 termen van L = 2. Hieruit is de derde term geïdentificeerd als degene die gerelateerd is met de epileptische activiteit. De positie in functie van het tijdsverloop is hier opnieuw iedere seconde bepaald op dezelfde manier als voorheen. De laatste methode gebruikt terug een BTD-(L, L, 1) ontbinding van een tensor geconstrueerd uit Hankel matrices. Deze tensor is gecomprimeerd tot een grootte van 100 in mode-2 en mode-3. De ontbinding is uitgevoerd in drie termen met L = 2, hiervan is de tweede geïdentificeerd als degene gerelateerd aan de epileptische activiteit. De resultaten hiervan zijn weergegeven in Figuur 3.9b waarbij de kleuren dezelfde betekenis hebben als in Figuur 3.8b. De overeenkomst tussen de gemiddelde lokalisatie van de drie methodes is hier beter dan in het vorige voorbeeld, zo is de maximale onderlinge afstand 0.25 relek . De methode die een variërende lokalisatie dient te detecteren, vindt in dit voorbeeld een groter verschil dan in de vorige twee gevallen. De maximale afstand tussen twee 61


b: Lokalisatie op basis van BTD-(L, L, 1)

a: SISCOM lokalisatie van het insult

1

1

0.5

0.5

0

LE

−0.5 −1

RE

−0.5

0

0.5

RE

0

1

−0.5 1

0.5

0

N

−0.5

−1

1

N

0.5


RE 0

LE

CWT BTD f=1 (mean)

−0.5

1

0.5

0

Hankel BTD

−1 −1

−0.5

Figuur 3.9: Links: lokalisatie van het epileptisch insult op basis van SISCOM data. Rechts: lokalisatie van het insult op basis van BTD-(L, L, 1) met drie verschillende methoden. De posities van de neus, het linker- en rechteroor zijn in de figuren respectievelijk aangeduid met N, LE en RE. De figuren tonen beide respectievelijk het achteraanzicht, het rechterzijaanzicht en het bovenaanzicht van de locatie.

geschatte lokalisaties van deze methode bedraagt 0.26 relek . Figuur 3.10 toont de veranderingen per coördinaat in functie van het tijdsverloop. Deze toont een licht stijgende trend, hoofdzakelijk in de x en z coördinaat. Figuur 3.9a geeft opnieuw de SISCOM opname ter referentie. We bemerken dat de overeenkomst met de BTD(L, L, 1) methoden redelijk goed is, doch minder consistent dan in het voorbeeld met MA 2. 0.6

0.4

Coordinate

0.2

0

−0.2

−0.4

X coordinate Y coordinate Z coordinate

−0.6

−0.8

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Relative position

Figuur 3.10: De evolutie in de verschillende coördinaten gedurende het verloop van het segment bepaald met een BTD-(L, L, 1) met de frequentie component rang-1.

62

Hoofdstuk 4

Conclusie Het onderzoek voor deze thesis is toegespitst op de voordelen van het gebruik van een BTD-(L, L, 1) tensorontbinding in het probleem van de lokalisatie van epileptische insulten op basis van EEG opnames. Hoofdstuk 1 heeft een inleiding tot de materie gegeven. Het BSS probleem is toegelicht in het kader van EEG lokalisatie. Het gebruik van tensor decompositie technieken als methode voor het oplossen van het BSS probleem is uitgelegd aan de hand van de uniciteitsvoorwaarden voor een ontbinding in rang-1 termen (CPD) alsook voor een ontbinding in termen van lage multilineaire rang (BTD-(L, L, 1)). De definities van beide ontbindingen zijn gegeven. Verder is er een beknopte inleiding gegeven tot het monitoren en analyseren van epileptische aanvallen aan de hand van EEG opnames. De resultaten van eerder onderzoek op basis van een CPD ontbinding van met CWT getensoriseerde EEG data zijn kort samengevat, alsook de tekortkomingen van een dergelijke methode. Hoofdstuk 2 behelst een uitgebreide simulatiestudie die de mogelijkheden van een BTD-(L, L, 1) ontbinding met betrekking tot lokalisatie en analyse van epileptische EEG data onderzoekt. Sectie 2.1 beschrijft de methodiek gebruikt voor het genereren van artificiële EEG data van epileptische insulten. Sectie 2.2 bespreekt de constructie van de derde orde tensor door een CWT van de EEG data waarbij de mexican hat wavelet gebruikt wordt. Dit is de methode die voor het grootste deel van de simulatiestudie onderzocht is. Een tweede methode, op basis van Hankel matrices, is in het laatste deel van de simulatiestudie onderzocht. In Sectie 2.3 zijn de verschillende types ruis beschreven die verder aan de simulaties toegevoegd worden als verstoring. Twee segmenten die respectievelijk musculaire en oculaire artefacten bevatten, zijn met CCA uit de klinische dataset geëxtraheerd. De MLSVD van deze segmenten toont aan dat de oculaire artefacten met een (L, L, 1) term beschreven kunnen worden. De musculaire artefacten bevatten minder structuur. In Sectie 2.4 wordt de analyse uitgevoerd met signalen die epileptische activiteit met een vaste lokalisatie bevatten. Eerst wordt het geval van epileptische activiteit met stationaire frequentie besproken. De belangrijkste conclusies hier zijn dat een CPD ontbinding in dit geval in staat is om de epileptische activiteit uit de data te scheiden en dit tot op zeer lage SNR. Het is wel belangrijk dat het aantal termen in de ontbinding (R) 63

4. Conclusie voldoende groot gekozen wordt, afhankelijk van het type ruis. Indien dit het geval is, zijn de resultaten voor de verschillende types ruis vergelijkbaar en kan de lokalisatie correct uitgevoerd worden tot op zeer lage SNR. De methode is vergeleken met SOBI als referentie BSS methode. De CPD analyse blijkt steeds beter in staat om de tijdreeks uit de data te scheiden, behalve als de ruis bestaat uit een beperkt aantal bronnen. De ruimtelijke component wordt tot op een vergelijkbare nauwkeurigheid uit de data gescheiden. We besluiten dat een CPD ontbinding van een met CWT geconstrueerde tensor die epileptische activiteit van stationaire frequentie en vaste lokalisatie bevat een goede BSS methode vormt. Ook indien er twee of meer bronnen met (verschillende) stationaire frequentie aanwezig zijn. Vervolgens wordt het geval van epileptische activiteit met veranderlijke frequenties besproken. Als voorbeeld is er een chirp signaal met afnemende frequentie gebruikt. Het blijkt dat de CWT van een dergelijk signaal een matrix met een rang groter dan 1 is. De temporele en frequentie component van de geconstrueerde tensor zijn in dit geval van lage rang, de ruimtelijke component blijft rang-1. Het wordt aangetoond dat een BTD-(L, L, 1) ontbinding in dit geval beter geschikt is om de epileptische activiteit uit de tensor te scheiden. Betreffende de ruimtelijke component is het verschil beperkt, maar de tijdfrequentie matrix van lage rang wordt beter geëxtraheerd met BTD-(L, L, 1). Een CPD ontbinding levert slechts een beste rang-1 benadering van de matrix op, ook wordt de epileptische activiteit over meerdere termen in de CPD gespreid indien R vergroot wordt. Aan de hand van een algoritme dat een BTD-(L, L, 1) ontbinding uitvoert op basis van een geclusterde CPD is aangetoond dat L = 2 optimaal is voor de epileptische activiteit. Hogere L laat de modellering van artefacten in de epileptische term toe. De lokalisatie op basis van het BTD-(L, L, 1) en CPD resultaat verschilt slechts minimaal. De BTD-(L, L, 1) ontbinding laat in dit geval dus hoofdzakelijk een betere analyse van het gedrag van het insult in functie van de tijd toe. Indien er meerdere bronnen met (verschillende) veranderlijke frequenties in de data aanwezig zijn, is er aangetoond dat een BTD-(L, L, 1) ontbinding beide bronnen uit de data kan scheiden indien de aanwezige ruis ook een (L, L, 1) structuur heeft. Indien er geen structuur in de ruis aanwezig is, ondervindt de BTD-(L, L, 1) methode moeilijkheden. In dit geval zijn beide bronnen wel uit de data en van elkaar gescheiden met behulp van een CPD ontbinding. Termen in de CPD die te wijten zijn aan eenzelfde bron kunnen, indien gewenst, wel handmatig geclusterd worden om zo de overeenkomst met de tijdfrequentie matrix te verhogen. In Sectie 2.5 worden twee voorbeelden besproken met een bron van epileptische activiteit die van locatie verandert gedurende het verloop van de aanval. Als eerste voorbeeld een bron die op hetzelfde moment ook nog van frequentie verandert. In dit geval zijn alle drie de componenten van de tensor geconstrueerd uit de zuivere epileptische activiteit van lage multilineaire rang. Een ontbinding in BTD-(L, M, N ) termen zou deze term dus mogelijk uit een ruizige tensor kunnen scheiden. We gebruiken echter een opeenvolging BTD-(L, L, 1) ontbindingen ieder uitgevoerd op opeenvolgende, overlappende segmenten uit de tijdreeks. Dit omdat de ruimtelijke 64

component rang-1 is indien de lokalisatie van de bron min of meer vast is gedurende het segment. Deze aanpak laat toe om tot op redelijke SNR de verandering in locatie van de bron te volgen. Het tweede voorbeeld behandelt een van locatie veranderend signaal met vaste frequentie. In dit geval heeft de tensor met zuivere epileptische activiteit een (L, L, 1) structuur waarbij de frequentie component rang-1 is. Met een BTD-(L, L, 1) ontbinding is het dus mogelijk om de veranderende lokalisatie in een enkele term te modelleren. Er wordt aangetoond dat de veranderingen in de locatie een subtiel effect hebben op de simulatie en hoofdzakelijk in de kanalen met laagste activiteit. Deze hebben een lagere SNR dan gemiddeld en bijgevolg is het slechts tot op redelijke SNR mogelijk om de veranderende lokalisatie in één term te modelleren. De resultaten tonen echter wel aan dat de veranderende locatie betrouwbaar bepaald wordt met BTD-(L, L, 1). In Sectie 2.6 wordt BTD-(L, L, 1) als BSS techniek toegepast op tensoren die uit EEG data geconstrueerd zijn door van ieder kanaal een Hankel matrix te vormen. In [16] is aangetoond dat dit tensoren met (L, L, 1) structuur oplevert voor exponentiële polynomen en dat de ontbinding uniek bepaald is. Het is geverifieerd dat de getensoriseerde ictale activiteit een (L, L, 1) structuur heeft, typisch met hogere L dan voor de CWT methode. Het is ook aangetoond dat ictale activiteit, zowel met stationaire als veranderlijke frequentie, tot op lage SNR betrouwbaar uit de ruizige tensor gescheiden kan worden indien de rang L goed gekozen is. Deze tensorisatie methode levert typisch tensoren met grotere dimensies op dan de CWT methode. Aan de hand van de meest significante singuliere vectoren uit de MLSVD kan de tensor echter gecomprimeerd worden, wat computationeel interessant is. Merk op dat de berekening van de MLSVD van een grote tensor computationeel ook veeleisend is. Het is aangetoond dat deze compressiemethode consistente resultaten oplevert, grotendeels onafhankelijk van hoe ver de compressie doorgedreven is. Uit Hoofdstuk 2 besluiten we dat BTD-(L, L, 1) in combinatie met CWT tensorisatie een geschikte BSS methode is om zowel de lokalisatie als analyse van een epileptische aanval uit te voeren indien de aanval van frequentie verandert, of een aanval met stationaire frequentie van locatie verandert. Indien zowel de locatie als de frequentie veranderen is een gesegmenteerde aanpak met BTD-(L, L, 1) geschikt. BTD-(L, L, 1) in combinatie met Hankel tensorisatie levert ook een goede BSS techniek aangezien deze (L, L, 1) structuur promoot. De resultaten met beide methodes zijn vergelijkbaar. In Hoofdstuk 3 worden drie voorbeelden van de analyse op klinische data besproken. Het eerste voorbeeld heeft, op enkele oculaire artefacten na, amper aanwezigheid van ruis. De vier gebruikte methodes, één op basis van CPD en CWT, twee op basis van BTD-(L, L, 1) en CWT en één op basis van BTD-(L, L, 1) en Hankel tensorisatie, vertonen een goede overeenkomst in de geschatte lokalisatie. Ze tonen ook aan dat de aanval van frequentie, noch van lokalisatie verandert. In het tweede klinische voorbeeld is er een frequentieverandering in het insult gedetecteerd aan de hand van BTD-(L, L, 1). De lokalisatie lijkt grotendeels onveranderlijk, maar er is wel een 65

4. Conclusie goede overeenstemming tussen de verschillende methodes onderling, alsook met de beschikbare SISCOM data die ter referentie gebruikt is. Het derde voorbeeld toont opnieuw een goede overeenkomst tussen de methodes onderling en met SISCOM. In dit voorbeeld is een subtiele verandering in lokalisatie opgemerkt aan de hand van BTD-(L, L, 1). We besluiten dat Hoofdstuk 3 de reeds gebruikte methodes gevalideerd heeft en de mogelijkheden van BTD-(L, L, 1) als BSS techniek voor EEG data geïllustreerd heeft. Mogelijkheden voor verder onderzoek over tensor decompositie technieken als BSS techniek voor EEG data zijn mogelijks het gebruik van een ontbinding in (L, M, N ) termen voor zowel een frequentie en lokalisatie verandering in een enkele term te beschrijven, het ontwikkelen van een methode die de lokalisatie real-time uitvoert en onderzoek naar alternatieve methoden om de tensor uit EEG data te construeren.

66

Bĳlagen

67

Bĳlage A

Gebruikte Software Een overzicht van de gebruikte software in dit werk: • Voor het opstellen van de simulaties en het uitvoeren van de bronlokalisatie is de FieldTrip toolbox [38] voor Matlab gebruikt. Meer bepaald de functies: – ft_dipolesimulation: voor het opstellen van simulatie aan de hand van een hoofdmodel, een dipoollocatie, een dipoolmoment en een gewenste tijdreeks. – ft_dipolefitting: voor het fitten van een dipool aan een tijdreeks. Gebruikt voor het bepalen van de locatie van het resultaat van een tensorontbinding. • Voor de constructie van de tensor volgens de eerste methode is er gebruik gemaakt van een CWT. Om na de tensorontbinding terug tot een EEG tijdreeks te komen, is de inverse CWT bepaald. Hiervoor zijn de functies cwtft en icwtft gebruikt, deze zijn beschikbaar in de Matlab Wavelet toolbox. Ze berekenen de wavelet transformatie en zijn inverse op basis van een fft algoritme. De functies cwt en icwtlin zijn ook beschikbaar voor het uitvoeren van de transformatie en zijn inverse. Deze tweede combinatie behoudt echter de amplitudeinformatie niet en heeft ook last van randeffecten, zoals geïllustreerd in Figuur A.1. Hierdoor is er steeds gebruik gemaakt van de eerste algoritmen. • De tensorontbindingen, zowel CPD als BTD-(L,L,1), zijn uitgevoerd met de Tensorlab toolbox [45] voor Matlab. De geclusterde BTD-(L,L,1) ontbinding is uitgevoerd aan de hand van een algoritme beschreven in [43]. • CCA is gebruikt om de musculaire en oculaire artefacten uit klinische data te scheiden. De gebruikte methode is beschreven in [5]. • SOBI is gebruikt als referentie BSS methode om het resultaat van een tensorontbinding mee te vergelijken. De gebruikte methode is beschreven in [2, 3, 8]. • MRICRON [39] is gebruikt voor het uitlezen en visualiseren van SISCOM data.

69

A. Gebruikte Software

40

Amplitude

20

0

−20

−40

Original signal cwt + icwtlin cwtft + icwtft

−60

−80 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Time (s)

Figuur A.1: Invloed van het effect van een CWT en de inverse transformatie op een sinus van 6 Hz met amplitude 1 (zwart). De keuze voor cwt/icwtlin (rood) geeft een verkeerd resultaat: de amplitude is niet behouden, er is een faseverschuiving van π en er zijn randeffecten. De combinatie cwtft/icwtft (groen) geeft een perfecte reconstructie. In beide gevallen is er geopteerd voor een mexican hat wavelet.

70

Bĳlage B

Artikel: Decomposition of a tensor constructed from ictal scalp EEG in rank-(Lr , Lr , 1) terms determines the evolution in location and frequency of the seizure Op de volgende bladzijden is het Engelstalig artikel toegevoegd dat een compendium vormt van de resultaten die eerder in deze thesis toegelicht zijn.

71

1

Decomposition of a tensor constructed from ictal scalp EEG in rank-(Lr , Lr , 1) terms determines the evolution in location and frequency of the seizure Daan Camps 1,∗ , Borb´ala Hunyadi 2 , Laurent Sorber 1 , Sabine Van Huffel 2 , Lieven De Lathauwer 2,3 1 2 3

KU Leuven, Department of Computer Science, Leuven, Belgium

KU Leuven, Department of Electrical Engineering, Leuven, Belgium

KU Leuven Campus Kortrijk, Subfaculty Science and Technology, Kortrijk Belgium

Abstract—In this paper a decomposition of a tensor in rank(Lr , Lr , 1) terms, or BTD-(Lr , Lr , 1), is used to extract ictal activity from noisy EEG data. Two methods to construct a thirdorder tensor from EEG data are discussed. The first method uses a Continuous Wavelet Transformation (CWT). The second method constructs Hankel matrices from every channel. Both methods are tested on artificial epileptic data. It is shown that when the seizure changes frequency or location, the decomposition is capable of extracting the changing behaviour. Finally, the method is also validated on a clinical example and is shown to be in agreement with SISCOM.

I. I NTRODUCTION Epilepsy is a neurological disease which is characterised by abnormal, synchronous discharges of neurons in the brain [1]. This definition is rather vague due to the large diversity of the phenomenon. Epileptic activity in the brain can be registered in EEG recordings. These provide a relatively inexpensive method to measure and analyse the activity in the brain. It has an excellent temporal resolution, determined by the sampling frequency, but a lower spatial resolution. Furthermore the method is sensitive to artefacts originating from nearby electrical equipment or, more importantly, from electrical activity of other physiological origins. Artefacts that are most common in EEG recorded during epileptic seizures are muscular and ocular artefacts. Muscular artefacts are due to uncontrollable contraction of muscles in the face, mostly the jaw, of the patient during the seizure [2]. Ocular artefacts are due to uncontrollable blinking of the eyelids. Both types of artefacts can impede the interpretation of the EEG by a physician. EEG recordings of patients suffering from epilepsy can be roughly divided in two parts: the ictal period, which is the period during which the seizure occurs and the interictal period in between seizures. The ictal period is characterised in EEG data as a period with strong rhythmical waves. The spectral signature of artefacts in ictal EEG typically overlaps with the frequency band of interest for the epileptic activity. Simple filtering techniques in the frequency domain are thus insufficient. This is where Blind Source Separation (BSS) techniques come into play to separate the ictal activity from ∗

Send offprint requests to: [email protected]

noisy data. BSS techniques try to determine the source signals from a (noisy) mixture of these signals. The problem can be denoted as follows, assume that A is the EEG data matrix which contains the different channels as its rows, then BSS aims to decompose this as: A = F ⋅ G + N,

(1)

in which F is an unknown mixing matrix, G is a matrix that contains the source signals and N represents additive noise. The goal of BSS is to estimate F and/or G, given only A. This factorization is in general however not unique and different methods make different assumptions to come to a solution. Canonical Correlation Analysis (CCA) [3] is an example of a BSS technique that solves the problem by constructing sources which are mutually uncorrelated but have a maximal autocorrelation. This technique has been successfully used to remove muscular artefacts from EEG [4], [5]. The approach followed in this paper is different. Where the factorization of a matrix (Eq. 1) is not unique without imposing additional constraints, the decomposition of a higher-order tensor can be unique under mild conditions. Transformation of ictal EEG data to a third order tensor with a continuous wavelet transform (CWT), followed by a decomposition of the tensor in rank-1 terms has previously shown to be capable of extracting and localizing the ictal activity [6], [7]. The decomposition in rank-1 terms is also known as a Canonical Polyadic Decomposition (CPD), and has its origin in [8]. It lasted some time until the concept came in use in psychometric research under the name of Candecomp [9] and Parafac [10]. The uniqueness result of Kruskal [11] was important for the current establishment of CPD as BSS technique. A recent overview of uniqueness results for CPD of a third-order tensor, together with new uniqueness results is given in [12], [13]. In this paper a more general decomposition of a tensor in rank-(Lr , Lr , 1) terms is used as BSS technique on epilepsy tensors. This decomposition is known as a Block Term Decomposition (BTD-(Lr , Lr , 1)) of a tensor and is

2

introduced in [14], [15], [16]. Two methods to construct a tensor from EEG data are used in combination with BTD-(Lr , Lr , 1). The first method uses a CWT of every channel as in [6], [7]. The influence of a change in location and frequency of the seizure on the rank of a CWT tensor with pure ictal activity has been studied. The use of BTD(Lr , Lr , 1) is justified in certain cases. The second method constructs the tensor with Hankel matrices. In [17] it has been shown that this leads to a tensor with an (L, L, 1) structure and BTD-(Lr , Lr , 1) can be used as BSS technique. The organization for the remainder of this paper is as follows: Section II discusses the notation, gives some basic definitions and introduces the tensor decomposition that is used. Section III describes the two methods that are used to form a third-order tensor from EEG data in detail. Section IV describes the results obtained from a realistic simulation study on artificial ictal EEG. Section V validates the method by discussing an example of the analysis on clinical data and compares the result with the available subtraction ictal SPECT co-registered with MRI (SISCOM) data. Section VI finally concludes the paper by summarizing the most important results.

Definition 3. The matrix representations A(1) , A(2) and A(3) of a third-order tensor A ∈ KI1 ×I2 ×I3 are defined as: (A(1) )i1 ,(i2 −1)I3 +i3 = (A)i1 i2 i3 ,

(A(2) )i2 ,(i3 −1)I1 +i1 = (A)i1 i2 i3 ,

(A(3) )i3 ,(i1 −1)I2 +i2 = (A)i1 i2 i3 ,

for all different values of the indices.

Definition 4. A third-order tensor A ∈ KI1 ×I2 ×I3 has multilinear rank (R1 , R2 , R3 ) iff rank(A(1) ) = R1 , rank(A(2) ) = R2 and rank(A(3) ) = R3 . The values R1 , R2 en R3 are respectively the mode-1, mode-2 and mode-3 rank. The definition of a BTD-(Lr , Lr , 1) of a tensor is given by [17]:

Definition 5. A decomposition of a tensor A ∈ KI1 ×I2 ×I3 in a sum of rank-(Lr ,Lr ,1) terms (1 ≤ r ≤ R) is a decomposition of the form R

A = ∑ (Dr ⋅ ETr ) ○ fr , r=1

in which the matrix Gr = Dr ⋅ ETr ∈ KI1 ×I2 has rank = Lr and the vector fr is nonzero. The number of terms R is assumed to be minimal.

II. N OTATION , DEFINITIONS AND TENSOR DECOMPOSITIONS

Vectors are denoted as a ∈ KI , matrices as A ∈ KI×J and third-order tensors as A ∈ KI×J×K . K stands for R or C. Two types of products are important in the use of the methods. These are defined below for N th-order tensors, the case of a third-order is straightforward. Definition 1. The mode-n product of a tensor A ∈ KI1 ×I2 ×⋯×IN with a matrix U ∈ KJ×In is denoted as A ×n U and is of size I1 × ⋯ × In−1 × J × In+1 × ⋯ × IN . The entries of the mode-n product are defined as: In

(A ×n U)i1 ⋯in−1 jin+1 ⋯iN = ∑ ai1 i2 ⋯in ⋯iN ujin . in =1

Definition 2. The outer product A ○ B of a tensor A ∈ KI1 ×⋯×IM and a tensor B ∈ KJ1 ×⋯×JN is the tensor defined by: (A ○ B)i1 ⋯iM j1 ⋯jN = ai1 ⋯iM bj1 ⋯jM ,

for all different values of the indices.

Generalization of the concept of rank of a matrix to a thirdorder tensor is not straightforward because the mode-n rank is in general not constant for all the modes. There are different ways to generalize rank. A possible generalization that is of interest, is called the multilinear rank of a tensor. In order to define this, the matrix representation of a tensor has to be defined first.

Fig. 1: BTD-(Lr , Lr , 1) of a tensor A.

Fig. 1 shows a graphical visualisation of the decomposition of a tensor in rank-(Lr , Lr , 1) terms. In [15] uniqueness results are derived for BTD-(Lr , Lr , 1). Conditions for uniqueness in case the underlying source signals are exponential polynomials are discussed in [17]. Note that CPD is special case of BTD-(Lr , Lr , 1) with Lr = 1 for all r.

Different types of algorithms for the calculation of both a CPD and a BTD-(Lr , Lr , 1) have been derived and discussed in the literature. First, there is the Alternating Least Squares (ALS) algorithm [9] for the calculation of the CPD in which the factor matrices are updated in an alternating manner. Later on new algorithms have been developed that have increased robustness and efficiency, an important category are the Nonlinear Least Squares (NLS) [18] algorithms. Both types of algorithms are available for both CPD and BTD-(Lr , Lr , 1) and have been used by means of their implementation in Tensorlab [19]. A last tensor concept that is used in this paper, is the multilinear SVD (MLSVD) of a tensor. This concept has its origin in [20] and the interpretation is given in [21]. It is a

3

generalisation of the SVD of a matrix. Here it is mostly used to interpret the structure of a tensor.

III. T ENSOR CONSTRUCTION A. CWT method For this method the EEG data is normalised by first subtracting the mean of each channel. Afterwards every channel is divided by its standard deviation. After the decomposition the amplitude information is restored by multiplying the spatial component with the standard deviations. In the next step the CWT of each channel is calculated. The mexican hat wavelet is used as wavelet, in line with [6], [22] where it has previously been used for the analysis of epileptic data. The CWT is calculated in 30 scales, which coincide linearly with the frequency range 1 − 30 Hz. This covers the frequency band in which epileptic activity is expected to occur. A CWT thus transforms the temporal data to a time-frequency plane representation. The different components in the tensor decomposition can be identified as the spatial component, the frequency component and the temporal component. This is the same tensorisation as used in [6], [7]. After the tensor decomposition the signal can reconstructed by an inverse CWT (ICWT) of the retrieved time-frequency plane.

B. Hankel method For this method the standard deviation is again extracted from the data before tensorisation. The mean is not subtracted because it can introduce an additional pole. It also is expected to be close to zero in every channel, so this has limited influence. Every channel, or row vector [a1 a2 ⋯ aN ], is then replaced by a (J × K) Hankel matrix with J + K − 1 = N . This matrix is structured as follows:

⎡ a1 a2 a3 ⋯ aK ⎤⎥ ⎢ ⎢a a3 ⋯ aK aK+1 ⎥⎥ ⎢ 2 ⎢ ⎥ ⋯ aK aK+1 aK+2 ⎥ . ⎢ a3 ⎢ ⎥ ⎢ ⋮ ⋰ ⋰ ⋰ ⋰ ⎥⎥ ⎢ ⎢aJ aJ+1 ⋯ aN −1 aN ⎥ ⎣ ⎦ It has constant entries along its anti-diagonals. If N is odd, a square matrix can be constructed with J = K = N2+1 . After the tensor decomposition the signal can be reconstructed by taking the mean along the anti-diagonals of the retrieved matrix. In [17] this method is introduced as a BSS technique. It is shown that the Hankel matrices have low rank causing the retrieved tensor to have an (L, L, 1) structure. The rank of the Hankel matrix depends on the number and nature of the exponential polynomials that are present in the data. A remark is that when N becomes large, J = K ≈ N2 also becomes large. The added mode in this tensorisation method is thus potentially of much larger size than for the

CWT method where it is chosen as 30. It is computationally more expensive to calculate the decomposition in this case. In order to overcome this problem, two solutions have been used. First, shorter segments of EEG data have been used for this method than for the CWT method. Second, the tensor A with large dimensions can be compressed to smaller size as in [17]. Therefore the multilinear rank of A is estimated by means of ˜ 2 and R ˜ 3 are upper bounds on its MLSVD. Assume that R respectively the mode-2 and mode-3 rank, where these two ˜ are the modes of the tensor that are large. Let U2 ∈ KJ×R2 en ˜ ˜ 2 mode-2 U3 ∈ KK×R3 be formed by respectively the first R ˜ singular vectors and the first R3 mode-3 singular vectors of ˜ ˜ A. Then the tensor T ∈ KI×R2 ×R3 of reduced size, given by: T = A ×2 UT2 ×3 UT3 ,

can be used in a BTD-(Lr , Lr , 1) decomposition. The computation time of the decomposition can be strongly reduced, but calculating the MLSVD of the large tensor is also computationally expensive. The use of smaller EEG segments is thus still meaningful. The approximation of the Hankel matrix of full size Gr ∈ KJ×K can be calculated from the result of the BTD(Lr , Lr , 1) decomposition of T as Gr = U2 ⋅ Hr ⋅ UT3 ,

with Hr ∈ KR2 ×R3 the matrix of low rank from the BTD(Lr , Lr , 1) decomposition of T . From Gr the original signal can again be estimated by taking the mean along its anti-diagonals. ˜

˜

IV. S IMULATION STUDY To generate artificial ictal activity a model is used that places a dipole in a three-shell spherical head model. The dipole is defined by its position in the head, its dipole moment and its waveform. The forward problem is solved in order to determine the intensities at the electrodes [23]. The head model is a volumeconductor consisting of three shells, which represent the brain, skull and scalp. These have conductivities 1 ⋅ 3.3 ⋅ 10−4 and 3.3 ⋅ 10−4 equal to respectively 3.3 ⋅ 10−4 , 16 /Ω/mm [7], [24]. The relative outer radii of the shells are respectively 0.88, 0.92 and 1.00 relek . In which relek is the outer radius of the scalp at which the electrodes are placed. The electrodes are placed at the same positions as in the clinical EEG recordings that are available. This is a placement of 21 electrodes along the 10-20 system [25] with additional T1 and T2 electrodes at the temporal lobe. The sampling frequency of the artificial data is chosen to be 250 Hz, which is consistent with the clinical data. In order to obtain a realistic simulation to test the methods, noise has to be added to the data. Two different types of noise are used: Gaussian random noise (G) and noise resembling muscular artefacts (MA). The second type of noise

4

is extracted with CCA from a heavily contaminated clinical EEG recording. The nine sources with lowest autocorrelation have been used to construct the noise matrix. Before adding the noise matrix N to the simulation matrix S, it is normalised to the same Frobenius norm, N=N⋅

∥S∥F . ∥N∥F

The observed matrix A is constructed as: A(λ) = S + λ ⋅ N.

Changing the parameter λ changes the amount of noise. This is quantified by the signal-to-noise ratio (SNR) as: SNR(λ) =

RMS(S) , RMS(λ ⋅ N)

in which RMS is the root mean square value. This is defined for a matrix M ∈ KK×L consisting of K channels, each with L datapoints, as: ¿ Á 1 K L−1 2 À RMS(M) = Á ∑ ∑ (M(k, l)) . K ⋅ L k=1 l=0

The last thing needed for performing the simulation study, is a method to estimate a seizure location from the result of the tensor decomposition. The localisation of the simulations is performed by a dipole fit on the extracted ictal EEG which is converted back to a matrix. The localisation of clinical seizures in the next section is based on solely the spatial component in the decomposition.

A. CWT method If the simulated ictal activity has a stationary frequency and fixed location, the MLSVD of the tensor containing pure ictal activity shows that it is rank-1 in all three modes. In this case a CPD on the noisy tensor is capable of reliably extracting the three rank-1 components of the ictal activity in one term, this is already confirmed in [7]. The method works best when the number of terms in the decomposition, R, is chosen well. If it is too low, part of the noise will be modelled in the rank-1 term containing ictal activity leading to a poorer fit. If it is too high, the data will be overmodelled. It becomes more interesting to use a BTD-(Lr , Lr , 1) decomposition if the seizure either changes in frequency or in location. The first case, which is observed in clinical seizures, is modelled by the use of a so-called chirp signal. This is a sinusoidal-like signal but changes from frequency during the course of the signal. Fig. 2a shows a 2s segment of a simulation with a chirp that has a frequency drop from 12 to 10 Hz. Fig. 2b shows the same simulation but with G+MA noise added. The MLSVD of the tensor with pure ictal activity reveals that in this case the spatial component stays rank-1, but the frequency and temporal component are of higher rank. They

a: Chirp: 12 − 10 Hz

b: with added noise

Original simulation


Simulation with noise added

7 uV

0

Fp2 F8 T4 T6 O2 F4 C4 P4 Fz Cz Pz Fp1 F7 T3 T5 O1 F3 C3 P3 T2 T1 0.5

1 Time (s)

1.5

2

28 uV

0

0.5

1 Time (s)

1.5

2

Fig. 2: Dipole simulation of ictal activity with fixed location and variable frequency. The SNR of the data on the right is 0.55.

are approximately rank-4, the fifth singular value is less than 1% of the first singular value. It is however observed that when using L > 2 in the term that describes the ictal activity, the result is more affected by noise than in the case L = 2. This result is confirmed by an algorithm that computes the BTD-(Lr , Lr , 1) from a clustering on a CPD and automatically determines the ranks of the different terms [26]. An explanation for this effect is that a higher L gives a higher degree of freedom such that part of the noise with similar spatial distribution will be modelled in the same term, leading to a poorer resemblance with the pure ictal activity. Fig. 3 gives the correlation between both the original and extracted time-frequency matrices and the original and extracted amplitude or spatial distribution and this for a description by a rank-(L, L, 1) term and a rank-1 term. It shows that in this case the use of BTD-(Lr , Lr , 1) leads to a better modelling of the temporal characteristics of the signal, the spatial distribution is modelled equally well. The localisation error is consequently comparable for the two cases. An important remark here is that the BTD-(Lr , Lr , 1) involves more rank-1 terms than the CPD, where only 2 terms are considered. If more than 2 terms are used in the CPD, multiple of these terms will be related to the ictal activity, each with a lower temporal correlation. The CPD retrieves in this case essentially a BTD-(Lr , Lr ,1) term. Fig. 4 shows an example of the extracted time-frequency matrices, both with CPD and BTD-(Lr , Lr , 1). The rank-1 CPD result shows no clear indication of a frequency change during the seizure, while this can be identified in the BTD(Lr , Lr , 1) result.

A second case in which a BTD-(Lr , Lr , 1) decomposition gives an advantage over CPD is when the seizure has a stationary frequency and a changing location. This a realistic situation since it is known that ictal activity tends to spread in the brain. The simulation is constructed by concatenating

5

80 simulations of 0.1s with 10 Hz ictal activity, each at a slightly different location such that the dipole travels a total distance of 1.17 relek along a straight line. This rather large distance is used because the effects of a changing location on the simulation are small and moreover largest in the channels with lowest activity. This is illustrated in Fig. 5a: the five channels with the largest change in activity are also these with lowest amplitude. Since the noise is distributed evenly over all the channels the average SNR in these five channels is lower than the overall SNR.

a: Correlation between time-frequency matrices 1

2 22

2 2

2 2 2 2

2 2 2 2

2 2

0.95

Correlation

1 0.9

0.85

1

a

0.8


0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Relative activity change / Relative amplitude (%)

MA − CPD (2 terms) 0.75

1

SNR

b: Correlation between amplitude distributions 1

Correlation

0.95

0.9

0.85

100 90 80 70 60 50 40 30 20 10 0


b

0.8

MA − BTD−(L,L,1)c MA − CPD (2 terms) 0.75 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.07

1


SNR 0.06

5

Frequency samples

Frequency samples

5

10

15

20

25

10

15

20

25

30

30 200

400

600

800

1000

1200

Time samples

1400

1600

1800

2000

200

400

600

800

1000

1200

1400

0.04

0.03

0.02

0.01

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Relative position

Fig. 5: Top: Statistics of the signal that changes in position. Black indicates the relative change in activity in a channel during the seizure, white indicates the maximal amplitude in each channel relative to the maximal amplitude over all channels. The arrow indicates the direction of the change in the channel. Bottom: Perpendicular distance between the path of the source and the estimated location from a dipole fit at 50 different time instances during the segment, retrieved from a single term. The change in location is modelled with term with L = 2.

b: BTD-(Lr , Lr , 1)

a: CPD

0.05

d (relek)

Fig. 3: Comparison between the result of a CPD with its terms clustered to a BTD-(Lr , Lr , 1) decomposition and a CPD decomposition in function of SNR for the description of a frequency changing signal with a single term. The rank of the ictal term, which is determined with the clustered algorithm, is indicated in the top figure. MA noise is added to the data. The used chirp signal lasts 8s or 2000 samples and is entirely used to form the tensor with the CWT method. It has a frequency drop from 12 to 6 Hz. The BTD-(Lr , Lr , 1) decomposition clusters 10 rank-1 terms together in 3 terms. From this result the ictal term is identified.

1600

1800

2000

Time samples

Fig. 4: Time-frequency matrices of a frequency changing signal modelled with respectively a term in CPD, and a term in BTD(Lr , Lr , 1). The same signal as in Fig. 3 is used.

The structure of the tensor with pure ictal activity is in this case also (L, L, 1) but here the frequency component is rank-1. A BTD-(Lr , Lr , 1) decomposition is thus capable of extracting the change in location in a single term. Fig. 5b shows this for 50 different locations in function of time at different SNR. Both the average SNR over all channels (TOT) and the SNR in the low activity channels (LA) are given. The method is capable of estimating the changing location, it is however noted that at lower SNR the tends to cluster towards the mean position of the dipole. This can again

6

B. Hankel method In this case the structure of the tensor with pure ictal activity is (L, L, 1) for both a stationary frequency and a variable frequency signal. L does however increase in the second case. The results of the method are shown in Fig. 6 for both a sine and chirp signal at three different SNR and different compression sizes. The sinusoidal signal has a frequency of 6 Hz, the chirp signal decreases in frequency from 12 to 6 Hz. Both signals contain 999 samples and are tensorised with J = K = 500. It is clear that this method reliably extracts the ictal activity from the noisy data, both the temporal and spatial characteristics are modelled well. Moreover there is less influence of noise for higher L in comparison with the CWT method and it gives consistent results for all compression sizes. a: Temporal correlation 1.01 1


0.99 0.98 0.97 0.96


0.95 0.94 0.93 0.92

50

100

150

200

250

300

350

400

Compressed size

b: Spatial correlation 1 0.998

Spatial correlation

0.996 0.994 0.992


0.99 0.988 0.986 0.984

50

100

150

200

250

300

350

400

Compressed size

Fig. 6: Results in function of compression at different SNR, both for a sine (L = 2) and a chirp signal (L = 18). The choice of L is based upon the MLSVD of the tensor with pure ictal activity.

V. C LINICAL VALIDATION The example discussed in this section is based upon an ictal EEG measurement categorised as having some presence of artefacts which have (limited) influence on the interpretation of the EEG. Fig. 7 shows the ictal rank-2 time-frequency

plane. This is extracted with BTD-(Lr , Lr , 1) in combination with the CWT method performed on the first 20s of the EEG after the start of the seizure. The result reveals a change in frequency during the seizure from about 5 Hz after 4s to about 8 Hz after 16s.

5

Frequency samples

be explained by the small influence the change in location has.

10

15

20

25

30 500

1000

1500


3500

4000

4500

5000

Fig. 7: Time-frequency matrix extracted from the clinical EEG with the CWT method and BTD-(Lr , Lr , 1) with the spatial component rank-1. The decomposition involved 4 rank-2 terms, the second term is shown and associated with the ictal activity.

Fig. 8 shows both the SISCOM location and the location determined with BTD-(Lr , Lr , 1). The results are in good agreement, which validates the method. The three different methods mutually differ, but the maximal distance between them is 0.39 relek , meaning that they are all still in agreement with SISCOM. The CWT method with f=1 is performed on the 20s segment and reveals no significant change in location during the seizure. The location is estimated every second from the ictal term, the maximal separation between these 20 estimated locations is 0.14 relek . VI. C ONCLUSION In this paper it has been shown that a BTD-(Lr , Lr , 1) decomposition in combination with ictal EEG tensors and the CWT method is capable to model a change in frequency of the seizure if the spatial component is rank-1. Furthermore it is also possible to model a change in location if the frequency component is rank-1. Construction of the tensor with Hankel method is also capable of extracting the ictal activity. This method promotes an (L, L, 1) structure for exponential polynomials and in particular oscillating signals. The methods have been tested and found to be reliable up to low SNR and thorough compression in case of the Hankel method. The methods are validated with clinical EEG data. They are in good agreement with SISCOM. BTD-(Lr , Lr , 1) has shown that the seizure in the clinical example undergoes a change in frequency, but has a fairly constant location. A possible starting point for further study could be the use of BTD-(Lr , Mr , Nr ) in combination with signals that change both in frequency and location. However in this case BTD-(Lr , Lr , 1) on a partitioned EEG is also capable to model both changes.

7

a: SISCOM location

b: BTD-(Lr , Lr , 1) location 1

1

0.5

0.5

0

LE

−0.5 −1

RE

−0.5

0

0.5

RE

0

1

−0.5 1

0.5

0

N

−0.5

−1

1

N

0.5


RE 0

LE

CWT BTD f=1 (mean)

−0.5

1

0.5

0

−0.5

Hankel BTD

−1 −1

Fig. 8: Top: Location of the seizure determined with SISCOM. Bottom: Location determined with BTD-(Lr , Lr , 1) in combination with the CWT method, both spatial rank-1 (s=1) and frequency rank1 (f=1), and the Hankel method. The ictal term from CWT s=1 is the same as in Fig. 7. The ictal term from CWT f=1 is determined as the first term in a BTD-(Lr , Lr , 1) involving 4 rank-2 terms. “sep” indicates the 20 different locations, which are close together. “mean” indicates the mean of these locations. The tensor with the Hankel method is formed from 999 samples starting 13s after seizure onset, J = K = 500. The decomposition involved 3 rank-2 terms and the ictal activity is associated with the third term. The positions of the nose, left- and right ear are respectively indicated as N, LE, RE.

R EFERENCES [1] W. Blume and et. al., “Glossary of descriptive terminology for ictal semiology: Report of the ilae task force on classification and terminology,” Epilepsia, 2001. [2] I. Gonchorova, D. McFarland, T. Vaughan, and J. Wolpaw, “EMG contamination of EEG: Spectral and topographical characteristics,” Clinical Neurophysiology, 2003. [3] M. Borga and H. Knutsson, “A canonical correlation approach to blind source separation,” IEEE PAMI, Tech. Rep., 2001. [4] W. De Clercq, A. Vergult, B. Vanrumste, W. Van Paesschen, and S. Van Huffel, “Canonical correlation analysis applied to remove muscle artifacts from the electroencephalogram.” IEEE Trans. Biomed. Engineering, no. 12, pp. 2583–2587. [5] A. Vergult, W. De Clercq, A. Palmini, B. Vanrumste, P. Dupont, S. Van Huffel, and W. Van Paesschen, “Improving the interpretation of ictal scalp EEG: BSS-CCA algorithm for muscle artifact removal,” Epilepsia, vol. 48, pp. 550–558, 2007. [6] E. Acar, C. Aykut-bingol, H. Bingol, R. Bro, and B. Yener, “Multiway analysis of epilepsy tensors,” Bioinformatics, 2007.

[7] M. De Vos, A. Vergult, L. De Lathauwer, W. De Clercq, S. Van Huffel, P. Dupont, A. Palmini, and W. Van Paesschen, “Canonical decomposition of ictal scalp EEG reliably detects the seizure onset zone,” NeuroImage, vol. 37, pp. 844–854, 2007. [8] F. Hitchcock, The Expression of a Tensor Or a Polyadic as a Sum of Products, ser. Contributions from the Department of Mathematics. sn., 1927. [9] J. D. Carroll and J.-J. Chang, “Analysis of individual differences in multidimensional scaling via an n-way generalization of “Eckart-Young” decomposition,” Psychometrika, vol. 35, no. 3, pp. 283–319, 1970. [10] R. Harshman, “Foundations of the parafac procedure: Models and conditions for an “explanatory” multi-modal factor analysis,” UCLA Working Papers in Phonetics, vol. 16, pp. 1–84, 1970. [11] J. B. Kruskal, “Three-way arrays: rank and uniqueness of trilinear decompositions, with applications to arithmetic complexity and statistics,” Psychometrika, vol. 18, pp. 95–138, 1977. [12] I. Domanov and L. De Lathauwer, “On the Uniqueness of the Canonical Polyadic Decomposition of Third-Order Tensors Part I: Basic Results and Uniqueness of One Factor Matrix,” SIAM J. Matrix Anal. Appl., to appear. [13] ——, “On the Uniqueness of the Canonical Polyadic Decomposition of Third-Order Tensors Part II: Uniqueness of the Overall Decomposition,” SIAM J. Matrix Anal. Appl., to appear. [14] L. De Lathauwer, “Decompositions of a higher-order tensor in block terms - part I: Lemmas for partitioned matrices,” Siam J. Matrix Anal. Appl., vol. 30, no. 3, pp. 1022–1033, 2008. [15] ——, “Decompositions of a higher-order tensor in block terms - part II: Definitions and uniqueness,” Siam J. Matrix Anal. Appl., vol. 30, no. 3, pp. 1033–1066, 2008. [16] L. De Lathauwer and D. Nion, “Decompositions of a higher-order tensor in block terms - part III: Alternating least squares algorithms,” Siam J. Matrix Anal. Appl., vol. 30, no. 3, pp. 1067–1083, 2008. [17] L. De Lathauwer, “Blind separation of exponential polynomials and the decomposition of a tensor in rank-(Lr , Lr , 1) terms,” Siam J. Matrix Anal. Appl., vol. 32, no. 4, pp. 1451–1474, 2011. [18] L. Sorber, M. Van Barel, and L. De Lathauwer, “Optimization-based algorithms for tensor decompositions: canonical polyadic decomposition, decomposition in rank-(Lr , Lr , 1) terms and a new generalization,” 2012. [19] ——, “Tensorlab v1.0, available online, url: http://esat.kuleuven.be/sista/tensorlab/,” February 2013. [20] L. R. Tucker, “Some mathematical notes on three-mode factor analysis,” Psychometrika, vol. 31, no. 3, pp. 279–311, 1966. [21] L. De Lathauwer, B. De Moor, and J. Vandewalle, “A multilinear singular value decomposition,” SIAM J. Matrix Anal. Appl., vol. 21, no. 4, pp. 1253–1278, 2000. [22] M. Latka, Z. Was, A. Kozik, and B. West, “Wavelet analysis of epileptic spikes,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys., vol. 67, pp. 105– 122, 2003. [23] R. Oostenveld, P. Fries, E. Maris, and J. Schoffelen, “Fieldtrip: Open source software for advanced analysis of meg, eeg, and invasive electrophysiological data.” Computational Intelligence and Neuroscience, vol. 2011, 2011. [24] T. Oostendorp, J. Delbeke, and S. D.F., “The conductivity of the human skull: Results of in vivo and in vitro measurements,” IEEE transactions on biomedical engineering, vol. 47, no. 11, 2000. [25] M. Nuwer, G. Comi, and E. R., “Ifcn standards for digital recording of clinical eeg,” Electroencephalography and Clinical Neurophysiology, 1998. [26] L. Sorber, M. Van Barel, and L. De Lathauwer, “Algoritmen voor de rang-(Lr ,Lr ,1) blok term decompositie,” Master’s thesis, KU Leuven, 2010.

Bĳlage C

Poster Op de volgende bladzijde is de poster ingevoegd die gebruikt is voor de masterproefbeurs van de master Wiskundige Ingenieurstechnieken.

79

C. Poster

Master Wiskundige ingenieurstechnieken

Masterproef Daan Camps

Promotoren Prof. Lieven De Lathauwer Prof. Sabine Van Huffel Academiejaar 2012-2013

Methode

• Klinische data: epileptische activiteit uit EEG data te scheiden. Lokalisatie in overeenstemming met SISCOM.

• Simulatiestudie: epileptische activiteit tot op zeer lage SNR uit de data te scheiden en een betrouwbare lokalisatie blijft mogelijk. → BTD-(L,L,1) nuttig bij veranderlijke frequenties

Resultaten

Lokalisatie van epileptische aanvallen met tensor decompositie technieken Situering

• Constructie van de tensor via CWT

• Identificatie en lokalisering epileptische term

• BTD-(L,L,1) / termen van lage ML rang

• Ontbinding van de tensor: • CPD / rang-1 termen

• Blind Source Separation (BSS) probleem: 𝑋 = 𝐴 ∗ 𝑆𝑇 • EEG metingen nemen een mengeling van verschillende bronnen gesitueerd in en rondom de hersenen waar. • Eigenschappen EEG: • Goedkoop • Uitstekende tijdsresolutie • Mindere ruimtelijke resolutie • Gevoelig aan artefacten • Epileptische insulten bestaan uit een ictale periode. In EEG waar te nemen als ritmische activiteit.

Doelstelling Een methode ontwikkelen om epileptische activiteit betrouwbaar uit EEG opnames te scheiden om zo de interpretatie en lokalisatie van de aanval te vergemakkelijken.

Figuur C.1: Poster voor de masterproefbeurs

80

Bibliografie [1] E. Acar, C. Aykut-bingol, H. Bingol, R. Bro, and B. Yener. Multiway Analysis of Epilepsy Tensors. Bioinformatics, 2007. [2] A. Belouchrani, K. Abed-meraim, J. F. Cardoso, and E. Moulines. Second Order Blind Separation of Temporally Correlated Sources., 1993. [3] A. Belouchrani and A. Cichocki. A Robust Whitening Procedure in Blind Source Separation Context. Electronics Letters, 36:2050–2051, 2000. [4] W. Blume and et. al. Glossary of Descriptive Terminology for Ictal Semiology: Report of the ILAE Task Force on Classification and Terminology. Epilepsia, 2001. [5] M. Borga and H. Knutsson. A Canonical Correlation Approach to Blind Source Separation. Technical report, IEEE PAMI, 2001. [6] R. Bro. PARAFAC. Tutorial and applications. Chemometrics and Intelligent Laboratory Systems, 38:149–171, 1997. [7] J. D. Carroll and J.-J. Chang. Analysis of individual differences in multidimensional scaling via an n-way generalization of “Eckart-Young” decomposition. Psychometrika, 35(3):283–319, 1970. [8] A. Cichocki and S. Amari. Adaptive Blind Signal and Image Processing: Learning Algorithms and Applications. John Wiley & Sons, Inc., 2002. [9] P. Comon and C. Jutten, editors. Handbook of Blind Source Separation, Independent Component Analysis and Applications. Academic Press, 2010. [10] I. Daubechies. Ten lectures on wavelets. Society for Industrial and Applied Mathematics, 1 edition, June 1992. [11] W. De Clercq. Advanced preprocessing techniques and nonlinear signal analysis applied to scalp electroencephalograms for the prediction of epileptic seizures. PhD thesis, 2005. [12] W. De Clercq, A. Vergult, B. Vanrumste, W. Van Paesschen, and S. Van Huffel. Canonical Correlation Analysis Applied to Remove Muscle Artifacts from the Electroencephalogram. IEEE Trans. Biomed. Engineering, (12):2583–2587. 81

Bibliografie [13] L. De Lathauwer. A link between the canonical decomposition in multilinear algebra and simultaneous matrix diagonalization. Siam J. Matrix Anal. Appl., 28:642–666, 2006. [14] L. De Lathauwer. Decompositions of a higher-order tensor in block terms - Part I: Lemmas for partitioned matrices. Siam J. Matrix Anal. Appl., 30(3):1022–1033, 2008. [15] L. De Lathauwer. Decompositions of a higher-order tensor in block terms - Part II: Definitions and Uniqueness. Siam J. Matrix Anal. Appl., 30(3):1033–1066, 2008. [16] L. De Lathauwer. Blind separation of exponential polynomials and the decomposition of a tensor in rank-(Lr , Lr , 1) terms. Siam J. Matrix Anal. Appl., 32(4):1451–1474, 2011. [17] L. De Lathauwer. Block Component Analysis, a new concept for blind source separation. In Proc. of 10th International Conference on Latent Variable Analysis and Signal Separation, pages 1–8, 2012. [18] L De Lathauwer, B De Moor, and J Vandewalle. A Multilinear Singular Value Decomposition. SIAM J. Matrix Anal. Appl., 21(4):1253–1278, 2000. [19] L. De Lathauwer and D. Nion. Decompositions of a higher-order tensor in block terms - Part III: Alternating Least Squares Algorithms. Siam J. Matrix Anal. Appl., 30(3):1067–1083, 2008. [20] M. De Vos. Decomposition methods with applications in neuroscience. PhD thesis, Faculty of Engineering, KU Leuven, 2009. [21] M. De Vos, L. De Lathauwer, and S. Van Huffel. Recent Advances in Biomedical Signal Processing, chapter Decomposition Techniques in Neuroscience, pages 1–25. In Gorriz et al. [25], 2009. [22] M. De Vos, A. Vergult, L. De Lathauwer, W. De Clercq, S. Van Huffel, P. Dupont, A. Palmini, and W. Van Paesschen. Canonical decomposition of ictal scalp EEG reliably detects the seizure onset zone. NeuroImage, 37:844–854, 2007. [23] C. Eckart and G. Young. The approximation of one matrix by another of lower rank. Psychometrika, 1(3):211–218, 1936. [24] I. Gonchorova, D. McFarland, T. Vaughan, and J. Wolpaw. EMG contamination of EEG: Spectral and topographical characteristics. Clinical Neurophysiology, 2003. [25] J.M. Gorriz, E.W. Lang, and J. Ramirez, editors. Recent Advances in Biomedical Signal Processing. Bentham Publishers, 2010. 82

Bibliografie [26] R. Harshman. Foundations of the PARAFAC procedure: Models and conditions for an “explanatory” multi-modal factor analysis. UCLA Working Papers in Phonetics, 16:1–84, 1970. [27] R. Harshman. Determination and proof of minimum uniqueness conditions for PARAFAC. UCLA Working Papers in Phonetics, 22:111–117, 1972. [28] F.L. Hitchcock. The Expression of a Tensor Or a Polyadic as a Sum of Products. Contributions from the Department of Mathematics. sn., 1927. [29] T.P. Jung, C. Humphries, and et. al. Extended ICA removes artifacts from EEG recordings. Advances in Neural information processing systems, 1998. [30] T. G. Kolda and B. W. Bader. Tensor Decompositions and Applications. Siam Review, 51(3):455–500, 2009. [31] J. B. Kruskal. Three-way arrays: rank and uniqueness of trilinear decompositions, with applications to arithmetic complexity and statistics. Psychometrika, 18:95– 138, 1977. [32] M. Latka, Z. Was, A. Kozik, and B.J. West. Wavelet analysis of epileptic spikes. Phys. Rev. E Stat. Nonlin. Soft Matter Phys., 67:105–122, 2003. [33] X. Liu and N. Sidiropoulos. Cramer-Rao lower bounds for low-rank decomposition of multidimensional arrays. IEEE Transactions on Signal Processing, 49:2074–2086, 2001. [34] M. Mørup, L. K. Hansen, C. S. Hermann, J. Parnas, and S. M. Arnfred. Parallel Factor Analysis as an exploratory tool for wavelet transformed event-related EEG. NeuroImage, 29(3):938–947, 2006. [35] H. Nam and et. al. Independent Component Analysis of ictal EEG in medial temporal lobe epilepsy. Epilepsia, 2002. [36] M.R. Nuwer, G. Comi, and Emerson R. IFCN standards for digital recording of clinical EEG. Electroencephalography and Clinical Neurophysiology, 1998. [37] T.F. Oostendorp, J. Delbeke, and Stegeman D.F. The Conductivity of the Human Skull: Results of In Vivo and In Vitro Measurements. IEEE transactions on biomedical engineering, 47(11), 2000. [38] R. Oostenveld, P. Fries, E. Maris, and JM Schoffelen. FieldTrip: Open Source Software for Advanced Analysis of MEG, EEG, and Invasive Electrophysiological Data. Computational Intelligence and Neuroscience, 2011, 2011. [39] C. Rorden, H.-O. Karnath, and L. Bonilha. Improving lesion-symptom mapping. Journal of Cognitive Neuroscience, 19:1081–1088, 2007. [40] P. Schönemann. An algebraic solution for a class of subjective metrics models. Psychometrika, 37(4):441–451, 1972. 83

Bibliografie [41] N.D. Sidiropoulos and R. Bro. On the Uniqueness of Multilinear Decomposition of n-way arrays. J. Chemometrics, 14:229–239, 2000. [42] A. Smilde, R. Bro, and Geladi P. Multi-way analysis with applications in the chemical sciences. John Wiley & Sons, 2004. [43] L. Sorber, M. Van Barel, and L. De Lathauwer. Algoritmen voor de rang(Lr ,Lr ,1) Blok Term Decompositie. Master’s thesis, KU Leuven, 2010. [44] L. Sorber, M. Van Barel, and L. De Lathauwer. Optimization-based algorithms for tensor decompositions: canonical polyadic decomposition, decomposition in rank-(Lr , Lr , 1) terms and a new generalization. 2012. [45] L. Sorber, M. Van Barel, and L. De Lathauwer. Tensorlab v1.0, Available online, URL: http://esat.kuleuven.be/sista/tensorlab/. February 2013. [46] J.M.F. Ten Berge and Sidiropoulos N. D. On uniqueness in CANDECOMP/PARAFAC. Psychometrika, 67:399–409, 2002. [47] L. R. Tucker. Some mathematical notes on three-mode factor analysis. Psychometrika, 31(3):279–311, 1966. [48] G.L. Wallstrom and et. al. Automatic correction of ocular artifacts in the EEG: a comparison of regression-based and component-based methods. Int J Psychophysiol, 2004.

84

KU Leuven Faculteit Ingenieurswetenschappen

2012 – 2013

Fiche masterproef Student: Daan Camps Titel: Lokalisatie van epileptische aanvallen met tensor decompositie technieken UDC : 621.3 Korte inhoud: Tensor decompositie technieken, met name een ontbinding in rang-1 termen of CPD, zijn reeds eerder succesvol toegepast om epileptische activiteit uit ruizige EEG data te scheiden en een lokalisatie van het insult uit te voeren. De tensor wordt hier typisch geconstrueerd door de CWT te berekenen van ieder kanaal. Tekortkomingen van deze methode komen naar voren als de epileptische activiteit van frequentie of lokalisatie verandert. Indien de frequentie veranderlijk is, zijn de temporele en frequentie component van de epileptische tensor verschillend van rang-1. Indien de lokalisatie verandert, zijn de ruimtelijke en temporele component verschillend van rang-1. In dit onderzoek is een ontbinding in termen van lage multilineaire rang of BTD-(L, L, 1) gebruikt om epileptische activiteit uit een ruizige EEG tensor, geconstrueerd met CWT, te scheiden. Afzonderlijk zijn een verandering in frequentie en een veranderlijke lokalisatie succesvol gemodelleerd met een BTD(L, L, 1) ontbinding en kunnen beide tot op lage SNR bepaald worden. De methode is uitvoerig getest in een uitgebreide simulatiestudie en vergeleken met SOBI als referentie BSS techniek. Nadien zijn er drie voorbeelden uit een klinische dataset besproken die een verdere validatie van de methode geven. Voor deze voorbeelden is de lokalisatie bepaald BTD-(L, L, 1) steeds in overeenstemming met de beschikbare SISCOM data. De resultaten van de BTD-(L, L, 1) ontbinding duiden in één van de drie voorbeelden op een frequentieverandering gedurende het insult. In een ander voorbeeld is er een kleine verandering in de lokalisatie van het insult opgemerkt met BTD-(L, L, 1). Naast de tensorisatie met CWT is er een tweede methode onderzocht die de tensor construeert door uit elk kanaal een Hankel matrix te vormen. Dit geeft aanleiding tot matrices van lage rang en bijgevolg een (L, L, 1) structuur in de tensor met zuivere epileptische activiteit. Deze methode laat in combinatie met BTD-(L, L, 1) opnieuw een betrouwbare scheiding toe van de epileptische activiteit uit de ruizige tensor tot op lage SNR. Bij de klinische validatie zijn de resultaten van deze methode in goede overeenstemming met de andere methode en de SISCOM data. Thesis voorgedragen tot het behalen van de graad van Master of Science in de ingenieurswetenschappen: wiskundige ingenieurstechnieken Promotoren: Prof. dr. ir. Lieven De Lathauwer Prof. dr. ir. Sabine Van Huffel Assessoren: Prof. dr. ir. Marc Van Barel Prof. dr. ir. Johan Suykens Begeleiders: ir. Borbála Hunyadi ir. Laurent Sorber

Lokalisatie van epileptische aanvallen met tensor decompositie technieken

Recommend Documents