Podsumowanie wyników i wnioski

W Tabeli 5.2. przedstawione zostały wyniki analizy połączeń funkcjonalnych 52 eksperymentów, oznaczonych identyfikatorem składającym się z daty oraz numeru skrawka badanego danego dnia. Nieciągłości w numeracji skrawka tego samego dnia wynikają z odrzucenia niektórych eksperymentów, co było spowodowane bardzo niską aktywnością nie pozwalającą na przeprowadzenie analizy lub zupełnym brakiem aktywności.

Obok identyfikatora eksperymentu przedstawione zostały: liczba zidentyfikowanych neuronów, gęstość połączeń pobudzających, gęstość połączeń hamujących, gęstość połączeń synchronicznych w zakresie rytmu gamma oraz połączeń poprzez wspólne źródło. Zaznaczona została też informacja, czy próbka pochodziła z mózgu myszy zdrowej (WT, ang. Wild Type), myszy z uszkodzonym jednym chromosomem X (Heterozygoty, HZ, ang. Heterozygous) lub myszy z uszkodzonymi dwoma chromosomami (Homozygota, oznaczenie ZZ).

125

Tabela 5.2. Gęstość połączeń w badanych kulturach organotypowych

Identyfikator _neuronów ^Liczba połączeń ^Gęstość pobudzającyc h [%] Gęstość połączeń hamujących [%] Gęstość połączeń „gamma” [%] Gęstość wspólnie pobudzanych [%] Syndrom Łamliwego Chromosomu X 2011-06-24-1 74 4.7845 0.6574 2.6297 0.2739 WT 2011-06-26-2 79 1.6183 0.6890 9.5978 0.5448 WT 2011-06-27-0 65 1.7515 0.5680 14.1775 0.4734 WT 2011-06-27-1 98 1.2495 0.1978 0.9371 0.1354 WT 2011-06-28-0 66 2.4564 0.5969 1.6758 0.3673 WT 2011-06-28-1 172 1.9707 0.8315 0.3008 0.1961 WT 2011-06-29-0 123 3.4371 0.5883 0.2313 0.2578 WT 2011-06-30-0 86 3.7047 1.1087 1.5955 0.2569 WT 2011-07-01-0 135 5.5034 1.4211 1.3169 0.4664 WT 2012-04-15-0 106 0.5785 0.1335 0.0445 0.1513 ZZ 2012-04-16-0 88 0.9556 0.2583 0.1550 0.1679 ZZ 2012-04-16-1 106 1.8779 0.6052 0 0.2403 WT 2012-04-16-2 23 1.5123 0.1890 0 0.1890 WT 2012-04-16-3 130 1.0533 0.3254 1.7574 0.1716 ZZ 2012-04-17-0 98 0.7289 0.3540 1.0829 0.0625 ZZ 2012-04-17-1 175 0.7967 0.1110 0.0620 0.1176 WT 2012-04-17-2 63 1.4361 0.5543 2.7463 0.3527 HZ 2012-04-17-3 236 0.4130 0.1867 0.6105 0.1077 WT 2012-04-18-0 20 1.2500 0.2500 0 0 HZ 2012-04-18-1 89 0.9216 0.2020 0.1010 0.1641 WT 2012-04-18-2 256 0.4211 0.0778 0 0.0717 WT 2012-04-18-3 100 1.4500 0.4300 4.0700 0.2500 WT 2012-04-19-0 107 0.8560 0.3930 0 0.1485 ZZ 2012-04-19-1 48 0.8247 0.3472 0.1302 0.0868 HZ 2012-04-19-2 143 0.8753 0.3276 0.0685 0.1223 WT 2012-05-21-0 165 0.9183 0.1726 0 0.1947 WT 2012-05-21-1 167 1.2227 0.3908 1.8789 0.1829 WT 2012-05-21-2 210 1.8662 0.3832 0 0.1859 WT 2012-05-22-0 174 1.0503 0.3072 1.7473 0.1817 HZ 2012-05-22-1 289 0.7890 0.1281 0 0.1030 HZ 2012-05-22-2 242 0.4286 0.0649 0.3227 0.0871 ZZ 2012-05-22-3 204 0.8026 0.1754 0.4638 0.1562 WT 2012-05-23-0 151 1.1929 0.3903 1.0263 0.1228 WT 2012-05-23-1 259 0.9630 0.2743 0.0328 0.1729 HZ 2012-05-23-2 293 1.1485 0.3529 0 0.2050 WT 2012-05-23-3 149 1.0180 0.3784 0.4324 0.1982 HZ 2012-05-24-0 94 2.7048 0.7130 0.1698 0.3282 ZZ 2012-05-24-1 355 0.7903 0.0968 0 0.1182 WT 2012-05-24-2 213 1.0690 0.3130 0 0.2138 HZ 2012-05-25-0 161 1.9405 0.2315 0.4784 0.2122 WT 2012-05-25-1 134 2.3335 0.7073 0 0.3342 WT 2012-05-25-2 118 2.9015 0.8403 0.9983 0.6751 WT 2012-05-25-3 227 2.2104 0.4561 0.1087 0.2523 ZZ 2012-05-26-0 95 0.1994 0.0111 0.0665 0.0443 HZ 2012-05-27-0 235 1.6931 0.3821 0 0.2318 WT 2012-05-27-1 163 1.2571 0.4931 2.7212 0.2484 WT 2012-05-27-2 121 2.4657 0.6147 0 0.2800 ZZ 2012-05-28-0 177 1.9758 0.3671 0 0.2713 WT 2012-05-28-1 133 1.3059 0.4410 0.6445 0.2827 WT 2012-05-28-2 185 1.5281 0.3682 0 0.2396 WT 2012-05-29-0 75 2.5067 0.3200 2.3822 0.3733 WT 2012-05-29-2 312 0.2322 0.0380 0 0.0544 WT

126

Wyniki z powyższej tabeli przedstawione zostały na wykresach, aby sprawdzić istnienie korelacji pomiędzy gęstościami połączeń danego typu. Na Rysunku 5-42 przedstawiony został wykres gęstości połączeń hamujących w funkcji gęstości połączeń pobudzających w skali liniowej. Dopasowanie prostej regresji wskazuje, że wartość oczekiwana współczynnika zależności liniowej wynosi 0.244. Oznacza to, że gęstość połączeń hamujących jest w przybliżeniu 4 razy niższa niż gęstość połączeń pobudzających. Wykres gęstości połączeń pobudzających i hamujących pozwala sądzić, iż w badanych tkankach występuje korelacja pomiędzy gęstością połączeń pobudzających i hamujących ze współczynnikiem korelacji liniowej Pearsona R = 0.842.

Nie są widoczne różnice pomiędzy zbiorami punktów reprezentujących myszy zdrowe a hetero i homozygot y Syndromu Łamliwego Chromosomu X. Gdyby Syndrom wpływał na gęstość połączeń (lub stosunek gęstości połączeń hamujących do pobudzających), zbiory punktów tworzyłyby rozdzielone klastry. Zauważyć można jedynie brak występowania gęstości połączeń pobudzających powyżej 2.8% oraz połączeń hamujących powyżej 0.75%. Być może cechą odróżniającą osobniki chore od próby kontrolnej byłby średni stosunek połączeń gęstości połączeń hamujących do pobudzających (czyli nachylenie krzywej z odpowiednio wysoką istotnością statystyczną), jednakże do tego celu potrzebna byłaby większa statystyka zbadanych przypadków.

Rysunek 5-42:. Wykres gęstości połączeń hamujących w funkcji połączeń pobudzających.

Gęstość połączeń gamma nie wykazuje korelacji liczbowej z gęstościami połączeń hamujących lub pobudzających. Również kultury pochodzące od myszy z Syndromem Łamliwego Chromosomu X nie wykazują znaczących różnic w porównaniu do kultur pochodzących od myszy zdrowych przy tym kryterium porównawczym.

Na ostatnich dwóch wykresach przedstawione zostały zależności gęstości połączeń zsynchronizowanych w rytmie gamma od gęstości połączeń pobudzających (Rysunek 5-43) oraz hamujących (Rysunek 5-44). Ze względu na bardzo duże różnice gęstości połączeń zsynchronizowanych w rytmie gamma (prawie 3 rzędy wielkości),

127

wykresy przedstawione zostały w skali półlogarytmicznej. Przedstawione wyniki nie pozwalają sądzić, aby gęstość połączeń funkcjonalnych w kulturach organotypowych mogła służyć jako wskaźnik obecności Syndromu Łamliwego Chromosomu X.

Rysunek 5-43:Wykres gęstości połączeń gamma w funkcji gęstości połączeń pobudzających.

128

P

W rozprawie doktorskiej zaprezentowano wyniki prac nad rozwojem systemu do równoczesnej rejestracji i stymulacji aktywności elektrycznej kultur organotypowych mózgu myszy oraz sporządzania map połączeń funkcjonalnych pomiędzy neuronami. Prace te podzielić można na część projektową oraz badawczą.

W części projektowej celem było opracowanie modułów prototypowego układu scalonego, łączącego w sobie funkcjonalność odczytu oraz stymulacji elektrycznej komórek nerwowych. Scalenie tych dwóch funkcji stanowi postęp w stosunku do poprzednich generacji systemów rozwijanych w Zakładzie Oddziaływań i Detekcji Cząstek WFiIS AGH we współpracy z Instytutem Fizyki Cząstek Uniwersytetu Kalifornijskiego w Santa Cruz. Pojedynczy układ scalony powinien umożliwić skonstruowanie zarówno systemu do badań in vitro, podobnego do opisanego w rozdziale 2, zawierającego matrycę kilkuset elektrod z kilkoma (lub kilkunastoma) układami scalonymi, jak również miniaturowego systemu do badań in vivo, składającego się z mniejszej liczby (1-2) układów scalonych oraz matrycy lub wszczepianej sondy mikroelektrodowej.

Pierwszym zadaniem autora rozprawy było opracowanie układu elektroniki odczytu, który umożliwiałby rejestrację pełnego spektrum sygnałów neurobiologicznych, czyli niskoczęstotliwościowych potencjałów polowych, potencjałów czynnościowych (~1 kHz) artefaktów stymulacyjnych, oraz superpozycji wymienionych sygnałów. Zadanie to zostało zrealizowane w postaci prototypowego układu scalonego, w którym pojedynczy kanał zawiera przedwzmacniacz oraz dwa układy filtrów pasmowo-przepustowych ze wzmacniaczami wyjściowymi. Pierwszy z tych układów o niższym wzmocnieniu i paśmie częstotliwości w zakresie 1 Hz - 300 Hz odpowiedzialny jest za rejestrację potencjałów polowych, drugi o wyższym wzmocnieniu i paśmie częstotliwości 500 Hz – 2 kHz służy do rejestracji potencjałów czynnościowych. Wzmocnienia w obu układach dobrane są tak, aby wykorzystać pełny zakres dynamiczny multipleksera analogowego, który minimalizuje liczbę linii wyjściowych. Multiplekser analogowy opiera się na nowatorskiej koncepcji sterowania odczytem kolejnych kanałów, dzięki której problem przesłuchu sygnału z rejestru przesuwnego w środku okresu multipleksowania został całkowicie wyeliminowany. Do próbkowania wyjścia multipleksera wykorzystany może być pojedynczy, dostępny komercyjnie układ ADC o częstotliwości próbkowania 2.5 MHz.

Drugim zadaniem w części projektowej było stworzenie przetwornika cyfrowo-analogowego o rozdzielczości 9 bitów do generacji sygnałów stymulacyjnych. Przetwornik zrealizowany został na bazie koncepcji sukcesywnego dzielnika ładunku. Projekt cechuje się niskim poborem mocy przy zachowaniu bardzo dobrej liniowości różniczkowej i całkowej. Przetwornik wyposażony jest w układ próbkująco-pamiętający oraz układ redukcji przesłuchów sygnałów cyfrowych, co w połączeniu z niską nieliniowością sprawia, że doskonale nadaje się do roli generatora sygnałów stymulacyjnych.

Część badawcza, wykonywana głównie w Instytucie Fizyki Cząstek Uniwersytetu Kalifornijskiego w Santa Cruz, koncentrowała się na wykorzystaniu istniejącego systemu do rejestracji sygnałów nerwowych dla opracowania metody analizy danych pozwalającej na automatyczną detekcję i klasyfikację połączeń funkcjonalnych między neuronami.

129

Uczestnictwo w eksperymentach rejestracji aktywności spontanicznej kultur organotypowych oraz poznanie technik ich preparatyki miało wysoką wartość poznawczą i pozwoliło na całościowe spojrzenie na proces ekstrakcji informacji na temat połączeń funkcjonalnych w żywych sieciach neuronowych. W analizie danych eksperymentalnych punktem wyjściowym była metoda korelacji wzajemnej, stosowana do poszukiwania połączeń funkcjonalnych pomiędzy neuronami. Przeprowadzona została analiza funkcji korelacji wzajemnej szeregów czasowych reprezentujących potencjały czynnościowe pod kątem rozróżnienia typu połączenia na podstawie kształtu funkcji. Zaproponowana została wydajna obliczeniowo metoda detekcji połączeń synaptycznych, analogiczna do stosowanej powszechnie metody drgań losowych, ale pozwalająca na wykonanie tej samej analizy w czasie nawet kilkaset razy krótszym. Analiza falkowa została wykorzystana do wykrywania aktywności par neuronów zsynchronizowanych, generujących potencjały czynnościowe z częstotliwością w zakresie fal mózgowych gamma (30-70 Hz). Łącznie sporządzone zostały mapy połączeń (nie wszystkie zaprezentowane w niniejszej rozprawie) pobudzających, hamujących oraz połączeń zsynchronizowanych w rytmie gamma dla danych pochodzących z 52 eksperymentów. Część z badanych kultur pochodziła od myszy z kontrolowanym defektem genetycznym Syndromu Łamliwego Chromosomu X, jednak nie udało się zarejestrować istotnych statystycznie różnic w stosunku do grupy kontrolnej myszy zdrowych.

Analiza funkcji korelacji wzajemnej neuronów połączonych synapsą pobudzającą, pozwoliła na obserwację 4 maksimów na wykresach ilorazu stałej czasowej narastania do stałej czasowej opadania funkcji korelacji. Prawdopodobnie maksima te reprezentują 4 rodzaje synaps o różnych przekaźnikach chemicznych i różnej dynamice synaptycznej. Jako kontynuację tych obserwacji, warto by przeprowadzić serię eksperymentów, w którym kultury organotypowe poddano by działaniu selektywnego blokera chemicznego działającego na jeden rodzaj synaps – brak jednego z maksimów pozwoliłby na 100% potwierdzenie postawionej tezy.

Rozwijane systemy do rejestracji i stymulacji sygnałów elektrycznych w żywych sieciach neuronowych, zwiększanie gęstości mikroelektrod, w połączeniu z metodami cyfrowego przetwarzania sygnałów pozwalają na coraz bardziej precyzyjną rejestrację aktywności pojedynczych neuronów. Niemniej jednak, nawet tak zaawansowane jak obecnie istniejące systemy borykają się z problemami związanymi z sortowaniem potencjałów czynnościowych, co ma swoje negatywne reperkusje na etapie tworzenia map połączeń funkcjonalnych. Ostatecznym celem rozwoju technologii powinna być możliwość przeprowadzenia eksperymentu na wycinku mózgu, w którym można będzie zarejestrować aktywność absolutnie wszystkich neuronów, sporządzić mapę połączeń funkcjonalnych, uwzględnić komunikację pojedynczymi impulsami oraz seriami impulsów. Porównanie tych informacji z obrazem anatomicznym wysokiej rozdzielczości pozwoli uzyskać pełny obraz działania sieci neuronowej i być może doprowadzić do pełnego jej opisu i zrozumienia.

130

D

A:S

Bezpośrednim wynikiem przeprowadzenia eksperymentu rejestracji aktywności elektrycznej komórek nerwowych jest zbiór danych – zarejestrowanych przebiegów napięcia, próbkowanych na 512 elektrodach z częstotliwością 20kHz (okres 50μs). Każdy przetwornik analogowo-cyfrowy pracuje z rozdzielczością 12 bitów, więc można przyjąć, że do zapisania pojedynczej próbki napięcia na dysku potrzeba 2 bajtów. Proste oszacowanie pozwala obliczyć ilość danych zapisywanych na dysku twardym w ciągu jednej sekundy:

^0-1 gdzie nb – liczba bajtów na jedną próbkowaną wartość napięcia, nch – liczba kanałów w systemie, fs – częstotliwość próbkowania (na 1 s. na kanał).

Prędkość rejestracji 20.48 MB/s pomnożona przez czas trwania eksperymentu 6 godzin daje wynik 442.368 GB, czyli prawie pół terabajta danych. Taka ilość informacji wymaga oczywiście zastosowania automatycznych procedur, przede wszystkim we wstępnej fazie wykrywania potencjałów czynnościowych.

Pierwszym krokiem jest wykrycie potencjałów czynnościowych generowanych przez pojedyncze neurony, interesujących z biologicznego punktu widzenia. Poważnym problemem jest ograniczona rozdzielczość przestrzenna elektrod. W praktyce, na każdej elektrodzie może zostać zarejestrowana aktywność wielu neuronów, jak również każdy neuron może pobudzać kilka elektrod. Dlatego kluczowe jest rozróżnienie pojedynczych neuronów na podstawie charakterystycznego kształtu potencjału czynnościowego.

Pierwszym krokiem jest identyfikacja momentów czasowych, w których mogły wystąpić potencjały czynnościowe. Znajdywane są punkty, w których amplituda przekroczyła, a następnie opadła poniżej pewnego progu napięcia. Typowo wartość progowa wynosi 60μV.

Przy małej ilości danych wystarczyłoby ocenić wzrokowo, czy punkt, w którym przekroczona została amplituda stanowi potencjał czynnościowy, czy tylko przypadkowy, chwilowy wysoki poziom szumu za pomocą oceny rozłożenia punktów w przestrzeni amplituda ⨯ szerokość potencjału. Wystąpienie grupy punktów (klastra) w takiej przestrzeni oznaczałoby wystąpienie potencjałów o podobnym kształcie, których źródłem jest najprawdopodobniej jeden neuron. Analiza wzrokowa wykresów amplitudy w funkcji szerokości potencjału byłaby czasochłonna i subiektywna. Co więcej, zazwyczaj klastry nakładają się na siebie, ponieważ dane z jednej elektrody zawierają aktywność kilku

komórek. Rozwiązaniem problemu separacji jest użycie danych pochodzących z kilku elektrod jednocześnie.

Rysunek A-1:Geometria elektrod – elektroda centralna, na której wykryte zostało wystąpienie nadprogowej amplitudy (1) oraz elektrody sąsiadujące (2-7).

131

Analiza automatyczna wykonywana za pomocą oprogramowania rozwijanego od kilku lat przez zespół z Santa Cruz korzysta z danych zarejestrowanych na wybranej elektrodzie, którą nazywać będziemy elektrodą centralną oraz 6-ciu elektrod ją otaczających. Geometrię elektrod przedstawia Rysunek A-1.

Pierwszym etapem analizy automatycznej jest stworzenie 182-wymiarowego wektora z surowych danych dla każdego zidentyfikowanego przekroczenia progu napięcia. Wektor składa się 26 wartości napięcia z okresu od 0.5ms przed do 0.8ms po przekroczeniu wartości progowej; z elektrody centralnej oraz sześciu sąsiadów (26 · ) P ięć ie ięciu ż c ebie ów ięci e e e wi R u e A-2.

Rysunek A-2: Przebiegi napięć na elektrodzie centralnej (1) i sąsiadujacych (2-7). Pokazane zostało 50 przebiegów [42].

Drugim etapem jest Analiza Głównych Składowych stosowana w celu zmniejszenia liczby wymiarów poprzez znajdywanie N najbardziej znaczących zmiennych (zmiennych o najwyższej wariancji). Zmienne te są liniowymi kombinacjami 182 punktów pomiarowych. Na przykład, dla N=5, każdy potencjał czynnościowy jest reprezentowany jako punkt w N-wymiarowej przestrzeni. Rysunek A-3 przedstawia wykres rozproszenia dwóch najbardziej znaczących zmiennych dla danych przedstawionych na Rysunku A-2.

Rysunek A-3:Wykres rozproszenia dwóch najbardziej znaczących zmiennych, określony metodą PCA. Dane w przestrzeni 5-wymiarowej (nie tylko 2-wymiarowej, jak pokazano na rysunku) zostały dopasowane do sumy ważonej 5-wymiarowej funkcji Gaussa. 4 elipsy przedstawiają dopasowane funkcje Gaussa, a dokładniej ich 2-wymiarowe rzuty [42].

W trzecim etapie, w przebiegach napięć rzutowanych na przestrzeń o zredukowanej liczbie wymiarów równej N, poszukuje się klastrów, które interpretuje się jako zbiory potencjałów czynnościowych pochodzące od pojedynczego neuronu. Wyróżnić można 4 grupy punktów, które potencjalnie stanowią aktywność pojedynczych neuronów. Uśrednione przebiegi napięcia na 7 elektrodach przedstawione zostały na Rysunku A-4.

132

Niektóre z klastrów mogą nie odzwierciedlać aktywności pojedynczego neuronu, ale nałożenie się aktywności dwóch lub więcej komórek. Dlatego czwartym etapem jest użycie funkcji autokorelacji, aby wyeliminować „zanieczyszczone” klastry. Funkcja autokorelacji to histogram różnic czasów pomiędzy potencjałami czynnościowymi dla wszystkich par znalezionych potencjałów czynnościowych. Z podstaw elektrofizjologii wiadomo, że pomiędzy dwoma potencjałami czynnościowymi występuje minimalny odstęp czasowy (okres refrakcji) równy ok. 1.5ms, podczas którego neuron nie może wygenerować potencjału czynnościowego po raz drugi. Dlatego jeżeli pomiędzy pewnymi parami potencjałów czynnościowych interwał czasowy jest krótszy niż 1.5ms, uznaje się taki klaster za „zanieczyszczony”.

Ostatnim etapem jest znajdywanie duplikatów. Wiele spośród zidentyfikowanych neuronów jest tak naprawdę aktywnością jednego neuronu zarejestrowaną na kilku sąsiadujących elektrodach. Identyfikacja duplikatów odbywa się na zasadzie poszukiwania koincydencji pomiędzy momentami czasowymi, w których wykryto aktywność poszczególnych neuronów.

Rysunek A-4: Przebiegi w poszczególnych wierszach przedstawiają uśrednione wartości zarejestrowanego napięcia dla odpowiadających im 4 grup sklasyfikowanych na Rysunku A-3.

133

D

B:W

Transformata falkowa stała się popularną metodą analizy zmienności mocy sygnału w szeregach czasowych. Poprzez rozkład szeregu czasowego do przestrzeni czasu i częstotliwości można określić dominujące składowe częstotliwościowe oraz ich zmienność w czasie. Chociaż istotą transformaty falkowej jest, podobnie jak transformaty Fouriera, rozkład widma częstotliwościowego, przewaga transformaty falkowej polega na dostarczeniu informacji, kiedy, a nie tylko czy wystąpiła składowa o danej częstotliwości harmonicznej. Niestety większość publikacji i opracowań na temat transformaty falkowej prezentuje ją jako metodę generowania atrakcyjnych, kolorowych ilustracji pomijając element opisu ilościowego oraz istotności statystycznej [108]. Niemniej jednak jest to metoda, na której oprzeć można algorytmy wykrywania kształtów na wykresach [106] jak również bardziej złożonych obiektów graficznych [109]. Własność ta została wykorzystana przez autora pracy właśnie do wykrywania kształtów wykresów funkcji związanych z synchronizacją aktywacji par neuronów.

Transformacie falkowej można poddać szereg czasowy xn, z równymi interwałami czasowymi t, gdzie n = 0, 1, 2 …N-1. W przedstawianym przypadku szeregiem czasowym xn będą wartości funkcji korelacji wzajemnej w przedziale [-80ms; 80ms]. Ponieważ transformata falkowa opiera się na operacji konwolucji, zwykle wybiera się kształt falki podobny do kształtu, który ma być automatycznie rozpoznawany. Spośród najczęściej stosowanych falek, falka Morleta jest najbardziej podobna do poszukiwanego. Falka Morleta opisana jest wzorem (B-1):

( ) e e B-1

gdzie – tzw. bezwymiarowy parametr czasowy ( – czas, s jest arbitralnie w b ą c wi ą wie ścią e u óbkowania), i – jednostka urojona,

0 - bezwymiarowy parametr częstotliwości falki Morleta.

Aby funkcja mogła być zastosowana jako falka, jej średnia musi być równa 0 oraz musi ona posiadać reprezentację zarówno w dziedzinie czasu, jak i częstotliwości. Parametr 0

może przyjmować tylko takie wartości, dla których powyższy warunek będzie spełniony. Analizując wzór (B-1) jest jasne, że falka Morleta jest funkcją zespoloną, iloczynem fali płaskiej i funkcji Gaussa. Przebieg falki w dziedzinie bezwymiarowego współczynnika czasu t/s przedstawiony został na Rysunku B-1.

Rysunek B-1:Falka Morleta, a) wykres wartości zespolonych, b) rzut wartości rzeczywistych, i c) rzut wartości urojonych

Jak widać, falka ma kształt spirali o zmiennym promieniu. Rzut wartości rzeczywistych jest symetryczny względem t/s, odpowiada składowej rzeczywistej funkcji cos(t/s)

134

pomnożonej przez funkcję Gaussa, natomiast rzut wartości urojonych odpowiada składowej urojonej funkcji fali płaskiej (sin(t/s)) pomnożonej przez funkcję Gaussa.

Ciągła transformata falkowa szeregu czasowego zdefiniowana jest jako konwolucja (splot) szeregu czasowego xn z przeskalowaną falką o odpowiednim przesunięciu:

( ) ∑ [^{( )} ]

B-2

gdzie (*) oznacza liczbę zespoloną sprzężoną, sumowanie po n = 0, …N-1 oznacza, że splot obliczany jest po całym szeregu czasowym xn.

Poprzez zmianę skali falki „s” i przesuwanie po indeksie czasu n, można zbudować obraz pokazujący amplitudę danej cechy w funkcji skali „s” jak również zmienność tej cechy w czasie. Brak indeksu „0” przy falce oznacza, że Falka jest znormalizowana. Mimo że jest możliwe obliczenie transformaty Falkowej w dziedzinie zespolonej, bardziej wydajne obliczeniowo jest przejście od dziedziny czasu do dziedziny częstotliwości.

Aby przybliżyć transformatę ciągłą, należy konwolucję przeprowadzić N razy dla każdej skali, gdzie N oznacza liczbę punktów w szeregu czasowym xn. Twierdzenie o konwolucji pozwala wykonać wszystkie N konwolucji równolegle w przestrzeni Fouriera przy pomocy dyskretnej transformaty Fouriera. Dyskretna transformata Fouriera szeregu xn wyraża się wzorem:

̂

N^{∑ e}

B-3

gdzie k=0, 1, … N–1 oznacza numerację wartości częstotliwości, dla których obliczana jest transformata.

W ograniczonym przedziale, transformata Fouriera falki Morleta ( ) dana jest przez funkcję ̂( ), opisaną wzorem (B-4).

̂ _{( ) ( )e} ( ) B-4

gdzie – częstość kołowa, H( ) jest funkcją Haevisidea (skok jednostkowy) równa 1 dla >0 i 0 w pozostałym przedziale.

Korzystając z twierdzenia o konwolucji, transformata falkowa wyrazi się jako odwrotna transformata Fouriera z iloczynu szeregu i falki w dziedzinie częstotliwości:

( ) ∑ ̂ ̂( )e

B-5

Częstość kątowa k zdefiniowana jest następująco:

{ N N N N ^B-6

Wykorzystując powyższe równania można obliczyć transformatę falkową dla danego czynnika skalującego „s” we wszystkich N punktach w jednym kroku.

Aby zachować pewność, że transformata falkowa (B-6) dla każdego czynnika skalującego jest porównywalna oraz że jest porównywalna z innymi szeregami czasowymi poddawanymi transformacie, konieczna jest normalizacja funkcji Falkowej tak ,aby miała ona energię równą jedności. Równanie normalizacji ma postać:

135 ̂( ) (

^{) ̂ ( ) B-7}

Wynik obliczenia transformaty Falkowej (wzór (B-7) po odwrotnej transformacie Fouriera) jest szeregiem liczb zespolonych. Można zatem wyróżnić część rzeczywistą Re( ( )) , część urojoną m( ( )) , amplitudę | ( )| , fazę c [ m{ ( )} Re{ ( )}] oraz najbardziej istotną wielkość : moc widmową falki | ( )| .

Aby łatwo porównywać różne wyniki mocy widmowych falki, pożądane jest stosowanie wspólnej normalizacji dla zakresu transformaty Falkowej. Stosując normalizację zgodnie z (B-7) i transformatę Falkową (B-5), wartość oczekiwana mocy widmowej | ( )| jest równa iloczynowi N razy wartość oczekiwana szeregu | ̂ | . Dla szeregu, który składałby się tylko z wartości szumu białego, wartość oczekiwana wynosi N, gdzie oznacza wariancję, N oznacza liczbę punktów w szeregu czasowym (w tym przypadku liczbę punktów funkcji korelacji wzajemnej). Podzielenie mocy widmowej transformaty Falkowej przez pozwala uzyskać wartość mocy widmowej względem poziomu szumu, przy czym obliczona została jako wartość oczekiwana szeregu czasowego po transformacie Fouriera | ̂ | .

136

B

1. Schuetze, S. M. The discovery of the action potential. Trends in neuroscience. 1983, Vol. 6, pp. 164-168. 2. Office of the Home Secretary, National Academy of Sciences. Biographical Memoirs. Wyd. 1. Washington D.C. : National Academies Press, 1982. pp. 178-211. Vol. 53. ISBN 0-309-03287-3.

3. Hodgkin, A. L. and Huxley, A. F. A Quantitative Description of Membrane Current and It's Application to Conductance and Excitation. Journal of Physiology. 1952, 117, pp. 500-544.

4. Thomas, Jr CA, et al. A miniature microelectrode array to monitor the bioelectric activity of cultured cells. Experimental Cell Research. 1972, Vol. 74(1), pp. 61-66.

5. Henze, D. A, et al. On the Origin of the Extracellular Action Potential Waveform: A Modeling Study.

Journal of Neurophysiology. 2006, Vol. 95, pp. 3113-3128.

6. Silicon-Substrate Intracortical Microelectrode Arrays. Kipke, D. R., et al. 2, 2003, IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, Vol. 11, pp. 151-155.

7. Mind Creators. [Online] maj 2012. http://www.mindcreators.com /NeuronBasics.htm.

8. Traczyk, W. Fizjologia człowieka w zarysie. wyd. 7. Warszawa : Wydawnictwo Lekarskie PZWL, 2002. 9. Matthews, G. G. Neurobiologia. Od cząsteczek i komórek do układów. Wyd. 1. Warszawa : Wydawnictwo Lekarskie PZWL, 1998. ISBN 83-200-2452-8.

10. Shepherd, G. M. The Synaptic Organization of the Brain. Wyd. 5. New York : Oxford University Press, 2004. 0-19-515955-1.

11. Beaulieu, C. and Colonnier, M. Number of neurons in individual laminae of areas 3B, 4γ, and 6aα of