Univerza v Ljubljani Fakulteta za elektrotehniko Janez Jamšek Spektri višjih redov kardiovaskularnih signalov na osnovi transformacije z valčki Doktorska disertacija Mentorica: doc. dr. Aneta Stefanovska Ljubljana, 2005 University df Ljubljana faculty df electrical engineering Janez Jamšek High°order Spectra of Cardiouascular Signals Based on IDavelet Transf orm DlSSERTATIDN Ljubljana, 2005 University of Ljubljana Faculty of Electrical Engineering Janez Jamšek High°Order Spectra of Cardiouascular Signals Based on Ulauelet Transf orm DlSSERTATIDN Mentor: ddc. dr. Aneta Stefandvska Ljubljana, 2DD5 ACKNDVVLEDGEMENTS I would like to thank my mentor, doc. dr. Aneta Stefanovska, for showing me the way into the world of cardiovascular dvnamics, and for ali the invaluable expert knowledge, support, perfection convergence, science intuition, and not to mention, for aH the Saturdays, Sundavs, and holidays spent discussing cardiovascular matter. This work would not be the same vvithout the influence of Lancaster University, vvhere professor Peter V.E. McClintok, hosted me in his Nonlinear Group. I am grateful to him for ali of his Cambridge perfection that has influenced me. His supervision was, by ali means, most pleasurable. I have met many great scientists, but above ali, I have met great people. Mitya D.G. Luchinsky showed me that science can be a real pleasure, and that social life must not be neglected, but rather accomplished in synchronization. "Start always from zero, rather than think twice about it, and there is always a way around it", is Mitya's inheritance. To Igor A. Khovanov, whom knowing was productive and problem-goal oriented. To Andriy Bandrivskyy, with whom we improvised a relaxation oscillatory kind of šport events after serious studying hours. Many thanks for ali his help and valuable advice whenever needed, and to Stefano Beri, whom it was a pleasure to study and socialize. To Alan Bernjak and Bojan Musizza, co-workers in the Group of Nonlinear Dvnamics and Svnergetic. Tam grateful to Alenka Flander, for ali her advices and help just when I needed it. For ali of their careful reading, I would like to thank Jožica Gračarin, Vesna Habinc, Bojan Ilich and Marko Mlakar. At last, I am grateful to my parents, Bernardka and Miroslav Jamšek for ali of their kind support and patience during my long studies. Razširjen« povzetek Spektri višjih redov kardiovaskularnih signalov na osnovi transformacije z valčki Razširjeni povzetek Ritmi so med najbolj očitnimi lastnostmi živih sistemov. Pojavljajo se na vseh nivojih bioloških organizacij, od enoceličnih do večceličnih organizmov, s periodami od dela sekunde do leta. Tako srčna in dihalna funkcija kot circadianov ritem1 v spanju in pri zavesti do ključnih periodičnih procesov za ohranjanje človeškega življenja. Kljub veliki navezavi s fiziologijo periodični pojavi niso omejeni le na žive sisteme. Primer so kemične reakcije, ki so jih prvi odkrili Bray leta 1921, Belousov in Zabotinski leta 1959 oziroma leta 1964 in drugi. Kardiovaskularni sistem Kardiovaskularni sistem je eden od osnovnih sistemov človeškega organizma. Vsem celicam organizma neprestano dovaja energijo in snovi, ki so potrebne za njihovo normalno delovanje, hkrati pa iz celic odnaša snovi, ki nastanejo z metabolizmom. Sestavljata ga srce in ožilje (arterije, kapilare in vene). Pretok, ki je enak celotnemu volumnu krvi (tj. 4 1 - 6 1 oziroma 7 % - 8 % telesne teže), sklene pot pri sproščenem, zdravem človeku po ožilju povprečno v eni minuti [55]. Tako dinamiko kardiovaskularnega sistema preučujemo na časovni osi okoli ene minute. Srce ima vlogo črpalke, ki poganja kri po sklenjenem krogu elastičnih žil. Pljuča lahko gledamo kot generator pritiska [55]. Krvni pretok, pritisk ter aktivnost pljuč in srca določajo dinamiko sistema krvnega obtoka. Raziskave so pokazale [8, 9, 10, 110-112, 117], da izmerjeni signali krvnega pretoka vsebujejo deterministično dinamiko, kar pomeni, da je sistem krvnega obtoka izid končnega števila podsistemov (avtonomnih oscilatorjev), med katerimi ima vsak svojo značilno frekvenco. Pri regulaciji krvnega pretoka in pritiska sodeluje pet podsistemov: srčni, respiratorni, miogeni, nevrogeni in metabolični sistem. Vsi ti sistemi so tudi pri zdravih ljudeh v mirovanju med seboj rahlo sklopljeni, zato njihove značilne frekvence niso stalne, temveč se spreminjajo s časom, njihove amplitude pa so modulirane [6, 7, 110, 111]. Med posameznimi oscilatorji lahko nastanejo fazne sklopitve in sinhronizacija, ki se pokažejo v L. circa - okoli; dies = dan. 'S J obliki povezav med njihovimi frekvencami in fazami [9, 79, 94, 97, 113, 114]. Fazna sklopitev je torej pojav določenih relacij med fazami medsebojno delujočih sistemov, medtem ko ni nujno, da med amplitudami obstaja korelacija. Sklopitve omogočajo izmenjavo informacije med procesi in so tako temelj za pravilno delovanje sistema krvnega obtoka. Frekvenca in amplituda vsake opazovane oscilacije nam priča o aktivnosti oscilatorja in učinku vseh sklopitev. Dobro je raziskana frekvenčna modulacija srčnega ritma v ritmu dihanja, znana pod imenom respiratorna sinusna aritmija [18, 34]. Bivariatna analiza v časovnem prostoru, ki je bila pred kratkim razvita za analizo sinhronizacije ali posplošene sinhronizacije pri kaotičnih in šumnih oscilatorjih, je pokazala, da obstaja sinhronizacija srčnega in respiratomega oscilatorja [3, 23, 112, 117]. Slika 1: Rešitev enačbe (1) v faznem prostoru za osnovni oscilator. Stabilni limitni cikel je r{ = ax s fazo $ = 27t#, če je $ = 0 pri t = 0. Sklopitve med podsistemi sistema krvnega obtoka obstajajo, narava njihovega delovanja pa je še nepojasnjena. Model kardiovaskularnega sistema lahko tako predstavimo s sistemom enačb petih sklopljenih podsistemov, od katerih lahko vsak avtonomno oscilira [115, 116]. Osnovna enota je preprost oscilator z limitnim ciklom, ki ga je opisal Poincare [128] dxi -± = alx'](al-rl)-27flxl2, (1) dt dxl —t = oclx'2{al-rj) + 27flx[, dt j -2 2 kjer je rt = ^x[ +x'2 . Oscilator vsebuje strukturno stabilnost in robustnost, ki ju določata fiziološko razumevanje in analiza izmerjenih signalov. Spremenljivki stanja xll'mx'2 opisujeta pretok in hitrost II /-2 2 kjer je rt = ^x[ +x'2 . Oscilator vsebuje strukturno stabilnost in robustnost, ki ju določata fiziološko razumevanje in analiza izmerjenih signalov. Spremenljivki stanja xll'mx'2 opisujeta pretok in hitrost Razširjeni povzetek pretoka z-tega oscilatorja. Vsak oscilator je določen s frekvenco / in amplitudo a\, konstanta a\ pa določa hitrost, s katero se vektor stanj približuje limitnemu ciklu, slika 1. Ker podsistemov ne moremo obravnavati ločeno, je proučevanje njihove narave in lastnosti posameznih medsebojnih sklopitev oscilatorjev zelo težko. Zanima nas narava in pomen faznih zvez med posameznimi avtonomnimi oscilatorji, ki lahko v primeru nelinearne sklopitve povzročijo novo odvisno komponento pri modulirani frekvenci, ki je vsota osnovnih frekvenc sklopljenih oscilatorjev. Sklopitve so torej ključnega pomena za razumevanje kardiovaskularnega in morda tudi celotnega človeškega sistema. Biološki signali so običajno pomešani s šumom. Za izločitev informacije o faznih, frekvenčnih in amplitudnih sklopitvah uporabljamo zapletene transformacije, ki poudarijo njihovo vsebino. Bispektralna analiza spada v skupino tehnik, zasnovanih na statistiki višjih redov, ki se lahko uporabljajo za analiziranje neGaussovih signalov, za pridobivanje fazne informacije, zmanjševanje Gaussovega šuma neznanih spektralnih oblik ter odkrivanje in določevanje nelinearnosti signalov [60, 68, 69]. V sledečem besedilu razširjamo bispektralno analizo za določanje uporabnih lastnosti iz nestacionarnih podatkov in demonstriramo spremenjeno tehniko z uporabo na testnih signalih sklopljenih oscilatorjev. Bispekter vsebuje statistiko tretjega reda. Ocenjevanje spektrov je osnovano na konvencionalnem Fourierevem tipu direktnega pristopa z izračunom momentov tretjega reda, ki so v primeru statistike tretjega reda enaki kumulantom tretjega reda [60, 66-69]. Klasično oceno bispektra dobimo kot povprečje ocenjenih momentov tretjega reda za vsak segment, na katere je razdeljen signal, in sicer zaradi zagotovitve statistične stabilnosti [67]. Moment tretjega reda je enak trojnemu produktu diskretne Fouriereve transformiranke pri diskretnih frekvencah k, l in k+l[69] (2) Bispekter B(k, t) je kompleksor, določata ga amplituda^4 in faza (kj,n) = k(n) + ,(n)-k+l(n). (5) Če sta dve frekvenčni komponenti k in / frekvenčno in fazno sklopljeni, 10 Hz. Za vsak primer je prikazan signal za prvih 15 s. (b) Močnostni spekter, (c) Bispekter \B\ izračunan iz K= 33 segmentov, s 66 % prekrivanjem in z uporabo Blackmanovega okna za zmanjšanje razlivanja in (d) njegov nivojni prikaz. Del bispektra nad/2 > 1,0 Hz je odrezan, ker trojica (1,1 Hz, 1,1 Hz, 1,1 Hz) povzroči visok vrh, kije fizično nepomemben, (e) Bifaza (/) in (f) biamplituda A za bifrekvenco (1,1 Hz, 0,24 Hz) z 0,3 s časovnim korakom in 100 s dolgim oknom za oceno diskretne Fouriereve transformiranke in uporabo okna Blackman. Frekvenčna modulacija. Z bispektrom želimo ugotoviti tudi parametrično frekvenčno modulacijo in jo razlikovati od kvadratične. Parametrična modulacija povzroči nastanek frekvenčnih komponent pri vsoti in razliki karakteristične frekvence prvega oscilatorja in frekvence modulacije drugega l)ll oscilatorja. Obe komponenti bi bili lahko tudi posledica nelinearne kvadratične sklopitve. Enačbo prvega oscilatorja modela (7) ustrezno spremenimo (12) Bispekter se razlikuje od tistega v primeru kvadratične sklopitve. Pri bifrekvenci (/1,^2) dobimo vrh, čeprav komponente drugega oscilatorja^ (ena komponenta trojice) ni prisotne v močnostnem spektru signala, je vrednost različna od nič zaradi šuma. Bifaza ne kaže razdobja konstantne bifaze. V primeru močnejše modulacije je ta manj spremenljiva in ni pogostih faznih 2n skokov, kar je samo dodatna indikacija, da gre za primer modulacije. SRČND-RESPIRATDRNA SKLDPITEV V raziskavi je sodelovalo šest moških starih od 25 - 27 let, brez evidence o srčni bolezni. Pred začetkom meritve je vsak ležal sproščeno 15 minut. En nabor meritev je bil izmerjen v normalnem sproščenem stanju pri spontanem dihanju ter nadaljnja dva do trije nabori meritev pri različnem enakomernem dihanju (počasnejšem/hitrejšem od spontanega). Meritve so trajale 20 minut pri spontanem in 12 minut pri enakomernem dihanju. Merili smo krvni pretok na štirih različnih mestih s podobnimi lastnostmi krvožilnega sistema: na obeh rokah (levo in desno zapestje) in obeh nogah (levem in desnem gležnju) z vzorčno frekvenco 40 Hz. Istočasno smo merili tudi električno aktivnost srca (EKG) in dihanja z vzorčno frekvenco 400 Hz. Pri vsaki meritvi je tako nastala podatkovna datoteka, kije vsebuje 7 signalov. Tehnika zajema podatkov je opisana v [112]. Vse skupaj je bilo zajetih 22 podatkovnih datotek. Signale krvnega pretoka smo predhodno obdelali. Odstranili smo zelo nizke in zelo visoke frekvence z uporabo oken z drsečim povprečjem; dolžine 200 s za trend in 0,2 s za visoke frekvence, hkrati pa smo jih prevzorčili na 10 Hz. Tako smo se izognili problemom prekrivanja [81]. Signale smo normalizirali med nič in ena in jim odstranili srednjo vrednost. Določili smo karakteristično srčno/1 in dihalno/2 frekvenco ter komponente pri njunih harmoničnih pozicijah. Slika 3 levo prikazuje časovni potek tako obdelanih signalov za primer enakomernega dihanja, počasnejšega kot v primeru spontanega dihanja, in desno detekcijo frekvenčnih komponent v močnostnem spektru. Sledil je izračun normaliziranega bispektra kot povprečje preko več segmentov, ki smo jim vsakokrat odšteli povprečno vrednost signala. Za ugotavljanje nelinearne kvadratične sklopitve smo za vsak signal obdelali 8 vrhov, kot so našteti v preglednici 1. 0 Time (s) 521 °>n °'98 l>96 f(Hz) Slika 3: Signal krvnega pretoka b(t), merjen istočasno na štirih različnih mestih. Vsakemu je odstranjen trend in visoke frekvence ter je prevzorčen, normaliziran in osrediščen. Signali so dolgi 521 s in prevzorčeni na vzorčno frekvenco/s= 10 Hz. (a) Signal, izmerjen na desnem zapestju ba(t) in njegov močnostni spekter; (b) levo zapestje bh(t) in njegov močnostni spekter; (c) desni gleženj bc(f) in njegov močnostni spekter; (d) levi gleženj bd(t) in njegov močnostni spekter. Za vsak vrh smo izračunali bifazo in biamplitudo. Frekvenčno ločljivost smo nastavili tako, daje bila vsaj 1/10 najnižje dihalne frekvence. Uporabili smo okno, dolgo 100 s za izračun bispektra, biamplitude in bifaze. Uporabljeno okno določa tudi časovno ločljivost. Za izključitev ugotovitve naključnih sklopitev smo se osredotočili na tiste, ki so trajale vsaj 10 period počasnejše - dihalne frekvence, to je približno 100 s ali krajše. Ker je hkratna frekvenčna in časovna ločljivost izključujoča po Heisenbergovem principu nedoločenosti [43], je izbor možnosti omejen in potreben je kompromis. Okno smo premikali vzdolž časovne vrste s korakom 0,1 s = (1//Š). Kritično vrednost za oceno biamplitude smo postavili v vseh primerih na 2, to je dvakrat več kot je povprečna vrednost bispektra v tako imenovanem notranjem trikotniku bifrekvenčne domene. IX Da bi lahko sklepali na kvadratično sklopitev, smo določili potrebne pogoje: (i) Konstantna bifaza vsaj 10 period počasnejše sklopitvene komponente; (ii) Istočasno prisoten plato bifaze za vseh šest (osem) vrhov; (iii) V času sklopitve ni nobenih faznih skokov in bifazne spremembe so znotraj intervala n radianov (šum); (iv) Biamplituda mora biti nad določeno kritično vrednostjo. Preglednica 1: Vrhovi in pripadajoče bifrekvence v bispektru kot posledica nelinearne kvadratične sklopitve med dvema oscilatorjema s karakterističnima frekvencama/i in^. Vrh Bifrekvenca 1 Vi,fi) 2 (fl-fljl) 3 (fr/2,2f2) 4 ((/i,2/2) 5 (fufrfi) 6 ifi+fiJrfi) 7 ifiJi) 8 Vufi) (a) f2(Hz) ^(Hz) 0,87 0,98 1,09 f (Hz) Slika 4: (a) Bispekter \B\ signala ba, izračunan iz K = 33 segmentov, 87 % prekrivanjem in Blackmanovim oknom za zmanjševanje razlivanja, (b) Nivojni prikaz dela bispektra/i,/^ < 1,4 Hz, ki nas zanima. Primer tipičnega bispektra za celotno frekvenčno področje prikazuje slika 4 (a). Razvidnih je več vrhov. Področje našega interesa je srčno-respiratorna sklopitev, področje okoli bifrekvence {f],fj), ki je podrobnejše prikazano na sliki 4 (b). Pri nižjih frekvencah so razvidne sklopitve, ki so lahko posledica sklopitev med miogenim in nevrogenim oscilatorjem. Teh v tem delu ne obravnavamo. Natančna analiza pokaže, da se v bispektru nahajajo vsi vrhovi, našteti v preglednici 1. Izračun bifaz in biamplitud za vse vrhove pa razkrije, da so v določenem časovnem intervalu Tqc izpolnjeni vsi pogoji za nelinearno kvadratično sklopitev. Primer vrhov, biamplitud in bifaz za vrhova ena in dva X Razširjeni povzetek prikazuje slika 5. Povzetek analize vseh signalov, v katerih je bila ugotovljena nelinearna sklopitev, pa je podan v preglednici 2. Čeprav so bili signali izmerjeni na šestih osebah, so v preglednici 2 podatki samo za pet oseb. Tudi pri šesti osebi smo ugotovili nelinearne sklopitve, vendar niso izpolnjevale zahtevanega časa trajanja. Samo za en primer spontanega dihanja so bili izpolnjeni vsi pogoji za nastop nelinearne sklopitve. Pri spontanem dihanju so epizode sklopitev kratke in fazni preskoki pogostejši. Slika 5: Analiza krvnega pretok signala bJJ), izračunanega iz K = 33 segmentov, 87 % prekrivanjem, z 0,1 s časovnim korakom in 100 s dolgim oknom za izračun diskretne Fouriereve transformiranke z uporabo Blackmanovega okna za zmanjševanje razlivanja za vrhova (a) 1 in (b) 2; levi stolpec, bispekter |Z?ba| s pripadajočim nivojnim prikazom; sredinski, biamplituda Aha; in desni, bifaza $,a. Vprašanje, ali kardiovaskularni sistem vsebuje deterministično dinamiko, je bilo že predmet številnih raziskav [6, 8, 109, 110, 112]. Številni rezultati potrjujejo, daje sistem, ki regulira krvni pretok, determinističen. Ali so rezultati bispektralne analize posledica deterministične ali stohatične komponente v signalih krvnega pretoka, preverimo z uporabo surogatov [102, 106]. V ta namen uporabimo metodo surogatov naključne faze [44, 102, 106, 123, 124]. Tako dobljeni signali imajo podobne spektralne lastnosti kot originalni signali krvnega pretoka, to je enako povprečno vrednost, enako varianco, enako avtokorelacijsko funkcijo in posledično enak močnostni spekter s to razliko, da ni faznih povezav, oziroma so rezultat linearnega Gaussovega procesa. Z bispektralno analizo ne ugotovimo nelinearnih sklopitev v surogatih signalov krvnih pretokov. Zaključimo lahko, da so fazne informacije vsebovane v kardiovaskularnem signalu krvnega pretoka deterministične narave. Aktivnost srca se izraža v vsaki krvni žili in je prisotna tudi v mikrocirkulaciji kapilarnega omrežja. Periferni krvni pretok regulirata zunanji (centralni) in notranji (lokalni) mehanizem in mora tako XI odražati aktivnost obeh [6, 7, 9, 108, 110]. Signali krvnega pretoka odražajo centralne in lokalne mehanizme regulacije v kardiovaskularnem sistemu. Signali, zajeti na različnih, med seboj precej oddaljenih mestih, so lahko zelo podobni. Čeprav odražajo pretok v kapilarni mreži, vsi vsebujejo informacijo o prostorsko invariantni periodični aktivnosti centralno generiranih srčnih in dihalnih signalov. Jakost periodičnih komponent v perifernem krvnem pretoku se spreminja s premerom žil in gostote omrežja, to je z lokalno upornostjo pretoka. Preglednica 2: Nelinearna kvadratična sklopitev zaznana v signalih krvnega pretoka. Za vsako meritev so bili istočasno izmerjeni krvni pretoki na štirih različnih mestih, kanali a-d. Tqc je interval, v katerem z bispektralno analizo ugotovimo, da sta srčni oscilator^ in respiratorni oscilator^es lahko nelinearno sklopljena. Produkt Tqc x fKs določa trajanje sklopitve. Za časa Tqc smo izračunali najvišjo biamplitudo za vrh 1, kije našega osnovnega zanimanja, največjo spremembo bifaze A$ njeno povprečno vrednost (j) in standardno deviacijo c%. Oseba Dihanje Kanal fbr (Hz) /res (Hz) T * qc (s) ■* qc-*7res -^lmax (arb. units) A^ (rad) (rad) cr (rad) 1 enakomerno a 1,08 0,11 105,7 11,6 190 1,11 8,92 0,20 1 enakomerno d 1,00 0,23 56,8 13,1 62 0,92 10,93 0,29 1 enakomerno b 0,97 0,34 18,9 6,4 50 0,84 0,47 0,28 2 spontano a 1,16 0,14 82,0 11,5 352 1,87 32,68 0,47 2 enakomerno c 1,05 0,10 89,5 9,0 122 1,48 4,05 0,34 2 enakomerno a 0,98 0,11 95,6 10,5 383 1,47 3,22 0,42 3 enakomerno d 1,08 0,13 56,5 7,3 334 1,29 2,21 0,48 3 enakomerno c 1,10 0,26 52,4 13,6 52 0,46 4,96 0,10 4 enakomerno d 1,01 0,10 105,6 10,6 407 2,47 0,58 0,18 4 enakomerno d 0,99 0,11 95,6 10,5 219 2,19 -6,51 0,76 5 enakomerno d 1,20 0,10 57,5 5,8 1009 2,05 5,88 0,67 Čeprav so bili izmerjeni signali krvnega pretoka zajeti na različnih, med seboj oddaljenih mestih (kanal a-d), vsi odražajo enake karakteristične srčne in respiratorne frekvence. Z bispektralno analizo signalov krvnega pretoka smo za vsako meritev dobili enake rezultate za vse signale, istočasno izmerjene na različnih mestih (kanal a-d) in tako potrdili, da se informacija o faznem razmerju ohranja, kar je v skladu s predhodnimi raziskavami. Bispekter, definiran kot (2), je poseben primer križnega spektra, ko so vsi trije signali enaki. Imenujemo ga tudi avto bispekter. Poleg signalov krvnega pretoka smo istočasno merili tudi signal EKG, signal dihanja in signal krvnega pritiska. Z uporabo križnega bi spektra preverimo srčno-respiratorno sklopitev z uporabo bivariatnih podatkov. Križni bispekter definiramo kot [69] XII Razširjeni povzetek BXYY(kJ) = X(k)Y(l)Y\k + l), (13) kjer sta X in Y diskretni Fourierevi transformiranki dveh različnih signalov x{t) in y{t) pri diskretnih frekvencah k, l in k+l. Izračunali smo križne spektre Z?cebb (kjer c pomeni križni, e signal e{i) in b signal b(t)), za primer, ko je x(t) signal EKG e{t) \ny(f) signal krvnega pretoka ba(t). Signal EKG nam primarno govori o srčni električni aktivnosti. Fazo prve, srčne komponente f\, v trojici {f\,fi,f\ +/2) dobimo tako neposredno iz EKG signala, dihalno komponento^ in komponento pri harmonski poziciji f\ +fi pa iz signala krvnega pretoka. Na podoben način kot (13) definiramo še križni bispekter Bxxx(k, J) in ga izračunamo za dva različna primera: (i) BcbTh, kjer je x(t) signal krvnega pretoka ba(t) in y(f) signal dihanja r{t). Signal r{t) najbolj direktno opisuje aktivnost respiratornega oscilatorja (ii) 5cprp, kjer je x{t) signal krvnega pritiska p(t) \ny{t) signal dihanja r{t). Z izračunom križnih bispektrov ugotovimo, daje informacija o sklopitvi med srčnim in respiratornim oscilatorjem neodvisna od signala, oziroma da je prisotna tudi v drugih kardiovaskularnih signalih. Križne spektre smo izračunali tudi za primer surogatov signalov e{i), r{i) in p{f) z naključno fazo. V tem primeru nismo ugotovili faznih sklopitev. Nelinearna sklopitev ali linearna sklopitev močno nelinearnih oscilatorjev. Naša študija je zasnovana na predpostavki, da sta srčni in respiratorni proces opisana kot šibki nelinearni oscilator in da so sklopitve med njima šibke [116]. Prikladno seje vprašati, kaj se zgodi, če predpostavke niso izpolnjene. Odgovor na to vprašanje smo poiskali na dva različna načina, z analitično aproksimacijo in digitalno simulacijo. V analizi v prilogi B obravnavamo generiranje harmoniko v paroma sklop lj enih šibko nelinearnih oscilatorjev. Ta potrjuje, da pri šibki sklopitvi nastopijo dodatne harmonične komponente pri 2o>2, 2oj\, \ ± o>i, 2&>i + 26)2, 3&>i ± 0)2, ki jih lahko povežemo s kvadratično sklopitvijo. Če gre za zadosti nelinearen oscilator in zadosti močno sklopitev, se lahko v principu pojavijo te in ostale kombinatorne frekvence kot posledica efekta drugega reda tudi v primeru linearne sklopitve. Vendar pa nastop teh kombinatornih frekvenc sam po sebi ni zadosten za izpolnitev pogojev za nastop nelinearne sklopitve v bispektru. Za to pri šibko do srednje močnih sklopitvah lahko vedno zanesljivo določimo, da gre za nelinearno sklopitev. Ko so nelinearnosti posebno močne ne moremo vedno pričakovati, da bo bispektralni pristop razkril zanesljivo informacijo o naravi sklopitve. Analiza je dopolnjena z digitalno simulacijo, s katero ugotavljamo področje ekstremnih pogojev, kjer pričakujemo, da bo bispektralni pristop neuspešen. Za generičen model izberemo van der Polov oscilator z dodatno nelinearnostjo, vsiljen z drugim relaksacijskim van der Polovim oscilatorjem v smislu aditivne sklopitve z dodanim Gaussovim šumom. Analiziramo podroben nabor parametrov za primera: (i) ko sta dva oscilatorja močno nelinearna, vendar linearno sklopljena; in (ii) ko sta nelinearna in nelinearno sklopljena. V najbolj ekstremnem primeru zelo močne linearne sklopitve in zelo močne dodane nelinearnosti ne moremo več razlikovati med močno nelinearnostjo oscilatorjev in močno nelinearno sklopitvijo. V tem primeru bispektralna metoda odpove. Kljub temu pa je več argumentov, ki podpirajo domnevo, da sta srčni in dihalni podsistem šibka nelinearna oscilatorja, ki sta šibko sklopljena. (i) Pri spontanem dihanju zdravih ljudi se pojavljajo le občasne in kratke epizode sinhronizacije [10, 99-101], kar nakazuje na relativno šibke sklopitve. (ii) Sinusna respiratorna aritmija je šibka pri spontanem dihanju in le malo močnejša pri zelo nizkih dihalnih frekvencah [23], kar ponovno podpira šibko sklopljen opis. (iii) Sklopitve lahko včasih popolnoma izginejo, kot je to pri komi [112]. Brez sklopitev se dinamika drastično poenostavi s popolno odsotnostjo sinhronizacije in modulacije. Dejstvo, da kljub majhni amplitudni spremenljivosti zaradi notranjega šuma ni opažene nobene spremembe naravnih frekvenc, nakazuje, da so sami oscilatorji kvečjemu šibko nelinearni, (iv) Če bi bili oscilatorji močno nelinearni in močne sklopitve (linearne), bi opazili veliko kombinatornih komponent okoli srčne frekvence. Analizirano pretirano močno sklopitveno področje je tako irelevantno za srčno-respiratorno sklopitev, ki jo ugotavljamo v tem delu. Razmerje do sinhronizacije. Dejstvo da lahko notranje sklopitve med oscilatorjema privedejo do sinhronizacije kot tudi do modulacije, je imelo za posledico veliko študij faznega razmerja med srčnim in respiratornim oscilatorjem [10, 24, 42, 46, 52, 71, 92, 95, 97, 100, 101, 113, 114, 118]. Prav možnost sinhronizacije nas je motivirala, da smo razvili nova orodja za nadaljnje raziskovanje sklopitev med sistemi: smer, jakost in še posebej naravo sklopitev. Informacijo o sklopitvah lahko dobimo s pomočjo bivariatnih podatkov (signal dihanja in EKG signal), z uporabo nedavno razvitih metod za analizo sinhronizacije, ali s posplošeno sinhronizacijo med kaotičnimi in/ali šumnimi sistemi (glej [72] in reference, ki so tam navedene). Tu nas zanima, ali se sinhronizacija pojavi ali ne, v pogojih, ko jasno zaznamo sklopitev. Z uporabo sinhrograma ugotovimo obstoj frekvenčne modulacije, ne pa sinhronizacije v primeru enakomernega nizkofrekvenčnega dihanja. Bispektralna analiza podaja drugačno informacijo kot jo dobimo iz sinhrograma. Razmerje do sinhronizacije v širšem smislu podrobneje podajamo v naslednjem poglavju. Vsiljen oscilator. Z uporabo novih razvitih tehnik za analizo smeri sklopitve [72, 93, 94, 103] je bilo pokazano [73, 117], da sta srčni in respiratorni sistem sklopljena obojesmerno. Vendar pa je vpliv dihanja prevladujoč (vodilni sistem) pri vseh dihalnih frekvencah, spontanih ali enakomernih [73, HIU Razširjeni povzetek 117]. Sklopitev med srčnim in respiratornim oscilatorjem lahko tako vidimo kot enosmerno: respiratorni sistem vodi srčnega. Poseben primer je enakomerno dihanje. Čeprav je med enakomernim dihanjem dihalna frekvenca konstantna, se primer razlikuje od primera vsiljenega oscilatorja (kjer je srčni oscilator vsiljen in respiratorni vodilo). Primer ponazorimo z generičnim modelom skoraj periodičnega Poincare oscilatorja, vsiljenega s periodično šibko zunanjo silo. Srčno-respiratorna sklopitev je bolj kompleksna kot vzeti primer, ki ne more povzročiti frekvenčnih komponent, ki jih opazimo pri srčno-respiratorni sklopitvi. Eksperiment enakomernega dihanja lahko razumemo kot sistem dveh sklopljenih oscilatorjev, čeprav je frekvenca enega vsiljena in konstantna (respiratorni sistem) sklopitev med obema spontana. Razmerje med bispektri in sinhronizacijo Sinhronizacija je osnovni pojav v fiziki, ki ga je v začetku moderne dobe znanosti prvič odkril Huvgens [37]. V klasičnem smislu pomeni sinhronizacija nastavitev frekvenc oscilatorjev zaradi šibkih medsebojnih vplivov [2, 4, 30]. Ne obstaja enotna definicija za sinhronizacijo. Najbolj osnovni sta frekvenčna in fazna sinhronizacija. Ti dve definiciji sta bili generalizirani na pojav sinhronizacije dveh ali več oscilatorjev, ki so periodični, šumi ali kaotičnih oscilatorji [79, 90, 91, 96]. V najbolj enostavnem primeru dveh periodičnih oscilatorjev je fazna sinhronizacija definirana kot sklopitev faz [79] (14) kjer sta n in m celi števili, ki opisujeta sklopitveno razmerje, i fazi oscilatorjev in £nek začetni fazni premik. Enačba (14) v ožjem smislu velja samo za kvaziperiodične oscilatorje. Za periodične oscilatorje je pogoj za fazno sklopitev ekvivalenten pogoju za frekvenčno sklopitev nf\ = mf2, kjer sta /i in^2 karakteristični frekvenci oscilatorjev. Kadar opazujemo sinhronizacijo v prisotnosti šuma, sinhronizacijo kaotičnih sistemov ali oscilatorjev z moduliranimi lastnimi frekvencami, fazna in frekvenčna sklopitev nista več ekvivalentni [101]. Kadar je šum močan, lahko pride do faznih preskokov in vidi se le težnja k sinhronizaciji. Fazno sinhronizacijo lahko razumemo kot pojav vrha v porazdelitvi ciklične relativne faze (15) HU in si jo razlagamo kot obstoj preferenčne stabilne vrednosti fazne razlike fa med dvema oscilatorjema. V takih primerih ne moremo enoumno odgovoriti na vprašanje o sinhronizaciji sistema, nanjo lahko gledamo le v statističnem smislu. Pri kardiovaskulamem sistemu s časovno spremenljivimi karakterističnimi frekvencami se lahko fazna sinhronizacija pojavi, za frekvence pa ni nujno, da so povezane. Držimo se zapisa (14) za fazno sinhronizacijo, v besedilu uporabljamo skrajšano sinhronizacija. Že odprto vprašanje razmerja bispektrov do sinhronizacije podrobnejše obdelamo na primeru signalov podgan med splošno anestezijo. Signali so bili že obdelani z metodami za analizo sinhronizacije [63]. Iz signalov smo izbrali in analizirali z bispektri dva signala, za katera dobimo v sinhrogramu zelo jasno izraženo epizodo sinhronizacije. Za podgane med splošno anestezijo smo merili električno aktivnost srca (EKG), dihanje, EEG in temperaturo. Izmerjenih je bilo 21 podgan, težkih 250 g, večina jih je bilo samcev. Prvih 11 smo uporabili za testiranje in umerjanje merilnih naprav ter za določanje kvalitete signalov. Z meritvijo smo začeli 4-7 minut po vbrizgu anestetika in končali, ko so se podgane začele spontano gibati. Meritve so trajale -70 min, vzorčna frekvenca je bila 1000 Hz. Med meritvijo so podgane ležale na trebuhu v Faradavevi kletki [63, 64]. 6 f(Hz) 8 0 0,5 1 1,5 2 2,5 3 f(Hz) Slika 6: 10 s detrendiranega, prevzorčenega, normaliziranega in osrediščenega signala (a) EKG e(t) in (c) dihanja r(f) signala rat20 med splošno anestezijo, -72 min dolg, vzorčen s/s = 50 Hz in njuna močnostna spektra (b) in (d). Vrh pri 6 Hz v močnostnem spektru e{t) nastopi, ko se podgana začne zbujati iz anestezije, strm prehod v trenutni srčni frekvenci okoli 40 minute, slika 7 (a). Za obe podgani smo izračunali križne bispektre Bceve, kjer pomeni c križni, e signal e{i) in r signal r{f). Signale EKG in dihanja smo predhodno obdelali. Odstranili smo zelo nizke in zelo visoke frekvence z uporabo oken z drsečim povprečjem; dolžine 200 s za trend in 0,04 s za visoke frekvence, hkrati pa hui Razširjeni povzetek smo jih prevzorčili na 50 Hz. Tako smo se izognili problemom prekrivanja [81]. Signale smo normalizirali med nič in ena in jim odstranili srednjo vrednost. Primer tako obdelanih signalov EKG e(i) in dihanja r{t) in njunih močnostnih spektrov za podgano rat20 je prikazan na sliki 6. Za oba signala smo najprej izračunali trenutno srčno in respiratorno frekvenco ter njuno razmerje in sinhrogram, ki so za podgano rat20 prikazani na sliki 7. Sledil je izračun križnih bispektrov in za prve štiri vrhove iz preglednice 1 še izračun biamplitude in bifaze. 1,83 10 20 30 40 Time (min) 63,17 Slika 7: (a) Trenutna srčna (b) dihalna frekvenca in (c) njuno frekvenčno razmerje za podgano rat20 med splošno anestezijo, (d) Srčno-respiratorni sinhrogram za podgano rat20. Frekvenčno ločljivost smo nastavili na 1/20 najnižje respiratorne frekvence, ki je bila okoli 1 Hz (potrebno je 20 s dolgo okno). Opazovali smo sklopitve, ki so trajale vsaj 10 period nižje sklopitvene frekvence, 10 krat (\//2) — 10 s. Glede na Heisenbergov princip nedoločenosti [43] smo zadovoljili potrebo po frekvenčni ločljivosti in izbrali 20 s dolgo okno za izračun bispektra, biamplitude in bifaze. Okno smo premikali vzdolž časovne vrste s korakom 0,1 s. Kritično vrednost smo za vse primere postavili na 2, to je dvakratno vrednost povprečne vrednosti bispektra na področju notranjega trikotnika. Križni bispektri razkrijejo nastanek in trajanje sinhronizacije, kot je prikazano na sliki 8. Ugotovimo tudi pojav nelinearnih sklopitev. Relacijo do sinhronizacije dopolnimo še z generičnim modelom. Pojav sinhronizacije lahko spremlja modulacija. V poglavju „Analiza sklopitev" smo pokazali primer modulacije brez sinhronizacije. Iz sinhrograma težko ugotovimo prisotnost modulacije, možno je le v primeru, daje ta zelo močna. Z generičnim modelom pa ugotavljamo zmožnosti bispektralne analize v XUII primeru, da istočasno nastopita sinhronizacija in modulacija. Uporabimo skoraj periodični vsiljen Poincare oscilator, ki ga periodično vodi šibka zunanja sila z dodano frekvenčno modulacijo. Njuno frekvenčno razmerje namenoma vzamemo za celoštevilsko; tako dobimo s sinhrogramom sinhronizacijo tudi, ko med njima ni sklopitve. Slika 8: Bifaza (j) in biamplituda A za vrh 1 (4,3 Hz, 1,05 Hz), izračunani z 0,1 časovnim korakom in 20 s dolgim oknom za oceno diskretne Fourierove transformiranke in uporabo Blackmanovega okna. Pri močni sinhronizaciji sta uspešna oba, sinhrogramom in bispekter. V bispektru opazujemo vrh 1 in obdobje konstantne bifaze, če je biamplituda nad kritično vrednostjo. Šibka sinhronizacija je težko zaznavna s sihrogramom, v bispektru pa je bifaza manj konstantna z več faznimi preskoki. Če sinhronizacije ni, potem ni vodoravnih črt v sinhrogramu, v bispektru pa vrha 1. Sinhrogram nas lahko zavede v primeru, da sklopitve med sistemoma ni, je pa njuno frekvenčno razmerje konstantno. V primeru sinhronizacije se lahko pojavi tudi nelinearna sklopitev, to pa ne velja nujno tudi v obratni smeri. Med njima ni enostavne povezave. Samo prisotnost modulacije lahko ugotovimo z bispektrom, zavede pa nas lahko hkratna prisotnost modulacije in linearne sklopitve. V tem primeru je potrebno analizirati več vrhov. Pomagamo si lahko z opazovanjem poteka faz frekvenčnih komponent v trojici. Skupni nastop močne nelinearne sklopitve in močne modulacije ni mogoče razločiti. Bispektralna analiza je bolj občutljiva na sklopitve med sistemi kot sinhrogram. Ugotovimo lahko tudi nastanek sinhronizacije, vendar v splošnem ne podaja enake informacije kot sinhrogram. Spektri viSjih redov na osnovi transformacije Z VALČKI Fouriereva transformacija je zasnovana na predpostavki (a) periodičnosti signala in (b) neskončno dolgih časovnih vrst [57, 58]. Ker nobena od predpostavk ni strogo izpolnjena za izmerjene signale, je določitev posameznih frekvenc v sistemu, ki vsebuje močne sklopitve, zelo zahtevna. Še težje je v nizkofrekvenčnem področju, ki nas še posebej zanima, saj je karakteristične frekvence še težje ločiti. Princip nedoločenosti Fouriereve transformacije omejuje zmožnosti ločevanja harmoničnih Razširjeni povzetek komponent v frekvenčnem področju bispektra [20, 69]. To lahko povzroči težave pri ugotavljanju nelinearne kvadratične sklopitve, ko je frekvenčni par blizu skupaj. Da lahko zagotovimo dobro ločljivost nizkih frekvenc, potrebujemo daljše odseke za izračun diskretnih Fourierevih transformirank. To hkrati zniža število odsekov in poslabša oceno bispektralnih mer. Daljši signali pa vodijo v nestacionarnost in varianca postane še večja [69]. Pomagamo si lahko s transformacijo z valčki (TV), ki je v nekem smislu posplošena Fouriereva transformacija [43] z dodano časovno ločljivostjo na bolj osnoven način, kot je to dovoljeno s KČFT [81]. TV je že bila uspešno uporabljena za obdelavo kardiovaskularnih signalov [3, 5]. Posplošitev bispektra na TV lahko omogoči ugotavljanje začasnih sprememb v faznih sklopitvah ali kratkotrajne sklopitve. Prav tako pričakujemo, da bo uspešna pri ločevanju širokih in sovpadajočih vrhovih na račun povečane časovne ločljivosti. VT preslika signal iz časovnega prostora v prostor čas-skala. Signal g{t) razstavi na družino v splošnem neortogonalnih funkcij %±, kijih dobimo s premikanjem in skaliranjem osnovnega valčka if{u). Valčna transformiranka Wg(s, t), ki jo je prvi vpeljal Morlet [28], je definirana kot [43] (16) Osnovni valček mora biti omejen tako v časovnem kot v frekvenčnem prostoru. Za analizo kardiovaskularnih signalov je najbolj primeren Morletov valček, kije Gaussova funkcija, modulirana s sinusnim valom [6, 8]. Poenostavljena oblika v časovnem prostoru je (17) Razmerje med skalo s in frekvenco /je f = fQ/s, kjer s jo določimo število period sinusa v oknu. Frekvenčna ločljivost se tako spreminja s skalo (frekvenco) - pri nizkih je boljša kot pri visokih, medtem ko je časovna ločljivost boljša pri visokih kot pri nizkih frekvencah. Valčni bispekter WB definiramo analogno v skladu s Fourierevo definicijo bispektra kot (18) kjer velja XIX Valčni bispekter meri vrednost fazne sklopitve, ki nastopi med valčnimi transformirankami pri skalah Si, s2 in s signala g(t) na intervalu T, tako daje pravilo vsote skal (frekvenc) (19) izpolnjeno. Valčni bispekter W2?je kompleksor, določata ga amplituda^ in faza (f> (20) Za vsako biskalo (s\t s2) ga lahko predstavimo kot točko v kompleksnem prostoru 9l[WB(s\, s2)] proti 3[WZ?(si, s2)], kar določa vektor. Njegova amplituda (dolžina) je biamplituda, bifaza pa je določena s kotom med vektorjem in pozitivno realno osjo. Trenutno bifazo izračunamo analogno kot v primeru (5) (21) Če sta dve skalni komponenti S\ in s2 skalno in fazno sklopljeni, $, = fa + 0s2, potem velja, daje bifaza enaka 0 (2n) radianov. V našem primeru je fazna sklopitev manj stroga, ker so lahko odvisne skalne komponente fazno zamaknjene. Za fazno sklopitev smatramo, če je bifaza konstantna (ni nujno enaka 0 radianom) za vsaj nekaj period višje skalne komponente biskale (sh s2). Istočasno opazujemo trenutno biamplitudo, ki lahko podaja relativno jakost sklopitve (22) V skladu s Fourierevo definicijo križnega bispektra lahko definiramo analogno križni valčni bispekter kot (23) Da lahko zadostimo frekvenčnemu pogoju vsote (19) pri visokih frekvencah, je potrebna višja frekvenčna ločljivost. Dosežemo jo tako, da Morletov valček spremenimo XX Razširjeni povzetek (24) Faktor d določa eksponentno upadanje Morletovega valčka. Večji kot je počasneje upada Gaussova funkcija in boljša je frekvenčna ločljivost ob hkratnem zmanjšanju časovne ločljivosti. Ustrezno frekvenčno ločljivost pri visokih frekvencah lahko zagotovimo tudi s faktorjem am, s katerim dosežemo, da širina okna ne upada hiperbolično, temveč počasneje. Za višje frekvence lahko uporabimo tudi konstantno širino valčka. S konstanto cn dosežemo normalizacijo, da transformacija ohranja energijo (moč). Za analizo kardiovaskularnih signalov izberemof0 = 1, daje razmerje med frekvenco in skalo/= l/s, in lažje interpretacije valčnih bispektrov. Na sliki 9 (a) in (b) je prikazan valčni bispekter za testni primer nelinearne kvadratične sklopitve med dvema Poincare oscilatorjema z dodanim Gaussovim šumom, signal Xu>, slika 2 (a). Slika 9: Rezultati za kvadratično sklopitev v prisotnosti aditivnega Gaussovega šuma, testni signal xm, dobljen z valčnim bispektrom. (a) Valčni bispekter \WB\, izračunan iz K = 33 segmentov, 66 % prekrivanjem, Tm2 = 8 s, Ge3 = 0,001 in z uporabo Morletovega valčka s konstantno dolžino 7hF = 40 s za izračun visokih frekvenc in (b) njegov nivojski prikaz. Del valčnega bispektra nad f2> 1,0 Hz je odrezan, ker trojica (1,1 Hz, 1,1 Hz, 1,1 Hz) povzroči visok vrh, ki je fizično nepomemben, (c) Bifaza (f> in (d) biamplituda A za bifrekvenco (1,1 Hz, 0,24 Hz), z 0,1 s časovnim korakom. Oba bispektra, valčni in Fourierev, podajata enako informacijo o sklopitvi. Razlika je vidna v obliki vrhov, ki so v primeru valčnega bispektra širši. To je pričakovati, saj je frekvenčna ločljivost za visoke frekvence nižja kot pri bispektru na osnovi Fouriereve transformacije. Hkrati je opazna večja časovna ločljivost valčnega bispektra predvsem iz časovnega poteka bifaze, slika 9 (c) in (d). 2 Dolžina Morletvega osnovnega valčka, s =1 s. 3 Vrednost Gaussove funkcije na robu Morletovega valčka. XXI Primerjava Fdurierevega in valčnega bispektra Vpeljano metodo valčnega bispektra uporabimo na signalih krvnega pretoka, ki smo jih že analizirali z uporabo Fourierevega bispektra na osnovi KČFT. Parametri valčne bispektralne transformacije so nastavljeni na tipske vrednosti, ugotovljene na podlagi testnih signalov. Valčni bispekter za celotno frekvenčno področje signala krvnega pretoka Z>a je prikazan sliki 10 (a). Razvidnih je več vrhov. Področje našega interesa je srčno-respiratoraa sklopitev, področje okoli bifrekvence (/j, f2), ki je podrobnejše prikazano na sliki 10 (b). Natančna analiza pokaže, da se v bispektru nahajajo vsi vrhovi, našteti v preglednici 1. Izračun bifaz in biamplitud za vse vrhove pa razkrije, da so v časovnem intervalu od 77,1 s do 170,4 s izpolnjeni vsi pogoji za nelinearno kvadratično sklopitev, katere dolžina je Tqc = 93,3. Primer vrhov, biamplitud in bifaz za vrhova 1 in 2 prikazuje slika 11. Slika 10: (a) Valčni bispekter | WB\ signala ba, izračunan iz K = 33 segmentov, 87 % prekrivanjem, Tm = 8 s, Ge ~ 0,001 in z uporabo Morletovega valčka s konstantno dolžino TH¥ = 80 s za izračun visokih frekvenc, (b) Nivojni prikaz dela bispektra fuf2 < 1,4 Hz, ki nas zanima. Primerjava rezultatov, dobljenih z valčnim bispektrom, s tistimi, dobljenimi s Fourierevim bispektrom, razkrije, da z obema metodama ugotovimo enako pozicijo vrhov v bispektru. Vrhovi, dobjeni z valčnim bispektrom, so širši kot v primeru Fourierevega bispektra. Časovni poteki biamplitud so zelo podobni. Pojavlja se enako število izstopajočih vrhov v enakih časovnih trenutkih. Izračun križno-korelacijskega koeficienta za biamplitudi prvega vrha obeh metod je enak 0,95. Podobno velja za časovni potek bifaz. Opazimo lahko, da so spremembe bifaze veliko bolj izrazite v primeru valjčnega bispektra. Čeprav si po obliki časovni poteki bifaz niso tako podobni, kot velja za biamplitude, pa se pojavljajo platoji s konstantno bifazo v enakih časovnih trenutkih. Z obema metodama ugotovimo XXII Razširjeni povzetek prisotnost kvadratične sklopitve. V primeru valčnega bispektra je za 2,3 s krajša, kar je 2,4 odstotna razlika. f,(Hz) 77,1 Time (s) Time (s) Slika 11: Analiza krvnega pretok signala ba(t), izračunanega iz K = 33 segmentov, 87 % prekrivanjem, z 0,1 s časovnim korakom, Tm = 8 s, Ge = 0,001 in z uporabo Morletovega valčka s konstantno dolžino Tw = 80 s za izračun visokih frekvenc, za vrhova (a) 1 in (b) 2; levi stolpec, bispekter |5ba| s pripadajočim ni vojnim prikazom; sredinski, biamplituda Aha; in desni, bifaza $,a. Če povzamemo, dobimo z obema metodama zelo podobne rezultate. Iz širine vrhov v bispektru lahko sklepamo, da je frekvenčna ločljivost Fourierevega bispektra višja, medtem ko časovna ločljivost valčnega bispektra višja, saj bolje zazna hitre bifazne spremembe. Valčna metoda ne da opazno boljših rezultatov. To velja za primer analize srčno-respiratorno sklopitve, in sicr pri določenih pogojih, po katerih sklepamo na nelinearno sklopitev. V tem primeru je izbrano okno za izračun Fourierevega bispektra ravno tisto, ki zadovolji potrebni frekvenčni ločljivosti in hkratni časovni ločljivosti. Posledično ne dobimo znatnih razlik med rezultati obeh metod. Razliko si oglejmo na primeru generičnega testnega signala dveh kvadratično sklopljenih Poincare oscilatorjev z dodanim Gaussovim šumom, i = 1,2, kjer sta/i = 1,1 Hz in f2 = 0,24 Hz: (25) XXIII Slika 12: Rezultati, dobljeni s Fourierevim (c)-(f) in valčnim (g)-(h) bispektrom za primer prekinjajoče se kvadratične sklopitve dveh Poincare oscilatorjev z dodanim Gaussovim šumom, (a) Testni signal xlh spremenljivka X\ prvega oscilatorja s karakteristično frekvenco/] = 1,1 Hz. Karakteristična frekvenca drugega oscilatorja je f2 ~ 0,24 Hz. Oscilatorja sta sklopljena enosmerno kvadratično z dvema različnima jakostima sklopitve: rj2 = 0,0 (1); in 0,2 (2). Sklopitev (2) je prisotna vsakih 20 s in traja 20 s. Signal je dolg 1200 s in vzorčen z vzorčno frekvenco/ = 10 Hz. Za vsak primer je prikazan signal za prvih 15 s. (b) Močnostni spekter, (c) Bifaza ^ in (d) biamplituda A za bifrekvenco (1,1 Hz, 0,24 Hz) izračunani s 100 s dolgim oknom za oceno diskretne Fouriereve transformiranke. (e) Bifaza in (d) biamplituda A za bifrekvenco (1,4 Hz, 1,4 Hz), izračunana iz K = 33 segmentov, z 81 % prekrivanjem, Ge = 0.001, Tm = 8 s, 0,1 s časovnim korakom in z uporabo Morletovega valčka s konstantno dolžino 7HF = 80 s za izračun visokih frekvenc. Pri izračunu valčnih bispektrov za visoke frekvence smo uporabili okno dolžine 20 s, za oceno bifaz in biamplitud za bifrekvenco, kjer se pojavi najvišji vrh v valčnem bispektru, pa dolžine 80 s. S tem smo zagotovili visokofrekvenčno ločljivost za visoke frekvence, kjer je kompleksen EEG signal MUHI Razširjeni povzetek koncentriran, okoli 1 Hz. Frekvenco smo povečevali s korakom 0,02 Hz, da smo ohranili frekvenčno pravilo vsote (19), okno pa premikali vzdolž serije s korakom 0,1 s. Kritična vrednost biamplitude je bila 2. Opazovali smo samo sklopitve med valovi delta iz EEG signala podgane. V posameznih odsekih zaznamo z valčnim bispektrom fazne sklopitve. Primer bispektra, biamplitude in bifaze za odsek c signala EEG podgane rat20 so prikazani na sliki 14. Če upoštevamo potrebne pogoje za nastop nelinearne sklopitve, vendar samo za vrh samosklopitve, dobimo fazne samosklopitve, zbrane v preglednici 3. Preglednica 3: Fazne sklopitve ugotovljene v EEG signalu podgane rat20. Tpc je časovni interval, za katerega z valčno bispektralno analizo ugotovimo, daje so valovi delta iz EEG signala fazno samosklopljeni pri bifrekvenci (/j /i). Produkt rpc x/j pove, koliko period valov delta iz EEG signala je trajala fazna sklopitev. V času trajanja Tpc smo izračunali še največjo biamplitudo vrha Amax ter največjo spremembo bifaze A$ njeno povprečno vrednost (f> in standardno deviacijo c^. Odsek (min) h (min) (Hz) T -* pc (s) TPcxf\ A -rlmax (del. enot) (rad) (rad) (rad) a 2,67 3,55 1,1 49 57 753 1,16 3,24 0,23 a 5,87 6,33 1,1 28 33 494 0,24 3,59 0,05 a 7,45 8,17 1,1 43 50 1224 0,81 -8,14 0,29 a 8,22 10,83 1,1 157 183 1894 2,53 -9,62 0,69 a 12,38 13,03 1,1 39 45 535 0,58 -6,05 0,12 b 18,83 19,50 1,3 40 67 922 1,31 -1,77 0,34 b 23,83 24,31 1,3 29 48 3249 1,33 -39,44 0,30 C 35,00 36,17 1,4 70 141 2931 2,43 8,45 0,62 C 36,42 36,75 1,4 20 40 839 0,99 3,64 0,23 C 36,92 37,45 1,4 32 65 617 1,09 7,92 0,34 d 40,67 41,53 1,6 52 130 558 1,02 4,50 0,31 Največ faznih sklopitev se nahaja v odseku a in c. Te sklopitve so najmočnejše in trajajo najdaljši čas Tpc. Če združimo vse biamplitude posameznih odsekov, potem izstopata dva vrhova. V času, ko nastopita vrhova je bifaza znotraj n intervala. Čas nastopa in trajanje faznih sklopitev sta prikazana na sliki 15. Ta dva dogodka zaznamo tudi v sinhrogramu srčno-respiratorne sklopitve, slika 13 (c), ko se pojavi in ko izgine sinhronizacija 4:1. V vmesnem času sinhrogram ne pokaže nobene sinhronizacije, z valčnim bispektrom pa zaznamo kratkotrajne fazne sklopitve, ki ne morejo sinhronizirati sklopljenih oscilatorjev, saj se bifaza spreminja uniformno. XXIX H* i 2,33 7,45 23,83 35,00 Time (min) 62,67 Slika 15: Fazne sklopitve rpc, ki izstopajo po amplitudi biamplitude (1000 del. enot in več) v preglednici 3. Prikazan je čas njihovega nastopa in trajanje. Prvi in zadnji Tpc sovpada z nastopom in izginotjem fazne sinhronizacije med srčnim in respiratomim oscilatorjem, slika 13 (c). Na podoben način kot smo analizirali srčno-respiratorno sklopitev bi lahko analizirali tudi sklopitev med srčnim oscilatorjem in valovi delta iz EEG signala in respiratomim oscilatorjem in valovi delta iz EEG signala, kar pa ni predmet tega dela. Valčni bispekter se je izkazal kot obetavno orodje za analizo EEG signalov med anestezijo. M Osnovna motivacija dela je bil razvoj orodja za proučevanje medsebojnega vpliva podsistemov kardiovaskulamega sistema. Z orodjem, ki bi bilo občutljivo na časovno spremenljivo vsebino, želimo povečati možnosti razkrivanja narave in pomena sklopitev. Sklopitve omogočajo izmenjavo informacije med procesi in so tako temelj za pravilno delovanje sistema krvnega obtoka. Frekvenca in amplituda vsake opazovane oscilacije nam pričata o aktivnosti oscilatorja in učinku vseh sklopitev. Sklopitve so torej ključnega pomena za razumevanje kardiovaskulamega in morda tudi celotnega človeškega sistema. Bolezensko stanje vodi v fiziološke spremembe, ki se odražajo v spremembah dinamičnih lastnosti kardiovaskulamega sistema. Poznavanje sklopitev pomeni možnost zaznavanja bolezenskega stanja v njegovem začetnem stadiju, ko je potreben minimalni zdravniški poseg. Razvito časovno občutljivo bispektralno metodo omogoča določitev narave sklopitev med sklopljenimi oscilatorjimi. Prednosti so: (i) možnost hkratnega opazovanja celotnega frekvenčnega področja; (ii) ugotavljanje sklopitev dveh ali več medsebojno sklopljenih oscilatorjev; (iii) določanje jakosti sklopitve; (iv) določanje narave sklopitve: linearna, nelinearna kvadratična ali parametrična v eni od frekvenc; (v) metoda je primerna za analizo šumnih signalov. Učinek sklopitev med srčnim in respiratomim oscilatorjem je prej epizodičen kot stalen in nespremenljiv. Frekvenčne in fazne sklopitve se izmenjujejo. Nelinearne sklopitve obstajajo tako med spontanim kot enakomernim dihanjem. Sklopitve med oscilatorji so šibke. Bispektralna in križna bispektralna analiza sta pokazali, daje informacija o sklopitvi med srčnim in respiratomim procesom lastna procesoma in prostorsko invariantna. Oba procesa sta centralnega izvora. Njuno fazno razmerje lahko opazujemo v signalu EKG, signalu krvnega pretoka in signalu krvnega pritiska, izmerjenih na različnih med seboj oddaljenih mestih. Nelinearna narava sklopitev med srčnim in respiratomim je vrojena in postane bolj izrazita, ko je frekvenca dihanja konstantna. Bispektralna analiza je torej sposobna določevati frekvenčne in fazne sklopitve med C V procesi in s tem primerna za razumevanje kardiovaskulamih signalov. V primerjavi s sinhrogramom je bolj občutljiva na sklopitve in manj občutljiva na šum. Zazna fazne sklopitve, vendar na drugačen način kot jih podaja sinhrogram. Med sklopitveno informacijo, dobljeno s sinhrogramom in z bispektrom, ne moremo določiti enostavne povezave. XXXI Valčna bispektralna metoda signalov krvnega pretoka ne da opazno boljših rezultatov kot Fourierova bispektralna metoda. To velja pri analizi srčno-respiratorno sklopitve pri določenih pogojih, po katerih sklepamo na nelinearno sklopitev. Za obdelavo signalov sta primernejša valčni in križno valčni bispekter kot bispekter na osnovi Fouriereve transformacije. Omogočata ugotavljanje kratkočasovnih faznih sklopitev, optimalno časovno in frekvenčno ločljivost, enostavno povezavo med skalo in frekvenco, neposredno interpretacijo, normalizacijo na energijo (moč), manjšo statistično napako, poljuben frekvenčni korak in sta računsko manj potratna. Iz dobljenih rezultatov valčne bispektralne analize signala EEG podgane med anestezijo lahko sklepamo, da so srčni in respiratorni kardiovaskularni oscilator in valovi delta iz EEG signala med anestezijo sklopljeni. Anestezija lahko vpliva na valove delta iz EEG signala tako, da ti vodijo srčni in respiratorni sistem v sinhronizacijo. Valčni bispekter lahko omogoči povezavo med teoretičnim kardiovaskulamim modelom in eksperimentalnimi meritvami. XXXII Izvirni prispevki k znanosti IZVIRNI PRISPEVKI K ZNANOSTI Pomembnejši izvirni prispevki k znanosti: 1. Ugotavljanje sklopljenih nelinearnih oscilatorjev z uporabo časovno občutljivih bispektralnih cenilk za biamplitudo in bifazo. Vpeljali smo časovno občutljivi bispektralni cenilki - biamplitudo (6) in bifazo (5) - za razkrivanje faznih sklopitev v univariatnih podatkih (poglavje 3, strani od 13 do 16). Pokazali smo, da je vpeljana metoda primerna za proučevanje sklopljenih nelinearnih oscilatorjev. Sposobna je meriti jakost fazne sklopitve z bispektralno cenilko - biamplitudo, katere vrednost je proporcionalna vrednosti koeficienta medsebojne sklopitve e (2.6) sklopljenih nelinearnih oscilatorjev (poglavje 2, strani od 7 do 12) in razkriti naravo sklopitve, t.j., ali je sklopitev dodana linearna, kvadratična ali parametrična v eni od frekvenc (poglavje 4, strani od 17 do 28). 2. Potrditev hipoteze o sklopitvi med srčnim in respiratornim oscilatorjem v človeškem kardiovaskularnem sistemu. Prvi smo uporabili vpeljano časovno bispektralno metodo na signalih krvnega pretoka za proučevanje narave srčno-respiratorne sklopitve. Kljub omejeni možnosti za razkrivanje metode med linearno in nelinearno sklopitvijo v ekstremnih razmerah je metoda uporabna, dokler sklopitve ne postanejo preveč zapletene, upoštevajoč fiziološko poznavanje kardiovaskularnega sistema (poglavje 5, strani od 29 do 56). Sklopitve med srčnim in respiratornim oscilatorjem so epizodične, prej kot stalne in nespremenljive. Frekvenčne in fazne sklopitve se izmenjujejo. Nelinearne sklopitve obstajajo tako med spontanim kot enakomernim dihanjem. Izrazitejše so pri enakomernem dihanju. Sklopitve med oscilatorji so šibke (poglavji 5 in 6, strani od 28 do 37 in od 57 do 76 ). Informacija o sklopitvi med srčnim in respiratornim procesom je lastna procesoma in prostorsko, znotraj CVS, invariantna. Oba procesa sta centralnega izvora. Njuno fazno razmerje lahko opazujemo v signalu EKG, signalu krvnega pretoka in signalu krvnega pritiska, izmerjenih na različnih med seboj oddaljenih mestih (poglavje 5, strani od 37 do 46). 3. Posplošitev bispektralnih cenilk na transformacijo z valčki Bispektralni cenilki - biamplitudo (22) in bifazo (21) - smo posplošili na transformacijo z valčki. Metoda je primernejša za ugotavljanje kratkočasovnih sklopitev. Omogoča za nas optimalno časovno in frekvenčno ločljivost. Valčni bispekter lahko omogoči povezavo med mm PRISPEVKI K ZNANOSTI teoretičnim kardiovaskularnim modelom in eksperimentalnimi meritvami (poglavji 7 in 8, strani od 77 do 106). 4. Ugotovitev sklopitve med srčnim in respiratornim sistemom ter delta valovi signala EEG z uporabo posplošenih bispektralnih značilk. Delta valovi iz EEG signala podgane med splošno anestezijo razkrivajo fiziološko razmerje med srčnim in respiratornim sistemom in delta valovi iz EEG signala. Srčno-respiratorna sinhronizacija je lahko posledica vodenja delta valov iz EEG signala podgane med splošno anestezijo (poglavje 9, strani od 107 do 114). Nekateri deli disertacije so objavljeni v naslednjih člankih: 1. J. Jamšek, A. Stefanovska, P.V.E. McClintock and I.A. Khovanov, Time-phase bispectral analvsis, Phys. Rev. E 68, 016201 (2003). 2. J. Jamšek, A. Stefanovska and P.V.E. McClintock, Nonlinear cardio-respiratorv interactions revealed by time-phase bispectral analvsis, Phvsics in Medicine and Biologv 49, 4407 (2004). Nekateri deli disertacije so bili predstavljeni na sledečih znanstvenih srečanjih: 1. J. Jamšek and A. Stefanovska, Bispectral analvsis of cardiovascular signals, Nonlinear Seminar, Department of nonlinear phvsics, Lancaster Universitv, United Kingdom (7.2. 2002). 2. J. Jamšek, A. Stefanovska and P.V.E. McClintock, Time-phase bispectral analvsis, basic theory and applications, Nonlinear Seminar, Department of nonlinear phvsics, Lancaster Universitv, United Kingdom (12.5. 2002). 3. J. Jamšek, A. Stefanovska and P.V.E. McClintock, Cardiovascular Svstem, Cardiovascular svstem, time-phase bispectral analysis, basic theory and application, 2nd Slovenia-Japan Seminar, Center for Applied Mathematics and Theoretical Physics University of Maribor, Slovenia (28.5.-5.6 2003). 4. J. Jamšek and A. Stefanovska, Quadratic cardio-respiratory coupling?, INTAS international Workshop, Department of Phvsics, University of Pisa, Italy (22-24.4. 2003). 5. J. Jamšek, A. Stefanovska and P.V.E. McClintock, Nonlinear cardio-respiratorv interaction, INTAS-ESF international Workshop, Ljubljana, Slovenia (10-13.11. 2003). 7 Bispectral analysis, a technique based on high-order statistics, is extended to encompass tirne dependence for the čase of coupled noisy nonlinear oscillators. It is applicable to univariate, as well as to multivariate, data obtained respectively from one or more of the oscillators. It is demonstrated for a generic model of interacting systems, whose basic units are Poincare oscillators. Their frequency and phase relationships are explored for different coupling strengths, both with and without Gaussian noise. The distinctions between additive linear, quadratic, and parametric (frequency modulated) interactions in presence of noise are illustrated. Bispectral analysis has been used to study the nature of the coupling between cardiac and respiratory activity. Univariate blood flow signals recorded simultaneously on both legs and arms were analysed. Coupling between cardiac and respiratory activity was also checked by use of bivariate data and computation of the cross-bispectrum betvveen ECG and respiratory signals and surrogate data of blood flow signals. Measurements were made on six healthy males, aged 25-27 years, during spontaneous breathing and during paced respiration, at frequencies, both lower and higher than that of spontaneous respiration. It was confirmed that the dvnamics of blood flow can be usefully considered in terms of coupled oscillators, and demonstrated that interactions betvveen the cardiac and respiratory processes are weak and time-varying, and that they can be nonlinear. Nonlinear coupling was revealed to exist during both spontaneous and paced respiration. Relation of bispectral analysis to svnchronization is outlined in the example of cardiovascular ECG and respiration signals of rats undergoing anaesthesia. Bispectrum proves to be more sensitive to interactions than the svnchrogram. It detects the phase svnchronization, and nevertheless, yields different information from that which can be resolved from a svnchrogram. Wavelet transform was incorporated into bispectrum and adopted to analyse cardiovascular signals using a Morlet wavelet as a mother wavelet. A tirne dependant biphase/biamplitude estimate, with higher frequency resolution at low frequencies, and higher tirne resolution at higher frequencies, was obtained. Its advantages, compared to Fourier based bispectrum, are discussed and demonstrated for a generic model of interacting systems, whose basic units are Poincare oscillators, in the application of CV blood flow signals, and in the application of EEG signal of rat undergoing anaesthesia. The wavelet bispectrum may provide a link between theoretical CVS models and experimental measurements. 1INTRODUCTION..............................................................................................................................1 2 CARDIOVASCULAR SYSTEM......................................................................................................5 2.1 BACKGROUND................................................................................................................................6 2.2 COUPLED OSCILLATORS.................................................................................................................7 3METHOD..........................................................................................................................................13 3.1 BISPECTRAL ANALYSIS................................................................................................................ 13 3.2 TIME-PHASE BISPECTRAL ANALYSIS............................................................................................ 14 3.3NORMALIZATION.......................................................................................................................... 15 4 ANALYSIS OF COUPLINGS........................................................................................................17 4.1 LlNEAR COUPLINGS...................................................................................................................... 18 4.2 LlNEAR COUPLINGS IN PRESENCE OF NOISE.................................................................................21 4.3 QUADRATIC COUPLINGS...............................................................................................................24 4.4 QUADRATIC COUPLINGS IN THE PRESENCE OF NOISE...................................................................26 4.5 FREQUENCY MODULATION IN THE PRESENCE OF NOISE...............................................................27 5 CARDIO-RESPIRATORY INTERACTIONS..............................................................................29 5.1 Data acojjisition......................................................................................................................29 5.2 Measurements..........................................................................................................................30 5.3Dataanalysis...........................................................................................................................32 5.4RESULTS.......................................................................................................................................33 5.5 Surrogates................................................................................................................................37 5.6 Global couplings.....................................................................................................................40 5.7cross-bispectrum.....................................................................................................................43 5.8discussion..................................................................................................................................45 5.8.1 Definition ofthe phase..........................................................................................................46 5.8.2 Nonlinear coupling, or linear coupling ofstrongly nonlinear oscillators?..........................47 5.8.3 Relationship to synchronization............................................................................................51 5.8.4 Synchronization, modulation andtype of coupling...............................................................52 5.8.5 Unidirectional or bidirectional coupling..............................................................................53 5.8.5.1 Forced oscillator.............................................................................................................53 6 BISPECTRAL RELATION TO SYNCHRONIZATION............................................................57 6.1 SYNCHRONIZATION DEFINITION...................................................................................................57 6.2 MEASUREMENTS..........................................................................................................................61 mm 6.3 Data analysis...........................................................................................................................62 6.4RESULTS.......................................................................................................................................62 6.4.1 Ratl6.....................................................................................................................................63 6.4.2 Rat20.....................................................................................................................................66 6.5 DISCUSSION..................................................................................................................................69 6.5.1 Synchronization and modulation...........................................................................................71 6.6 CONCLUSIONS..............................................................................................................................74 7 HIGH-ORDER SPECTRA BASED ON WAVELET TRANSFORM.........................................77 7.1 WAVELET TRANSFORM................................................................................................................77 7.1.1 Discretization........................................................................................................................80 7.1.2 JVavelet transform adoptedto CVsignals.............................................................................80 7.2 WAVELET BISPECTRUM DEFINITION.............................................................................................82 7.2.1 JVavelet bispectrum transform adoptedto CV signals..........................................................84 7.3 WAVELET BISPECTRUM EXAMPLE OF TEST SIGNAL.....................................................................86 7.4 DISCUSSION..................................................................................................................................88 8 FOURIER AND WAVELET BISPECTRUM COMPARISON..................................................91 8.1 WAVELET BISPECTRUM OF CV BLOOD FLOW SIGNALS................................................................91 8.1.1 Results...................................................................................................................................91 8.1.2 Fourier and wavelet bispectrum results comparison............................................................94 8.1.3 Results interpretation............................................................................................................95 8.2 FOURIER AND WAVELET BISPECTRUM ADVANTAGES AND WEAKNESS........................................99 8.3 OTHER POSSIBLE METHODSFOR BISPECTRUM ESTIMATION.......................................................102 8.4 DISCUSSION................................................................................................................................104 9 OTHER POSSIBLE APPLICATIONS OF WAVELET BISPECTRUM METHOD.............107 9.1 Measurements........................................................................................................................109 9.2 Data analysis.........................................................................................................................109 9.3 Results.....................................................................................................................................110 9.4 Discussion................................................................................................................................113 10SUMMARY...................................................................................................................................115 11 CONCLUSIONS...........................................................................................................................123 REFERENCES..................................................................................................................................125 INDEX................................................................................................................................................137 APPENDIX........................................................................................................................................143 A. Variance of the bispectrum estimate.................................................................................143 B. Generation of harmonics......................................................................................................145 Ndmenclature oscillator's amplitude multiplication factor for additional Morlet mother vvavelet stretching biamplitude blood flow signal discrete bispectrum adapted discrete bispectrum cross-bispectrum complex conjugation normalization constant reconstruction constant exponential decav of Gaussian function of Morlet wavelet noise intensitv oscillator's characteristic (cyclic) frequency sampling frequency frequency step forcing amplitude Gaussian function of Morlet vvavelet edge value coupling term unit imaginary discrete frequency (in bins) triplet number of data segments discrete frequency (in bins) L • L area for bispectrum frequency averaging number of samples for tirne averaging discrete Fourier transform window size third-order moment discrete tirne number of periods entering the windowed signal number of samples in the interval T: {T0 - T/2 < t < T0 + 772} percentage of segments overlapping blood pressure signal discrete power spectrum XXXIX r limit cycle radius r(t) respiration signal s(t) surrogate data t continuous tirne A/ tirne step Thf fixed Morlet wave length for high frequencies Tm Morlet mother vvavelet length T 1 qc Time duration of quadratic coupling w(n) window W(s, t) wavelets transform WB wavelet bispectrum WBC cross-wavelet bispectrum X state vector x{n) discrete signal X(k) discrete Fourier transform a oscillator's limit cycle approaching rate £ coupling coefficient M coupling coefficient rj coupling coefficient TJm frequency modulation strength coeficient ¥ phase difference betvveen interacting oscillators M zero-mean white Gaussian noise 71 constant = 3.141592... * biphase A initial biphase (0 angular frequency; cd = 2nf Q observed frequency of interacting oscillator Notation convention. For conciseness the notation B(f\,f2) is often used in this thesis to denote the discrete bispectrum at the discrete bifrequency (k, t) which corresponds most closely to the normalised bifrequency (fufi)- If the frequencies/i,72 fall exactly on discrete frequency bins, then (k = Mfi, l = Mf2). If they fall in between the discrete bins then (k&Mf], I« Mf2). XL ACRDNVMS AND AB B REVI ATI O N S BIS CV CVS DFT EEG ECG FT FB HOS IT RSA SDFT STFT WB WT Bispectral Index CardioVascular CardioVascular System Discrete Fourier Transform ElectroEncephaloGram ElectroCardioGram Fourier Transform Fourier Bispectrum High-Order Statistics Inner Triangle Respiratory Sinus Arrhythmia Selective Discrete Fourier Transform Short tirne Fourier Transform Wavelet Bispectum Wavelet Transform XLI Most real systems are nonlinear and complex. In general, they may be regarded as a set of interacting subsystems; given their nonlinearity, the interactions can also be expected to be nonlinear. Phase relationships between a pair of interacting oscillators can be obtained from bivariate data (i.e., where the coordinate of each oscillator can be measured separately) by use of the methods recently developed for analysis of synchronization, or generalized synchronization, between chaotic and/or noisy systems. Not only can the interactions be detected [79], but their strength and direction can also be determined [72, 93, 94, 103]. The next logical step in studying the interactions among coupled oscillators must be to determine the nature of the couplings; the methods developed for svnchronization analysis do not provide us with the means to answer this question. Studies of higher order spectra, or polyspectra, offer a promising way forvvard in digital signal processing. The approach is applicable to interacting systems, quite generally, regardless of vvhether or not they are mutually svnchronized. Following the pioneering work of Brillinger and Rosenblatt [11], increasing applications of polyspectra have appeared in a diversity of fields, such as: telecommunications, radar, sonar, speech, biomedical, geophysics, imaging systems, surface gravity waves, acoustics, econometrics, seismology, nondestructive testing, oceanography, plasma physics and seismology. An extensive overview can be found in [120]. The use of bispectrum as a means of investigating the presence of second-order nonlinearity in interacting harmonic oscillators has been of particular interest during the last few years [20, 25, 48, 68, 74, 85]. Systems are usually taken to be stationary. For real systems, however, the mutual interaction among subsystems often results in time-variability of their characteristic frequencies. Frequency and phase couplings can occur temporally, and the strength of coupling betvveen pairs of individual oscillators may change with tirne. In studying such systems, bispectral analysis for stationary signals, based on tirne averages, is no longer sufficient. Rather, the tirne evolution of bispectral estimates is required. Priestley and Gabr [80] were probably the first to introduce the time-dependent bispectrum for harmonic oscillators. Most of the subsequent work has been related to the time-frequency representation and is based on high-order cumulants [5, 26]. The parametric approach has been used to obtain approximate expressions for the evolutionary bispectrum [87]. Furthermore, Perry and Amin have proposed a recursion method for estimating the time-dependent bispectrum [76]. Dandavvate and ■■■■i I 1 Giannakis have defined estimators for cvclic and time-varying moments, and cumulants of cyclostationary signals [16]. Schack et al. [98] have recently introduced a tirne-varying spectral method for estimating the bispectrum and bicoherence: the estimates are obtained by filtering in the frequency domain and then obtaining a complex time-frequency signal by inverse Fourier transform. Their assumption is, hovvever, that the interacting oscillators are harmonic. Millingen et al. [61, 62] introduced the wavelet bicoherence, and were the flrst to demonstrate the use of bispectra for studying interactions among nonlinear oscillators. They used the method to detect periodic and chaotic interactions betvveen two coupled van der Pol oscillators, but vvithout concentrating on time-phase relationships in particular. It has long been known that the heart of a healthy human subject in repose does not beat regularly. The rhvthmic variation in the heart rate occurring at the frequency of respiration is known as respiratory sinus arrhythmia (RSA), which can be seen in [3, 18, 23, 34, 117] and references therein; it is not the only arrhvthmia [112]. In fact, at least^zve characteristic frequencies can be seen in blood flow signals, at ~1 Hz, 0.3 Hz, 0.1 Hz, 0.04 Hz and 0.01 Hz. The first two components correspond to the cardiac and respiratory oscillators respectively. The component at -0.1 Hz is often attributed to intrinsic mvogenic activity. The other two correspond respectively to neural and endothelial related metabolic activity. The wavelet transforms of such signals have been discussed in detail [6, 7, 111]. The cardiac and respiratory systems can be perceived from the nonlinear dvnamics point of view as coupled autonomous oscillators, each with its own characteristic frequency [112, 115, 116]. It is respiratory sinus arrhvthmia - the rhvthmic fluctuations of electrocardiographic R-R intervals, or the rhvthmic modulation of the instantaneous cardiac frequency - that provides the most obvious manifestation of their coupling. Although the interaction betvveen the cardiac and respiratory rhvthms has been known to exist since the early works by Hales [29] and Ludwig [56], the underlying physiological mechanisms are not completely understood. In his recent revievv, Eckberg [23] discusses several possible mechanisms for respiratory gating, of both central and peripheral origin: central, secondary to efferent respiratory motoneurone activity; and peripheral, secondary to afferent neural activity from pulmonary and thoracic stretch receptors. He presents a wide range of evidence favouring the influence of respiration on R-R interval fluctuations (as well as on the fluctuations in systolic blood pressure that are strongly correlated to the R-R fluctuations), rather than the influence of peripheral baroreceptor physiology as an origin of modulation. As the physiological mechanisms of cardio-respiratory coupling are not fully understood, even less is knovvn about the nature of this coupling, e.g., vvhether it is linear, quadratic, or of an even higher-order. p 1 INTRDDUCTIDN In addition to the modulation, a mutual adjustment of the cardiac and respiratory rhythms may occur, leading to their synchronization: in a conscious healthy subject at rest, the cardiac and respiratory systems have been shown to synchronize for short periods of tirne [10, 99-101]. The state of the system is characterized by the interactions and couplings between the oscillatory physiological processes. For instance, in anaesthesia, the cardiac and respiratory systems synchronize for more extended periods of tirne [114]. It has also been shown that during anaesthesia, RS A is reduced [1, 83]. Not only C V signals are relevant for studying CVS. Neural C V subsystem coupling information is incorporated in the brain waves. Bispectral index (BIS) is a processed electroencephalogram (EEG) parameter that purports to measure the level of hypnosis in anaesthetized patients. The BIS was formulated retrospectively using a large database of EEG recordings and clinical correlative data. It incorporates parameters derived by high-order spectra, and is one of the most widely applied cases of high-order spectra [89] use. The EEG measures electrical activity in the brain, i.e., brainwaves of different frequencies and short-lived evoked potentials that occur when the brain responds to sensory input. Quantification of nonlinear quadratic phase-coupling between EEG signal components has been established since G. Dumermuth's pioneered investigations using bispectral analysis in 1969 [51]. A number of EEG studies have been published using the mathematical tools of high-order spectra analysis EEG [12, 27, 53, 65]. Recently the depth of anaesthesia vvas found to be related to svnchronization states between the cardiac and respiratory oscillators [63, 64]. The interactions can be detected by analysis of the recorded tirne series, and their strength and direction can also be determined [72, 73, 93, 94, 103]. The next logical step in studying interactions among the coupled oscillators must be to determine the nature of the couplings from the tirne series. A long-term aim is therefore to develop a coupled oscillator model that can provide a description of the system, quantify the couplings and relate their values to its different states of health or disease. We may thus aim for improved techniques of early diagnosis, better assessment of the efficacy of treatment for a range of cardiovascular diseases, and perhaps quantification of depth of anaesthesia. The thesis is organised as follows: Chapters 2 provides a brief overview of the human cardiovascular system, viewed from a point as a dvnamical system, and outlines cardio-respiratory interaction; Chapter 3 briefly describes bispectral analysis and introduces our development of a new approach [40] that introduces tirne dependence to the bispectral analysis of univariate data, while focusing on the time-phase relationships between two (or more) interacting systems. Chapter 4 presents our demonstration/testing of the aforementioned technique on a well-characterized simple model. 3 Examples of different kinds of interaction among the subsystems, e.g., additive linear or quadratic, or parametric frequency modulation, both with, and without, the consideration of zero-mean white Gaussian noise. Chapter 5 gives application of the new technique to univariate cardiovascular (CV) blood flow signals, which reflect the activities of both the local and central mechanisms of cardiovascular regulation [41]. We summarise how measurements are made, discuss how the resultant data are analysed, and present the results. In 5.8 Discussion, it is shown that the cardiac and respiratory processes can be nonlinearly phase coupled. Chapter 6 establishes the relation between bispectral analysis and svnchronization. Chapter 7 introduces the time-frequency variant resolution in the bispectrum, using vvavelet transform. Chapter 8 compares bispectrum based on Fourier transform, and bispectrum based on vvavelet transform. Chapter 9 displays these new potentials of usage, using the čase of an electroencephalogram signal. Chapter 10 provides an overview of the main results of the work, and outlines areas where more work is required. Chapter 11 presents drawn conclusions. Appendix A details the normalization techniques used for comparison of the different measurements. Appendix B provides an analysis of harmonic generation by a pair of weakly-coupled weakly-nonlinear oscillators. 4 Z Cardidvasdular system 2 CARDIDVASCULAR SYSTEM 2.1 Backgroimd 2.2 Coupled oscillators Rhythms are among the most conspicuous properties of living systems. They occur at ali levels of biological organization, from unicellular to multicellular organisms, with periods ranging from fractions of a second to years, Tab. 2.1. In humans, the cardiac and respiratory functions and the circadian rhvthms of a sleep and wakefulness point to the key role of periodic processes in the maintainance of life. In spite of their close association with physiology, however, periodic phenomena are by no means restricted to living systems. Oscillatory chemical reaction was discovered by Bray in 1921, the reaction of Bray, Belousov and Zhabotinsky reaction in 1959, respectively 1964 and others. Tab. 2.1: A list of the main biological rhythms, classified according to increasing period. Hovvever, oscillatory behaviour does not always possess a simple periodic nature. Thus, both in chemistry and biology, oscillations sometimes present complex patterns of bursting, in which successive trains of high-frequency spikes are separated at regular intervals by phases of quiescence. In the phase space, sustained oscillations correspond to the evolution towards a closed curve called a limit cycle by Poincare of it's uniqueness and independence from initial conditions (Minorskv, 1962). Time taken to travel along the closed curve represents the period of oscillations. VVhen a single limit cycle exists, the system always evolves tovvards the same closed curve characterized by a fixed amplitude and period, for a given set of parameter values, regardless of the initial state of the system. Moreover, sustained oscillations of limit cycle type can be viewed as a temporal dissipative structure (Prigogine, 1968). 2A Background The human cardiovascular or circulatory system is one of the basic systems that plays an essential role in the maintenance of a constant internal body environment [3, 55]. It distributes matter and energy to the cells and removes by-products of their metabolism. The cells extract matter and energy from the blood which is pumped by the heart into the network of vessels. The lungs, where the blood becomes oxygenated, are also part of the cardiovascular network. The circulatory system can be divided into two parts: the pulmonary circulation, which moves blood trough the lungs for exchange of oxygen and carbon dioxide and systemic circulation, which supplies ali other tissues. Both, the pulmonary and systemic circulation have a pump, an arterial system, capillaries and a venous system. Arteries and arterioles function as a distribution system. Capillaries serve to exchange diffusible substances between blood and tissue. Venules and veins serve as collection and storage vessels that return the blood to the heart. The heart of a relaxed, healthy subject, pumps an amount equivalent to the total amount of blood in the body (i.e., 4-6 litres or 7-8 % of the body weight) in approximately one minute [55]. Thus, in cardiovascular dvnamics we consider the dynamics of blood distribution trough the cardiovascular netvvork on a tirne scale of approximately one minute. It can be characterised by the dvnamics of the blood flow and the blood pressure in the system, and the activity of the lungs and the heart pump. The research has shown [8, 9, 10, 110-112, 117], that cardiovascular blood flovv signals posses a deterministic dvnamic, meaning, that the cardiovascular system is a result of a final number of subsystems (oscillators), each with its own characteristic frequency. Five subsystems take part in blood flovv regulation: cardiac, respiratory, myogenic, neural and endothelial metabolic system. They ali regulate blood flow that they have in common. The function of the heart is manifested as electric potentials spread across the heart muscle, and as a mechanical pump that rhythmically expels the blood into the arterial netvvork approximately once per 6 2 Cardiovascular system second (1 Hz4) [55]. However, the period of the heart cycle is not constant but, rather, varies in tirne. The lungs can be seen as a pressure generator [55]. The frequency of respiration also varies betvveen 0.15 and 0.3 Hz consequently, the flow and the pressure change in an oscillatory fashion with tirne, and do so on several different tirne scales. Peripheral blood flow depends from the systems that control the blood vein diameter - the blood flow resistance. The smooth-muscle cells in the vessel walls respond continually to the changes in intravascular pressure, which is known as the myogenic response. This intrinsic rhvthmic activity of the vessels, caused by the pacemaker cells in the smooth muscles of their walls, is called mvogenic activity. The vessels contractions respond to the changes in the blood flow and blood pressure and are related to the characteristic frequency at approximately 0.1 Hz [50]. Beside the mvogenic activity at least two more systems influence the veins resistance. The first one is the neural control, provided by the autonomous nerve innervations of vessel. Having it's origin in some centres in the brainstem that are connected to other parts of the central nervous system, and sensors throughout the whole network of vessels, it provides svnchronization of the function of the entire system. It's characteristic frequency is approximately 0.04 Hz [45]. The second one corresponds to metabolic activity. A number of substances, which are required for cellular metabolism or are produced as metabolites, have an effect on the state of contraction of vascular musculature. The rhvthmic regulation of vessel resistance to the blood flow, initiated by concentrations of metabolic substances, can be related to characteristic frequency at approximately 0.01 Hz [3, 8, 110]. 2.2 Coupled oscillators Cardiovascular subsystems do interact among them selves even at healthy human in rest, therefore their characteristic frequencies are not constant but, rather, varie in tirne and their amplitudes are modulated [3, 18, 23, 34, 117]. Each svstem can be regarded as an autonomous nonlinear oscillator with its own characteristic frequency that can be detected in the power spectrum of the cardiovascular blood flow signal. Phase couplings and svnchronization can onset among separate oscillators. It has been shown recently [9, 79, 94, 97, 113, 114] that cardiac and respiratory oscillations are svnchronized. Svnchronization or adjusting in tirne onsets vvhen two or more nonlinear oscillators are coupled. It can be seen in form of their frequency or phase interactions. After Huvgens the svnchronization is frequency adjustment of autonomous harmonic oscillators5 due to their weak 4 Characteristic frequencies differ from human to human. Given values are result of estimated characteristic frequencies obtained by frequency analysis of blood flow signals for čase of young, healthy male humans. 5 Oscillators can be chaotic or noisy. interaction. Phase coupling is thus onset of some phase relations among interacting systems while it doesn't necessary exist any amplitude correlation. The coupling effect on the interacting oscillator's behaviour depends on coupling strength. While weak couplings result in oscillator's characteristic frequency varying, strong couplings can cause qualitative changes in system behaviour named as phase transitions [13]. If there were no couplings among cardiovascular oscillators sharp peaks would occur in power spectrum and the characteristic frequency ratio would be proportional what is not the čase for healthy human in rest where the ratio is at most rational. Appearance of a resonance is a phenomenon, where the system breaks down and the death occurs. Couplings enable information exchange among processes and are basic for normal cardiovascular system activity. Understanding of physiological nature of these couplings is of essential meaning for the understanding of the whole system. Each detected oscillations frequency and amplitude describes the oscillator's activity and impact of ali couplings. It is not possible to measure separate oscillator's activity, therefore various methods are used to analyse phase, frequency and amplitude couplings. A subsystem of the cardiovascular system may be represented as an oscillator, described by a state vector x which satisfies [115, 116]: (2.1) where x is the nonlinear rate function and s denotes the parameter of the oscillator. A simple limit cycle oscillator, as proposed by Poincare [128], can be used to describe a basic unit of the system. It possesses the structural stability and robustness necessitated by physiological understanding and the analysis of measured time series. The state variables x[ and x\ describe the flow and velocity of the flow contributed by the /-th oscillator. 8 2 Cardidvascular system Five oscillators are assumed to contribute to the blood flow trough the cardiovascular system: cardiac, respiratory, myogenic, neural and endothelial related metabolic activity. Each of them is characterized by a frequency / and amplitude ax. The constant oj determines the rate at which the state vector approaches the limit cycle. In polar coordinates, with radius rx and angle oo the state of the oscillator on the limit cycle is essentially described by the angle , the phase. The periodic solution travels around with period Tj = \lf{, Fig. 2.1. Fig. 2.1: The phase plane solutions of differential equations (2.2) for abasic oscillator. The asymptotically stable limit cycle is r\ = ax with phase ^ = 2nf\t, if $ = 0 is taken at t = 0. However, the characteristic frequencies of the cardiovascular system vary in time [6, 7, 110, 111]. Therefore, besides the autonomous part, we suppose that there also exists a component resulting from mutual interactions. Accordinglv, we add a coupling term to differential equations (2.2) Hi .(x{9xJ2) , (2.6) where et is a coupling coefficient. Hi , (x(, xJ2) represents ali possible influences from the rest of the system on the z-th oscillator. In čase of cardio-respiratory interaction, the coupling terms are not symmetrical, i.e., impact of respiration on the heart differs from the impact of the heart on respiration. If there is no coupling betvveen the oscillators (e, = 0), each will oscillate at its own frequency and the state vector in the four dimensional phase space approaches the so-called attracting invariant torus. Analogously to the one-dimensional system discussed above, in the limit as / —► qo, the original system of four differential equations can be reduced to a two dimensional system describing the flow on a two dimensional torus. The amplitudes of both oscillators define the torus, while the flow on a torus can be described entirely in terms of the rate of change of phases of the first ($) and second (^) oscillator. In the uncoupled čase, and so the phase difference increases at a constant rate, determinated by the differences between the natural frequencies of both oscillators (/i and^)- If two oscillators are loosely coupled {et« 0), so that each has only small effect on the other, the invariant torus does not vanish, but is only slightly different. There states are close, \xl(t) - x2(t)\ ~ 0, but remain different. Different types of svnchronization may be expected, depending on the type of coupling. Svnchronization is defined as phase locking or frequency entrainment [79]. In čase of cardiovascular system, with time-varying characteristic frequencies, phase synchronization may onset, while the frequencies may or may not be entrained, we use a weaker condition for phase locking [79] \n\ - m(j>i - d\< const, (2.8) where n and m are integers, (jh, fc are the phases of two oscillators and d is some phase shift. In čase of real systems, measured data contain some noise. It can be instrumental, numerical (resulting from 10 quantization of analogue signals) or physiological (the effect of interactions with the rest of the system on the measured quantity). For weak noise the phase difference ^m = nfr - m(/>2 would be expected to fluctuate in a random way around a constant value. In čase of strong noise, phase slips would occur. h Bispectral analysis belongs to a group of techniques, based on high-order statistics (HOS) that may be used to analyse non-Gaussian signals, to obtain phase information, to suppress Gaussian noise of unknown spectral form, and to detect and characterize signal nonlinearities [60, 68, 69]. In what follows we extend bispectral analysis to extract useful features from nonstationary data, and we demonstrate the modified technique by application to test signals generated from coupled oscillators. 3.1 Bispectral analysis The bispectrum involves third-order statistics. Spectral estimation is based on the conventional Fourier type direct approach, through computation of the 3rd-order moments which, in the čase of 3rd order statistics, are equivalent to 3rd-order cumulants [60, 66-69]. The classical bispectrum estimate is obtained as an average of estimated 3rd-order moments (cumulants) (3.1) where the 3rd-order moment estimate is performed by a triple product of Discrete Fourier Transforms (DFTs) at discrete frequencies k, l and k+l [69]: (3.2) with i = 1,..., K segments into which the signal is divided to try to obtain statistical stability of the estimates [67]. Algorithm for bispectrum estimation is described in detail in [39]. The bispectrum B(k, /) is a complex quantity, defined by magnitude A and phase (3.3) Consequently, for each {k, /), its value can be represented as a point in a complex space, %l[B(k, /)] versus 3[B(k, /)], thus defining a vector. Its magnitude (length) is known as the biamplitude. The phase, which for the bispectrum is called the biphase, is determined by the angle between the vector and the positive real axis. As discussed in detail in [39], the bispectrum quantifies the relationships among the underlying oscillatory components of the observed signals. Specifically, bispectral analysis examines the relationships between the oscillations at two basic frequencies, k and / and a harmonic component at the frequency k +/. This set of three frequencies is known as a triplet (k, /, k+ /). The bispectrum B(k, /), a quantity incorporating both phase and power information, can be calculated for each triplet. A high bispectrum value at bifrequency (k, l) indicates that there is at least frequency coupling within the triplet of frequencies k, /, and k ± /. Strong coupling implies that the oscillatory components at k and / may have a common generator. Such components may svnthesize a new component at the combinatorial frequency, k± /, if a quadratic non-linearity is present. 3.2 Time-phase bispectral analysis The classical bispectral method is adequate for studving stationary signals whose frequency content is preserved over tirne. We now wish to encompass tirne dependence within the bispectral analysis. In analogy with the Short Time Fourier transform (STFT), we accomplish this by moving a tirne window w(ri) of length Macross the signal x(n), calculating the DFT at each window position [81] (3.4) Here, k is the discrete frequency, n the discrete tirne and rthe tirne shift. The choice of window length M is a compromise betvveen achieving optimal frequency resolution and optimal detection of the time-variability. The instantaneous biphase is then calculated: from Eqs. (3.2) and (3.3), it is (3.5) If the two frequency components k and / are frequency and phase coupled, ^+1 = 0.88 Hz and below/J < 0.3 Hz are cut, because the triplets (1.1 Hz, 1.1 Hz, 1.1 Hz) and (0.24 Hz, 0.24 Hz, 0.24 Hz) produce high peaks that are physically meaningless. (c) Adapted biphase $, and (d) biamplitude Aa for bifrequency (1.1 Hz, 0.24 Hz), using a 0.3 s tirne step and a 100 s long Blackman window for estimating the DFT. The peak (1.1 Hz, 0.24 Hz) indicates that oscillations at those pairs of frequencies are at least linearly frequency-coupled. Frequency coupling alone is sufficient for a peak in the bispectrum to occur. Although the situation can in principle arise by coincidence, frequency and phase coupling together strongly imply the existence of nonlinearities. To be able to distinguish between different possible couplings we calculate the adapted biphase Fig. 4.2 (c). During the first 400 s of test signal X\A, where no coupling is present, the adapted biphase changes continuously between 0 and 2tt radians. For the same tirne of observation it can be seen that the adapted biamplitude is 0, Fig. 4.2 (d). During the second and third 400 s of the signal X\A, a constant adapted biphase can be observed indicating the presence of linear coupling. The value of the adapted biamplitude is higher in the čase of stronger coupling. The coupling constant 772 can be obtained by normalization, and we are thus able to define the relative strengths of different couplings. 20 When the oscillators are coupled bidirectional the frequency content of each of them changes and components 2/1 and 2f2 are generated. Both of these characteristic frequencies can be observed in the tirne series of each oscillator. Two combinatorial components are also present in their spectra, 2/1 -f2 and/1 - 2/2, assuming that/1 > f2. In analysing bidirectional coupling the procedure described above can be extended and two combinatorial components should be analysed in the same way. Making use of the calculated instantaneous phases of both oscillatory components we also construct a synchrogram (Fig 4.1 (c)), as proposed by Schafer et al. (see Ref. [79] and the references therein), and can immediately establish whether or not the coupling also results in synchronization. The instantaneous phases can also be used to calculate the direction and strength of coupling, using the methods recently introduced by Schreiber, Rosenblum et al. and Paluš et al. [72, 93, 94, 103]. 4.2 Linear couplings in presence ofnoise We now test the method for the čase where noise is added to the variable x\ of the first oscillator: (4.5) Here, %j) is zero-mean white Gaussian noise, (4(0) = 0> (4(0» £$)) = D$J), and D = 0.08 is the noise intensity. In this way we obtain a test signal X\B,(t), Fig. 4.3 (a). For nonzero coupling strength /72, the component at frequency position/3 can stili be seen in the power spectrum despite the noise Fig. 4.3 (b). The adapted biphase Fig. 4.3 (f) can clearly distinguish betvveen the presence and absence of coupling. When coupling is weaker, the adapted biamplitude [Fig. 4.3 (g)] is lower and the adapted biphase is less constant. The bispectrum for the signal jcib, shown in Fig. 4.4 (a), differs from that in the čase of zero noise, Fig. 4.1 (d). Noise raises two additional peaks positioned at (1.1 Hz, 0.24 Hz) and (0.86 Hz, 0.24 Hz). The former could be the result of interaction; the latter is due to the method: the sum of the frequencies in this bifrequency pair gives the frequency of the first oscillator. Close inspection of the (1.1 Hz, 0.24 Hz) peak by calculation of the biphase gives Fig. 4.4 (c). When coupling is present the characteristic frequency of the second oscillator appears in the power spectrum Fig. 4.3 (b). Two frequencies of high amplitude result in a small peak even if no harmonics are present at the sum and/or difference frequencies. The second peak is not of interest to us. It can easily be checked whether a phase coupling exists among the bifrequencies from the tirne evolution of the biphase. Fig. 4.3: Results in the presence of additive Gaussian noise. (a) Test signal x1B, variable Xi of the first oscillator with characteristic frequency/j = 1.1 Hz. The characteristic frequency of the second oscillator is 72 = 0-24 Hz. The oscillators are unidirectionally and linearh/ coupled with three different coupling strengths; r]2 = 0.0 (1); 0.1 (2); and 0.2 (3). Each coupling lasts for 400 s at a sampling frequency^ = 10 Hz. Only first 15 s are shown in each čase. (b) Its povver spectrum and (c) svnchrogram. (d) Adapted bispectrum \Ba\ using K = 33 segments, 66 % overlapping and the Blackman window and (e) its contour view. The parts of the \Ba\ above^ > 0.79 Hz and below/i < 0.3 Hz are omitted because the triplets (1.1 Hz, 1.1 Hz, 1.1 Hz) and (0.24 Hz, 0.24 Hz, 0.24 Hz) produce a high peak that is physically meaningless. (f) Adapted biphase $, and (g) adapted biamplitude Aa for bifrequency (1.1 Hz, 0.24 Hz), using 0.3 s tirne step and 100 s long window for estimating the DFT using the Blackman window. In general, besides estimating bispectral values, one can also observe the tirne dependences of the phase and amplitude for each frequency component, and their phase relationships. This applies particularly to frequencies that form a bifrequency giving a high peak in the bispectrum or adapted bispectrum. Synchrograms, Fig. 4.1 (c) and 4.3 (c), are obtained by first calculating the instantaneous phase of each oscillator, and then their phase difference [79]. The phase difference in this čase is 22 4 Analysis df cduplings ; ''':* ;':1: between two fixed frequencies. We do not calculate their instantaneous frequencies, although it is possible to follow the frequency variation by calculating the phase difference at neighbouring bifrequencies around the observed one and showing them simultaneously on the same plot. Examples of the phase difference *F= fa - fo between the phases of the first (f>\ and the second & interacting oscillators are shown in Fig. 4.4 (e) and (f). Bispectral relation to svnchronization is discussed in detail in Chapter 6. Fig. 4.4: Bispectrum \B\, calculated from the signal xlB presented in Fig. 4 (a), using K = 33 segments, 66 % overlapping and Blackman window to reduce leakage and (b) its contour view. (c) Biphase and (d) biamplitude A for bifrequency (1.1 Hz, 0.24 Hz), using a 0.3 s tirne step and a 100 s long window for estimating the DFTs using a Blackman window. (e) Phase difference 5P between $\ of the characteristic frequency component/! of the first oscillator and ^ of the characteristic frequency component^ of the second oscillator, for tirne step \lfs and (f) at each period of lowest frequency \lf2 in the bifrequency pair (1.1 Hz, 0.24 Hz), using interpolation and 100 s long window for estimating DFTs using Blackman window. 23 4.3 Quadratic couplings We now assume that two Poincare oscillators can interact with each other nonlinearly. A quadratic nonlinear interaction generates higher harmonic components in addition to the characteristic frequencies [69]. In order to study an example where the first/1 = 1.1 Hz and second^S - 0.24 Hz oscillators are quadratically coupled, we change the coupling terms in the model (3.1) to quadratic ones Fig. 4.5: Results for quadratic coupling in the absence of noise. (a) The test signal xiC, variable x\ of the first oscillator with characteristic frequency/i = 1.1 Hz. The characteristic frequency of the second oscillator is^ = 0.24 Hz. Oscillators are unidirectionally and quadratically coupled with three different couplings strengths: t]2 = 0.0 (1); 0.05 (2); and 0.1 (3). Each coupling lasts for 400 s at sampling frequency/s = 10 Hz. Only the first 15 s are shown in each čase. (b) Its power spectrum. (c) The bispectrum \B\ using K = 33 segments, 66 % overlapping and the Blackman window to reduce leakage and (d) its contour view. The part of the bispectrum abovef2 > 1,0 Hz is cut, because triplet (1.1 Hz, 1.1 Hz, 1.1 Hz) produces a high peak that is not physically significant. Clearly, the test signal jcic presented in Fig. 4.5 (a) for three different coupling strengths (no coupling rj2 = 0 (1); and weak couplings rj2 = 0.05 (2), rj2 = 0.1 (3)) has a richer harmonic structure. In addition to the characteristic frequencies, it contains components with frequencies 2/1, 2f2,f\ +f2 and/i -f2 Fig. 4.5 (b). Eq. (4.6) also indicates that, as well as having a particular harmonic structure, the components of the signal X\c also have related phases, 2$, 2^2, \- 1.0 Hz is cut, because the triplet (1.1 Hz, 1.1 Hz, 1.1 Hz) produce a high peak that is physically meaningless. (e) The biphase ^ and (f) biamplitude A for bifrequency (1.1 Hz, 0.24 Hz), with a 0.3 s tirne step and a 100 s long window for estimating DFTs using the Blackman window. 26 4 Analysis of cduplings 4.5 Frequency modu lati on in the presence ofnoise We are also interested of being able to detect parametric frequency modulation and to distinguish it from quadratic coupling. Parametric modulation produces frequency components at the sum and difference of the characteristic frequency and the modulation frequency, i.e., the same two frequency components that can also result from quadratic coupling. Let us now consider an example where the first oscillator/i = 1.1 Hz is frequency modulated by the second onQ/2 = 0.24 Hz. For this purpose the equations of the first oscillator become (4.7) The model parameters a^2, a\,i and the noise intensity D are chosen to be the same as in the previous examples. We thus obtain a test signal ;t1E. It is the tirne evolution of the variable X\ of the first oscillator, presented in Fig. 4.8 (a) with the corresponding power spectrum 4.8 (b) for three different parametric frequency modulation strengths: no modulation rjm = 0 (1); and modulation rjm = 0.1 (2), r/m = 0.2 (3). The bispectrum of the test signal jciE, Fig. 4.8 (c), exhibits several high peaks. The highest are at bifrequencies (1.1 Hz, 0.86 Hz), (0.86 Hz, 0.24 Hz) and (1.1 Hz, 0.24 Hz), in addition to the (1.1 Hz, 1.1 Hz) peak. They also appear in the čase of quadratic coupling. In general, hovvever, the other peaks that appear for quadratic coupling are absent. The reason is that although the component of the second oscillator f2 (one component of the triplet) is not present in the power spectrum, its value is not exactly zero. Observing the biphase, no epochs of constant biphase can be observed, although for strong frequency modulation the biphase is less variable. In the power spectrum, Fig. 4.8 (b), no component rises above the noise level at frequency f2, of the bifrequency pair, where the bispectrum peaks. This is an indication that there is parametric coupling between the oscillators as there is a high value of biamplitude. The biphase changes runs between 0 and 2tc, and is modulated in the absence of noise. There are also no rapid 2n phase slips of the kind that are normal if no modulation is present. In the absence of couplings and modulation, but with noise present, there would be no such peaks in the power spectrum and bispectrum. 21 Fig. 4.8: Results for parametric frequency modulation in the presence of additive Gaussian noise. (a) The test signal x1E, of variable xx of the first oscillator with characteristic frequency/i = 1.1 Hz frequency modulated by the second oscillator/> = 0.24 Hz with three different frequency modulation strengths; rjm = 0.0 (1), 0.1 (2) and 0.2 (3). Each frequency modulation lasts for 400 s, at sampling frequency/s =10 Hz. Only the first 15 s are shovvn in each čase. (b) Its power spectrum. (c) The bispectrum \B\ calculated with K = 33 segments, 66 % overlapping and using the Blackman window to reduce leakage and (d) its contour view. The part of the bispectrum above^ > 1.0 Hz is cut, because the triplet (1.1 Hz, 1.1 Hz, 2.2 Hz) produces a high peak that is physically meaningless. (e) The biphase (j) and (f) biamplitude A for bifrequency (1.1 Hz, 0.24 Hz), with a 0.3 s tirne step and a 100 s long window for estimating the DFTs using the Blackman window. 28 S Cardid-respiratdry interactidns 5 CARDID-RESPIRATDRY INTERACTIDNS 5.1 Data acauisition 5.2 Measurements 5.3 Data analysis 5.4 Results 5.5 Surrogates 5.6 Global couplings 5.7 Cross-bispectrum 5.8 Discussion 5.8.1 Definition ofthe phase 5.8.2 Nonlinear coupling, or linear coupling ofstrongly nonlinear oscillators? 5.8.3 Relationship to synchronization 5.8.4 Synchronization, modulation andtype of coupling 5.8.5 Unidirectional or bidirectional coupling 5.8.5.1 Forced oscillator 5.1 Data acquisition The interaction between two harmonic components can in practice contribute to the power at their sum and/or difference frequencies. We assume that the cardiac and respiratory oscillators are weakly coupled and can interact with each other nonlinearly. The coupling is assumed to be weak in part because of the transient/episodic character of cardio-respiratory synchronization in healthy subjects; the assumption of weak nonlinearity is on account of several factors including the lack of combinatorial components near the cardiac frequency. We return to these questions and discuss them in more detail at the end of Sec. 5.8.2. A quadratic interaction will give rise to higher harmonic components vvith frequencies 2/i, 2f2,f\ +f2 and/i -f2, in addition to the characteristic frequencies [ 68, 69]. As well as having a particular harmonic structure, the components also have phases that are related, 2, and five others. 32 To investigate the cardio-respiratory coupling eight peaks were analysed for each signal, as shown in Tab. 5.2. To be able to compare results, a normalization procedure was performed. The maximum biphase and biamplitude were calculated for each peak. The frequency resolution was set to be 1/10 of the lowest respiration frequency or better. The slowest-paced breathing, f2, was approximately 0.1 Hz, so that a window of 100 s or longer was necessary for estimation of the bispectrum, biphase and biamplitude. Short plateaus in the estimated biphase occur frequently. To exclude coincidence interactions we focused on those that lasted for at least about 10 periods of the lower coupling frequency/2- The length of the window also determines the time resolution. In the čase of the slowest-paced breathing, where^ was 0.1 Hz we were seeking approximately 10 times (1/ f2) = 100 s long epochs of constant biphase. Therefore a window length of 100 s or less was necessary to meet the criterion for time resolution. Due to the Heisenberg uncertainty principle [43], the scope for choice of window length is limited, and compromise is needed between time and frequency resolution. The window was moved along the time series with a minimum time step of \/fs = 0.1 s, vvhere^Š is the sampling frequency. The critical value for the biamplitude estimate to be considered valid was set in ali cases to 2, i.e., tvvice the average value of the bispectrum within its so-called inner triangle (IT). To be able to conclude that quadratic coupling exists, we require several condition to be fulfilled: (i) A constant biphase during at least 10 periods of the lower interacting component; (ii) Biphases for ali six (eight) peaks must be present at the same time as the biphase plateau; (iii) No phase slips must occur during the coupling, and the biphase variations must stay within n rad interval (the biphase being expected to be more or less constant, depending on the coupling strength and noise intensity: phase slips are frequent when the interaction is extremely weak; they are mostly due to noise, but sometimes caused by modulation; strong modulation is expected to result in a biphase with fewer phase slips); (iv) The biamplitude must be above the specified critical value, i.e., be more than tvvice the average bispectrum value vvithin the IT. 5.4 Results Examples of detrended, resampled, blood flow signals are presented in the left-hand column of Fig. 5.1. These signals correspond to the čase of paced respiration slovver than the natural frequency. Their calculated frequency content is presented in the right hand column. The peak at ~0.98 Hz belongs to cardiac activity,/i; that at ~0.11 Hz to respiratory activity,/2, which was also obtained directly as a check by use of a piezo sensor. Although the characteristic frequencies differ from person to person, they ali lie within defmed frequency bounds. Stefanovska proposed in [112] that the respiration 33 S CARDID-RESPIRATaRY INTERACTIDNS frequency interval should be defined as 0.14-0.6 Hz with a median frequency 0.3 Hz. The spontaneous respiratory frequency for person 2 was 0.14 Hz, i.e., it fell at the lower limit of this interval. The slowest/?ace» 5 Cardid-respiratdry interactions attributable to cardio-respiratory coupling; at (0.87 Hz, 0.11 Hz) which we assume to be coupling betvveen the respiratory component f2 and the difference f\ - f2, that could be due to a nonlinear coupling mechanism; and two peaks attributable to interaction with lower cardiovascular characteristic components. The latter interactions (with the intrinsic myogenic and neural oscillators) are not of interest in the present context. Also, other lower frequency peaks can be seen in the bispectrum. Their positions can be seen in the bispectrum contour view shown in Fig. 5.2 (c). Fig. 5.2: (a) The bispectrum \B\ for a signal ba, calculated with K = 33 segments, 87 % overlapping and using the Blackman window to reduce leakage. (b) Part of the bispectrum fhf2 < 1.4 Hz that is of our interest and (c) its contour view. The bispectrum is sensitive to tirne-variations of the frequency components, yielding in the bispectrum a characteristic diagonal elongation of peaks. The cardiac frequency/i spans 0.93-1.02 Hz. Although the respiratory frequency^ extends from 0.09 Hz to 0.12 Hz, this large range is actually the result of a single deep breath: the respiratory frequency (being paced) is constant for most of the tirne, leading to a high bispectrum. The cardio-respiratory bifrequency coupling consequently has a wide frequency range resulting mainly from variation of the cardiac frequency (in Fig. 5.2 it is elongated along the/i axes). 35 fj(Hz) 1.16 768 172.4 Time (s) 471 76.8 172.4 Time (s) 47' Fig. 5.3: Blood flow analyses for signal ba(t), calculated with K = 33 segments, 87 % overlapping, using a 0.1 s time step and a 100 s long window for estimating the DFT, with a Blackman window to reduce leakage, for peaks (a) 1, (b) 2, (c) 3, (d) 4, (e) 5 and (f) 6; len column, the bispectrum |Z?ba| with its corresponding contour plots; middle, the biamplitude A\,a; and right, the biphase $,„. 36 5 Cardid-respiratdry interactions The peaks corresponding to cardiac activitv are lower, mainly due to tirne frequency variations. In the presence of quadratic nonlinear coupling, peaks should be present at ali six of the bifrequencies summarized in Tab. 5.2. Significant values of the bispectrum (i.e., exceeding twice the average bispectrum \vithin the IT) were obtained near ali these bifrequencies (see the left-hand column of Fig. 5.3), showing that the peaks are indeed present. Time (s) Fig. 5.4: (a) Determination of quadratic nonlinear coupling interval Tqc and (b) maximum biphase variation A^ for signal ba for peak 4. Once the peaks at the defined bifrequencies had been confirmed, the tirne biphase and biamplitude were calculated at the bifrequency peaks. The time interval Tqc during which quadratic coupling persisted was determined; if ali 6 peaks fulfilled our conditions, then the Tqc interval was calculated for ali peaks and the boundaries were defined such that the biamplitude for ali the peaks in Tqc interval would be above the condition Fig. 5.4 (a). For each Tqc interval, the maximum variation of the biphase A$ was determined as shown on Fig. 5.4 (b). It can be seen that the biamplitude during the time interval from 76.8 s to 172.4 s meets our criterion of being more than tvvice as large as the average bispectrum in the IT: see middle column of Fig. 5.3 (a) to (f). The biphase in this time interval, 95.6 s long (shaded area), remains constant within a 1.47 rad interval, i.e., there are no phase slips. The biphases at bifrequencies 1, 2, 3 and 5 are very constant; those at 4 and 6 are less so, but they stili remain within the n rad interval. 5.5 Surrogates In many areas of signal processing, an important problem exists in determining whether an observed time series is deterministic, contains a deterministic component, or is purely stochastic. The surrogate 31 RACTIONS data method provides a rigorous way to determine if an observed tirne series has a statistically significant deterministic component [102, 106]. The question of whether or not the cardiovascular system possesses a deterministic dynamic was already subject to many research activities by a variety of analyses [6, 8, 109, 110, 112]. These provided notable evidences that the system that regulates blood flow is deterministic. In this work we use the surrogate data method to validate bispectral analysis to determine whether the obtained biphases from the cardiovascular blood flovv signals are a result of the deterministic dvnamic of blood flovv, or as a result of its stochastic component. For this purpose, we use the phase randomization method of generating the surrogates [44, 102, 106, 123, 124]. The discrete Fourier transform is computed from the original data, which consists of amplitude and phase at each frequency, and then shuffled to destroy ali correlations. The surrogate is the inverse discrete Fourier transform of the shuffled data. Besides phase shuffling, one can also perform phase randomization or data shuffling. Each amplitude is replaced by the amplitude of the same frequency as in the original data. After the inverse Fourier transform, the amplitude of the surrogate data is adjusted by applying a nonlinear transform to give the surrogate data the same distribution as in the original data. In this way we obtain surrogate data which have similar spectral properties as that of the original data. The surrogate data sequence has the same mean, the same variance, the same autocorrelation function, and therefore, the same power spectrum as the original sequence, but the phase relations are destroved. The generated surrogate data are output of a linear Gaussian process. We posit a null hvpothesis: H0: Quadratic nonlinear coupling is present. We determine the null hvpothesis on surrogate data Sba(0 generated from blood flow signal ba{t). The obtained surrogate signal Sba(t), Fig. 5.5 (a), has almost the same power spectrum. Fig. 5.5 (b), as it is the povver spectrum Pba of the original signal ba(t), Fig 5.1 (a). As expected, we also obtained a similar bispectrum for the whole frequency domain, Fig. 5.5 (c). The peaks are lower, and less evident, compared to bispectrum of ba(t), Fig. 5.2. (b). 38 5 Cardid-respiratdry interactidns Fig. 5.5: (a) Surrogate of blood flow signal from the right wrist sha(t), its power spectrum (b), and (c) its bispectrum \B\, calculated with K = 33 segments, 87 % overlapping and using the Blackman window to reduce leakage. Fig. 5.6: Blood flow analvses for surrogate signal sha(t), calculated with K = 33 segments, 87 % overlapping, using a 0.1 s tirne step, a 100 s long window for estimating the DFT, and using a Blackman window to reduce leakage, for peaks (a) 1, (b) 2 and (c) 3; left column, the bispectrum |5S| with its corresponding contour plots; middle, the biamplitude As; and right, the biphase ^s. By calculating the biphases for ali of the six bifrequencies summarized in Tab. 5.2, we conclude that the necessarv conditions for quadratic nonlinear coupling are not fulfilled. Fig 5.6 shows peaks, 39 NTERACTIONS biamplitudes and biphases for the first three bifrequencies of bispectum Bs. It can be seen that ali the biphases, , bc and b^ measured on channels c-d. As excepted the maximum amplitude of the peak is positioned at the same bifrequency (0.98 Hz, 0.11 Hz) as already seen for peak 6 of the signal ba measured on channel a, Fig. 5.8, midddle 40 5 Cardid-respiratdry interactidns column. The correlation of the biamplitude and biphase for signals ba and b\>-b& for ali peaks 1-6 is very high. For peak 6, for signal ba and bd, is it is 0.85 for both, the correlation of the biamplitude and biphase, as can also be seen from their tirne evolution presented in the right-hand and middle column of Fig. 5.8. Fig. 5.7: Left: the bispectrum |B| calculated with K = 33 segments, 87 % overlapping, and using the Blackman window to reduce leakage for the test signal (a) bh, (b) bc and (c) bd for peak 1 with its contour view. Middle: the biamplitude A, and right the biphase (/>; a 0.1 s tirne step and 100 s long window were used for computation of the DFT using a Blackman window. The biamplitude meets our amplitude criterion vvithin the same tirne interval from 76.8 s to 172.4 s, and the biphase is also constant during this interval. We obtain the same coupling information at ali four measuring sites for ali the peaks (1-8). The results obtained from time-bispectral analvses of the measured signals are summarized in Tab. 5.3. Inspecting the data in Tabs. 5.1 and 5.3 we see no obvious correlation between the average tidal volume and the onset of nonlinear cardio-respiratorv interaction. II RACTICJNS Fig. 5.8: Left: the bispectrum |B| calculated with K = 33 segments, 87 % overlapping and using the Blackman window to reduce leakage for the test signal (a) bh, (b) bc and (c) bA for peak 6 with its contour view. Middle: the biamplitude A, and right the biphase (jr, a 0.1 s tirne step and 100 s long window were used for computation of the DFT using a Blackman window. 12 Tab. 5.3: Quadratic nonlinear couplings detected in blood flow signals marked by * in Tab. 5.1. For each measurement, four blood flow signals were measured simultaneously at different sites, channels a-d. Tqc is the tirne interval during which the bispectral analysis showed that the heart oscillator f^ and the respiratory oscillator fies might be nonlinear coupled. The product of Tqc -frt% tells us over how many respiratory periods the interaction persisted. During Tqc the maximum biamplitude is calculated for the peak 1 that is of our primary interest. In addition, the maximum variation of the biphase A0, its average value (/) , and its standard deviation c^ were calculated during Tqc. 5.7 Cross-bispectrum The bispectrum as defined in [40] can be seen as a special čase of the cross-bispectrum when the three signals are the same. In addition to the blood flow signals, the ECG e(t), respiration r{f) and blood pressure p{i) were also simultaneously recorded. This gave us the possibility of globally checking the coupling between cardiac and respiratory activity using bivariate data. Let us define the cross-bispectrum as [69] (5.1) where X and Y are discrete Fourier transforms of two different signals x(t) and y{i) at discrete frequencies k, l and k + L We calculated the cross-bispectrum Bcehb (where c stands for cross, e for signal e{i) and b for signal b(t)), for the čase where x(t) is the ECG signal e(t) and y{t) is the blood flow signal ba(t). The ECG signal tells us primarily about the cardiac electrical activity. The phase of 43 iNTERACTICJNS the first, cardiac component/ls in the triplet (fufiifi+fi) is thus directly extracted from the ECG signal. The respiratory component f2 and the component at the harmonically related position f\ + f2 are extracted from the blood flow signal. We also define the cross-bispectrum as (5.2) and calculate it for 2 different cases: (i) i?cbrb, where x(t) is the blood flow signal bjf) and y(t) is the respiration signal r(t). The signal r(i) most directly describes the activity of the respiratory oscillator. Therefore the phase of the second component in the triplet (f\,fi,f\ +fi) is directly extracted from the respiratory signal, (ii) 5cprp, where x(t) is the blood pressure signal p{t) and y(t) is the respiration signal r(f). By calculating the latter cross-bispectrum we are interested in establishing whether the information about coupling between the heart and the respiratory oscillators is signal-independent, i.e., whether it is also present in other CV signals. Fig. 5.9: Results for the cross-bispectrum. Row (a) shows the first 25 s of the signals: left ECG signal e(f), middle respiration signal r(t), and right differentiated blood pressure signal p(t); whereas row (b) shows their power spectra. The sampling frequency was fs = 400 Hz. In the lower rows (c) the biphase (f> and (d) the biamplitude A are shovra; a 0.1 s tirne step and 100 s long window were used for computation of the DFT using a Blackman window. The biphase and the biamplitude were calculated using the cross-bispectrum (left) Behh, (middle) Bhlh and (right) 5prp. m We proceeded as discussed above in Sec. 5.3. for each of the three different cross-bispectrum cases. The time evolution of the signals e(t), r{t) and p{t) and their power spectra are presented in Fig. 5.9 (a) and (b). For the cross-bispectrum Z?cebb ali the peaks from 1 to 8 are present at the same bifrequencies as in the čase of auto-bispectrum of the blood flow signal. Since the power spectrum of the respiratorv signal Prexhibits only components of the respiratorv oscillator we cannot expect peaks 5, 6 and 8 to appear in the cross-bispectrum BcbTb and Bcprp. Ali the rest of the peaks are present at the same bifrequencies as in čase of the auto-bispectrum. The biamplitude meets our amplitude criterion vvithin the same time interval from 76.8 s to 172.4 s, for ali peaks; moreover the biphase is constant during this time interval. Examples of the biamplitude and biphase time evolution for peak 1 for the three different cross-bispectrum are presented in Fig. 5.9 (c) and (d). Cross-bispectrum were also calculated with surrogate data, where the phases of the frequency components of the signals e(t), r(t) and p(i) were randomized. No phase couplings were detected in this čase. The coupling information among cardiac and respiratory process seems to be signal independent. 5.8 Discussion The signals were measured on six persons, whereas in the first column of Tab. 5.3 data are only provided for five; the sixth person also showed evidence of nonlinear couplings, but over time, these intervals were too short to fulfil our required conditions. Four blood flow signals, simultaneously measured in different places (channels a, b, c, and d), were available for each recording. We usually first analysed the one with the most distinctive characteristic frequencies in its power spectrum and, if our criteria were fulfilled, checked the other three signals as well. The time interval, Tqc, during which quadratic coupling persisted, was determined; if ali 6 peaks fulfilled our conditions, then the Tqc interval was calculated for ali peaks, and the boundaries were defined such that the biamplitude for ali the peaks in Tqc interval vvould be above the condition. Also shown in Tab. 5.3 are three cases where the couplings Tqc lasted less than 10 • 1//S, where \/f2 is the longest respiratory period, since they could be detected very distinctly and clearly. Column Aimax is the maximum biamplitude for peak 1 during the Tqc interval. The strength of the coupling is, in general, not correlated with its duration. For each Tqc interval, the maximum variation A of the biphase, its average value, and its standard deviation, were calculated. 45 ::......mSttttttm Fig. 5.10: (a) Example of blood flow signal during spontaneous breathing for person 5, channel a and (b) its power spectrum. (c) Bispectrum |B| and its contour view calculated with K = 33 segments, 87 % overlapping, and using the Blackman window to reduce leakage. There was only a single čase of spontaneous breathing for which the coupling lasted long enough (82 s without phase slips) to fulfil the criteria. During the spontaneous respiration phase slips are relatively frequent, and epochs with constant biphase are short. Fig. 5.10 shows an example of data for spontaneous breathing (a) and its power spectrum (b). Signals of this kind are, in general, very difficult to analyse due to their large time-variable frequency content, which results in broadened and coalescing peaks, as seen in Fig. 5.10 (c). Peaks 1-6 cannot easily be resolved due to time-frequency resolution restrictions. 5.8.1 Definition of the phase As already mentioned above in Sec. 2.2, the cardiac and respiratory systems can be perceived as coupled autonomous oscillators [6, 7, 9, 13, 59, 75, 99-101, 110, 111, 115]. Using a bispectrum based on Fourier transform, which is a decomposition of the signal in terms of complex exponential (sinusoidal) components, each component can be represented as a point in a complex space, 9?[X(A:)] versus 3[X(&)], thus defining a vector, where X is DFT and A: is a discrete frequency. Its magnitude represents a power, whereas the phase is determined by the angle between the vector and the positive real axis. The phase of an oscillator is defined as the phase of the sinusoidal component that lies closest to the characteristic frequency of the oscillator, with a corresponding spectral peak. Thus our phase definition differs from that in Kuramoto phase reduction, [54]. In this way it is possible to study the phase relations and resolve the nature of the couplings. 46 5 Cardid-respiratciry interactidns 5.8.2 Nonlinear coupling, or linear coupling of strongly nonlinear oscillators? Our study is based on the assumptions that the cardiac and respiration processes can be described as weakly nonlinear oscillators and that the interaction betvveen them is also weak [116]. It is pertinent to investigate what happens when these assumptions are not fulfilled. We have addressed the question in two different ways based on analytic approximation and digital simulation, respectively. The analysis in Appendix B considers harmonic generation by a pair of coupled, weakly nonlinear, oscillators. It confirms that, for weak coupling, the appearance of additional harmonics at la^, 2a>\, 2a>\ ± 2o>2, a>\ ± a>2, 3coi ± a>2 can confidently be associated with the presence of quadratic coupling. For a sufficiently nonlinear oscillator and sufficiently strong coupling, these and other combinatorial harmonics can in principle be generated, as a second-order effect, even for linear coupling. As will be illustrated below, however, the appearance of these combinatorial harmonics does not in itself fulfil the necessary conditions to support a conclusion that there is nonlinear coupling when tested by bispectral analysis: the high biamplitude and constant biphase may be absent. In any čase, the bispectral approach cannot be expected to yield reliable information about the nature of the coupling when the nonlinearities are extremely strong. We have therefore complemented the analysis of Appendix B vvith a digital simulation, exploring the range of extreme conditions where the bispectral approach is expected to fail. We have chosen to simulate a generic model (5.3) of the van der Pol relaxation oscillator, vvith an additional nonlinear term, linearly driven by another relaxation van der Pol oscillator via an additive coupling [125] (5.3) The activity of the oscillators is described by the two state variables X{ and_yi, vvhere % and a>\ = 2nf[ are constants, and i = 1,2 denote respectively the driven or the driving oscillator. 771 is a constant that sets the strength of the additional nonlinear term in the driven oscillator, and ju\ is the strength of the coupling. Here, Č[i) is zero-mean white Gaussian noise, (%j)) = 0, {č{i), ^(0)) = DS(t) and D = 0.8 is the noise intensity. Following the pioneering work of van der Pol and van der Mark [125], the parameters are set to S\ = 70, and £2 = 3. A detailed parameter space analysis has been completed, showing that a situation indeed exists for which the bispectral technique fails to distinguish betvveen the two situations when: (i) two oscillators are strongly nonlinear, but linearly coupled; or (ii) when they are nonlinear and nonlinearly coupled. As an illustration, bispectral analysis was performed for two coupled van der Pol oscillators, with and without added Gaussian noise, for different sets of parameters. In the first čase the strength of the additional nonlinearity was changed, while jU\ was kept constant: (a) ju\ = const. = 1; i)\ was 0, 0.5, 1, 2, 5, 10, 12 and 15, and from 20 to 90 was varied with step 10 and 93 was also included. (b) ju\ = const. = 25; J]\ was varied from 1 to 10 vvith step 1, values 12, 15, 18 and 25 were also considered, and again values from 20 to 60 with step 10. In the second čase the strength of the coupling was changed vvhile 771 was kept constant: (a) r\\ - const. = 1; ju\ was varied from 0.1 to 1 vvith step 0.1, from 1 to 10 vvith step 1, then values 2.5, 3.5, 11, 15, 20, 25, 30, 35, 120 and 200 vvere considered, and again values from 40 to 100 vvith step 10. (b) 771 = const. = 5; /4 = 25. (c) /71 = const. = 15: /4 was varied from 1 to 10 vvith step 1, values 12, 15, 17, 20, 24 and 25 vvere considered, and again from 30 to 60 vvith step 10. The test signal X\7(t) is the variable X\ of the driven oscillator, recorded as a continuous time series. For the first 400 s, the strength of the additional nonlinear term, i.e., 771 = 15 was very strong and the coupling, i.e., ju\ = 1 was relatively strong; ju\ vvas then substantially increased to 25, vvhereas the strength of the additional nonlinear term vvas decreased to 1. After a further 400 s, the strength of the additional nonlinear term vvas increased back to 15. The first 5 s and corresponding povver spectrum for each coupling strength are shovvn in Fig. 5.11 (a) and (b). The peak 1.1 Hz in the absence of coupling, labelled as/i, represents the driven cardiac oscillator; and^ = 0.24 Hz represents the driving respiratory oscillator. These frequencies are deliberately chosen such that their ratio is not an integer. The povver spectra, Fig. 5.11 (b), for ali the three different cases of the strengths of the linear coupling and additional nonlinear term exhibit rich frequency content. As the coupling gets stronger, and/or the strength of the additional nonlinear term increases, the frequency content of the signal X\? becomes richer. The povver spectra clearly exhibit components at the harmonically related positions/i +f2 and/i -fi- The principal domain of the bispectrum for the test signal jciF, Fig. 5.11 (c) and (d), shovvs a peak at the bifrequency (0.96 Hz, 0.24 Hz) that is of our primary interest. A vvindovv length of 100 s vvas chosen to calculate the instantaneous biphase and biamplitude, Fig. 5.11 (e) and (f), and vvas moved across the signal in 0.1 s steps. The vvhole signal is analysed as a single entity, but transients caused by the changes in coupling and/or in the strength of the additional nonlinear term are removed prior to processing. 18 5 Cardid-respiratdry interactidns Fig. 5.11: Digital simulation illustrating a situation where the bispectral method fails. The simulation models two unidirectionallv, very stronglv coupled, relaxation van der Pol oscillators with additional very strong nonlinear terms, in the presence of additive Gaussian noise. (a) The test signal XiF of variable X\ of the forced van der Pol oscillator with characteristic frequency f{ = 1.1 Hz periodicallv forced at frequency^/2 = 0.24 Hz for three different coupling strengths fix ~ 1 (1), 25 (2), and 25 (3), with strengths of the additional nonlinear term >ft = 15 (1), 1 (2), and 15 (3). Each coupling lasts for 400 s, at sampling frequency^ = 50 Hz. Only the first 5 s are shown in each čase. (b) The power spectrum of x1F. (c) Its bispectrum \B\ calculated with K = 34 segments, 67 % overlapping and using the Blackman window to reduce leakage and (d) its contour view. (e) The biphase (j) and (f) the biamplitude A for bifrequency (0.96 Hz, 0.24 Hz), with a 0.1 s time step and a 100 s long window for estimating the DFTs using the Blackman window, for cases (1), (2) and (3). Note that only for (3) are both conditions (i.e., high enough biamplitude, and constant biphase) for the (incorrect) inference of nonlinear coupling satisfied. During the period of relatively weak coupling ju\ = 1 and strong nonlinearities, no peak is present in the bispectrum as can be seen from the biamplitude, Fig. 5.11 (f), which remains far below unity (0.012). Moreover, at value t/i = 0.7 the frequency component at modulation position/i +f2 appears in the power spectrum for the first time. The modulation components f\ ±fi become large and almost equal in size in the power spectrum, but not until T]\ = 15. Hovvever, even then, not ali the necessary 19 ISTERACTiaNS peaks (also peak at bifrequency (/i, 2f2)) in the bispectrum are present and the method correctly resolves the absence of nonlinear coupling, even though the biphase is constant, Fig. 5.11 (e). By grossly exaggerating the strength of the coupling, with ju\ = 25, the frequency components that might also arise from nonlinear coupling become large. This results in an substantial increases of the biamplitude at bifrequency {f\,fi) in the bispectrum; the biphase is then non-constant, however, and increases continuously. Again the required conditions for the identiflcation of nonlinear coupling are not fulfilled. In the most extreme example shown, with very strong coupling and very strong additional nonlinearity, i.e., with rj\ = 15, fi\ =25, we are unable to distinguish betvveen strong nonlinearity of the oscillators and strong nonlinear coupling. In a signal coming from a "black box", the observed frequency components could mistakenly be attributed to nonlinear coupling: the bispectrum, Fig. 5.11 (c), contains ali the necessary peaks (only three of them are visible, the rest of them being much smaller, although fulfilling the necessary amplitude condition). There are other frequency components that could result from nonlinear coupling, and the biphase remains constant. In this čase the method clearly fails. There are, however, compelling arguments suggesting that the cardiac and respiratory subsystems should be in fact treated as weakly nonlinear oscillators that are weakly coupled. (i) In healthy subjects, breathing spontaneously, only occasional and ^r/^episodes of synchronization are seen [10, 99-101], indicative of relatively weak coupling. (ii) Sinus arrhvthmia is small at spontaneous breathing frequencies and only slightly larger at very low breathing frequencies [23], again supporting a weak-coupling description. (iii) The couplings can sometimes decrease almost to vanishing point, e.g., in coma [112]. Without couplings, the dvnamics becomes drastically simplified - with complete absence of svnchronization or modulation. The fact that virtually no variability is seen in any of the natural frequencies, despite small amplitude variations attributable to internal noise, suggests that the oscillators themselves are at most weakly nonlinear. (iv) If there were strong oscillator nonlinearity, and strong (but linear) coupling, we would observe many combinatorial components around the cardiac frequency, which is not the čase (see Fig. 5.1 (a) - (d)). The excessively strong-coupling regime explored in the above simulations would appear, therefore, to be largely irrelevant to the cardio-respiratory interaction that we study in this work: our bispectral technique [40] should be applicable as we have assumed. 50 5 Cardid-respiratdry interactidns 5.8.3 Relationship to synchronization The fact that inter-oscillator interaction can give rise to synchronization, as well as to modulation, has excited much attention in studies of the phase relationships between the cardiac and the respiratory oscillators [10, 24, 42, 46, 52, 71, 92, 95, 97, 100, 101, 113, 114, 118]. Indeed, it was the possibility of synchronization that motivated us to develop new techniques to investigate further the interactions between the systems: the direction, strength and, in particular, the nature of the couplings. They can be obtained from bivariate data (respiration and ECG signal) by use of the methods recently developed for analysis of synchronization, or generalized synchronization, between chaotic and/or noisy systems (see [79] and references therein). We now consider whether or not synchronization is onsets under conditions where we have clear evidence of interaction. Fig. 5.12 (d) shows a cardio-respiratory synchrogram based on the stroboscopic technique [79] for person 2 during paced respiration, supplementing bispectral results for signal ba. It is constructed by plotting the normalized relative phase of a heartbeat within m respiratory cycles [79] (5.4) where ?* is the time of &-th heart beat and (^ is the instantaneous phase of respiration. In perfect n\m phase locking, ^(/k) attains exactly the same n different values within each m adjacent respiratory cycles and the synchrogram consists of n horizontal strips. The instantaneous phase of the cardiac activity was obtained by characteristic or marker events - the R peaks in the ECG signal. A 2n increase of phase is attributed to the interval between subsequent R peaks. The instantaneous phase of respiration was obtained in a similar way, using zero-crossing as the marker event. For m = 1 we cannot see any horizontal structure that would resolve n:\ phase locking during 77-172 s. By differentiating the instantaneous phases we obtain the instantaneous (a) cardiac and (b) respiratory frequencies and (c) their ratio as shown in Fig. 5.12. In the histogram of frequency ratio (not shown) two peaks appear, one tX/j/fi = 9 and the other atfz/fi = 10, as can also be seen from the Fig. 5.12 (c). The synchrogram presented in Fig. 5.12 (e) shows the case where m = 9. The vertical inclined lines suggest that the frequency ratio is almost constant, whereas phase drift is present most of the time. No synchronization is evident in either synchrogram. It looks as though the cardiac oscillation has a tendency to synchronize with the respiratory one, but cannot tune due to the slow-paced respiration frequency. In other words, it appears that, for most of the time, the cardiac frequency 51 ITERACTIONS is just modulated by the respiratory rhythm, i.e., what is commonly referred to as respiratory sinus arrhythmia. Fig. 5.12: (a) Instantaneous cardiac frequency /j, (b) instantaneous respiratory frequency f2 and (c) their frequency ratio for person 2 during paced respiration, supplementing bispectral results for blood flow signal ba. (d) Synchrogram based on stroboscopic technique calculated for m = 1 and (e) m =9. Using the synchrogram technique we thus detect the existence of frequency modulation but no synchronization for the čase of paced, low frequency, respiration. We conclude, therefore, that bispectral analysis yields different information from that which can be resolved from a synchrogram. The relation to synchronization in more general way is going to be discussed in detail in the subsequent chapter. 5.8.4 Synchronization, modulation and type of coupling Just as frequency modulation does not necessarily coexist with synchronization between the cardiac and respiratory systems, there is not a one-to-one correspondence between quadratic coupling and modulation. As seen in Fig. 5.12 (a), the cardiac frequency is continuously modulated by the 52 5 Cardid-respiratdry interactidns respiratory frequency throughout the recording, whereas quadratic coupling is reliably detected only in the interval between 77 s and 172 s. Using bispectral analysis we can obtain additional information about the possible presence of modulation. Inspection of the tirne evolution of the amplitudes and phases for ali the frequencies that constitute the triplets for bifrequencies summarized in Tab. 5.3 has shown that the observed nonlinear coupling does not correspond directly to frequency modulation. It remains an open question whether the observed interchange and coexistence betvveen modulation, svnchronization and other manifestations of coupling arise through complexity of the cardio-respiratory interaction itself, or through indirect interactions involving the other oscillatory processes (see [110] and the references therein) known to exist in the cardiovascular system. 5.8.5 Unidirectional or bidirectional coupling We have studied the combinatorial frequencies, that arise from the influence of the respiratory on the cardiac system. This naturally begs the question of whether the coupling is unidirectional or bidirectional. Using the same method, one could equally well analyse the combinatorial frequencies responsible for the influence of cardiac on the respiratory system. In fact, using newly developed algorithms for analvsis of the direction of coupling [72, 93, 94, 103] has already shown [73, 117], that the two systems are bidirectionally coupled. 5.8.5.1 Forced oscillator The effect of respiratory system is, however, dominant (i.e., is the driving system) at ali respiratory frequencies, whether paced or spontaneous [73, 117]. The interaction betvveen the cardiac and respiratory oscillators can be seen as unidirectional: the respiratory system drives the cardiac one. A particular čase is the čase of the paced respiration. Although, during paced respiration, the respiration frequency is kept constant, the situation differs from that of a forced oscillator (with the cardiac oscillator being driven, and the respiration oscillator being the drive). Paced respiration experiments can in fact be perceived as a state of the system of two coupled oscillators, where, although the frequency of one of them (respiration) is forced and kept constant, the interaction betvveen the tvvo oscillators remains spontaneous. To illustrate how this happens, we use a generic model (5.5) of an almost periodic, Poincare oscillator periodically, driven by a weak external force J J INTERACTIONS (5.5) The activity of the oscillator is described by the two state variables x andy, a, a and co\ are constants and F is the forcing amplitude with frequency a>i and initial force phase $>. Here, Č{t) is zero-mean white Gaussian noise, <^(/)> = 0, {Č{t\ ^(0)> = Dd(f), and D = 0.1 is the noise intensity. The parameters of the model are set to a= 1, and a = 0.5. Fig. 5.13: Results for a forced Poincare oscillator in the presence of additive Gaussian noise. (a) The test signal *io, of variable x of the forced oscillator with characteristic frequency f\ = 1.1 Hz periodically forced at frequency^ =0.24 Hz with three different forcing amplitudes; F = 0.0 (1), 0.1 (2) and 0.2 (3). Each forcing lasts for 400 s, at sampling frequency^ = 10 Hz. Only the first 15 s are shown in each čase. (b) Its power spectrum. (c) The bispectrum \B\ calculated with K = 33 segments, 66 % overlapping and using the Blackman window to reduce leakage and (d) its contour view. The part of the bispectrum above^ > 1.0 Hz is cut, because the triplet (1.1 Hz, 1.1 Hz, 2.2 Hz) produces a high peak that is physically meaningless. (e) The biphase \ is constant, but only until 33 minutes when 4:1 svnchronization is detected in the svnchrogram. At this point, biphase (/h is no longer constant, and biamplitude A\ begins to decrease until it reaches 0, at which point, the svnchronization stops. Peak 3 is not present during this tirne. The data of Rat20 was specifically chosen, as there is constant frequency ratio from approximately 15 minutes to 30 minutes, Fig. 6.8 (c), but no svnchronization can be seen in the svnchrogram, Fig. 6.8 (d). Bispectrum resolves more information about the coupling. From the Fig. 6.11 (a), (b) and (d) there is high biamplitude during this tirne vvhereas (c) is low. The biphase for peak 1 is constant from 25 10 minutes onwards. This strongly suggests that modulation takes plače. In the following chapter this will be discussed in detail. Tab. 6.1: Quadratic nonlinear couplings detected in ratl6 and rat20 signals. Tqc is the time interval during which the bispectral analysis showed that the heart oscillator, f^, and the respiratory oscillator, fres, might be nonlinear coupled. The product of Tqc x fKS provides us with the amount of respiratory periods over which the interaction persisted. During Tqc, the maximum biamplitude is calculated for peak 1, which is of primary interest to us. In addition, the maximum variation of biphase A$ its average value (j) , and its standard deviation a^, were calculated during Tqc. It is not possible to discern from synchrogram whether the horizontal strips are due to synchronization or modulation. Moreover, both phenomena can overlap. Nevertheless, if the modulation is very strong, then the horizontal strips are not equidistant, and the modulation can be detected. 6.5.1 Synchronization and modulation Synchronization analysis, based on the mutual prediction approach [93, 94] and on information-theoretic functional [70, 73] for rat 16 and rat20, showed that during the synchronization episodes, the respiratory system is dominant (i.e., is the driving system) and drives the cardiac system [63]. In Sec. 5.8.5.1, we have demonstrated the čase of a forced oscillator illustrating the unidirectional interaction between the cardiac and respiratory oscillators, where the respiration frequency is kept constant. In this example, cardiac and respiratory oscillators are not synchronized. In Sec. 4.5, we showed a numerical example of frequency modulation. Again, there was no synchronization of the interacting oscillators. Besides the modulation, when two oscillators are interacting, whether unidirectionally or bidirectionally, synchronization can also onset. These phenomena are distinct, although they can overlap. There can be modulation vvithout svnchronization, synchronization vvithout modulation, or a combination of both effects. To illustrate how efficient bispectral analysis is in the latter two examples, we use a generic model (6.7) of an almost periodic, Poincare oscillator, periodically driven by a weak external force, as in case of the model (5.5) with additional frequency modulation (6.7) The activity of the oscillator is described by the two state variables, x andy. a, a and a>\ are constants. F is the forcing amplitude with frequency a^ and initial force phase fa and i]m is the strength of modulation by the forcing oscillator. Here, <£{f) is zero-mean white Gaussian noise, <^(/)> = 0, (£{i), <^(0)) — DS(J), and D = 0.2 is the noise intensity. The parameters of the model are set to a= 1, and a = 0.5. The test signal xm(t) is the variable x of the driven oscillator, recorded as a continuous time series. For the first 400 s, there was no forcing, i.e., the amplitude was set to zero. It was then raised to a small constant value 0.1, without frequency modulation. After a further 400 s, the forcing was increased to 0.2 and the modulation strength was set to 0.2 (moderate). The corresponding power spectrum for the first 15 s and for each forcing strength are shown in Fig. 6.12 (a) and (b), in order to demonstrate the changes in spectral content and behaviour caused by the forcing. The peak labelled as f\ = 1 Hz represents the driven cardiac oscillator, and the peak labelled f2 = 0.2 Hz represents the driving respiratory oscillator. These frequencies are deliberately chosen to have an integer ratio 5:1. In case of n:\ locking, the effect of the forcing can be twofold. It causes modulation on the period of the oscillator that occurs with the period of the forcing, and the force adjusts the average period of oscillations, i.e., synchronization. Synchrogram for the test signal xm, Fig. 6.13, exhibits 5:1 synchronization for the whole signal duration. In the first 400 s, there is no interaction, forcing or modulation on the oscillator. The frequencies and the force of the oscillator are constant and in integer ratio, which is the reason for the synchronization appearing in the synchrogram. In general, one should be cautious when interpreting synchrogram, as it can be misleading. For this case, the bispectrum will be completely flat, without the appearance of any peaks, as can be seen from Fig. 6.12 (f). From 400 s to 800 s, weak forcing is present. External force tries to change the amplitude as well as the phase of the oscillation. The amplitude is stable, whereas the phase is neutral (it is neither stable nor unstable). Weak force influences the oscillator phase that results in synchronization. A similar case without synchronization was already discussed in Sec. 5.8.5.1 Adapted bispectrum resolves that linear interaction takes place. If biphase is constant for the bifrequency (f\,f-i), then we can conclude that I 12 S BlSPECTRAL RELATIDN TD SYNCHRONIZATIDN they are phase coupled. This results in synchrogram in horizontal strips - synchronization. If the frequency ratio was rational, then horizontal strips would not appear in the synchrogram. Thus, the bispectrum yields the correct information about the coupling. X! -P Time (s) Time (s) 1150 Fig. 6.12: Results for a forced and frequency modulated Poincare oscillator in the presence of additive Gaussian noise. (a) The test signal jtiH, of variable x of the forced oscillator, with characteristic frequency /i = 1 Hz periodicalh/ forced at frequency^ = 0.2 Hz, with three different forcing amplitudes; F= 0.0 (1), 0.1 (2) and 0.2 (3) and three different modulations ; rjm = 0.0 (1), 0.0 (2) and 0.2 (3). Each forcing lasts for 400 s, at sampling frequency^ = 10 Hz. Only the first 15 s are shown in each čase. (b) Its povver spectrum. (c) The bispectrum \B\ calculated with K = 33 segments, 66 % overlapping and using the Blackman window to reduce leakage, and (d) its contour view. (e) The biphase (p, (f) biamplitude A for bifrequency (1,0 Hz, 0,2 Hz) peak 1, (g) biphase $ and (h) biamplitude A for bifrequency (1.2 Hz, 0.8 Hz) peak 6, with a 0.1 s tirne step, and a 100 s long window, for estimating the DFTs using the Blackman window. In the last 400 s of the test signal xm, moderate forcing and moderate modulation take plače. The synchronization is preserved, as can be seen from the synchrogram, Fig. 6.13. Combination of forcing and modulation could be misleading in detecting nonlinear interaction using the bispectrum as ali the frequency components, except for the appearance of 2/1 (the 2/1 -f2 is present) in the power spectrum, Fig. 6.12 (b) (3). It is not as evident as in the čase presented in chapter 4.7, vvhere only the modulation 13 SYNCHRQNIZATiaN Fig. 6.13: Synchrogram for test signal xm. takes plače, but one can determine it by calculating ali the necessary peaks for the nonlinear coupling. Biamplitude for the peak of primary interest is very high and the biphase is constant, Fig. 6.12 (e) and (f), whereas for peak 6, it is not, Fig. 6.12 (g). Also significant for the modulation, is the appearance of very high peaks (1 and 2), compared to the others (3-6). From the constant biphase for peak 1, we can conclude that the interacting oscillators are in phase and synchronized. If the phase would not be that constant, then it would be difficult to say that they are synchronized, as it was seen in the čase of rat20, where only modulation was detected, since the biphase of peak 2, Fig. 6.11 (b), is not similarly tirne dependant at that tirne, as is the čase for peak 1. 6.6 Conclusions For the čase of synchronization, we can conclude as follows: • Strong synchronization. Synchrogram and bispectrum (cross-bispectrum) provide us with the same results. We only inspect the peak at bifrequency of our primary interest. When there is a strong synchronization betvveen two interacting oscillators, then clear horizontal strips appear in synchrogram indicating synchronization. In bispectrum, the synchronization is indicated with a high biamplitude value and constant biphase at bifrequency of our primary interest (/], h\ • Weak synchronization. Horizontal strips in synchrogram can hardly be detected or they cannot be detected. Bispectrum results in moderate biamplitude and less constant biphase, with more phase slips at bifrequency of our primary interest (fufi)- No synchronization. The synchrogram contains no horizontal strips. Bispectrum results in zero biamplitude at bifrequency of our primary interest (fufi). Although there is no coupling, synchronization can onset in synchrogram, due to constant frequency ratio. 11 S BlSPECTRAL RELATIDN TO SYNCHRDNIZATIDN As synchronization can take plače simultaneously with frequency modulation/forcing/nonlinear coupling, we can further conclude: • • • Synchronization and nonlinear coupling. Nonlinear coupling can appear while synchronization takes plače, whereas while there is nonlinear coupling, synchronization does not necessarily onset. There is no obvious link between the two phenomena. Analysis of rats undergoing anaesthesia shows that nonlinear coupling occurs during synchronization, while analysis of CV blood flow signals of humans in resting shows nonlinear coupling when synchronization is not present. Frequency modulation. Frequency modulation can be detected using bispectrum. The peak of primary interest (/j, f2), and the second peak at bifrequency (/j - f2, f2), are high in comparison to other peaks (that may not be present at ali) that appear in the čase of nonlinear interaction. Their biphase is constant if the frequency modulation is strong. Frequency modulation and forcing. Instantaneous presence of both phenomena can be misleading in detecting the nonlinear coupling. It is necessary to check ali the peaks, 1-6. Biphase for peaks 4 and 6 is not constant. It is recommended to analyse the peak of primary interest (fi,f2) for adopted bispectrum. Biphase should be constant. In this way, it is possible to resolve this kind of interaction between two oscillators. Nevertheless, one should be careful in interpreting the bispectrum results when strong frequency modulation and strong forcing in the presence of strongly noisy data take plače. Observing the phases of each frequency component in the triplet can be helpful. Frequency modulation and nonlinear coupling. When strong frequency modulation and nonlinear coupling take plače simultaneously, it is not possible to detect modulation. When the frequency modulation is weak, this can be seen as undulated biphase at peak of primary interest (/i, f2). It is difficult to be sure, as it could also be the čase of a weak nonlinear coupling. We conclude, therefore, that bispectral analysis is more sensitive to interactions and is more noise robust than the svnchrogram. It detects the phase svnchronization, and nevertheless, yields different information from that which can be resolved from a svnchrogram. Frequency modulation interaction can be detected, vvhereas it is not always possible to resolve it if it simultaneously occurs with other types of interactions. • 15 V HlBH-ORDER SPECTRA BASED DN VVAVELET TRANSFDRM V HlGH-ORDER SPECTRA BASED DN VVAVELET TRANSFDRM The Fourier transform is based on presumptions (a) of the periodicity of the signal and (b) of infinitely long signal series [57, 58]. Because neither assumption is strictly true for the measured signals, the determination of separate frequencies in a system that possesses strong couplings is very demanding. The difficulty is even greater in the low frequency range, which is of our particular interest, where the characteristic frequencies are close to each other, and are therefore even harder to separate. The uncertainty principle of the Fourier transform limits its ability to separate harmonic components in the frequency domain of the bispectrum [20, 43]. This might cause problems for detection of quadratic phase couplings in the čase of frequency pairs that are close together. To ensure good resolution of low frequencies, we need longer sections for calculation of the discrete Fourier transform. This immediately decreases the number of sections possible and weakens the bispectrum estimation. Hovvever, we cannot use longer signals, because they lead to nonstationarity, and the variance consequently becomes even larger [69]. 7.1 Wavelet Transform Wavelet analysis can be seen as a generalization of the Fourier analysis [43] by adding tirne resolution - in a more fundamental way than is permitted by the Short-Time Fourier Transform (STFT) [81]. Wavelet analysis has been applied with considerable success to cardiovascular data [6, 8]. The generalization of bispectrum to vvavelet analysis may be expected to be able to detect temporal variation in phase coupling or short-lived couplings, and čope with broadened and coalescing peaks TI IASED ON VVAVELET TRANSFORM that cannot be resolved due to time-frequency resolution restrictions by the bispectrum based on STFT. Morlet first introduced wavelet analysis [28]. Within this transform, the window length is adjusted to the frequency currently being analysed. It is a scale independent method. Window function is called a mother wavelet or basic wavelet yAji). It can be any function ifAii) that satisfies the wavelet admissibility condition [43] (7.1) This function introduces a scale s (its width) into the analyses. Commitment to any particular scale is avoided by using all possible scaling of ifiu). The mother wavelet is also translated along the signal to achieve time localization. Thus, a family of generally non-orthogonal basis function is obtained [43] The wavelet transform Wg(s, i) is a mapping of the function g{t) onto the time scale plane. Not every function can be used as the mother wavelet. Only those that enable us to reconstruct the original function g(t) from its wavelet transform Wg(s, t) are admissible. The inverse continuous wavelet transform is defined as [43]: (7.4) where the constant C is determined by the shape of the mother wavelet [43] 18 V HlGH-DRDER SPECTRA BA5ED DN VVAVELET TRANSFORM The function can be interpreted as the energy density of the signal in the tirne scale plane, also called a scalogram. Applying Parseval identity {e, a) = (e, a) into Eq. (7.4) we obtain (7.9) (7.10) (7.11) The wavelet transform provides a multiplying constant and phase shift Q]27lft information about g inside the window that is determined by instantaneous scale and shape of the mother wavelet. 19 A BASED ON VVAVELET TRANSFORM 7.1.1 Discretization In numerical applications, scale s and time / are restricted to discrete values only. The natural discretization of the scaling parameter is sm = cf, where m e Z, and the step is a positive number 0.8, the value of the second term in (7.14) is so small that it can be ignored in practice, and a simplified expression for the Morlet wavelet in the time domain is [43] (7.15) The corresponding wavelet family consists of Gaussians, centred at a time / with standard deviation s. In the frequency domain, we have Gaussians with a central frequency f=fjs and a standard deviation of Mljlm. Therefore, the wavelet transform at a given scale s can also be interpreted as band-pass filtering, giving an estimation of the contribution of the frequencies in this band. The relation between the scale and the central frequency for the Morlet wavelet is [43] (7.16) The frequency resolution changes with frequency; at low frequencies (large scales), the resolution is better than at the high frequencies (small scales). Accordingly, the time resolution is better for high frequency than it is for low frequency components. In order for peaks to be detected at/i and/ (fi > /), they must be separated by at least one half of the standard deviation of the peak at the higher frequency, namely f\ - f2> fi/4nf0. The choice of/ determines the current frequency resolution. By choosing/ = 1, a simple relation between scale and frequency was obtained/= lis. To obtain the energy density in the time-scale plane, an approximation of the continuous wavelet transform was calculated using the Morlet mother wavelet discretized with a = 1.05 and r = 1 s. However, to make the three-dimensional plots of the transform clearer, time was not discretized as tn = ncfr, but tn =«rwas used instead. In this way, the transform is over sampled in time for large scales. Slow events are examined with a long window, whilst a shorter window is used for faster events, Fig. 7.1. The Morlet wavelet [28], a Gaussian window, i.e., a Gaussian function modulated by a sin wave, is used. Thus, for our purpose, the best time-frequency localization within the limits of the uncertainty principle can be achieved. For details see [43]. For the Morlet mother wavelet, the value C, Eq. (7.5), equals C= 1.0132 [8]. mmmm 81 IC3H-DRDER SPECTRA BASED ON VVAVELET TRANSFDRM Fig. 7.1: (a) Real (black line) and imaginarv (grey line) part of the Morlet mother wavelet for scale s = 1 and (c) for scale, 5 = 5, (b) and (d) its Fourier transforms. In both cases,y^ = 1. 7.2 Wavelet bispectrum definition The definitions are completely analogous to the definitions used in Fourier analysis [67]. The Wavelet Bispectrum (WB) is given by where (7.17) (7.18) The WB measures the amount of phase coupling in the interval T that occurs between wavelet components of scale lengths s\ and s2 and s of signal g(t), such that the frequency sum-rule is satisfied (7.18). It is a complex quantity, defined by magnitude^ and phase (f) (7.19) 82 7 HlGH-DRDER SPECTRA BASED ON VVAVELET TRANSFDRM Consequently, for each ($1, s2), its value can be represented as a point in a complex space, 9{[WB(s\, s2)] versus 3[Wi?(si, s2)], thus defining a vector. We define its magnitude (length) as biamplitude, and the phase, which is determined by the angle between the vector and the positive real axis, as biphase. The instantaneous biphase is then calculated: from Eqs. (7.17) and (7.19), it is (7.20) If two scale components S\ and s2 are scale and phase coupled, such that the frequency sum-rule is satisfied (7.18). For ease of interpretation, the WB is plotted in the (/i,/2)-plane, rather than in the ($i, s2)-plane. It has the same symmetries in frequency domain as in the case of Fourier based Bispectrum (FB). The non-redundant region is the principal domain of the wavelet bispectrum. Similarly, the principal domain can be divided into two triangular regions in which the wavelet bispectrum has different properties: the inner triangle (IT), and the outer one. The IT is of our interest. IGH-DRDER SPECTRA BASED ON VVAVELET TRAN5FDRM 7.2.1 Wavelet bispectrum transform adopted to CV signals Relation between frequency (scale) and the width of the window that is used for calculation of the wavelet transform, is hyperbolic. Logarithmic-logarithmic scale is natural for its presentation, whereas, to be able to comply with the (7.18) frequency (scale) sum-rule, we need to achieve better frequency (scale) resolution for high frequencies (low scales) (according to the CV frequency bands proposed in [112]), as can be achieved using the Morlet wavelet when/0 is chosen to be 1 for the reason of the interpretation of the wavelet bispectrum. Otherwise, nearby peaks at high frequencies cannot be resolved. We introduce a parameter d into the Morlet wavelet that determines the exponential decay of the Gaussian (7.23) This also decays the Morlet wavelet, and thus permits suitable combination of time and frequency (scale) resolution to be selected. The time resolution is At = sd, given by the decay of the exponential part of the wavelet. As d increases (d > 1), frequency (scale) resolution improves, whereas time resolution deteriorates. We do not impose the condition that the wavelets must be orthogonal, as we wish to choose the frequencies in the analysis procedure freely, and not restricted to s e {2n}. This implies a certain redundancy in the wavelet transform coefficients, which must be taken into account upon interpreting the results. The parameter d is calculated so that the Gauss function decays to 0.001 for each scale (d is between 2.5 and 2.6). A high value of d causes a non-zero value of Morlet window at its edges that results in side lobes in wavelet bispectrum. If d would be infinite, than Morlet window would become a unit window, and wavelet transform would become Selective Discrete Fourier Transform (SDFT) [47]. Parameters am and cn are discussed in the following text. Frequency resolution for high frequencies is yet insufficient. It is necessary to increase the length of the Morlet wavelet for high frequencies. This can be obtained in different ways. Fig. 7.2 shows the hyperbolicous decay of the Morlet wavelet length with the increasing frequency (solid line). The wavelet length can be multiplied by a factor am, that is for the one with the lowest frequency of interest, and then increases with the increasing frequency 84 V HlGH-ORDER SPECTRA BASED DN VVAVELET TRANSFDRM am = 2/max /mm , (7.24) where/= l/s is the frequency of observation, and/min andfmax define the frequency range of interest. In this way, we obtain a dotted line on Fig. 7.2. As the wavelet length is prolonged for high frequencies, the frequency resolution increases, whereas the tirne resolution deteriorates. Other ways to obtain the necessary frequency resolution is by using a fixed wavelet length for ali high frequencies - Fig. 7.2, dash-dot line. We propose to use a 20-80 s long vvavelet in the čase of analysing the human cardiovascular blood flow signal of a normal, healthy subject at rest. Fig. 7.2: Length of the Morlet wavelet 7VW depended from frequency (scale). Morlet wavelet (solid line), adapted Morlet wavelet (dotted line), and flxed vvavelet length for high frequencies (dash-dot line). WB estimation using the proposed Morlet vvavelet as a mother vvavelet encounters a normalization problem. For each scale, a vvindovv of different length is used. In the čase of a signal composed of different frequency components, but equal Fourier povvers, this vvould result in different vvavelet spectral energies for separate frequencies. Tvvo couplings among different frequencies with the same Fourier povvers and the same nature of coupling, vvould result in different coupling strength in vvavelet I 1—1/2 bispectrum. In (7.2), a factor \s\ is used to ensure energy preservation. We choose to use a factor l/iVw instead, vvhere 7VW is the Morlet vvindovv length. Constant cn, Eq. (7.23), equals to 3.9487- n~ In this way, vve can compare results obtained by FB and WB, since both preserve energy. Normalization of the WB is applied in the same way as on the FB, discussed in Sec. 3.3. The normalized WB indicates the average level of quadratic nonlinear phase coupling and, in a way, serves as an indicator of hovv non-Gaussian the signal is [31]. The critical values for the WB and biamplitude estimates vvere normalized to 1. If the estimated value is higher than the average value of WB in the IT, then it is taken as valid. By critical value, it is meant that a value exceeds the noisy background (other than Gaussian), rounding, and estimation errors. 85 iASED ON VVAVELET TRANSFDRM 7.3 Wavelet bispectrum example of test signal Results of WB are illustrated on a numerically generated test signal of two Poincare oscillators quadratically coupled in the presence of additive Gaussian noise - Eq. (4.6), Sec. 4.4, test signal xm, presented in Fig. 4.7 (a). Test signal jciD was already analysed using Fourier based bispectrum in Sec. 4.4. Results obtained using WB, depend on the choice of parameters set for WB calculation. In the first step, the WB is calculated for the whole bifrequencies domain. Similarly, as in the čase of FB, one has to set the parameters: K - number of segments into which the signal X\d is divided to try to obtain statistical stability of the estimates; O - percentage of segments overlapping; and L - L ■ L area for bispectrum frequency averaging. These parameters have already been discussed in detail in [39]. In čase of FB, one chooses tapering window (Hamming, Hanning, Blackman or other), vvhereas WB ušes the Morlet mother vvavelet. Parameters that can further be chosen are: ■ Tm - Morlet mother vvavelet length. By choosing^o = 1, a simple relation betvveen scale and frequency is obtained: / = \/s. Mother vvavelet for/= 1Hz is then stretched and compressed. The prevailing choice of the Tm is from 8 s (i.e., ± 4 s) to 12 s [6, 7]. ■ d - exponential decay of Gaussian function of Morlet vvavelet. Rather than setting the parameter d, we set Gaussian function of Morlet vvavelet edge value to Ge, so that it decays to some small value. The prevailing choice is from 0.01 to 0.0001 (see Sec. 7.2.1 for details). ■ A/- frequency (scale) step. It can be chosen arbitrarily, vvhereas, to be able to comply vvith the (7.18) frequency (scale) sum-rule, the prevailing choice is to be at least 1/10 of the slovvest frequency in the bifrequency pair of our interest. ■ am - multiplication factor for additional Morlet wavelet stretching (see Sec. 7.2.1 for details). One can set basis, and the constant in the povver. Morlet vvavelet length is multiplied by a factor am. The factor equals the one of lovvest frequency of interest and then increases vvith the increasing frequency. Either one chooses to use factor am vvhose prevailing choice for basis is 2, and the constant in the povver is 1.8, or one chooses to fix the Morlet mother vvave length for high frequencies Tw. The prevailing choice of the fixed vvindovv length for higher frequencies is from 20 s to 80 s. One should start vvith a 40 s long vvindovv. If the peaks at the bifrequencies of our interest are distinct than shorter fixed window can be used, othervvise longer fixed vvindovv must be used. 86 7 HlGH-DRDER SPEGTRA BASED ON VVAVELET TRAN! Once WB is obtained, and longer lasting (bispectral averaging over K segments eliminates short lasting couplings) phase and/or frequency coupling are detected, biphase and biamplitude tirne evolution are estimated for bifrequencies of our interest (see Tab. 5.2 for details). Parameters K, O and L do not influence the estimation. Additional parameters can be chosen: ■ At - tirne step. Minimum tirne step is defined with sampling frequency fs and equals Mnin=l//^ It should be set to such a value that epochs of constant biphase of approximately 10 times the slovvest period (1/ f2) of the bifrequency pair (fufi) can De detected, i.e., at least 1/10 of the epoch length. ■ Lt - number of samples for tirne averaging. In the čase of signals with non-Gaussian noise, one can reduce the noise by averaging in the tirne domain over Zt samples, before and after the tirne of observation. WB for the test signal x\D is presented in Fig. 7.3 (a). Parameters are set according to prevailing choice for CV signal analysis as described above (Tm = 8 s, Ge = 0.001, Af= 0.01 Hz, TH¥ = 40 s and At = 0.1 s). Averaging is neither used in frequency, nor in tirne domain (L = 0, Zt = 0). Peaks appear at bifrequencies (1.1 Hz, 0.24 Hz), (0.86 Hz, 0.24 Hz), (0.62 Hz, 0.48 Hz), (0.86 Hz, 0.48 Hz), (1.1 Hz, 0.48 Hz), (1.1 Hz, 0.86 Hz) and (1.34 Hz, 0.86 Hz). Fig. 7.3: Results for quadratic couplings in the presence of additive Gaussian noise, test signal xw, obtained with the wavelet bispectrum for comparison with the Fourier bispectrum, See. 4.4. (a) The wavelet bispectrum \WB\ calculated with K = 33 segments, 66 % overlapping, Tm = 8 s, Ge = 0.001 and using fixed Morlet wavelet length of 7hf = 40 s for high frequencies calculation and (b) its contour view. The part of the wavelet bispectrum above fi > 1.0 Hz is removed, because the triplet (1.1 Hz, 1.1 Hz, 1.1 Hz) produces a high peak that is physically meaningless. (c) The biphase (j) and (d) biamplitude A for bifrequency (1.1 Hz, 0.24 Hz), with a 0.1 s tirne step. 81 As before, the self-coupling peaks at (1.1 Hz, 1.1 Hz) and (0.24 Hz, 0.24 Hz) are of no interest, so they are removed from the wavelet bispectrum. We obtain the same information as with the FB, Fig. 4.7 (c). The most obvious difference is in the shape of the peaks. They are much wider than in the čase of FB. This is expected, since frequency resolution for high frequencies is lower than in the čase of FB. Fig. 7.3 (c) shovvs biphase and (d) biamplitude for the peak of our principal interest at bifrequency (1.1 Hz, 0.24 Hz). We obtain the same information as with FB, Fig. 4.7 (e) and (f). The biphase is constant in the presence of quadratic coupling. Coupling strength can be determined from the biamplitude by normalization. The WB possesses the FB conceraing the noise robustness. The results for non-zero coupling are quite different from those where coupling is absent, Fig. 7.3 (d). From the biphase tirne dependence, it can be seen that the WB is better at detecting biphase changes, since its tirne resolution is higher than in the čase of FB. 7.4 Discussion In the čase of bispectral analysis of cardiovascular interaction, the time-dependent biphase/biamplitude estimate was estimated with an STFT, using a window of constant length. The optimal window length depends, hovvever, on the frequency being studied. The effective length of the window used for each frequency can be varied by applying the wavelet transform. If natural frequencies of the oscillators lie within a relatively narrow frequency interval, then STFT is sufficient for good tirne and phase/frequency localization. With broader frequency content, however, the vvavelet transform, or selective discrete Fourier transform, needs to be applied. The vvavelet and cross-wavelet bispectrum was defined analogous to the deflnitions used in Fourier based bispectrum and cross-bispectrum. By doing this time-dependant biphase/biamplitude estimate with higher frequency resolution at low frequencies, higher tirne resolution at higher frequencies was obtained. The vvavelet bispectral analysis was adopted for analysing cardiovascular signals. For a mother vvavelet modulated Gauss function, the Morlet mother vvavelet vvas used. The vvavelet bispectral analysis vvas illustrated on a test signal. Since the tirne resolution of vvavelet bispectrum is higher, and the frequency resolution is poorer at high frequencies compared to FB, it is necessary to ensure sufficient frequency resolution before interpretation of the results. Poor frequency resolution vvould result in poor/incorrect localization of characteristic frequencies. Too high of a tirne resolution could result in extremely high sensitiveness to noise and statistical error, that vvould result in phase slips and incorrect oscillator coupling determination. It vvas necessary to raise frequency resolution for high frequencies, as vvell as to preserve the scale (frequency) sum condition necessary 88 V HlGH-DRDER SPECTRA BASED ON VVAVELET TRANSFDRM for bispectrum estimation. Wavelet bispectrum results are parameter set dependant. Parameter impact on wavelet bispectrum estimation and detailed comparison with Fourier based bispectrum are discussed in the subsequent chapter. 89 B Fdurier and vvavelet bispectrum cdmparisdi B Fdurier and vvavelet bispectrum cdmparisdn In the follovving section, wavelet bispectrum is compared in detail with Fourier based bispectrum. First, vvavelet bispectrum is applied to CV blood flovv signals that were alreadv used for studving cardio-respiratorv interactions using the Fourier based bispectrum method in Sec. 5. Its benefits and vveakness over the Fourier based bispectrum method are then compared in detail and discussed. 8.1 VVavelet bispectrum of CV blood flovv signals For the comparison of WB and FB, we illustrate results obtained for the blood flovv signal ba(t), shovved on Fig. 5.1 (a), left column in Sec. 5, used for cardio-respiratorv interactions analvsis. Data analvsis vvas preformed in the same manner as in Sec. 5, vvhereas instead of using FB, WB vvas used. Parameters for WB evaluation vvere set according to the prevailing choice recommended in section 7.3. 8.1.1 Results WB for the vvhole frequency domain for signal ba(t) is presented in Fig. 8.1 (a). A very high peak located at bifrequency (0.11 Hz, 0.11 Hz), belonging to the respiratorv self-coupling, can be seen in the \WBha\, Fig. 8.1 (b). At least three other peaks are clearlv evident: at (0.98 Hz, 0.11 Hz) attributable to cardio-respiratorv coupling; at (0.87 Hz, 0.11 Hz), vvhich vve presume to be coupling betvveen the 91 BISPECTRUM CDMPARISDN respiratory components and the difference f\ - f2; and peak attributable to interaction with lovver cardiovascular characteristic components. Their positions can be seen in the | WBha\ contour view shown in Fig. 8.1 (c). Close inspection of \lVBha\, Fig. 8.1 (b), resolves that ali the necessarv peaks, according to Tab. 5.2, arise as a possible result of a nonlinear interaction betvveen the two oscillators f\ and f2 are present. Characteristic frequency at ~0.98 Hz belongs to cardiac activitv, f\\ that at ~0.11 Hz to respiratory activity, f2, Peaks 1 to 6 are presented in Fig. 8.2, left column. Time evolution for biamplitude and biphase for ali the peaks are shown in Fig. 8.2, mid and right columns. Fig. 8.1: (a) The wavelet bispectrum \WBha\ for signal bjj), calculated vvith K= 33 segments, 87 % overlapping, Tm = 8 s, Ge = 0.001, and using a THF = 80 s long fixed Morlet wavelet for estimating high frequencies. (b) Part of the wavelet bispectrum/j, ^ < 1.4 Hz that is of our interest, and (c) its contour view. The tirne interval Tqc, during which quadratic coupling persisted, was determined. If ali 6 peaks fulfilled our conditions (see Sec. 5.3 for details), then the Tqc interval was calculated for ali peaks, and the boundaries were defined such that the biamplitude for ali the peaks in Tqc interval would be above the condition. It can be seen that the biamplitude during the tirne interval from 77.1 s to 170.4 s meets our criterion of being more than tvvice as large as the average wavelet bispectrum in the IT domain; 92 B FDURIER AND VVAVELET BISPECTRUM CDMPARISDi Fig. 8.2: Wavelet bispectrum results for blood flow signal bt(t), calculated with K = 33 segments, 87 % overlapping, using a 0.1 s tirne step, Tm = 8 s, Ge = 0.001, and THV = 80 s long fixed Morlet wavelet for estimating high frequencies for peaks (a) 1, (b) 2, (c) 3, (d) 4, (e) 5 and (f) 6; left column, the wavelet bispectrum | WBba| with its corresponding contour plots; middle, the biamplitude Aba; and right, the biphase = Z)<5(/), and D = 0.08 is the noise intensity. The parameters of the model are set to ct\ = 1, ct\ = 0.5 and a2, a2 = 1. In this way, we obtain a test signal Xu(t), presented in Fig. 8.4 (a), with the corresponding power spectrum for two different coupling strengths, which are interchanging every 20 s: no coupling r/2 = 0; and weak coupling rj2 = 0.2, Fig. 8.4(b),(l)and(2). 100 120 200 Time (s) 300 100 120 200 Time (s) 300 Fig. 8.4: Results for tirne intermittent quadratic couplings in the presence of additive Gaussian noise analysed using Fourier bispectrum. (a) The test signal jch, variable xx of the flrst oscillator with characteristic frequency/i = 1.1 Hz. The characteristic frequency of the second oscillator is f2 = 0.24 Hz. The oscillators are unidirectionally and quadratically coupled with two different coupling strengths: r]2 = 0.0 (1); and 0.2 (2). The coupling (2) is present every 20 s and lasts for 20 s. The signal is 1200 s long and sampled with sampling frequency^ = 10 Hz. Only the first 15 s are shown in each čase. (b) Its power spectrum. (c) The bispectrum \B\ calculated with K = 33 segments, 66 % overlapping and using the Blackman window to reduce leakage and (d) its contour view. The part of the bispectrum above^ > 1.0 Hz is cut, because the triplet (1.1 Hz, 1.1 Hz, 1.1 Hz) 91 T BISPECTRUM CDMPARISDN produces a high peak that is physically meaningless. (e) The biphase and (f) biamplitude A for bifrequency (1.1 Hz, 0.24 Hz), calculated with a 100 s long window for estimating DFTs and (g) biphase ^ and (h) biamplitude A, calculated with a 130 s long window for estimating DFTs. In both cases, a 0.3 s time step and the Blackman window was used. Fig. 8.5: Results for time intermittent quadratic couplings in the presence of additive Gaussian noise analysed using vvavelet bispectrum. (a) The wavelet bispectrum |W5| calculated with K = 33 segments and 66 % overlapping. The part of the bispectrum abovef2 > 1.0 Hz is removed, because the triplet (1.1 Hz, 1.1 Hz, 1.1 Hz) produces a high peak that is physically meaningless. (c) The biphase and (d) biamplitude A for bifrequency (1.1 Hz, 0.24 Hz), calculated with a Ge = 0.01 and (e) biphase and (f) biamplitude A, calculated with a Ge = 0.0001. In both cases, a 0.1 s time step, Tm = 8 s and Tn¥ = 20 s long fixed Morlet wavelet for estimating high frequencies was used. Fig. 8.5 (a) shows wavelet bispectrum, and Fig. 8.5 (b) shows its contour view, obtained for signal Jtn. It gives the same information about the peak's relative amplitude and bifrequency position. Wavelet bispectrum parameter set, (THF = 100 s, Tm = 8 s and Ge = 0.01), was deliberately chosen in such a way that time resolution for high frequencies was increased. Results of biamplitude and biphase estimates are presented in Fig. 8.5 (c) and (d). Biamplitude clearly exhibits episodes when coupling is present and when it is not, i.e., every 20 s. From biphase time evolution, episodes of constant biphase can be seen. Episodes of constant biphase, i.e., without phase slips, last longer than the coupling between the 98 a Fourier and wavelet bispectrum cdmparisdi oscillators, i.e., more than 20 s. By taking into account the biamplitude condition for quadratical coupling occurrence, (Sec. 5.3), it is possible to determine the correct quadratic coupling time persistence. Increasing time resolution of wavelet bispectrum even higher, (rHF " 100 s, Tm = 8 s and Ge = 0.0001), we obtain an even more realistic result. Whereas biamplitude, Fig. 8.5 (f), is estimated to be poorer due to lower frequency resolution, but still preserving true information about the time of quadratic coupling occurrence. The biphase, Fig. 8.5 (e), correctly exhibits 20 s of constant biphase episodes, and then 20 s of constantly growing biphase episodes. 8.2 Fourier and wavelet bispectrum advantages and weakness Time and frequency resolution. To observe a given frequency, the signal must be observed over at least one period of this frequency what excludes the time localization. Due to the Heisenberg uncertainty principle [43] sharp localization in time and frequency are mutually exclusive (8.2) where At is time interval and Af frequency band. The equality holds if and only if the window is Gaussian. They are defined as [43] where w is in general some window function and ||w||2 is its norm. For STFT the representation of function g(u) in time-frequency plain G(t, j) has not sharp time and frequency parameters, but represents an interval around centre time t or frequency/. Time-frequency window is (8.4) For wavelet transform the representation of Wg(s, t) in time-frequency plain is (8.5) 99 where to and To are mother vvavelet centres of gravitv in tirne and frequency plain and corresponding deviations of At0 and A/o- Note that the centre of the tirne window depends only on parameter /, vvhereas the centre of the frequency window depends only on parameter s. On the contrary to the STFT the vvavelet transform's frequency resolution changes vvith frequency (lovv frequencies have better frequency resolution) and so does the tirne resolution (high frequencies have better tirne resolution). The ratio betvveen centre frequency f(s) = fo /s and bandvvidth AJ{s) = Afo/s is equal to fo /Afo and is scaling independent. In general WB detects intermittent phase couplings vvhereas FB averages out most of the tirne relevant information. Triplet (/i, fi,fs) results in a high peak in bispectrum if the coupling condition f$ —f\ +f2, is satisfied. Nevertheless, the coupling condition needs to be satisfied only vvithin the frequency resolution. This condition is less strict in vvavelet bispectral analysis. For example if there is a mismatch in coupling frequency A/ such that f3 = f\ + f2 + Af, and A/ is larger than the frequency resolution of the Fourier bispectrum but smaller than the vvavelet bispectrum frequency resolution corresponding to a high bispectrum value for the triplet (/i, f2, /3), than the vvavelet bispectrum will peak vvhereas the Fourier bispectrum vvill not. Increasing the frequency resolution of vvavelet bispectrum by increasing the length of the Morlet vvavelet for high frequencies results in gradually approximate results as obtained vvith the Fourier bispectrum. On the other hand if there is a short lasting coupling present in the signal the Fourier bispectrum cannot detect the coupling due to large tirne vvindovv used, vvhereas the vvavelet bispectrum vvill detect the coupling if assuming the coupling has a certain minimum duration. Wavelet bispectrum allovvs intermittent couplings to be detected. Applying Fourier bispectrum to real data vve have to ensure the necessary frequency resolution to be able to distinguish separate frequency components and at the same tirne achieving sufficiently tirne resolution to be able to detect the onset of the couplings among CVS oscillators. The scope for choice of vvindovv length is limited due to the Heisenberg uncertainty principle [43], and compromise is needed betvveen tirne and frequency resolution. Wavelet bispectrum based on Morlet mother vvavelet in contrast to the Fourier bispectrum enables us to gain the optimum tirne and frequency resolution at the same tirne vvhat is an advantage compared to the Fourier bispectrum. Since the tirne resolution of the vvavelet bispectrum is higher and the frequency resolution is poorer at high frequencies compared to Fourier bispectrum it is necessary to ensure sufficient frequency S Fdurier and vvavelet bispectrum cdmpariso resolution before interpretation of the results. Poor frequency resolution would result in poor/incorrect localization of characteristic frequencies. To high tirne resolution could result in to high sensitiveness to noise and statistical error, what would result in phase slips and incorrect oscillators coupling onset determination. Setting the tirne and frequency resolution so that episodes of approximately 10 periods of the lower coupling frequency are detectable and its characteristic frequency can be estimated to at least one tenth or less should be considered. Frequency step. Once the window length is chosen the frequency resolution is set and fixed for the Fourier based bispectrum. This is not the čase when using vvavelet transform. Since the vvavelet transform is continuous vve can choose the frequency step arbitrary. In this way, the transform can be over sampled in tirne for large scales but vve are not concerned for the inverse transform. Energy preservation. CV signals are power signals [81]. The Fourier bispectrum is based on the DFT vvhich gives the signals energy (power). In čase of vvavelet bispectrum the normalization is necessary to obtain the signals energy (povver). Statistical error. Integrating over finite tirne series in order to calculate the vvavelet bispectrum causes noise contribution to it's estimation. It is called statistical noise level since it is the value of vvavelet bispectrum that vvould be attained by a vvhite noise input signal, and is caused by finite statistics (i.e., using a limited number of values in the integrating or averaging process). Beside the noise contribution there is al so error estimate, vvhich is the product of uncertainties in the determination of the individual vvavelet bispectrum coefficients [43, 61, 62, 69]. To calculate the vvavelet bispectrum Eq. (7.17) the vvavelet coefficients are determinated for each JVw = Tfs samples in the interval T: {T0 - T/2 < t < T0 + 772} and averaged Eq. (7.3). Let us assume that ali the estimates of the vvavelet bispectrum are independent, than the averaged vvavelet bispectrum suffers a statistical error due to summation over 7VW values. Similarly in čase of Fourier bispectrum the summation is carried out over NI M ensembles, where M is the number of points in each statistically independent ensemble for vvhich M-points Fourier transform is calculated. The statistical error in the Fourier bispectrum decays as IN , and a factor of M more points are needed to obtain the same statistical error as vvith the vvavelet bispectrum. From this point of view, the vvavelet bispectrum represents a significant improvement in the tirne resolution of the bispectrum. Hovvever, the vvavelet coefficients are not ali statistically independent, since the chosen vvavelet family is not orthogonal. Each coefflcient is calculated by evaluating Eq. (7.3) integrating over the range -oo < t < +co. Due to the periodicity s of the vvavelets of scale s (Fig. 7.1), tvvo statistically independent BISPECTRUM CDMPARISON estimates of the wavelet coefficients are separated by a tirne s/2, or a number of points M(s) = sfjl, where fs is the sampling frequency. Thus the summation done in the evaluation of the wavelet bispectmm is not really carried out over Nw points, but only over 7VW /max(M(s)), where the maximum is taken over the values of s that come into account for the evaluation of a specific value of the wavelet bispectrum. An estimate for the statistical noise le vel in WB(f\,f?) is (8.6) The statistical error in the determination of WB(f],f2) can be deduced from Eq. (7.17). Each factor in this equation that is obtained by integrating over T suffers an error w , so that the error is estimated by (8.7) Eqs. (8.6) and (8.7) imply that vvavelet bispectral analysis is able to detect coherent signals in extremely noisy data, provided the coherency remains constant during sufficiently long times, since the noise contribution falls off rapidly with increasing N. Bispectrum interpretation. By choosing/0 = 1 a simple relation betvveen scale and frequency can be obtained /= l/s. In this čase the interpretation of the wavelet bispectrum is the same as for the FB othenvise it is not straightfonvard. Computation. Default wavelet bispectrum window length drops hyperbolically, whereas Fourier bispectrum used fixed window length. Wavelet bispectrum is therefore computational less demanding and much faster. Also relatively short data sequences are sufficient to perform an analysis, in contrast to the Fourier bispectrum that needs long tirne series to obtain both sufficient frequency resolution and statistics. 8.3 Other possible methods for bispectrum estimation A hvbrid betvveen Fourier and vvavelet transform is Selective Disctrete Fourier Transform (SDFT) which can also be used to perform the bispectrum calculation. It is a modified STFT first introduced by Keselbrener and Akselrod [47]. Like STFT it is tirne dependent FT. The time-frequency sensitivity 102 B Fdurier and vvavelet bispectrum cdmparisdn is obtained by windowing with a window of specific length around the analysed data point for estimating each spectral component. Low frequencies are expected to vary slowly on the other hand high frequencies are expected to show fast or sudden changes. For each frequency of interest, a DFT calculation is preformed, while the tirne window around the considered data point is selected inversely proportional to the frequency of interest. This is similar to wavelet transform's stretching and compressing of the mother wavelet. Therefore narrow windows are used for estimating high frequencies and wide ones for low frequencies, what implies estimating low frequencies with good frequency resolution and high frequencies with good tirne resolution. For each tirne of interest spectral components are calculated using different length of window. The window duration is chosen so that T = NJf. Parameter Np, Np e Z, is the number of entire periods entering the windowed signal. A high value of Np will lead to poor tirne resolution (wide window), while on the contrary, a small value could lead to a less reliable estimation of the spectral components in čase of noisy signals. The value is determined experimentally to gain best results, usually betvveen the range of 3 and 7 [38]. Leakage may appear in the spectrum, if the signal entering the rectangular window is not periodic, or at least if the amplitude of the end points is not equal. In order to remove such leakage, the data is usually convolved with some kind of smoothing window, such as Hamming, Hanning or Blackman windows. Their role is to taper the windowed data in order to make the two end points amplitudes smoothly equal. Besides the leakage removal, this tapering windows also improve the tirne resolution of the time-dependent spectral analysis. SDFT and WT provide similar results. Both transforms are using a specific window length to estimate each spectral component. SDFT ušes convolution with Blackman, Hanning, Hamming, or other taper window whereas WT ušes different mother vvavelets such as Morlet or other. Both methods enable a choice betvveen a good tirne and a good frequency resolution. We can change frequency and tirne resolution by changing parameters, but we cannot gain both of them simultaneously, according to the Heisenberg uncertainty principle. The WT obtained by Morlet wavelet enables an optimal time-frequency resolution, while using SDFT it can be approached by an appropriate choice of parameters. They can both be normalised to energy. The main difference betvveen the transforms is that the WT is continuous whereas the SDFT is not. IU J 8.4 Discussion Wavelet bispectrum was applied to CV blood flow signals. Parameters were set to the prevailing choice, Sec. 7.3. The WB method is suitable for studving cardio-respiratorv interaction from the CV blood flow signals. Results obtained with WB analysis are the same as the ones obtained with the Fourier based bispectrum method. There are no obvious advantages of WB over FB when detecting cardio-respiratory interaction that fulflls the conditions defined in Sec. 5.3. Our motivation was to develop a method that will be able to provide insight into the nature of the CV subsystem couplings. Dvnamics of CV blood flow can be considered in terms of coupled oscillators. There are at least five subsystems that take part in blood flow regulation: cardiac, respiratorv, mvogenic, neural and metabolic svstem [8, 9, 10, 110-112, 117]. Analysing interaction among cardiac and respiratory system is thus the first step taken. The effect of respiration on heart rate has been the most intensively studied. The question, what do the revealed nonlinear cardio-respiratory couplings mean, and how they arise, are yet to be resolved. One possibility is that they result from nonlinearity of the carotid baroreceptor-cardiac reflex [22]; another is that they are attributable to the active involvement of the peripheral vessels during cardiac and respiratory wave propagation in the netvvork, or they are due to modulation of the cardiac filling pressure during respiratory movements [126]. In order to be able to analyse ali other CVS interactions like cardio-myogenic, cardio- neural, cardio-metabolic, respiratory-myogenic, respiratory-neural, respiratorv-metabolic, myogenic-neural, we need to use the wavelet bispectral method. Namely an important feature of CV signals is that they are nonlinear, time-varying, and subject to fluctuation [3, 18, 23, 34, 117]. In low frequency range, which is of our particular interest, the characteristic frequencies are close to each other, and are therefore even harder to separate. The uncertainty principle of the FT, limits its ability to separate harmonic components in the frequency domain of the bispectrum [20, 69]. This might cause problems for detection of the quadratic phase couplings in the čase of frequency pairs that are close together. To ensure good resolution of low frequencies, we need longer sections for calculation of the discrete Fourier transform. This immediately decreases the number of sections possible and weakens the bispectrum estimation. However, we cannot use longer signals, because they lead to nonstationarity, and the variance consequently becomes even larger [69]. Moreover, determining short-lasting couplings, shorter than 10 times the lovver period in a bifrequency pair, makes the Fourier based bispectrum incapable of coping with the necessary time-frequency resolution, as was also clearly demonstrated using a model of coupled oscillators (8.1). 104 B Fourier and vvavelet bispectrum comparisoi Tab. 8.2: Summary of Fourier bispectrum (FB), using STFT and wavelet bispectrum (WB), tirne and frequency discretized using adapted Morlet mother wavelet comparison for CV signals analysis. * -denotation in the table means that no method is in advantage. Aside from WB being able to trace fast changes of high frequency components, and being able to locate slov/ frequency components at the same tirne, there are many other advantages of vvavelet bispectrum (adopted for CV signals) over the Fourier bispectrum. An overvievv is shown in table 8.2. Since WB is continuous, and FB is not, it allows an arbitrary frequency step to be chosen, and thus a better frequency component location. It allows intermittent phase couplings to be detected, vvhereas Fourier bispectrum averages out most of the tirne relevant information. The Heisenberg uncertainty principle, [43], limits simultaneous tirne and frequency resolution. Using the wavelet bispectrum, the optimum tirne and frequency resolution can be achieved; there is a simple relation between scale and frequency; and it has smaller statistical error; and is computationally less demanding. The only drawback of WB, compared to FB, is that it has to be normalized to obtain signal energy, and it is not orthogonal. Normalization can be preformed, vvhereas we are not concerned with the inverse vvavelet transform. A hybrid betvveen Fourier and vvavelet transform, the selective discrete Fourier transform, can also be used to perform the bispectrum calculation. It is a modified STFT, that v/as first introduced by 7 Bispectrum for signal xw(t) (see Sec. 4.4) was computed for the whole IT of the principal domain, Fig. 3.1. Fourier bispectrum \B\ was computed with K = 34 segments, 67 % overlapping and using the Blackman window to reduce leakage. Wavelet bispectrum | WB\ was computed with K = 34 segments, 67 % overlapping, Tm = 8 s, Ge = 0.001 and using fixed Morlet vvavelet length of THF = 40 s for high frequencies calculation. 100 % is the number of computations preformed for FB. 81 % less computations are necessary for WB estimation. If we use hyperbolically decreasing Morlet vvavelet length then computation of WB is approximately 40 times faster than computation of FB (only 2.6 % of FB computations are necessary). 105 :t bispectrum cdmparisdn Keselbrener and Akselrod [47], but the wavelet transform is more adequate, since it is continuous, whereas the SDFT is not. 106 9 DTHER POSSIBLE APPLICATIDNS OF VVAVELET BISPECTRUM METHDI 9 DTHER POSSIBLE APPLICATIONS OF VVAVELET BISPECTRUM METHDD C V signals are not the only ones relevant for studying CVS. Neural CV subsystem coupling information is incorporated in the brain waves. Synchronization among separate brain centres indicates an interaction between them. In most cases, it is associated with oscillatory behaviour in specific structures, frequencies and behaviour states. Electroencephalogram (EEG) measures electrical activity in the brain, i.e., brainwaves of different frequencies, and short-lived evoked potentials that occur when the brain responds to sensory input. Generally, low frequency oscillations originate from larger structures than do high-frequency oscillations. In certain conditions, such as general anaesthesia, svnchronization can be seen in EEG measurements as organized, distinguishable patterns. These patterns depend on the anaesthetic agent and the level of anaesthesia [107]. Delta waves are the slowest oscillating waves (0-4 cvcles per second). They are associated with a deep dreamless sleep, trance state, lucid dreaming, increased immune functions and hvpnosis, and are thus expected to occur during anaesthesia. In some cases, svnchronization can be observed over large distances in the consistent tirne lags between signals, using cross-correlation techniques. However, the cross-correlation techniques might encounter problems when compared to signals that are not stationary, or to oscillations that are weakly related. In the concept of phase svnchronization of chaotic oscillators [90], the solution is approached with the consideration of two tirne series, originating from two coupled oscillators. The amount of coupling can be quantified from the phase difference of the signals. Spatial heterogeneity in EEG during anaesthesia is often studied by means of amplitude and spectral estimate-based methods [21, 86, 122]. The quantification of quadratic phase-coupling betvveen EEG signal components has been established since G. DumermuuVs pioneered investigations using bispectral analvsis in 1969 [51]. A number of 101 E APPLICATIONS DF VVAVELET BISPECTRUM METHDD EEG studies have been published using the mathematical tools of high-order spectra analysis EEG [12, 27, 53, 65]. Additionallv, the so-called bispectral index (BIS) [89] is a frequently used parameter for the quantification of anaesthesia and sedation depth [84, 88, 98]. The BIS is a statistically based, empirically derived complex parameter that is composed of a combination of time domain, frequency domain and high-order spectral sub parameters. By means of bispectral analysis, non-linear interactions of EEG signal brainwaves can be quantified. The stationarity of the signal is one important prerequisite for consistent bispectrum estimation. Generally, the mathematical property of stationarity cannot be obtained from real EEG signals. Therefore, methods for spectral analysis of non-stationary signals have been introduced. Time-frequency distributions, vvavelet transform, and time-variant autoregressive moving average (ARMA) modelling are the most prominent classes of such approaches with sufficient time and frequency resolution. By means of time-frequency analysis, only transient linear relations of, and betvveen, signal components can be captured. Transient nonlinear interactions are undetectable. Therefore, approaches for time-variant bispectral analysis have been developed [5, 69, 98]. While these approaches concentrate on time-frequency domain (second-order spectra) or on shape in (third-order spectra) frequency-frequency domain, we extract time information related to coupling from the frequency-frequency domain of bispectrum, i.e., the biamplitude and the biphase. Recently, the svnchronization index technique was applied to signals of rats undergoing anaesthesia [63, 64]. EEG signals contain several time-varying frequency components. The most dominant ones are in the delta frequency range. Similar signal pattern was observed for ali rats analysed while undergoing anaesthesia. At the beginning, there is one dominant, slightly-varying, frequency component around the central frequency of 2 Hz. In its surrounding, there are higher frequency components that are not distinctly at the beginning. The predominant frequency component in the EEG signal vanishes when rats started to move and breath spontaneous. Synchronization indexes have been calculated for the čase of delta waves of EEG and ECG, delta waves of EEG and respiratory, and ECG and respiratory signals. Synchronization was distinctive only in the latter čase. There is one general pattern that occurs in ali cases: 2:1 or 3:1 svnchronization at the beginning, which eventually transits to 4:1 or 5:1, and then later returns back to 3:1 or 2:1. At the end of the signal there is no svnchronization what is in connection with rat transition from deeper to less deep anaesthesia (wakening). Furthermore, direction and strength of the coupling was studied. While there can be seen that in the first part the respiratory oscillator drives the oscillations contributed to delta waves of EEG signal, there cannot be made any conclusions about cardiac and delta waves of EEG signal direction and strength of the interaction. 108 9 DTHER PDSSIBLE APPLICATIONS O F VVAVELET BISPECTRUM METHD The depth of anaesthesia is related to synchronization states between cardiac and respiratory oscillator [63, 64]. Anaesthesia deepness can be extracted from the EEG signal therefore relation betvveen the cardio-respiratory synchronization and the bispectrum of EEG signal is expected. In this chapter we apply the wavelet bispectrum to EEG signal. For the demonstration of its applicability we use the EEG signal from the rat20 as already analysed in Chapter 6. By estimating the biamplitude and the biphase we wish to test whether the wavelet bispectrum can extract the same information from the univariate EEG signal as can be obtained from svnchrogram of bivariate ECG and respiration signals. 9.1 Measurements Measurements have already been discussed in the Sec. 6.1. In this analysis we use only the EEG signal measured on rat20 undergoing the anaesthesia. 9.2 Data analysis The EEG signal was first pre-processed. Both very low and very high frequencies were removed by use of moving average windows: drift with a 200 s long window; and high frequencies with a 0.04 s window while, and at the same tirne, the signal was resampled to 50 Hz. By using the moving average before resampling, we avoid problems of aliasing. The signal has been further normalized betvveen zero and one and its mean value was subtracted. For clearer interpretation of the results we divided -63 minutes long EEG signal of rat20 to four parts a, b, c and d, each of them containing only one phenomenon, i.e., svnchronization or no svnchronization. First we calculated the vvavelet bispectrum for each separate EEG part for the whole frequency domain. We used fixed, 20 s long, vvavelet for calculation of high frequencies as discussed in Sec. 7.2.1. Than we estimated the biphase and the biamplitude for the highest peak appearing in the vvavelet bispectrum. For this calculation 80 s long fix vvavelet vvas used to increase frequency resolution for high frequencies since the signal EEG is highly complex and the signal povver is concentrated at approximately 1 Hz. Frequency step vvas equidistant and set to 0.02 Hz to preserve the Eq. (7.18) condition. Morlet window vvas moved along the tirne series vvith a tirne step of 0.1 s. The critical value for the biamplitude estimate to be considered valid vvas set in ali cases to 2, i.e., tvvice the average value of the WB vvithin its so-called inner triangle (IT). »Ig—i_m 109 PLJCATIONS OF VVAVELET BISPECTRUM METHDD 9.3 Results Example of detrended, resampled and it's mean value subtracted EEG signal for rat20 undergoing anaesthesia is presented on Fig. 9.1 (a) and its power spectrum (b). As it can be seen from the power spectrum that it's power is concentrated in range 0.5-5 Hz with maximum between 1 Hz and 2 Hz. The synchrogram between the rat20 ECG and respiration signal is presented in Fig. 9.1 (c). Four distinctive parts can be seen; a and c where synchronization 4:1 takes plače and b and d where no synchronization is evident. The EEG signal was divided according to marked parts and analysed separately for the sake of clarity. Figure 9.2 (a) - (d) shows WB for each part separately with its contour view. In ali the WB-s there is a dominant peak. In the part a it is located around bifrequency (1.1 Hz, 1.1 Hz), in part b around bifrequency (1.3 Hz, 1.3 Hz), in part c around bifrequency (1.4 Hz, 1.4 Hz) and in the last part around bifrequency (1.6 Hz, 1.6 Hz). For each part biamplitude and biphase were calculated for the dominant, the highest peak appearing in the obtained WT. A longer fixed Morlet vvavelet was used, Jrf = 80 s, for estimating high frequencies as in the čase of \JVBa.d\ calculation, where fixed Morlet wavelet of 20 s was used for estimating high frequencies. The reason is that \WBa.d\ exhibit reach bispectral contents, therefore higher frequency resolution is necessary. In ali the cases the highest peak appears at the so-called self-coupling bifrequencies. We observe self-phase couplings of delta waves of EEG signal. Fig. 9.1: (a) 10 s of detrended, and removed zero mean signal EEG(t) and its power spectrum (b) for the čase of rat20 undergoing anaesthesia, -63 minutes long at sampling frequency fs = 50 Hz. (c) Cardio-respiratory synchrogram for rat20 divided into four parts a-d. Biamplitude Az\, Fig. 9.3 left (a), shows some coupling activity over the whole tirne vvith two higher peaks from approximately 7 minutes until approximately 10 minutes. During this tirne the biphase fa, Fig. 9.3 right (a), tends to be within n interval. We detect phase coupling Tpc. During part b the 9 DTHER PDSSIBLE APPLICATIONS OF VVAVELET BISPECTRUM METHD biamplitude Ahi, Fig. 9.3 left (b), shows three short lasting peaks while the biphase $,i, Fig. 9.3 right (b), tends to decrease ali the tirne of observation. In part 3 the biamplitude AcU Fig. 9.3 left (c), shows a very distinct and high peak lasting from 35 minutes to 36.17 minutes, when the synchronization is the strongest. In this tirne the biphase fa, Fig. 9.3 right (c), changes for 2.43 rad what is within the n interval and can be treated as phase coupling. The last part d shows no coupling, biamplitude Adu Fig. 9.3 left (d), except at the very beginning where the biphase $u, Fig. 9.3 right (d), is also constant whereas othervvise is not. Fig. 9.2: Results for the rat20 undergoing anaesthesia. (a) The wavelet bispectrum \WBa\ for part a calculated with 74 % overlapping and (b) its contour view, (c) | WBh\ for part b calculated with 65 % overlapping and (d) its contour view, (e) \WBC\ for part c calculated with 81 % overlapping and (f) its contour view and (g) \WBd\ for part d calculated with 49 % overlapping and (h) its contour view. In ali cases of WB K was set to 33 segments, whereas Ge = 0.001, Tm = 8 s and TH¥ = 20 s long fixed Morlet wavelet for estimating high frequencies was used. By applying the necessary conditions for the nonlinear quadratic coupling to be present to analysis of rat20's EEG signal, we consider the conditions only for the peak of observation - the self-coupling peak, Sec. 5.3. Ali detected phase couplings are summarized in Tab. 9.1. ATIONS CJF VVAVELET BISPECTRUM METHOD so Time (min) 62.67 50 Time (min) 62-67 Fig. 9.3: Left column the biphase , and its standard deviation c^ were calculated during Tpc. 112 9 DTHER PD5SIBLE APPLICATIONS OF VVAVELET 9.4 Discussion From Tab. 9.1, it can be seen that phase couplings onset mostlv in part a and part c. These couplings are the strongest and have the longest lasting tirne, Tpc. If we put ali biamplitudes from Aa\ to Adi together as one tirne evolution of biamplitude, and do this similarlv with biphases from fa to fat then two (or three) peaks stand out by their biamplitude. At the times where the peaks appear, the biphase tends to be constant (within the n interval). The tirne of onset and phase coupling duration is shown in Fig. 9.4. These two events can also be detected, as shown in the svnchrogram Fig. 9.1 (c), when the svnchronization 4:1 onsets and when it disappears. In the tirne betvveen these events, the svnchrogram does not show svnchronization to be present, whereas from the WB, there are visible, short-time phase-coupling events that cannot svnchronize the interacting oscillations as the biphase decavs uniformly. "S S o & H 2.33 7.45 23.83 35.00 Time (min) 62.67 Fig. 9.4: In Tab. 9.1, phase couplings, Tpc, that stand out by their biamplitude (1000 and more). Their onset times and duration are shown. The first and the last Tpc coincide with the onset and disappearance of the phase svnchronization between the cardiac and respiratory oscillators, Fig 6.8 (d). We can conclude that cardiac and respiratory CV systems and delta waves of EEG signal are coupled during the anaesthesia. The question arises vvhether the delta waves of EEG signal drive the cardiac and respiratory systems in svnchronization when rat20 is undergoing anaesthesia, or slow respiration svnchronizes with the cardiac system that further influences the delta waves of EEG signal. There is one applausive hvpothesis of anaesthesia reducing the RSA [83, 84]. Anaesthesia may stimulate inhibitory glvcine and GABAergic synapses in the NTS-NA axis, whose projections then inhibit higher brain centres, such as the limbic system. It also modulates the level of endogenous hvpothalamic peptides, responsible for the natural control of brain metabolism that are known to affect vagal control of cardiac rhvthm [82], and which would be affected if anaesthesia artificially reduces brain metabolism [1]. That is, anaesthesia may affect the delta waves of EEG signal in such a manner that they drive the cardiac and respiratory system to svnchronization. Similarly, as we analysed the cardio-respiratory interaction, we could also analyse cardio-delta waves of EEG signal and respiratory-delta waves of EEG signal interaction to study their interaction and the nature of couplings, whereas this is not the intent of this work. Wavelet bispectrum proved to be a ■mmm IIJ »PLICATIONS DF VVAVELET BISPECTRUM METHDD promising tool to analyse EEG signals during anaesthesia. It can detect the phase synchronization onset and its disappearance. The bispectrum is the first high-order spectrum after ordinary spectral analysis, and is useful for the investigation of non-linear interactions of the lowest-order (i.e., quadratic interactions). Powerful noise reduction is an integral part of the standard technique, as non-coherent contributions are averaged out, and weak coherent signals can be detected in very noisy data [68, 69]. The bispectral method has been extended to encompass tirne dependence, and has demonstrated the potential of the extended technique to determine the type of couplings among interacting nonlinear oscillators [40]. Time-phase couplings can be observed by calculating the bispectrum and adapted bispectrum, and obtaining the time-dependent biphase and biamplitude. The method has the advantage that it allows an arbitrary number of interacting oscillatory processes to be studied. Recently introduced methods for synchronization analysis among chaotic and noisy oscillations (see [79] and references therein) have stimulated applications to a variety of different systems. Methods for quantifying the strength and identifying the direction of couplings, based on nonlinear dynamic or information theory approaches, have recently been proposed [72, 93, 94, 103]. In this work, the question of the type of coupling that may result in synchronization was addressed, and a method was proposed for its analysis. It is applicable to both univariate data (a single signal from the coupled system) and multivariate data (a separate signal from each oscillator). Millingen et al. [61, 62] have analysed multivariate data using a combined wavelet and bispectral method, and have discussed its application in the field of chaos analysis. Here we have concentrated on univariate data and illustrated the potential of the time-phase bispectral method for the detection of higher-order couplings in the presence of noise. The possibility of using univariate data is of particular importance when dealing with real signals, as in practice, we often cannot observe and measure the separate subsystems directly, but only their combination, which is intrinsically difficult. Most of the methods proposed so far for svnchronization analysis and detection of the direction of couplings are based on bivariate or multivariate data [72, 79, 93, 94, 103]. In conjunction with frequency or time-frequency filtering [113, 121] or mode decomposition [36] to obtain two or more "separate" signals, these methods can be used for univariate data as well. Svnchronization can also be detected in univariate data through an analysis of angles and radii [42] in return tirne maps [119]. The time-phase bispectral method proposed in this work is not only applicable to the synchronization analysis of univariate data, but also, at the same tirne, allovvs one to determine the nature of the couplings among the interacting nonlinear oscillators. Its benefits include: (i) the possibility of iij observing the whole frequency domain simultaneously; (ii) detecting that two or more subsystems are interacting with each other; (iii) quantification of the strength of the interaction; and (iv) determination of whether the coupling is additive linear or quadratic, or parametric in one of the frequencies. We have shown the method to be suitable for the analysis of noisy signals. Although it was shown that the technique works effectively on a well-characterized simple model, there are some difficulties to be faced and overcome in applying it to real problems, e.g., to data from the cardiovascular system. Understanding the content of the bispectrum, and identification of the peaks of interest, are not always straightfonvard. To appreciate which peaks are the ones to focus on, one has to be aware of the basic properties of the system and its fundamental frequencies. Distinguishing a quadratic interaction from parametric frequency modulation may be easy when the coupling/modulation is relatively strong, but becomes more difficult in the čase of relatively weak coupling/modulation. In the latter čase, observing each phase in the triplet separately can be helpful. Also, it is not always an easy task to distinguish between quadratic interaction and parametric frequency modulation in cases when both of them occur simultaneously. Furthermore, where the possible basic frequencies are relatively close, it will be hard to detect them separately. This could cause particular problems in the detection of quadratic phase couplings where frequency pairs are close together. Although it is possible in principle to study an arbitrary number of interacting oscillators, it is advisable in practice to study them in pairs: knovvledge of the basic frequency of each is necessary. The blood flow signal contains a great deal of information and is exceptionally challenging in relation to processing. It possesses components whose amplitudes and frequencies vary in tirne. Moreover, the interactions among its characteristic oscillations also vary in tirne, and their nature (frequency, phase, linear and/or quadratic couplings) also changes, giving rise to the observed complexity of cardiovascular dvnamics. Bispectral analysis has provided insight into the nature of the couplings. Results in this work support the inference that the dvnamics of blood flovv can usefully be considered in terms of coupled oscillators. Application to the cardio-respiratory interaction has shown for the first tirne that nonlinear coupling is present [41]. Although evidence for couplings bevond second-order has not been sought, higher-order coupling may also exist. In this work, it was shown that the effect of the coupling between the cardiac and respiratory oscillations is episodic, rather than fixed and permanent. Moreover, an interchange between frequency and phase couplings is also present, as demonstrated by the evolution of their time-biphase. 116 1 D SUMMARY Nonlinear coupling was revealed, and shown to exist during spontaneous, as well as during paced, respiration. Episodes with nonlinear coupling were detected in 11 out of the 22 recordings, and lasted between 19 s, in the čase of high respiratory frequency (f2p = 0.34 Hz), and 106 s, in the čase of low paced frequency of respiration (/->p = 0.11 Hz). The episodic nature of the cardio-respiratory interaction in a healthy human during spontaneous and paced respiration had already been demonstrated using quite different techniques of analysis [10, 46, 95, 97, 101, 104]. It allows us to infer that the inter-oscillator coupling is probably relatively weak. There are, however, compelling arguments suggesting that the cardiac and respiratory subsystems should be, in fact, treated as weakly nonlinear oscillators that are weakly coupled. (i) In healthy subjects, breathing spontaneously, only occasional and brief episodes of synchronization are seen [10, 99-101], indicative of relatively weak coupling. (ii) Sinus arrhythmia is small at spontaneous breathing frequencies and only slightly larger at very low breathing frequencies [23], again supporting a weak-coupling description. (iii) The couplings can sometimes decrease almost to a vanishing point, e.g., in coma [112]. Without couplings, the dynamics become drastically simplified - with complete absence of synchronization or modulation. The fact that virtually no variability is seen in any of the natural frequencies, despite small amplitude variations attributable to internal noise, suggests that the oscillators themselves are, at most, weakly nonlinear. (iv) If there were strong oscillator nonlinearity, and strong (but linear) coupling, we would observe many combinatorial components around the cardiac frequency, which is not the čase. Using bispectral and cross-bispectral analysis, it was also shown that the coupling information among cardiac and respiratory processes is inherit from the processes and spatially invariant. Both processes are of central origin, and their phase relationships can be observed in ECG, blood flow and blood pressure signals derived from widely separated sites. It would appear that the information is incorporated within the wave motion of the blood propagating through the vessels. It is of interest to compare results of this experiment with that of [19]. They had also studied physiological tirne series while respiration was being paced at a constant rate; in addition, they also guided the ventilatory amplitude (tidal volume) so as to produce a sinusoidal modulation envelope with a period of 60 s. It resulted in oscillations of the same period in several physiological quantities, including the R-R intervals, blood pressure, and the cardiac stroke volume and output. In this study, the ventilatory amplitude was left to the subjecfs spontaneous choice. The slow amplitude modulation of [19] was selected to mimic the pattern of Cheyne-Stokes respiration, which is often associated with heart failure. The question addressed in this work was to examine whether or not the relationship between cardiac and respiratory oscillations can be nonlinear, without any amplitude modulation. It is difficult to decide whether the nonlinear interactions, shown to occur episodically in the present study, would or vvould not have occurred if the amplitude of respiratory oscillations had also been controlled, m | as in the study by [19]. It was shown, however, that quadratic coupling exists even when both the frequency and the amplitude of respiration are spontaneous. We may therefore conclude that the nonlinear nature of the interaction betvveen the cardiac and respiratory oscillations is inherent, and that it becomes more pronounced when the frequency of respiration is kept constant. The question of what the nonlinear couplings mean, and how they arise, are yet to be resolved. One possibility is that they result from nonlinearity of the carotid baroreceptor-cardiac reflex [22]; another is that they are attributable to the active involvement of the peripheral vessels during cardiac and respiratory wave propagation in the netvvork. A full understanding of these couplings is essential to gain insight into the physiology and pathophysiology of cardiovascular dynamics, as well as for the construction of mathematical models that offer novel possibilities for obtaining clinically relevant physiological information. It can be concluded that bispectral analysis provides a promising tool for the determination of frequency and phase couplings in the processing and understanding of cardiovascular signals. One of the coupling phenomena that can onset among interacting oscillators is phase svnchronization. The most natural way to study the phase relations among interacting oscillators is by using svnchronization techniques. As bispectral analyses provides information about frequency and phase coupling, and moreover, the nature of the coupling oscillators, we draw its relation to the svnchronization. In the čase of strong phase synchronization, the svnchrogram and bispectrum provide the same results. The phase svnchronization is indicated with a high biamplitude value and constant biphase at bifrequency of our primary interest (f\,fi). Weak phase svnchronization usually does not result in obvious horizontal strips in svnchrogram, vvhereas bispectrum results in moderate biamplitude and less constant biphase, with more phase slips at the bifrequency of our primary interest. No phase svnchronization results in zero biamplitude at the bifrequency of our primary interest. Bispectral analysis is more sensitive to interactions than the svnchrogram. It detects the phase svnchronization, and nevertheless, yields different information from that which can be resolved from a synchrogram. Synchronization can take plače simultaneously with other type of interaction, such as frequency modulation, forcing and/or nonlinear coupling. Nonlinear coupling can appear while phase svnchronization occurs, whereas svnchronization does not necessarily appear while there is nonlinear coupling. There is no obvious link between the two phenomena. Analysis of rats undergoing anaesthesia shows that nonlinear coupling occurs during svnchronization while analysis of CV blood flow signals of humans in resting shows nonlinear coupling while no svnchronization is present. 118 1 D SUMMARY Frequency modulation alone can be detected using bispectrum. Peak of primary interest (f\,fi) and the second peak at bifrequency if\-fi,fi) are high in comparison to other peaks (which may not be present at ali) that appear in the čase of nonlinear interaction, and their biphase is constant (strong frequency modulation). Instantaneous presence of frequency modulation and forcing can be misleading in detecting the nonlinear coupling. Thus, it is necessary to check ali the 1-6 peaks, and analyse the peak of primary interest (/i, f2) for modified bispectrum, and observe the phases of each frequency component in the triplet. When strong frequency modulation and nonlinear coupling are occur simultaneously, detection of the frequency modulation is not possible. In the čase of bispectral analysis of cardiovascular interaction, the time-dependent biphase/biamplitude estimate was estimated with a short tirne Fourier transform, using a window of constant length. The optimal window length depends, however, on the frequency being studied. The effective length of the window used for each frequency can be varied by applying the wavelet transform, or the selective discrete Fourier transform. For demonstration purposes above, the natural frequencies of the oscillators were chosen to lie within a relatively narrow frequency interval. An STFT was therefore sufficient for good tirne and phase/frequency localization. With broader frequency content, however, the wavelet transform or selective discrete Fourier transform needs to be applied. Wavelet transform to higher-order spectra was presented in this work. The wavelet and cross-wavelet bispectrum was defined analogous to the definitions used in Fourier based bispectrum and cross-bispectrum. By doing this time-dependant biphase/biamplitude estimate with higher frequency resolution at low frequencies, higher tirne resolution at higher frequencies was obtained. The vvavelet bispectral analysis was adopted to analyse cardiovascular signals. For a mother wavelet modulated Gauss function, the Morlet mother wavelet was used. Wavelet bispectral analysis was illustrated using a test signal. Since the tirne resolution of vvavelet bispectrum is higher, and the frequency resolution is poorer at high frequencies compared to Fourier based bispectrum, it is necessary to ensure sufficient frequency resolution before interpretation of the results. Poor frequency resolution vvould result in poor/incorrect localization of characteristic frequencies. Too high of a tirne resolution could result in extremely high sensitiveness to noise and statistical error, which vvould result in phase slips and incorrect oscillator coupling determination. It was necessary to raise frequency resolution for high frequencies, and also to preserve the scale (frequency) sum condition necessary for bispectrum estimation. Wavelet bispectrum results are parameter set dependant. Wavelet bispectrum was applied to CV blood flow signals. Parameters were set to the prevailing choice. The vvavelet bispectral method is suitable for studying cardio-respiratory interaction from the 113 CV blood flow signals. Parameter impact on wavelet bispectrum estimation and detailed comparison with Fourier bispectrum was preformed. Results obtained with WB analysis resemble the ones obtained with Fourier bispectrum method. There are no obvious advantages of WB over FB when detecting cardio-respiratory interaction that fulfils the defined necessary conditions for quadratical coupling occurrence. Our motivation was to develop a method that will be able to provide insight into the nature of the CV subsystems couplings. Analysing interaction among cardiac and respiratory systems is thus the first step taken. In order to be able to analyse ali other CVS interactions, such as cardio-myogenic, cardio-neural, cardio-metabolic, respiratory-myogenic, respiratory-neural, respiratory-metabolic, mvogenic-neural, and others, we need to use the wavelet bispectral method. Namely, an important feature of CV signals is that they are nonlinear, time-varying and subject to fluctuation [3, 18, 23, 34, 117]. In the low frequency range, which is of our particular interest, the characteristic frequencies are close to each other, and are therefore even harder to separate. The uncertainty principle of the Fourier transform limits its ability to separate harmonic components in the frequency domain of the bispectrum [20, 69]. This might cause problems for detection of the quadratic phase couplings in the čase of frequency pairs that are close together. To ensure good resolution of low frequencies, we need longer sections for calculation of the discrete Fourier transform. This immediately decreases the number of sections possible and weakens the bispectrum estimation. However, we cannot use longer signals, because they lead to nonstationarity, and the variance consequently becomes even larger [69]. Moreover, determining short-lasting couplings, shorter than 10 times the lower period in a bifrequency pair, makes the Fourier based bispectrum incapable of coping with the necessary time-frequency resolution, as was also clearly demonstrated using a model of nonlinearly coupled Poincare oscillators. Aside from WB being able to trace fast changes of high frequency components, and being able to locate slow frequency components at the same tirne, there are many other advantages of wavelet bispectrum (adopted for CV signals) over the Fourier bispectrum. Since WB is continuous, and the FB is not, it allovvs an arbitrary frequency step to be chosen, and thus a better frequency component location. It allovvs intermittent phase couplings to be detected, whereas Fourier bispectrum averages out most of the tirne relevant information. The Heisenberg uncertainty principle, [43], limits simultaneous tirne and frequency resolution. Using the vvavelet bispectrum, the optimum tirne and frequency resolution can be achieved: there is a simple relation betvveen scale and frequency; it has smaller statistical error; and is computationally less demanding. The only dravvback of WB compared to FB, is that it has to be normalized to obtain signal energy, and it is not orthogonal. Normalization can be preformed, whereas we are not concerned with the inverse vvavelet transform. 120 1 D SUMMARY A hybrid betvveen Fourier and wavelet transform, the Selective Discrete Fourier Transform (SDFT), can also be used to perform the bispectrum calculation. It is a modified STFT, that was first introduced by Keselbrener and Akselrod [47], but the wavelet transform is more adequate, since it is continuous, whereas the SDFT is not. Not only CV signals are relevant for studying CVS. Neural CV subsystem coupling information is incorporated in the brain waves. Svnchronization among separate brain centres indicates an interaction between them. In certain conditions, such as general anaesthesia, svnchronization can be seen in EEG signals [107] as the onset of delta braimvaves (0-4 cvcles per second). Generallv, EEG signals do not satisfy the stationarity condition, which is generally necessary for high-order spectra estimation. Therefore, methods for spectral analysis of non-stationary signals have been introduced. Time-frequency distributions, wavelet transforms and time-variant autoregressive moving average (ARMA) modelling are the most prominent classes of such approaches with a sufficient tirne and frequency resolution [12, 27, 51, 53, 65]. Spatial heterogeneity in EEG during anaesthesia is often studied by means of amplitude and spectral estimate-based methods [21, 86, 122]. These approaches concentrate on time-frequency space (second-order spectra) or on shape in (third-order spectra) frequency-frequency space. Using the introduced tirne dependant wavelet bispectrum, one can extract tirne information related to coupling from the frequency-frequency space of bispectrum, i.e., the biamplitude and the biphase. Recently, the svnchronization index technique was applied to signals of rats undergoing anaesthesia [63, 64]. EEG is a highly complex signal, which contains several time-varying frequency components. The signal power is concentrated at delta brainwave range, at approximately 1 Hz - 2 Hz. The predominant frequency component in the EEG signal vanishes when rats started to move and breath spontaneous. Svnchronization indexes have been calculated for the čase of delta waves of EEG and ECG, delta waves of EEG and respiratory, and ECG and respiratory signals. Svnchronization was distinctive only in the latter čase. The depth of anaesthesia is related to svnchronization states betvveen the cardiac and respiratory oscillators. We applied the vvavelet bispectrum to the EEG signal of rat20 undergoing anaesthesia, in order to demonstrate its applicability. The analysis showed that the cardio-respiratory synchronization can be detected in the EEG signal. It can be distinctively seen from the evolution of the biamplitude and the biphase, when the svnchronization 4:1 onsets and when it disappears. We can conclude that cardiac and respiratory C V systems and delta waves of EEG signal are coupled during the anaesthesia. The physiological question arises vvhether cardio-respiratory synchronization might be a consequence of, or result in, delta waves of EEG signal when rats are undergoing anaesthesia. One possible hvpothesis explains that anaesthesia reduces RS A [83, 84]. 121 : Anaesthesia may stimulate inhibitory glycine and GABAergic synapses in the NTS-NA axis, whose projections then inhibit higher brain centres such as the limbic system. It also modulates the level of endogenous hypothalamic peptides responsible for the natural control of brain metabolism that are known to affect vagal control of cardiac rhythm [82], and which would be affected if anaesthesia artificially reduces brain metabolism [1]. That is, anaesthesia may affect the delta waves of EEG signal in such a manner that they drive the cardiac and respiratory systems to synchronization. Similarly, as we analysed the cardio-respiratory interaction, we could also analyse cardio-delta waves of EEG signal and respiratory-delta waves of EEG signal interaction to study their interaction and the nature of couplings, whereas this is not the intent of this work. Wavelet bispectrum proved to be a promising tool to analyse EEG signals during anaesthesia. It can detect the delta waves of EEG signal phase coupling states. A long-term aim is therefore to develop a coupled oscillator model that can provide a description of the system, quantify the couplings and relate their values to its different states of health or disease. The vvavelet bispectrum may provide a link betvveen theoretical CVS models and experimental measurements. Higher order spectral methods can be used to study arbitrary interactions among coupled oscillators: of quadratic, cubic, or even higher order. In this work we have concentrated on the lovvest interaction, using the third-order spectrum or bispectrum. It has been suggested to proceed to the calculation of even higher-order spectra than the third-order bispectrum. For higher orders, the volume of the calculations rises substantially, and the method becomes increasingly demanding, numerically. However, the difficulty is not mathematical, since the generalization of the bispectrum to higher order is straightforward, but practical. The representation and interpretation of such high-order spectra become increasingly difficult vvhich is the limiting factor. We aspect new, powerful computer visualization tools to open up this direction of development in the coming future. 122 The time-phase bispectral method allows us to determine the nature of the couplings among the interacting nonlinear oscillators. Its benefits include: (i) the possibility of observing the whole frequency domain simultaneously; (ii) detecting that two or more subsystems are interacting with each other; (iii) quantification of the strength of the interaction; (iv) determination of whether the coupling is additive linear or quadratic, or parametric in one of the frequencies; and (v) the method is suitable for the analysis of noisy signals. The effect of the coupling between the cardiac and respiratory osciliations is episodic, rather than fixed and permanent. Frequency and phase couplings interchange. Nonlinear coupling exists during spontaneous as well as during paced respiration. The inter-oscillator coupling is relatively weak. Bispectral and cross-bispectral analysis showed that the coupling information among cardiac and respiratory processes is inherit from the processes, and is spatially invariant. Both processes are of central origin, and their phase relationships can be observed in ECG, blood flow and blood pressure signals derived from widely separated sites. The nonlinear nature of the interaction between the cardiac and respiratory osciliations is inherent, and it becomes more pronounced when the frequency of respiration is kept constant. Bispectral analysis is capable of determination of frequency and phase couplings in the processing and understanding of cardiovascular signals. Bispectral analysis is more sensitive to interactions and is more noise robust than the synchrogram. It detects the phase synchronization, and nevertheless, yields different information from that which can be resolved from a synchrogram. A simple relation between the synchrogram and the bispectrum revealed information cannot be drawn. Results of CV blood flow signals analysed using wavelet bispectral method gave the same results as in the čase of using the Fourier based bispectrum method. There are no obvious advantages of wavelet bispectral method over the Fourier bispectral method, when detecting cardio-respiratory interaction that fulfils the conditions defined for quadratical coupling onset. 123 It is recommended that the wavelet and cross-wavelet bispectrum are applied to cardiovascular signals, rather than to the Fourier based bispectrum. They allow intermittent phase couplings to be detected, optimum tirne and frequency resolution, simple relation between scale and frequency, direct interpretation, normalization to signal energy, smaller statistical error, arbitrary frequency step and are computationally less demanding. We can conclude that cardiac and respiratory CV systems and delta waves of EEG signal are coupled during the anaesthesia. Anaesthesia may stimulate inhibitory glycine and GABAergic synapses in the NTS-NA axis, whose projections then inhibit higher brain centres, such as the limbic system. It also modulates the level of endogenous hypothalamic peptides, responsible for the natural control of brain metabolism that are known to affect vagal control of cardiac rhvthm, and which would be affected if anaesthesia artificially reduces brain metabolism. Anaesthesia may affect the delta waves of EEG signal in such a manner that it drives the cardiac and respiratory system to svnchronization. A long-term aim is therefore to develop a coupled oscillator model that can provide a description of the system, quantify the couplings and relate their values to its different states of health or disease. The wavelet bispectrum may provide a link betvveen theoretical CVS models and experimental measurements. 124 References References [I] M.T. Alkire, C.D.J. Pomfrett, Toward a monitor of depth: Bispectral index (BIS) and respiratory sinus arrhythmia (RSA) both monitor cerebral metabolic reduction during isoflurane anaesthesia, Anaesthesiology 87, A421 (1997). [2] A. Andronov, A. Vitt and S. Khaykin, Theory of oscillations (Pergamon Press, Oxford, 1966). [3] R.M. Berne, M.N. Levy (Eds.), Physiology (Mosby, St. Louis, Missouri, 1998). [4] I. Blekhman, Synchronization of dynamical systems (Nauka, Moscow, 1971). [5] B. Boashash, P.J. 0'Shera, Polynomial Wigner-Ville distributions and their relationship to time-varying higher order spectra, IEEE Trans. Signal Processing 42, 216 (1994). [6] M. Bračič and A. Stefanovska, Wavelet-based analysis of human blood-flow d>namics, Buli. Math.Biol. 60,919(1998). [7] M. Bračič and A. Stefanovska, Wavelet analysis in studying the dynamics of blood circulation, Nonlinear Phenom. Complex Syst. 2, 68 (1999). [8] M. Bračič Lotrič, Couplings among Subsystems that regulate Blood Flow, Ph.D. Thesis, University of Ljubljana, 1999. [9] M. Bračič, P.V.E. McClintock and A. Stefanovska, Stohastic and Chaotic Dynamics in the Lakes (Melville, Nevv York: American Institute of Phvsics, 2000). [10] M. Bračič Lotrič and A. Stefanovska, Svnchronization and modulation in human cardiorespiratory system, Physica A 283, 451 (2000). [II] D.R. Brillinger and M. Rosenblatt, Spectral Analysis ofTime Series (Nevv York, Willey, 1967). [12] T.H. Bullock, J.Z. Achimovvicz, R.B. Duckrovv, S.S. Spencer, V.J. Iragui-Madoz, Bicoherence of intracranial EEG activitv, Electroenceph. Clin. Neurophysiology 103, 661 (1997). 125 [13] A.C. Burton and R.M. Taylor, A study of the adjustment of peripheral vascular tone to the requirements of the regulation of body temperature, Am. J. Phvsiol. 129, 565 (1940). [14] V. Chandran, Two-dimensional Bispectral Analysis and Leakage effects on the Statistics of the Bispectrum, Ph.D. Thesis, Washington State University, 1990. [15] V. Chandran and S.L. Elgar, Mean and variance of estimates of the bispectrum of a harmonic random process - an analysis including leakage effects, IEEE Transactions on Signal Processing 39,2640(1991). [16] A.V. Dandavvate and G.B. Giannakis, Asvmptotic theory of mixed tirne averages and kth-order cvclic-moment and cumulant statistics, IEEE Transaction on Information Theory 41, 216 (1995). [17] I. Daubechies, Ten Lectures on JVavelets, (Philadelphia, SIAM, 1992). [18] C.T.M. Davies and J.M.M. Neilson, Sinus arrhvthmia in man at rest, J. Appl. Phvsiol. 22, 947 (1967). [19] L.C. Davies, D.P. Francis, A. Crisafulli, A. Concu, A.J.S. Coats and M. Piepoli, Oscillations in stroke volume and cardiac output arising from oscillatory ventilation in humans, Exp. Physiol. 85, 857 (2000). [20] V. DeBrunner, M. Ozavdin, and T. Przebinda, Resolution in time-frequency, IEEE Transactions on Signal Processing 47, 783 (1999). [21] J.C. Drummond, CA. Brann, D.E. Perkins, D.E. Wolfe, A comparison of median frequency, spectral edge frequency, a frequency band power ratio, total power, and dominance shift in the determination of depth of anesthesia, Acta. Anaesthesiology Scand. 35, 693 (1991). [22] D.L. Eckberg, Nonlinearities of the human carotid baroreceptor-cardiac reflex, Circ. Res. 47, 208 (1980). [23] D.L. Eckberg, The human respiratory gate, J. Phvsiol. 548, 339 (2003). 126 References [24] P. Engel, G. Hildebrandt, and H.-G. Scholz, Die Messung der Phasenkopplung zwischen Herzschlag und Atmung beim Menschen mit einem neuen Koinzi-denzmeBgerat, Pflugers Arch. 298,258(1968). [25] J.W.A. Fackrell, Bispectral Analysis of Speech Signals, Ph.D. Thesis, University of Edinburgh, 1996. [26] J.R. Fonollosa, C.L. Nikias, Wigner higher order moment spectra: Defmition, properties, computation and application to transient signal analysis, IEEE Trans. Signal Processing 41, 245 (1993). [27] RJ. Gajraj, M. Doi, H. Mantzaridis, G.N.C. Kenny, Analysis of the EEG bispectrum, auditory evoked potentials and the EEG power spectrum during repeated transitions from consciousness to unconsciousness, Br. J. Anaesthesy 80,46 (1998). [28] A. Grossmann and J. Morlet, Decomposition of Hardy functions into square integrable wavelets of constant shape, SIAM J. Math. Anal. 15, 723 (1984). [29] S. Hales, J. Physiol, Statistical Essays II, Haemastatisticks (London, Innings Manby, 1773). [30] C. Hayashi, Nonlinear oscillations inphysical systems (McGraw-Hill, New York, 1964). [31] M.J. Hinch, Testing for Gaussianity and linearity of a stationary tirne series, J. Time Ser. Anal. 3, 169(1982). [32] M.J. Hinch, Detecting a transient signal by bispectral analysis, IEEE Transactions on Acoustics. Speech. and Signal Processing 38, 1277 (1990). [33] M.J. Hinch, On the principal domain of the discrete bispectrum of a stationary signal, IEEE Transactions on Signal Processing 43, 2130 (1995). [34] J.A. Hirsch and B. Bishop, Respiratory sinus arrhythmia in humans: How breathing pattern modulates heart rate, Am. J. Physiol. 241, H620 (1981). [35] U. Hoffman, A. Yanar, U.K. Franzeck, J.M. Edvards, A. Bollinger, The frequency histogram: A new method for the evaluation of laser Doppler flux motion, Microvascular Res. 40, 293 (1990). 121 [36] NE. Huang, Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N. Yen, C.C. Tung, and H.H. Liu, The empirical mode decomposition method and the Hilbert spectrum for non-stationary tirne series analvsis, Proč. R. Soc. Lond. A 454, 903 (1998). [37] C.H. Hugenii (Huvgens), Horologium Oscillatorium (Apud F. Muguet, Parisiis, France, 1673). [38] J. Jamšek, A. Stefanovska, Selektivna diskretna Fourierova analiza, Zbornik sedme Elektrotehniške in računalniške konference ERK '98 B, 339 (1998). [39] J. Jamšek, Bispectral Analysis ofCardiovascularSignals, M.Se. Thesis, Universitv of Ljubljana, 2000. [40] J. Jamšek, A. Stefanovska, P.V.E. McClintock and I.A. Khovanov, Time-phase bispectral analvsis, Phys. Rev. E 68, 016201 (2003). [41] J. Jamšek, A. Stefanovska and P.V.E. McClintock, nonlinear cardio-respiratorv interactions revealed by time-phase bispectral analvsis, Phvsics in Medicine and Biologv 49, 4407 (2004). [42] N.B. Janson, A.G. Balanov, V.S. Anishchenko and P.V.E. McClintock, Phase relationship betvveen two or more interaeting processes from one-dimensional tirne series. I. Basic theorv, Phys. Rev. E 65, 036211 (2002). [43] G. Kaiser, A Friendly Guide to JVavelets (Boston, Birkhauser, 1994). [44] D. Kaplan and L. Glass, Understanding Nonlinear Dynamics (New York, Springer, 1995). [45] J. Kastrup, J. Bhlow and N.A. Lassen, Vasomotion in human skin before and after local heating recorded with laser Doppler flowmetry. A method for introduetion of vasomotion, Int. J. Microcirc. 8, 205(1989). [46] T. Kenner, H. Passenhofer and G. Schwaberger, Method for the analysis of the entraiment betvveen heart rate and ventilation rate, Pflugers Archiv. 363, 263 (1976). [47] L. Keselbrener, S. Akselrod, Selective Discrete Fourier Transform Algorithm for Time-Frequency Analvsis, IEEE Transactions on Biomedical Engineering 43, 789 (1996). 128 References [48] Y.C. Kim, J.M. Beall, E.J. Powers, and R.W. Miksad, Bispectrum and nonlinear wave coupling, Phys. Fluids 32, 258 (1980). [49] Y.C. Kin and E.J. Powers, Digital bispectral analysis and its applications to nonlinear wave interactions, IEEE Transactions on Plasma Science PS-7, 120 (1979). [50] R.I. Kittney, O. Rompelhan, Analysis of the interaction of the human blood pressure and thermal system, In: J. Perkins (Ed.), Biomedical Computing, Pitman Medical, London, 49 (1977). [51] B. Kleiner, P. J. Huber, G. Dummermuth, Analysis of the interrelations betvveen frequency bands of the EEG by means of the bispectrum, Electroenceph clin. Neurophysiology 27, 693 (1969). [52] H. Koepchen, Rhythms in Physiological Systems (Springer-Verlag, Berlin, 1991). [53] M. Koskinen, T. Seppanen, J. Tuukkanen, A. Yli-Hankala, V. Jantti, Propofol anesthesia induces phase synchronization changes in EEG, Clinical Neurophysiology 112, 386 (2001). [54] Y. Kuramoto, Chemical Oscillations, JVaves, and Turbulence (Berlin, Springer, 1984). [55] J.R. Levick, An Introduction to Cardiovascular Physiology (London, Arnold, 2000). [56] C. Ludwig, Beitrage zur Kenntnis des Einflusses des Respirationbewegungen of den Blutumlaf im Aortensystem, Arch. Anat. Physiol. Wiss. Med. 13, 242 (1847). [57] V.K. Madiseti and D.B. Williams, The Digital Signal Processing Handbook (Florida, CRC Press, 1998). [58] R. Malek-Madani, Advanced Engineering Mathematics, (Reading, Addison Wesley Longman, 1998). [59] A. Malliani, M. Pagani, F. Lombardi and S. Cerutti, Cardiovascular neural regulation explored in the frequency domain, Circulation 84, 482 (1991). [60] J.M. Mendel, Tutorial on higher-order statistics (spectra) in signal processing and system theory: Theoretical results and some applications, Proceedings of the IEEE 79, 278 (1991). 129 | [61] B.Ph. van Milligen, C. Hidalgo, and E. Sanchez, Nonlinear phenomena and intermittency in plasma turbulence, Physical Rev. Lett. 74, 395 (1995). [62] B.Ph. van Milligen, E. Sanchez, T. Estrada et al., Wavelet bicoherence: a new turbulence analysis tool, Phys. Plasmas 2, 3017 (1995). [63] B. Musizza, Vzorčne povezave med biološkimi sistemi: Pristop k ugotavljanju globine anestezije, B.Sc. Thesis, University of Ljubljana, 2000. [64] B. Musizza, F. Bajrovič, P.V.E McClintock, M. Paluš, J. Petrovčič, and A. Stefanovska, Cardio-respiratory and neural interactions in anaesthesia, Nature, submited. [65] J. Muthuswamy, D.L. Sherman, N.V. Thakor, Higher-order spectral analysis of burst patterns in EEG, IEEE Trans. Biomed. Eng. 46, 92 (1999). [66] A.K. Nadi, Robust estimation of third-order cumulants in applications of higher-order statistics, IEE Proceedings-F 140, 380 (1993). [67] A.K. Nadi, Higher-order statistics in signal processing (Cambridge, Cambridge University Press, 1998). [68] C.L. Nikias and J.M. Mendel, Signal processing with higher-order spectra, IEEE Signal Processing Magazine 7, 10 (1993). [69] C.L. Nikias and A.P. Petropulu, Higher-order spectra analysis:A nonlinear signal processing framework (Englewood Cliffs, Prentice-Hall, 1993). [70] M. Paluš, Kolmorogov entropy from tirne series using information-theoretic functionals, Neural Netvvork World 7, 269 (1997). [71] M. Paluš and D. Hoyer, Surrogate data in detecting nonlinearity and phase svnchronization, IEEE Engineering in Medicine and Biology 17, 40 (1998). [72] M. Paluš, V. Komarek, Z. Hrnčif and K. Štebrova, Svnchronization as adjustment of information rates: detection from bivariate tirne series, Phys. Rev. E 63, 046211 (2001). 130 REFERENCES [73] M. Paluš and A. Stefanovska, Direction of coupling from phases of interacting oscillators: an information-theoretic approach, Phys. Rev. E 67, 055201 (R) (2003). [74] H. Parthasarathv, S. Prasad and S.D. Joshi, An ESPRIT-Like method for quadratic phase coupling estimation, IEEE Transactions on Signal Processing 43, 2346 (1995). [75] J. Penaz, Maver VVaves: Historv and methodologv, Automedica 2, 135 (1978). [76] R.J. Perry and M.G. Amin, On computing and implementing the running bispectra, IEEE Transaction on Signal Processing 43, 1017 (1995). [77] L.A. Pflug, G.E. loup, and J.W. loup, Sampling requirements and aliasing for higher-order correlations, J. Acoust. Soc. Am. 94, 2159 (1993). [78] L.A. Pflug, G.E. loup, and J.W. loup, Sampling requirements for nth-order correlations, J. Acoust. Soc. Am. 95, 2762 (1994). [79] A.S. Pikovskv, M.G. Rosenblum, and J. Kurths, Svnchronization; A universal concept in nonlinear sciences (Cambridge, Cambridge Universitv Press, 2001). [80] M.B. Priestlev and M.M. Gabr, Multivariate Analysis: Future Directions (North Holland, 1993). [81] J.G. Proakis and D.G. Manolakis, Digital Signal Processing (New Jersev, Prentice-Hall, 1996). [82] V.M. Pokrovskv, O.E. Osadchiv, Regulatorv peptides as modulators of vagal influence on cardiac rhvthm, Can. J. Phvsiol. Pharmacol. 73, 1235 (1995). [83] C.J.D. Pomfrett, J.R. Sneyd, M. Beech, T.E.J. Healy, Variation in respiratory sinus arrhvthmia may reflect levels of anaesthesia, British Journal of Anaesthesia 67, 6216 (1991). [84] C.J.D. Pomfrett, Heart rate variability, BIS and the depth of anaesthesia, British Journal of Anaesthesia 82, 659(1999). [85] M.R. Raghuveer, Time-domain approaches to quadratic phase coupling estimation, IEEE Transactions on Automatic Control 35, 48 (1990). 131 [86] I.J. Rampil, A primer for EEG signal processing in anesthesia, Anesthesiologv 89, 815 (1998). [87] T.S. Rao and K.C. Indukumar, Spectral and wavelet methods for the analvsis of nonlinear and nonstationar tirne series, J. Franklin Inst. 33, 425 (1996). [88] M. Renna, R. Venturi, Bispectral index and anaesthesia in the elderv, Minerva Anaesthesiologv 66, 398 (2000). [89] C. Rosow, P.J. Manberg, Bispectral Index Monitoring, Anesthesiol Clin North America 19, 947 (2001). [90] M.G. Rosenblum, A.S. Pikovskv, J. Kurths, Phase svnchronization of chaotic oscillators, Phys. Rev. Letters 76, 1804(1996). [91] M.G. Rosenblum, A.S. Pikovskv, J. Kurths, From phase to lag svnchronization in coupled chaotic oscillators, Phys. Rev. Letters 78, 4193 (1997). [92] M.G. Rosenblum, A.S. Pikovskv, C. Schafer, P. Tass, and J. Kurths, Handbook of Biological Physics (Elsevier, 2000). [93] M.G. Rosenblum, and A.S. Pikovsky, Detecting direction of coupling in interacting oscillators, Phys. Rev. E. 64, 045202 (2001). [94] M.G. Rosenblum, L. Cimponeriu, A. Bezerianos, A. Patzak, and R. Mrovvka, Identification of coupling direction: Application to cardirespiratory interaction, Phys. Rev. E 65, 041909 (2002). [95] F. Rsachke, Temporal Disorder in Human Oscillatory Systems (Springer-Verlag, Berlin, 1987). [96] N.F. Rulkov, M.M. Sushchik, L.S. Tsimring, and H.D.I. Abarbanel, Generalized svnchronization of chaos in directionally coupled chaotic svstems, Phys. Rev. E. 51, 980 (1995). [97] S. Rzeczinski, N.B. Janson, A.G. Balanov and P.V.E McClintock, Regions of cardiorespiratory svnchronization in humans under paced respiration, Phys. Rev. E 66, 051909 (2002). IJLi References [98] B. Schack, H. Witte, M. Helbig, Ch. Schelenz, M. Specht, Time-variant non linear phase coupling analysis of EEG burst patterns in sedated patiens during electroencephalic burst suppression period, Clinical neurophysiology 112, 1388 (2001). [99] C. Schafer, Analysis of synchronization in complex systems: Application to physiological data, Ph.D. Thesis, University of Potsdam, 1998. [100] C. Schafer, M.G. Rosenblum, J. Kurths, H.H. Abel, Heartbeat svnchronized with ventilation, Nature 293, 239 (1998). [101] C. Schafer, M. G. Rosenblum, H.H. Abel and J. Kurths, Svnchronization in the human cardiorespiratory svstem, Phys. Rev. E 60, 857 (1999). [102] T. Schreiber and A. Schmitz, Improved surrogate data for nonlinearity tests, Phys. Rev. Lett. 77, 635(1996). [103] T. Schreiber, Phys. Measuring information transfer, Rev. Lett. 85, 461 (2000). [104] H. Seidel and H. Herzel, Analysing Entrainment of Heartbeat and Respiration with Surrogates, IEEE Eng. Med. Biol. Mag. 17, 54 (1998). [105] I. Sharfer and H. Messer, The bispectrum of sampled data: Part I-detection of the sampling jitter, IEEE Transactions on Signal Processing 41, 296 (1993). [106] M. Small and K. CA. Tse, Detecting determinism in tirne series: The method of surrogate data, IEEE Trans, on Circuits and Sys. 50, 663 (2003). [107] T.B. Šolan, Anaesthetic effect on electrophysiologic recordings, J. Clin. Neurophysiology 15, 217(1998). [108] T. Soderstrom, A. Stefanovska, M. Veber and H. Svenson, Involvement of svmpathetic nerve activity in skin blood flow oscillations in humans, Am. J. Phvsiol. 284, H1638 (2003). [109] A. Stefanovska and P. Krošelj, Correlation integral and frequency anlysis of cardiovascular functions, Open Sys. & Infomation Dyn. 4, 457 (1997). 133 [110] A. Stefanovska and M. Bračič, Physics of the human cardiovascular systems, Contemporary Phys. 40,31(1999). [111] A. Stefanovska, M. Bračič and H.D. Kvernmo, Wavelet analvsis of oscillations in the peripheral blood circulation measured by laser Doppler technique, IEEE Trans. Biol. Med. Eng. 46, 1230 (1999). [112] A. Stefanovska and M. Bračič, Reconstructing cardiovascular dvnamics, Control Eng. Pract. 7, 161 (1999). [113] A. Stefanovska and M. Hožič, Spatial synchronization in the human cardiovascular system, Prog. of Theor. Phys. Suppl. 139, 270 (2000). [114] A. Stefanovska, H. Haken, P.V.E. McClintock, M. Hožič, F. Bajrovič and S. Ribarič, Reversible transitions betvveen svnchronization states of the cardiorespiratory system, Phys. Rev. Lett. 85, 4831 (2000). [115] A. Stefanovska, M. Bračič Lotrič, S. Strle and H. Haken, The cardiovascular system as coupled scillators?, Phvsiol. Meas. 22, 535 (2001). [116] A. Stefanovska, D.G. Luchinsky and P.V.E. McClintock, Modelling couplings among the oscillators of the cardiovascular system, Physiol. Meas. 22, 551 (2001). [117] A. Stefanovska, Cardiorespiratory interactions, Nonlinear Phenom. Complex Syst. 5, 462 (2002). [118] K. Stutte and G. Hildebrandt, Untersuchungen iiber die Koordination von Herzschlag und Atmung, Pflugers Arch. 289, R47 (1966). [119] K. Suder, F.R. Drepper, M. Schiek and H.H. Abel, One-dimensional, nonlinear determinism characterizes heart rate pattern during paced respiration, Am. J. Physiol. 275, HI092 (1998). [120] A. Swami, G.B. Giannakis and G. Zhou, Bibliography on high-order statistics, Signal Processing 60,65(1997). I3H [121] P. Tass, M.G. Rosenblum, J. Weule, J. Kurths, A. Pikovsky, J. Volkmann, A. Schnitzler, H.-J. Freund, Detection of n:m phase locking from noisy data: Application to magnetoencephalography, Phys. Rev. Lett. 81, 3291 (1998). [122] J. H. Thinker, F.W, Sharbrough, J.D. Michenfelder, Anterior shift of the dominant EEG rhythm during anesthesia in the Java monkey, Anesthesiology 46, 252 (1977). [123] J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, J.D. Farmer, Testing for nonlinearity in tirne series: The method of surrogate data. Physica D 58, 77 (1992). [124] J. Theiler, D. Prichard, Using "Surrogate Surrogate Data' to calibrate the actual rate of false positives in tests for nonlinearity in tirne series, Fields Inst. Comm. 11, 99 (1997). [125] B. van der Pol and J. van der Mark, The heartbeat considered as a relaxation oscillation, and an electrical model of the heart, Phil. Mag. 7, 763 (1928). [126] M.B. Visscher, A. Rupp, and F.H. Scott, Respiratory wave in arterial blood pressure, Am. J. Phvsiol. 70, 586 (1924). [127] G. Zhou, and G.B. Giannakis, Retrieval of self-coupled harmonics, IEEE Transactions on Signal Processing 43, 1173(1995). [128] V.S. Zykov, G. Bordiougov, H. Brandtstadter, I. Gerdes, and H. Engel, Periodic forcing and feedback control of nonlinear lumped oscillators and meandering spiral waves, Phys. Rev. E 68,016214(2003). ■■■■ 135 131 —E— ECG...............................See electrocardiogram EEG.........................See electroencephalogram electrocardiogram..........................................30 electroencephalogram..............................3, 107 —F— Fourier transform...........................................79 frequency locking.......................See frequency synchronization frequency modulation........................52, 71, 75 frequency resolution...................See resolution frequency step..............................................101 frequency sum-rule........................................82 —G— Gaussian noise....................................See noise —H— harmonics......................................................47 heart cycle........................................................7 Heisenberg uncertainty principle.............33, 99 high-order statistics........................................13 HOS............................See high-order statistics Huygens...........................................................7 instantaneous phase.......................................51 —K— Kuramoto.......................................................59 —L— laser Doppler technique.................................30 limit cycle....................................................5, 8 Poincare........................................................8 —M— magnitude................................................14, 82 mechanism extrinsic......................................................40 intrinsic......................................................40 modulation frequency....................................................27 moments........................................................13 Morlet.................................See mother wavelet mother wavelet..............................................78 Morlet........................................................80 mutual prediction approach..........................See synchronization mutual synchronization.....See synchronization —N— noise........................................................10, 60 Gaussian..............................................21, 97 nonidentical systems.....................................60 normalization................................................16 null hypothesis..............................................38 —O— oscillator forced.........................................................53 Poincare.....................................................17 relaxation...................................................47 van derPol.................................................47 —P— paced breathing.............................................30 Parseval identity............................................79 periodicity.....................................................77 phase.......................................................14, 82 definition...................................................46 randomization method...............................38 slips............................................................60 phase locking..........See phase synchronization principal domain...........................................15 inner triangle..............................................15 outer triangle..............................................15 —R— relative phase.................................................51 resampling.....................................................32 resolution frequency...........................................81, 100 tirne....................................................84, 100 138 lNDEX respiration paced....................................................30, 53 spontaneous................................................30 respiratory sinus arrhythmia......................2, 52 RS A...............See respiratory sinus arrhythmia —S— scale...............................................................78 scalogram......................................................79 SDFTSee Selective Discrete Fourier Transform selective discrete Fourier transform............102 Short Time Fourier Transform......................14 STFT..........See Short-Time Fourier Transform stroboscopic technique..................................51 subsvstems.......................................................8 cardiac..........................................................6 endothelial....................................................6 mvogenic......................................................6 neural...........................................................6 respiratorv....................................................6 surrogate data................................................38 svnchrogram..................................................51 svnchronization...........................10, 51, 57, 74 complete.....................................................60 frequency...................................................58 generalized.................................................60 global.........................................................59 high-order..................................................59 lag..............................................................59 mutual........................................................58 mutual prediction approach.......................71 phase..........................................................58 phase svnchronization...............................10 region.........................................................58 strong.........................................................74 weak..........................................................74 —T— torus..............................................................10 —U— uncertaintv principle.................See Heisenberg uncertaintv principle unidirectional................................................53 uniformlv....................................................113 univariate data......................................See data —V— van der Pol...................................See oscillator —W— wavelet bispectrum.......................................82 adopted Morlet vvavelet.............................84 biamplitude................................................83 biphase.......................................................83 cross-bispectrum.......................................83 definition...................................................82 normalization.............................................85 parameters.................................................86 scale sum-rule............................................82 vvavelet transform...................................77-80 definition...................................................78 discretization.............................................80 Morlet vvavelet..........................................80 mother vvavelet..........................................78 scale...........................................................81 weakly nonlinear.....................................47, 50 CDNTRIBUTIDNS TD SCIENCE CDNTRIBUTIDNS TO SCIENCE Original contributions to science: 1. The study of interacting nonlinear oscillators, using time-phase bispectral estimators of biamplitude and biphase. We have introduced, for the first tirne, time-phase bispectral estimators of biamplitude (3.6) and biphase (3.5) for unveiling phase coupling information from univariate data (Chapter 3, pages 13 to 16). We have shown that the introduced method is suitable for studving interacting nonlinear oscillators, and that it is capable of quantifying the strength of the interaction by bispectral estimate - biamplitude, whose value is proportional to the coupling coefficient value £•(2.6) of coupled nonlinear oscillators (Chapter 2, pages 7 to 12), and revealing the nature of the coupling, i.e., whether the coupling is additive linear or quadratic, or parametric in one of the frequencies (Chapter 4, pages 17 to 28). 2. Cardio-respiratory coupling in human cardiovascular svstem hvpothesis conflrmation. We have applied, for the first tirne, the time-phase bispectral method to cardiovascular blood flow signals, in order to studv the nature of cardio-respiratorv interactions. Despite the limitations of the method used to resolve between linear and nonlinear couplings in extreme conditions, the method is applicable before the couplings become too complex, considering the phvsiological knowledge of the svstem (Chapter 5, pages 29 to 56). Couplings betvveen cardiac and respiratorv oscillations are episodic, rather than fixed and permanent. Frequency and phase couplings interchange. Nonlinear coupling exists during spontaneous, as well as during paced, respiration. It becomes more pronounced vvhen the frequency of respiration is kept constant. The inter-oscillator coupling is relatively weak (Chapter 5 and 6, pages 29 to 37 and 57 to 76). The coupling information among cardiac and respiratory processes is inherent from the processes and is spatially invariant. Both processes are of central origin, and their phase relationships can be observed in ECG, blood flow, and blood pressure signals derived from widely separated sites (Chapter 5, pages 37 to 46). 3. Bispectum estimates generalization to wavelets. We have generalized bispectum estimates - biphase (7.20) and biamplitude (7.21) - to wavelets. The method is suitable for intermittent phase coupling detection, while providing optimum tirne TO THE SCIENCE and frequency resolution. It may provide a link bervveen theoretical CVS models and experimental measurements (Chapter 7 and 8, pages 77 to 106). 4. Coupling between cardiac and respiratorv svstems, and delta waves of EEG signals determination using generalized bispectral estimates. Delta waves of EEG signals of rats undergoing anaesthesia, reveal a phvsiological relation between cardiac and respiratorv svstems and delta waves of EEG signals. The cardio-respiratorv svnchronization might be a consequence of delta waves of EEG signals of rats undergoing anaesthesia (Chapter 9, pages 107 to 114). Portions of this thesis were published in the following papers: 1. J. Jamšek, A. Stefanovska, P.V.E. McClintock and I.A. Khovanov, Time-phase bispectral analvsis, Phys. Rev. E 68, 016201 (2003). 2. J. Jamšek, A. Stefanovska and P.V.E. McClintock, Nonlinear cardio-respiratory interactions revealed by time-phase bispectral analvsis, Phvsics in Medicine and Biology 49, 4407 (2004). Portions of this thesis were presented in the follovving scientific meetings: 1. J. Jamšek and A. Stefanovska, Bispectral analysis of cardiovascular signals, Nonlinear Seminar, Department of nonlinear physics, Lancaster University, United Kingdom (7.2. 2002). 2. J. Jamšek, A. Stefanovska and P.V.E. McClintock, Time-phase bispectral analysis, basic theory and applications, Nonlinear Seminar, Department of nonlinear physics, Lancaster University, United Kingdom (12.5. 2002). 3. J. Jamšek, A. Stefanovska and P.V.E. McClintock, Cardiovascular System, Cardiovascular system, time-phase bispectral analysis, basic theory and application, 2nd Slovenia-Japan Seminar, Center for Applied Mathematics and Theoretical Phvsics University of Maribor, Slovenia (28.5.-5.6 2003). 4. J. Jamšek and A. Stefanovska, Quadratic cardio-respiratory coupling?, INTAS international Workshop, Department of Phvsics, University of Pisa, Italy (22-24.4. 2003). 5. J. Jamšek, A. Stefanovska and P.V.E. McClintock, Nonlinear cardio-respiratory interaction, INTAS-ESF international Workshop, Ljubljana, Slovenia (10-13.11. 2003). B |42 Appendix APPENDIX A. Variance ofthe bispectrum estimate B. Generation ofharmonics A. Variance ofthe bispectrum estimate In order to interpret bispectral values from a finite length tirne series, the statistics of bispectrum estimates must be known [49]. To achieve statistical stabilitv, the tirne series is divided into K segments for averaging. Phases of different segments are independent of each other and random variables over [0, 2n). When there is a large number of segments, the estimate gains statistical stabilitv at the expense of power spectral and bispectral resolution. For a real signal, with a finite number of points, the compromise between bispectral resolution and statistical stabilitv may be expected at K approximately 30. Estimates are subject to statistical error, such as bias and variance. An estimate must be consistent, that is the statistical error must approach zero in the mean-square sense as the number of realizations becomes infinite. Here we neglect the effects of finite tirne series length, we assume that they are sufficiently long. Let us consider the bias and the variance of the bispectrum estimate B(k,l). The expected value of B(kJ) will be (Al) as K becomes infinite, X\ is the DFT ofthe z-th segment. Thus, B(kJ) can be taken as an unbiased estimate [14, 15, 25]. Its variance will be m Note that the variance is inverselv proportional to K. From a mathematical statistics point of view, it is a nontrivial task to compute the quantity in the bracket in terms of low-order spectra, but one may write a good approximation [14, 15, 25] in which čase the variance will be Note that it is a consistent estimate in the sense that the variance approaches zero as K becomes infinite. The variance is proportional to the product of the powers (P(k) = E[X(£)Z*(A:)]) at the frequencies k, l and k + I. Consequently, a larger statistical variability is introduced in estimating larger values in the bispectrum. Finally, the variance is proportional to [1 - b2(k, /)], where the bicoherence b2 is a normalized bispectrum, b2(k, l) = E[B(k,l)]21 [P(k)P(l)P(k + /)]. That is, when the oscillations at k, l and k + l are nonlinearly coupled (b2 - 1), the variance approaches zero, and when the components are statistically independent (b2 ~ 0), the variance is proportional to the power at each spectral component [14, 15, 25]. Brillinger and Rosenblatt [11] have investigated the asvmptotic mean and variance of Fourier-type estimates of high-order spectra and proved that under certain assumptions the £-th-order spectral estimate is asymptotically unbiased and Gaussian distributed and that estimates of different-order are asymptotically independent. The variances of the real and imaginary parts of the bispectrum are asymptotically (i.e., for large K) Gaussian and are equal, var{Re[#(&,/)]} = var{lm[ B(k,l)]}. For a perfect phase-coupled triplet, the variances of the real and imaginary parts are equal to zero. In the čase of no coupling, there is an identical contribution to the variances from the real and imaginary parts of the estimate of the bispectrum. The total variance is a sum of individual (i = 1,..., K) contributions, because different triplets are mutually statistically uncorrelated in the absence of phase coupling. Partial coupling can be expected to result in a combination of perfectly phase-coupled oscillations, and oscillations with randomly changing phases. m Appendix B. Generation of harmonics In this appendix we show which harmonics appear in the spectrum of a weakly driven weakly nonlinear oscillator, and in particular we establish which harmonics correspond to quadratic coupling. The analysis that follows is for a Poincare oscillator, but a similar result also follows for, e.g., a van der Pol or other oscillators of similar type. Consider an oscillator of the form x, =-xlrl -oxyx + Q(x2,xl), 9\ =~y\r\ +®i*p (B1) n =aQxf +y2y -o), where the term Q(x2, x\) corresponds to a coupling of the main oscillator (x\, y\) to another one (x2, j^). If the coupling is non-zero, the system cannot be solved exactly analytically. An approximate solution for small coupling and weak nonlinearity can, however, stili be obtained. The amplitude and phase in this čase vary only slightly and they can be expanded about (A0, fo) as 145 For small /?, # a the Eq. (B3) corresponding to those for fi and / can be solved approximately. For the simplest linear coupling of form (B6) one obtains Then, in the spectrum of the variable (B8) one observes the following harmonics: a>\, oh, 2co\ ± o>i. In the case of quadratic coupling (B9) there appear additional harmonics: 2a>i, 2cq\, 2g>\ ± 2a>i, (0\ ± cfy, 3co\ ± a^. In the limit under consideration, with small nonlinearity and weak coupling, the appearance of these additional combinational harmonics can confidently be associated with the presence of a nonlinear coupling. It is of course the case that, for a nonlinear oscillator, all sorts of combinational harmonics can in principle appear even for linear coupling. However, the generation of these harmonics is a second-order effect which becomes significant only for large nonlinearity and large coupling coefficients. Under the latter circumstances, just the appearance of particular combinational harmonics cannot necessarily be related to a given type of coupling and some further analysis is then required. 146 I hereby declare that the research described in this thesis, and the thesis itself, is the original, and šole work, performed by the author, under the guidance of doc. dr. Aneta Stefanovska, the mentor. The results which were done in collaboration with other colleagues are published in presented articles. Other assistance from colleagues is stated in the acknovvledgments. The published results of other authors are presented in the references.