Detekce korové aktivity vyvolané míšním neurostimulátorem u pacientů trpících nezvladatelnou neuropatickou bolestí Lukáš Svoboda1 , Pavel Sovka1 , Andrej Stančák2 1
České vysoké učení technické v Praze, Fakulta elektrotechnická 2 Univerzita Karlova v Praze, 3. lékařská fakulta
[email protected],
[email protected],
[email protected]
Abstrakt: Spektrální analýza pomocí výkonové spektrální hustoty (PSD) může být úspěšně aplikována při detekci změn v EEG. Užití této metody je demonstrováno na EEG získaném od 10 pacientů, trpících nezvladatelnou neuropatickou bolestí v oblasti beder a dolních končetin. Tato bolest je tlumena pomocí elektrické stimulace dorzální páteřní míchy (SCS). S využitím spektrální analýzy bylo zjištěno, že SCS vybudí oscilace v mozkové kůře. U osmi pacientů byl při SCS zaznamenán signifikantní nárůst PSD na stimulačním kmitočtu.
1.
Úvod
Pro vyhodnocení změn signálů a vzájemných vazeb mezi nimi v současnosti existuje řada metod. Své stálé místo si udržují postupy založené na Fourierově transformaci, k nimž patří i výkonová spektrální hustota (PSD) použitá s několika rozšířeními k detekci změn výkonu EEG signálu na konkrétních kmitočtech. Data byla získána v rámci spolupráce s Ústavem normální, patologické a klinické fyziologie 3. lékařské fakulty Univerzity Karlovy. Jedná se o EEG záznamy 10 pacientů trpících failed back surgery syndromem (FBSS), což je onemocnění způsobující nezvladatelnou neuropatickou bolest v oblasti beder a dolních končetin, která často odolává léčbě farmaky [4]. Jedním z nefarmakologických postupů, které pomáhají tuto bolest mírnit, je elektrická stimulace dorzální páteřní míchy (spinal cord stimulation – SCS), při níž dochází k významnému potlačení intenzity vnímané bolesti. Korové projevy míšní stimulace již byly u omezeného počtu pacientů zkoumány pomocí fMR [2, 3], podstata elektrofyziologických změn korové aktivity při míšní stimulaci však dosud objasněna nebyla. V příspěvku je popsán postup, jímž bylo prostřednictvím PSD prokázáno, že SCS vyvolává v mozkové kůře oscilační aktivitu na stejném kmitočtu, jako je stimulační, a jímž byly lokalizovány oblasti, kde k tomu dochází.
2.
Detekce statisticky významných změn ve spektru EEG
Jako první bylo třeba zjistit, zda působením SCS v mozkové kůře oscilační aktivita vůbec vzniká. Předpokladem bylo, že vyvolané rytmy budou mít kmitočet shodný se stimulačním a mohou se vyskytovat v kterékoli oblasti kůry. K dispozici jsme měli EEG záznamy pořízené pomocí 111 povrchových elektrod, a to dva od každého pacienta – jeden při zapnuté stimulaci, druhý bez ní. 2.1.
Odhad PSD
Vyskytne-li se rozdíl mezi PSD párových EEG záznamů pořízených při zapnuté, respektive vypnuté míšní stimulaci, znamená to, že stimulace ovlivňuje elektrickou aktivitu mozku. Pokud se výkonový nárůst (ve srovnání se záznamem pořízeným při vypnuté stimulaci) projeví na stimulačním kmitočtu, je pravděpodobné, že SCS evokuje korové rytmy o stejném kmitočtu, jako je ten její. Výpočet odhadu PSD z reálných dat je odvozen v řadě publikací [1, 5], proto jsou uvedeny pouze výsledné vztahy v diskretizovaném tvaru, jak byly aplikovány při analýze našich EEG záznamů. Pro konzistentní odhad výkonové spektrální hustoty bylo použito vyhlazování spektra Welchovou metodou, jejíž jednotlivé kroky a mezivýsledky přibližuje schéma na obrázku 1.
Obrázek 1: Princip vyhlazování spektra Welchovou metodou Jedná se o přímou metodu využívající diskrétní Fourierovy transformace váhovaných segmentů vstupního signálu. Vstupní posloupnost x[n] délky N rozdělíme na K navazujících úseků xk [l], každý o délce L vzorků – platí tedy N = KL (existuje i varianta, kdy se sousední úseky překrývají, tedy neplatí N = KL). Každý segment je poté přenásoben
zvoleným typem okna w[l] shodné délky (je nutné použít některé z hladkých oken, například Hammingovo) a pomocí DFT transformován na dílčí modifikovaný periodogram Sˆk [m], jak popisuje vztah ¯ ¯2 L−1 ¯X ¯ j2πml ¯ 1 T N ¯ Sˆk [m] = x[l + kL] · w[l] · e− L ¯ , k ∈ h0, K − 1i, K = , (1) ¯ ¯ WL¯ L l=0
kde
L−1
W =
1X 2 w [l]. L l=0
(2)
Tak dostaneme celkem K dílčích modifikovaných periodogramů, jejichž zprůměrováním získáme požadovaný vyhlazený odhad: K−1 1 Xˆ ˆ Sx [m] = Sk [m]. K k=0
(3)
Pomocí průměrování dosáhneme K-násobného snížení celkového rozptylu odhadu ve srovnání s periodogramem počítaným z celé vstupní posloupnosti. Výkonová spektrální hustota spočtená podle vztahu (3) byla transformována do decibelové stupnice: SˆxdB [m] = 10 log Sˆx [m] [dB]. (4) Průběh dílčího periodogramu a vyhlazeného výkonového spektra jsou na obrázku 1 uvedeny v tomto logaritmickém měřítku. 2.2.
Intervaly spolehlivosti
Posloupnost vzorků EEG signálu x[n] lze považovat za normálně rozdělenou, proto také její Fourierova transformace (reálná i imaginární část) je normálně rozdělena, neboť se jedná o transformaci lineární [5]. Rozdělíme-li vstupní posloupnost x[n] na K úseků, sčítáme při průměrování K náhodných veličin s χ22 rozdělením. Vyhlazený odhad má proto χ2d rozdělení o d = 2K stupních volnosti. Interval spolehlivosti pro vyhlazený odhad výkonové spektrální hustoty je definován vztahem ( ) dSˆx [m] dSˆx [m] Pr ≤ Sx [m] ≤ 2 = 1 − α, (5) χ2d; 1− α χd; α 2
2
kde Pr{·} označuje pravděpodobnost jevu ve složené závorce, d = 2K je počet stupňů volnosti a α je hladina významnosti určující, s jakou pravděpodobností se skutečná hodnota Sx [m] bude vyskytovat mimo spočtený interval (obvyklé hodnoty α jsou 0,01, nebo 0,05). Symboly χ2d; α a χ2d; 1− α představují kritické hodnoty χ2 rozdělení s d stupni volnosti 2 2 na hladině významnosti α2 , které jsou tabelovány. Vztah (5) vyjadřuje interval spolehlivosti pro jednu frekvenční čáru“. Pro výkonové ” spektrum vyjádřené v decibelech (logaritmickém měřítku) podle rovnice (4) je interval pro všechny kmitočty shodný. Je popsán rovnicí [5] ) ( d d ≤ SxdB [m] ≤ SˆxdB [m] + 10 log 2 = 1 − α. (6) Pr SˆxdB [m] + 10 log 2 χd; 1− α χd; α 2
2
2.3.
Parametry výpočtu
Vlastnosti odhadu výkonové spektrální hustoty výrazně závisí na volbě parametrů jako je délka segmentu a typ váhovacího okna, překryv mezi segmenty či počet diskrétních frekvencí, pro které se PSD počítá. Hodnoty parametrů byly stanoveny tak, aby odhad PSD vyhověl několika protichůdným požadavkům: • vysoké frekvenční rozlišení PSD, • co nejmenší pronikání energie mezi sousedními spektrálními čarami“, ” • potlačení složek, jejichž frekvence není celočíselná (jako základní jednotku uvažujeme 1 Hz), • úzký interval spolehlivosti. Parametry použité při naší analýze shrnuje tabulka 1: Vstup sig se stim. bez stim.
nfft 1024 1024
Fs 1024 1024
Parametry win hamming(1024) hamming(1024)
over 0 0
p – 0,95
Výstupy spc conf f 1×513 – 1×513 1×513 2×513 1×513
Tabulka 1: Vstupní parametry a výstupní veličiny při výpočtu PSD Význam symbolů je následující: • nfft – počet diskrétních frekvencí, v nichž je výkonová spektrální hustota počítána. Číslo je vhodné volit jako mocninu dvojky, nejlépe shodně s délkou segmentu, aby při výpočtu FFT nedocházelo k jeho doplňování nulami. Pro reálný vstupní signál, jímž EEG záznam je, má pak výsledné spektrum délku nF2F T + 1. • Fs – vzorkovací frekvence vstupního signálu v hertzech. Na výpočet PSD nemá žádný vliv, slouží však k přepočtu normované frekvenční osy na kmitočty odpovídající skutečnému frekvenčnímu rozsahu signálu. • win – váhovací okno, jímž je segment signálu násoben. Důvodem je potlačení pronikání energie jedné spektrální složky do složek sousedních, které by snížilo kontrast“ ” mezi jednotlivými spektrálními čarami“. Váhování tento nežádoucí jev zmírní, ” ovšem za cenu snížení rozlišovací schopnosti (míra poklesu frekvenčního rozlišení závisí na zvoleném typu okna). • over – překryv mezi sousedními segmenty. Musí být menší než délka okna. Částečně eliminuje vliv použití váhovacích oken a zvyšuje počet průměrovaných dílčích modifikovaných periodogramů. Je-li nenulový, není možné přesně spočítat interval spolehlivosti. • p – představuje požadovanou pravděpodobnost, se kterou se skutečné hodnoty výkonové spektrální hustoty náhodného signálu mají vyskytovat ve spočteném intervalu spolehlivosti conf. Platí p = 1 − α, kde α je hladina významnosti (typicky α = 0,05 nebo α = 0,01).
3.
Výsledky
Výkonové spektrum je z hlediska vlivu míšní stimulace na aktivitu mozku nejdůležitějším ukazatelem. Podává mám informace o tom, do kterých elektrod se stimulací vybuzená aktivita šíří a kde se projevuje nejsilněji. Na frekvenčních průbězích vidíme, zda při SCS vznikají také vyšší harmonické složky a zda se skutečná hodnota stimulačního kmitočtu shoduje s deklarovanou, případně jak se liší. Tabulka 2 uvádí pro každého pacienta hodnotu použitého stimulačního kmitočtu, frekvenci zjištěnou ze spočtených výkonových spekter a číslo elektrody s maximálním rozdílem mezi PSD se stimulací a PSD bez stimulace na zjištěném kmitočtu. Číslo pacienta Stim. kmitočet deklarovaný [Hz] Stim. kmitočet zjištěný [Hz] Číslo elektrody
A01
A02
A03
A04
A05
A06
A07
A08
A09
A10
60
45
100
60
100
85
85
85
85
65
62
45
89
29
–
83
84
84
–
64
50
50
37
50
38
51
50
37
50
77
Tabulka 2: Deklarované a zjištěné hodnoty stimulačního kmitočtu a elektrody s nejvyšším výkonovým rozdílem na této frekvenci Nejvíce byl nárůst výkonu při SCS patrný na výkonových spektrech pacientů A01, A02, A07 a A08. Naopak žádný významný rozdíl se nepodařilo nalézt na PSD záznamů pořízených od pacientů A05 a A09, kde buď míšní stimulace korové oscilace nevybudila, nebo se tato aktivita v EEG záznamech neprojevila natolik, aby ji bylo možné považovat za statisticky významnou. Problematické bylo také hodnocení výkonového spektra u pacienta A03, neboť deklarovaná stimulační frekvence 100 Hz může interferovat s druhou harmonickou složkou síťového rušení. Totéž může nastat u již zmíněného pacienta A05. Na obrázku 2 jsou zachyceny výkonové spektrální hustoty prostorově filtrovaných EEG signálů z elektrod uvedených v posledním řádku tabulky 2 pro interval frekvencí 10–150 Hz a prvních pět pacientů. Modře je vykreslena výkonová spektrální hustota EEG bez přítomnosti stimulace, zeleně její interval spolehlivosti a červeně výkonové spektrum při zapnuté SCS. Deklarovaný stimulační kmitočet je naznačen černou tečkovanou čárou. Frekvenci, na níž došlo k maximálnímu nárůstu výkonu při SCS, označuje červená šipka. U pacientů A03, A06 a A10 byly jako významné zvoleny kmitočty blízké stimulačnímu, na nichž při SCS došlo ke zvýšení výkonu nad interval spolehlivosti PSD bez stimulace. Rozdíl PSD se stimulací a bez ní, třebaže je z tohoto pohledu statisticky významný, však v těchto případech dosahuje jen velmi malých hodnot. Nelze proto s jistotou tvrdit, že je důsledkem SCS, jelikož mohl vzniknout například celkovým vzájemným posunem spekter či vlivem šumu. Příčinou mohl být také vyšší rozptyl hodnot PSD se stimulací, které pak snadno přesáhnou interval spolehlivosti spočtený z PSD bez stimulace (jejíž rozptyl může být menší). 3.1.
Povrchové mapy PSD
Na obrázku 3 jsou povrchové mapy spočtené průměrováním map rozdílové výkonové spektrální hustoty všech deseti pacientů na zjištěném stimulačním kmitočtu fstim a nejbližších okolních frekvencích fstim − 1, fstim + 1. Jsou důkazem toho, že korové oscilace vznikající
Obrázek 2: Výkonová spektra s nejvyšším rozdílem výkonu na zjištěném stimulačním kmitočtu – pacienti A01 až A05
jako důsledek SCS jsou u jednotlivých pacientů lokalizovány shodně v oblasti pod elektrodami 37, 50, 51 a 63. Hodnoty rozdílu nacházející se v intervalu spolehlivosti, tedy statisticky nevýznamné, byly vynulovány. Jediným vážnějším nedostatkem PSD je neschopnost potlačit širokopásmové typy rušení, které zasahují všechny frekvence a výkonový rozdíl vyvolaný SCS (zvláště při vyšších stimulačních kmitočtech, které mají menší energii) v nich může snadno zaniknout.
Obrázek 3: Průměrné povrchové mapy rozdílových výkonových spekter pro zjištěný stimulační kmitočet a dvě nejbližší okolní frekvence – průměr všech pacientů
4.
Statistická významnost rozdílů mezi PSD se stimulací a bez ní ve skupině všech pacientů
K ověření statistické významnosti předpokladu, že míšní stimulace vyvolá u většiny pacientů na stimulačním kmitočtu zvýšenou elektrickou aktivitu některých částí mozkové kůry, byl použit párový t-test. Vstupní data jsou uvedena v tabulce 3. Číslo pacienta A01 A02 A03 A04 A05 A06 A07 A08 A09 A10
Frekvence Číslo fstim [Hz] elektrody 62 50 45 50 89 37 29 50 – 38 83 51 84 50 84 37 – 50 64 77
Výkonová spektrální hustota [dB] Ss∗ Sn∗ Ss∗ − Sn∗ 6,213519 18,169980 −6,691207 7,867134 −0,932345 0,749472 −1,291417 2,540490 −3,035791 1,050898
−4,001903 2,738629 −8,311514 4,006663 −2,290253 −1,280231 −4,073170 −2,944647 −2,664648 0,190503
10,215421 15,431351 1,620307 3,860471 1,357909 2,029703 2,781754 5,485138 −0,371143 0,860395
Tabulka 3: Maximální výkonové rozdíly na stimulačním kmitočtu pro jednotlivé pacienty Vysvětlivky k tabulce: Frekvence fstim je kmitočet, na kterém byl (na uvedené elektrodě) u každého z pacientů zjištěn největší výkonový rozdíl Ss∗ − Sn∗ mezi spektrálním výkonem se stimulací Ss∗ a bez stimulace Sn∗ .
Normální rozdělení hodnot náhodných výběrů Ss∗ a Sn∗ bylo ověřeno pomocí Lillieforsova testu normality na hladině významnosti 1 %. Nulovou hypotézou bylo tvrzení, že ve výkonech EEG při stimulaci a bez ní není rozdíl, tedy H0 : µ∆ = 0, alternativní hypotézou tvrzení, že při stimulaci výkon narůstá, tedy H1 : µ∆ > 0, kde µ∆ = µs −µn je rozdíl (skutečných) středních hodnot veličin reprezentovaných náhodnými výběry ve čtvrtém a pátem sloupci tabulky 3. Hladina významnosti testu byla zvolena α = 0,05, počet stupňů volnosti je d = 9 (počet pacientů minus jedna).
Výběrový průměr výkonového rozdílu přes všechny pacienty (hodnoty v posledním 2 sloupci tabulky 3) je M∆ = 4,32713 [dB]. Výběrový rozptyl těchto hodnot S∆ = 24,1596 [dB2 ]. Testová statistika má hodnotu T = 2,78391. Kritická hodnota kvantilu Studentova t-rozdělení je tc = 1,83311. Protože T > 1,83311, zamítáme nulovou hypotézu ve prospěch hypotézy alternativní na hladině významnosti 5 %. Lze tedy tvrdit, že nárůst výkonu EEG při stimulaci je v rámci celé skupiny pacientů účastnících se studie statisticky významný.
5.
Závěry
Dosažené výsledky a jejich praktický přínos lze shrnout do několika bodů: • Bylo zjištěno, že EEG signály nasnímané při SCS vykazují na stimulačním kmitočtu u osmi z deseti pacientů signifikantní nárůst výkonu na elektrodách v centrální mediální oblasti skalpu. • Statistická významnost výkonového rozdílu v rámci skupiny všech pacientů byla potvrzena pomocí t-testu na hladině významnosti 5 %. • Výsledky prokázaly, že vznik indukovaných rytmů v primární senzomotorické kůře odpovídající dolním končetinám má přímou souvislost s elektrickou míšní stimulací. • Získané poznatky lze využít k optimalizaci hodnoty stimulačního kmitočtu a k zpřesnění polohy pro implantaci stimulační elektrody tak, aby došlo k vybuzení oscilací pouze v těch centrech, která somatotopicky souvisí s částmi těla postiženými bolestí.
Poděkování Tento výzkum byl podporován z výzkumného záměru MŠMT MSM6840770012 Tran” sdisciplinární výzkum v biomedicínském inženýrství 2“ a z grantu GAČR 102/03/H085 Modelování biologických a řečových signálů“, IGA MZ ČR NR 8287-3/2005. ”
Reference [1] Bendat, J. S.; Piersol, A. G. Random Data: Analysis and Measurement Procedures. John Wiley & Sons, New York, 1971. [2] Hautvast, R. W.; et al. Relative changes in regional cerebral blood flow during spinal cord stimulation in patients with refractory angina pectoris. The European Journal of Neuroscience (1997), 1178–1183. [3] Kiriakopoulos, E. T.; et al. Functional magnetic resonance imaging: A potential tool for the evaluation of spinal cord stimulation: Technical case report. Neurosurgery (1997), 501–504. [4] Kozák, J.; et al. Metodické pokyny pro farmakoterapii akutní a chronické nenádorové bolesti. Bolest – časopis pro studium a léčbu bolesti (2004), 9–18. [5] Sovka, P.; Uhlíř, J. Číslicové zpracování signálů. Vydavatelství ČVUT, Praha, 2002.
