UFDC Home  UF Theses & Dissertations   Help 
Material Information
Thesis/Dissertation Information
Subjects
Notes
Record Information

Full Text 
PAGE 1 1 OPTIMIZATION AND DATA MINING IN HEALTHCARE: PATIENTS CLASSIFICATION AND EPILEPTIC BRAIN STATE TRANSITION STUDY USING DYNAMIC MEASURES, PATTERN RECOGNITION AND NETWORK MODELING By JICONG ZHANG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2011 PAGE 2 2 2011 J icong Z hang PAGE 3 3 To my p arents in China PAGE 4 4 ACKNOWLEDGMENTS I would like to thank Dr. Panagote M. Pardalos Dr. Basim M. Uthman, Dr. Deng Shan Shiau, Dr. J. Chris Sackellares and Dr. J. Cole Smith for their guidance and support. During the past f our and a half years their energy and enthusiasm have motivated me to become a better researcher. I would also like to thank Dr. William Hager, Dr. Joseph P. Geunes and Dr. Jean Philippe Richard for serving on my supervisory committee and for their insightful comments and suggestions. I have been very lucky to have the chance to work with many excellent people in the department of Industrial and Systems Engineering ( I.S.E ) and Center for Applied Optimization (C.A.O) I would like to thank Dr. Chang Chia Liu, Dr. Pando G. Georgiev, Dr. Yongpei Guan Petros Xanthopoulos, Dr. Vladimir Boginski, Dr. Stas Busygin, Dr. Stanislav Uryasev, Dr. Vitaliy Yatsenko, Dr. Ravindra K. Ahuja, Dr. Joseph C. Hartman, Dr. Petar Momcilovic, Dr. Ashwin Arulselvan, Dr. Onur Seref, Dr. Nikita Boyko, Dr. Altannar Chinchuluun, Dr. Michael Bewernitz f rom B.D.L. Dr. Oleg Shylo, Dr. Alla Kammerdiner, and others in C.A.O for their valuable friendship and kindly assistance. I thank the staff in the Department of Neurology at the Gainvesville VA Medical Center including Scott Bearden and David Juras for th eir co llaboration in collecting research data I thank the principal research scientist Dr. Wangcai Liao in Cyberonics Inc. for his sharing of his experience in research and product development of implantable medical devices I thank Joy Mitchell from North Florida Foundation for Research and Education, Inc. for providing travel support in my research and conference. I acknowledge the financial support provided by College of Engineering and I.S.E department throughout my Ph.D. study. I acknowledge the collaboration of Optima Neuroscience Inc. by providing research data. I acknowledge the financial support PAGE 5 5 provided by Cyberonics Inc. during my research work there. I also acknowledge the collaboration and financial support provide d by VA hospital in Tampa. There are many colleagues I have worked with and many friends I hav e made at University of Florida, VA hospital or the collaborating companies, who gave me their help and support but were not mentioned here, I would like to take this opportunity to thank all of them for being with me throughout this wonderful journey. Lastly, I wish to thank my lovely family in my home country China my father Gang Zhang, my mother Dongmei Li, and my grandfather Zhitang Zhang for their selfless sac rifice and constant love. PAGE 6 6 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 9 LIST OF FIGURES ................................ ................................ ................................ ........ 10 LIST OF ABBREVIATIONS ................................ ................................ ........................... 12 ABSTRACT ................................ ................................ ................................ ................... 14 CHAPTER 1 INTRODUCTION OF EPILEPSY ................................ ................................ ............ 16 Definition of Neurological Disorder and Epilepsy ................................ .................... 16 Categorization of Epilepsy ................................ ................................ ...................... 16 2 CAUSES, DIAGNOSIS AND TREATMENT OF EPILEPSY ................................ .... 19 Neurobiological Causes of Epilepsy ................................ ................................ ....... 19 Genetic Disorders ................................ ................................ ............................. 19 Developmental Disorders ................................ ................................ ................. 20 Acquired Epilepsies ................................ ................................ .......................... 21 Diagnosi s and Treatment of Epilepsy ................................ ................................ ..... 22 3 PATIENTS CLASSIFICATION BETWEEN GENERALIZED NCSE AND TME USING EEG QUANTIFICATION FOR DIAGNOSIS OF STATUS EPILEPS TICUS ................................ ................................ ................................ ...... 25 Significance ................................ ................................ ................................ ............ 25 Approach ................................ ................................ ................................ ................ 26 Results ................................ ................................ ................................ .................... 31 Discussion and Conclusion ................................ ................................ ..................... 37 4 MINIMUM PREDICTION ERROR MODELS ................................ ........................... 40 LSE Model ................................ ................................ ................................ .............. 41 Basic LSE Model ................................ ................................ .............................. 41 Another A pproach of U nderstanding the LSE Model ................................ ........ 45 MMSE Model Prediction ................................ ................................ ......................... 49 Other Prediction Methods ................................ ................................ ....................... 51 5 CAUSAL RELATIONS BETWEEN TWO OR MULTIPLE TIME SERIES ................ 53 Two Variable (Time Series) Models ................................ ................................ ........ 53 PAGE 7 7 Three Variable (Time Series) Models ................................ ................................ ..... 59 Multiple Variable Models ................................ ................................ ......................... 61 Predict from the P ast of I tself ................................ ................................ ........... 63 Predict from the P ast of I tself and the P ast of O thers ................................ ....... 64 Predict from the P ast of I tself and the P resent S tate and the P ast of O thers ... 65 Predict from the P ast of I tself and the T emporally C omplete I nformation of O thers ................................ ................................ ................................ ........... 67 6 EPILEPTIC BRAIN NETWORK MODELING AND BRAIN STATE TRANSITION STUDY ................................ ................................ ................................ .................... 70 Network M odeling of B rains ................................ ................................ .................... 70 Neocortical Epilepsy and Its Clinical Challenges ................................ .................... 71 Existing Solutions in Clinical Applications ................................ ............................... 72 Epileptic Brain Network Modeling for Neocortical Epilepsy ................................ ..... 73 EEG Analysis on Networks in Epileptic Brains ................................ ........................ 76 7 BRAIN STATE TRANSITION STUDY FOR SEIZURE DETECTION AND PREDICTION, USING MULTIPLE FEATURE EXTRACTION, PRINCIPAL COMPONENT ANALYSIS AND CLASSIFICATION ................................ ............... 79 Purpose and Scope ................................ ................................ ................................ 79 Procedure Design and Performance Evaluation Criteria ................................ ......... 79 Brain State Classification of Diff erent Epileptic Stages ................................ ........... 82 Conjecture ................................ ................................ ................................ ........ 82 Method s ................................ ................................ ................................ ............ 82 Results and Conclusion ................................ ................................ .................... 83 Events D etection U sing S ub B ands D ivision, M ultiple F eature E xtraction, and Classification in High Dimensional Space on A mbulatory EEG ........................... 90 Event Type : Eye Close ................................ ................................ ..................... 90 Event Type: Eye Flutter ................................ ................................ .................... 91 Event Type: Eat & Chewing ................................ ................................ .............. 92 Patient Specific D etection of Intractable Seizures U sing S ub B ands D ivision, F eature E xtraction, and Classification in High Dimensional Space on S calp EEG R ec orded from P ediatric S ubjects ................................ ............................... 93 Database and Background ................................ ................................ ............... 93 Results of C ase 1 as A n E xample ................................ ................................ .... 94 8 R WAVE/BEAT DETECTION FROM SURFACE ECG WITH MUSCLE ARTIFACTS FOR CARDIAC BASED SEIZURE DETECTION ............................... 97 Background and Purposes ................................ ................................ ...................... 97 Design Features of the Proposed Beat Detection Algorithm ................................ ... 98 Performance and Standard Testing Results of the Proposed Algorithm ............... 100 Test: Beat Detection Accuracy with Human ECG ................................ ........... 105 Test : Beat Detection Accuracy with Noise Stressed Human ECG ................ 106 PAGE 8 8 Test: Beat Detection on Range and Accuracy of Heart Rate with Simulated ECG ................................ ................................ ................................ ............ 108 Test: Beat Detection on Range of R Wave Amplitude and QRS Duration with Simulated ECG ................................ ................................ .................... 110 Test: Beat Detection on Tall T Wave Rejection Capability with Simulated ECG ................................ ................................ ................................ ............ 111 LIST OF REFERENCES ................................ ................................ ............................. 113 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 120 PAGE 9 9 LIST OF TABLES Table page 3 1 List of EEG recordings of generalized NCSE and TME patients in this study .... 30 3 2 Summary the numerical values of the centroids and the mean distance to centroid of all 11 subjects o f the study. ................................ ............................... 35 8 1 Results of the p roposed b eat d etection a ccuracy with h uman ECG from MIT BIH Arrhythmia Database ................................ ................................ ................. 105 8 2 Results of the proposed beat detection algorithm (20 ms match window) ........ 106 8 3 Results of the proposed beat detection algorithm (5 ms match window) .......... 107 8 4 Results of IMEC beat detection algorithm (5 ms match window) ...................... 107 8 5 Results of the t est of the p roposed b eat d etection algorithm on r ange and a ccuracy of h eart r ate with s imulated ECG ................................ ...................... 109 8 6 Result s of the t est of IMEC b eat d etection algorithm on r ange and a ccuracy of h eart r ate with s imulated ECG ................................ ................................ ...... 109 8 7 Results of the t est of the p roposed b eat d etection algorithm on r ange of R wave a mplitude and r ange of QRS d uration with s imulated ECG ..................... 110 8 8 Re sults of the t est of IMEC b eat d etection algorithm on r ange of R wave a mplitude and r ange of QRS d uration with s imulated ECG .............................. 111 8 9 Results of the t est of the p roposed b eat d etection algorithm on t all T wave r ejection c apability with s imulated ECG ................................ ............................ 112 8 10 Results of the t est of IMEC b eat d etection algorithm on t all T wave r ejection c apability with s imulated ECG ................................ ................................ .......... 112 PAGE 10 10 LIST OF FIGURES Figure page 3 1 10 second EEG recording of generalized NCSE and TME. ............................... 29 3 2 Scatter plot of the three nonlinear measures (STL max, Phase, ApEn) for two different patients. ................................ ................................ ................................ 31 3 3 Plot of nonlinear dynamical measures (STL max, Phase, ApEn) for all 11 patients together. ................................ ................................ ................................ 32 3 4 Centroids of each patient. ................................ ................................ ................... 33 3 5 Plot of the evolution curve of the mean distances, from EEG samples to the centroid, for 11 patients. ................................ ................................ ..................... 34 3 6 Pairwise comparison using nonlinear dynamic measures of all 6 NCSE (Patient id 1~6) and 5 TME (Patient id 7~11) patients ................................ ........ 36 4 1 The LSE pr ediction ................................ ................................ ............................. 48 7 1 Plot of the two largest components of patient 48 from interictal to ictal period ... 84 7 2 Plot of the three largest components of patient 48 from interictal to ictal period ................................ ................................ ................................ ................. 84 7 3 Plot of the two largest components of patient 51 from interictal to ictal period ... 85 7 4 Plot of the three largest components of patient 51 from interictal to ictal period ................................ ................................ ................................ ................. 85 7 5 Plot of the two largest components of patient 52 from interictal to ictal period ... 86 7 6 Plot o f the three largest components of patient 52 from interictal to ictal period ................................ ................................ ................................ ................. 86 7 7 Plot of kernel classification effect between i nterictal period and preictal period, using polynomial kernel, of Patient 48. ................................ ................... 87 7 8 Plot of kernel classification effect between interictal period and preictal period, using polynomial kernel, of Patient 51 ................................ .................... 87 7 9 Plot of kernel classification effect between interictal period and preictal period, using polynomial kernel, of Patient 52 ................................ .................... 88 7 10 Plot of kernel classification effect betwee n interictal period and preictal period, using quadratic kernel, of Patient 48 ................................ ...................... 88 PAGE 11 11 7 11 Plot of kernel classification effect between interictal period and preictal period, using quadratic kernel, of Patient 51 ................................ ...................... 89 7 12 Plot of kernel classification effect between interictal period and preictal period, using quadratic kernel, of Patient 52 ................................ ...................... 89 7 13 Eye Close Detection Results of Subjec t 1 ................................ .......................... 90 7 14 Eye Close Detection Results of Subject 2 ................................ .......................... 91 7 15 Eye Flutter Detection Results of Subject 1 ................................ ......................... 91 7 16 Eye Flutter Detection Results of Subject 2 ................................ ......................... 92 7 17 Eat & Chewing Detection Results of Subject 1 ................................ ................... 92 7 18 Eat & Chewing Detection Results of Subject 2 ................................ ................... 93 7 19 Detection Result of Intractable Epileptic Seizure in chb01_04.edf ...................... 95 7 20 Detection Result of Intractable Epileptic Seizure in chb01_16.edf ...................... 95 7 21 Detection Result of Intractable Epileptic Seizure in chb01_18.edf ...................... 96 7 22 Detection Result of Intractable Epileptic Seizure in chb01_26.edf ...................... 96 8 1 Flow chart of multiple features extraction in time domain, frequency domain and other transformed domains and feature infusion ................................ ....... 101 8 2 P lot of the selection process that individualized neighborhood widths are partly determined by the likelihood score of each candidate heartbeat. .......... 102 8 3 Plot of dynamical scaling of the neighborhood widths according to previous heart rate statistics, in a given time window ................................ ..................... 103 8 4 Plot of calculating previous heart rate statistics using exponentially decreasing weighting (an example) for a given time window. .......................... 104 8 5 Plot of a new feature built into the proposed beat detection algorithm that can help distinguish the true heart beats from a class of false ones ....................... 104 PAGE 12 12 LIST OF ABBREVIATION S AED Antiepileptic Drug(s) ApEn Approximate Entropy AR Auto Regression ARIMA Auto Regressive Integrated Moving Average ARMA Auto Regressive Moving Average BZD Benzodiazepines CNS Central Nervous System CT X ray Computed Tomography ECG Electrocardiography ECoG Electrocorticography EEG Electroencephalogram FAR False Alarm Rate fMRI Functional Magnetic Resonance Imaging FMX Frequency Modulation Extended Range GC Granger Causality LSE Least Square Estimation (model) MMSE Minimum Mean Square Error (model) MRI Magnetic Resonance Imaging NCSE Nonconvulsive Status Epilepticus PCA Principal Component Analysis P D P robability of Detection PMRS Pattern Match Regularity Score PPV Positive Predictivity Value S TD Standard Deviation PAGE 13 13 STLmax the Maximum Short Term Lyapunov exponent S TM Standard deviation Minimum S TX Standard deviation Maximum SVM Support Vector Machine T EG Teager Energy TLSE Total Least Square Estimation (model) TME Toxic Metabolic Encephalopathy PAGE 14 14 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy OPTIMIZATION AND DATA MINING IN HEALTHCARE: PATIENTS CLASSIFICATION AND EPILEPTIC BRAIN STATE TRANSITION STUDY USING DYNAMIC MEASURES, PATTERN RECOGNITION AND NETWORK MODELING By J icong Z hang August 2011 Chair: Panos M Pardalos Major: Industrial and Systems Engineering Epilepsy is the second most common chronic neurological disorder that affected over 50 million people in the world, of which one f ifth were never treated Recurrent unprovo ked seizures, which are used to characterize epilepsy, are transient symptoms of abnormal, excessive or synchronous neuronal activities in the brain. Seldom cured with medication or surgery, epilepsy is generally being contro lled clinically Patients with refractory epilepsy may feel stressed or live in depression with the uncertainty of seizure occurring time in daily lives, since most epileptic s eizures come all of a sudden. A seizure detection and/or prediction method can help identify the occurrence of seizures and alleviate the uncertainties, hen ce improving the quality of life of epileptic patients. In the dissertation the epilepti c brain state transition and epileptic brain network we re investigated with dynamic measures, granger causality, network modeling and optimization With the electroenc ephalogram ( EEG ) recordings epileptic brain network was constructed using dynamic features an d dependency or causal functions The relationship between the change of epileptic brain network and the evolution of epileptic seizure was further investigated a nd s ome indications were found for PAGE 15 15 identifying the pre ictal and ictal period s Furthermore, alternative methods using m ultiple feature s extraction and classification were applied to investigat e the epileptic brain state transition using principle component analysis, support vector machine, and skeleton methods In addition, t he electrocardiogram ( ECG ) based seizure detection methods were also investigated Last but not least, a patient classification method for the fast diagnosis of nonc onvulsive status epilepticus from toxic metabolic encephalopathies was presented PAGE 16 16 CHAPTER 1 IN TRODUCTION OF EPILEP SY To study an epileptic brain, it is necessary to understand what is epilepsy is and what are the major types of epilepsy. Definition of Neurological Disorder and Epilepsy Epilepsy is the second most commo n chronic neurological disorder that affected over 50 million people or population where neurological disorder i s the disturbance in structure or function of the central nervous system (CNS) resulting from developmental abnormality disease, injury or toxin About o ne fifth of patients wit h epilepsy were never treated Epilepsy can occur at any age of human, but is more likely to occur in young children or people over the age of 65 years. Chronic seizure conditions or r ecurrent unprovoked seizures, which are used to characterize epilepsy, are transient symptoms of abnormal, excessive or synchronous neuronal activities in the brain. However, epilepsy itself should not be understood as a single disorder, but rather as syndrome with vastly divergent symptoms (all involving episodic abnormal el ect rical activities) in the brain. Categorization of Epilepsy E pilepsy refers to a diverse set of abnormalities, with different causes different symptoms and clinical expressions and different treatment s There are many types of epilepsies : s ome types ar e subtle that can be easily missed by an observer ; and some are dramatic as to be frightening (Engel 1989 ; Engel & Pedley 1998). The e pilepsies can be categorized according to the extent of causes into : (1) i diopathic e pilepsy (2) s ymp tomatic e pilepsy and (3) c ryptogenic e pilepsy where i diopathic epilepsy has no apparent cause and may have genetic causes; s ymptomatic PAGE 17 17 epilepsy has a cause that has been identified; and c ryptogenic epilepsy has a likely cause that has not been identified. Usually e pilepsies are also categorized by clinicians on the basis of the characterization behavior and associated pattern of brain electrical discharge recorded by the EEG Alt hough t his categorization scheme ignores the underlying bases (molecula r and cellular) for the epileptic seizures, it results in a useful and complex system for diagnosis and treatment ( Commission on Classification and Terminology of the International League Against Epilepsy 1981, 1989 ) According to the areas of the brain of abnormal electr ical discharge, epileptic seizures can be categorized into partial seizures and generalized seizures W h ile the abnormal electrical discharge in partial seizure s is limited to a local area of the brain, many parts of the whole brain are involved at once in generalized seizures I t is worth mentioning that some seizures are caused by non epileptic medical conditions known as non epileptic events For partial seizure s, it can further be categorized into complex partial seizures and simple partial seizures : t he former type affects consciousness and generally is characterized by automatisms ; whereas the latter type do not affect consciousness and may cause sudden, jerking motion s of the body and affects vision or hearing In the category of partial e pilepsy benign focal ep ilepsy of childhood is idiopathic, and temporal lobe epilepsy and frontal lobe epilepsy are symptomatic or cryptogenic For generalized seizures, there are major types: absence seizures, tonic clonic seizures, myoclonic seizures, tonic seizures and atonic seizures The vast majority of generalized epilepsies are idiopathic. However, s ome generalized epilepsies, for example west syndrome and Lennox Gastaut syndrome, are symptomatic or PAGE 18 18 cryptogenic. In addition, t here is a type calle d secondarily generalized seizures, which starts as a partial seizure in one limited area of the brain, then spreads throughout the brain and becomes generalized. Sometimes the partial seizure phase is so brief that it is hardly noticed. The generalized co nvulsive phase of th e se seizures usually lasts no more than a few minutes, the same as primary generalized seizures. Besides the self limited seizure types, there are also continuous seizure types, called status epilepticus (SE). Status epilepticus is usua lly a life threatening condition in which the brain is in a state of persistent seizure. It can be subcategorized into generalized status epilepticus and partial status epilepticus, including generalized tonic clonic status epilepticus, clonic status epile pticus, absence status epilepticus, tonic status epilepticus, myoclonic status epilepticus, and epilepsia partialis continua of Kojevnikov, Aura continua, limbic status epilepticus, hemiconvulsive status with hemiparesis, respectively Furthermore, some seizures are non epileptic. An example is psychogenic seizures which are psychological in nature and can happen in people diagnosed with or without epilepsy. S ometimes psychogenic seizures look like true epileptic seizure s but no abnormal electrical acti vities occur There are also some medical conditions that have symptoms similar to epileptic seizures such as narcolepsy, heart stroke, cardiac arrhythmia, and low blood sugar. Epileptic seizures and non epileptic events can be distinguished by diagnostic tests, typical examples as continuous EEG monitoring functional MRI and MRI. PAGE 19 19 CHAPTER 2 C AUSES DIAGNOSIS AND TREATM ENT OF EPILEPSY Neurobiological Causes of Epilepsy Because s eizures can occur in normal brain reflecting normal responses to injury or stress occurren ce of an isolated single seizure does not necessarily reflect brain path ology. However, repetitive seizures can indicate some chronic pathological condition of the central nervous system During the time between the seizures, or called interictal period, the brain functions in an apparently normal fashion. For most types of epileps y, the occurrence of seizure is episodic and apparently unpredictable, which makes it particularly difficult to understand the abnormality of pathological condition Rather than categorizing epilepsies on the basis of the characterization behavior during ictal period and the pattern of brain electrical discharges recorded by EEG, i t is also possible to categorize many epileptic disorders on the basis of the neurobiological cause with advances in the understanding of mechanisms of seizures and of the processes by which a chronic seizure state is acquired A general classification scheme, that covers most types of epilepsy includes three categories: genetic d isorders; developmental disorders; and acquired epilepsies. Genetic Disorders It is believed that the mutations of single genes are the causes of a good number of epilepsies (Berkovic & Scheffer 1999). As is known, many genes code for ion channels or for neurotransmitter receptor proteins, which affect brain excitability A bnormalities in th o se proteins can part ially explain some seizure discharge s In addition, s tructural brain abnormalities, in lissencephaly or tub erous sclerosis for example (Spreafico et al. 1999) are causes of o ther epileps ies which are associated PAGE 20 20 with single gene mutations. A t present, the understanding of a bnormalities in the gene products that result in developmental disorders is still at an early stage and there are functions of the i mplicated genes that need to be elucidated. More commonly, some epileps ies such as c hildhood absence epilepsy (CAE) and juvenile myoclonic epilepsy (JME), are disorders associated with multiple genes. This provides the possibility of applying genetics in the diagnosis of epilepsy specifically determining seizure susceptibility in the brain It is worth mentioning that, c haracterizing and identifying the genes involved in multigenic seizure disorders is a major challenge in this field because it is critical to understand which genes play a significant role in determining seizure predisposition. Further gaining an insightful understand ing of the cellular pathways through which the genie abnormalities lead to neurophysiologic dysfunction is important Overall, i dentifying epilepsy genes is a critical start in developing new drugs or other treatments with specific molecular targets though ther e is still a long way to go to completely understand the aberrant central nervous system function Developmental Disorders The genetic and molecular signals in a normal brain ha ve been investigated and understood in the last few decades benefited from the development in neuroscience research (Sanes & Donoghue, 2000) Neurological disorders, specifically epilepsy, may be caused by any abnormalities in the process of cell proliferation, migration, differentiation and establishment of connections (Schwartzkroin et al. 1995). T he incidence of neurological dysfunction is low g iven the potential for developmental mistakes Especially, the percentage of seizure occu rrence in the adult population is con siderably lower than that in the children partly due to abnormalities in the brain development (Spreafico et al. 1999), and partly due to the seizure prone properties of PAGE 21 21 the immature brain (Nehlig et al. 199 2 ). It is observed that t he epilepsy incidence in the children with other developmental disorders is much higher than that in normal chi ldren This reflect s some underlying brain patholog y Considering (1) seizure incidence among the children is high er than that among the adults ; (2) epilepsies may be caused by abnormalities of genetics or brain development ; and (3) the human brain is vulnerable to perinatal trauma ; epilepsy is considered primarily as a set of disorders of childhood by a large number of epileptologists Some phenomena are observed surely not by coincidence : e pilepsies are often induced by insults during the young ages, and seizure syndromes are usually highly correlated with developmental brain abnormalities A commonly accepted explanation to th ose phenomena is the highly plastic nature of t he immature CNS. Overall speaking t he immature CNS is more resilient compared with the adult CNS where brain trauma may lead to cell death or loss of function s T he i mmature neurons are not as sensitive as the mature neurons to injurious stimuli. Specifically, t he immature brain are capable of produc ing new neurons to a much greater extent than the adult CNS c an and the i mmature nerve cells make new connections more readily during the disruption of the normal organization. Alt hough the plasticity of the immature brain helps preserve a large amount of behavioral function s t he new circuits may be established by injured neural elements and may not be entirely normal. Thus, neuron al injury and in correct connections are considered two major contributors in the development of epileptic electrical discharges during the seizures Acquired Epilepsies Compared with genetic and development al disorders, the acquired epilepsies are solely induced by some particular insult s T he brain was o riginally normal and would othe rwise be normal if the insults had n o t happened. Such insults could be medical PAGE 22 22 condition s such as an infection or a tumor, and/or a particular traumatic injury. The hypothes i s is that : the injured brain may lose some critical elements which control the electrical excitability of the CNS, and may create inappropriate new connections Howe ver, whether or not some insults will lead to an epileptic state is individual based. A stimulus could be a trigger for epil eptogenesis on one subject and not on another. S ubtle factors differences genetic or environmental may exist concerning seizure susceptibility and abnormalities of brain development All those provide a background determining which insults plays out on a subject. For instance, a preexisting injury though subtle, may make the chance higher that an individual with a serious infection develops epilepsy at a later time. This implies that epilepsy is a progressive disorder (Heinemann et al., 1996), developin g with active neuronal processes. And this concept is often discussed in the study of pediatric epilepsies Thus, an y initial ly subtle insult has the possibility of kindl ing into a serious and troubling condition over a period of time Diagnosis and Treatment of Epilepsy Accurate diagnosis is a critical first step for the treatment of epilepsy. There are many disorders that have chang es in behavior and clinical symptoms similar to those of epilepsy which makes it difficult to determine if a person is having an epileptic seizure and to diagnos e the type of seizure or epilepsy syndrome One of the most important information that is need ed for the diagnosis is: what happened during the seizure occurrence or the ictal period Since how the patients, the caregiver s and other witnesses give the information to the doctor and the healthcare professionals is important Yet even with accurate description s of events, other records and/or tests are also needed to learn more about the brain: what is causing the events underlyingly and where exactly the problem is located. T ypical PAGE 23 23 records and tools include but are not limited to medical history, blood tests, continuous video EEG monitoring, ECoG mo nitoring and brain imaging tests such as functional MRI CT and MRI. Th e EEG and ECoG recordings give information about the electrical activity of the brain and the imaging tests can tell in some sense what the brain looks like. Different sources of information including the individual s feeling description and the models how seizures may affect the way the brain works are put together for diagnosis of epilepsy After a diagnosis of seizure or epilepsy of sufficient conf idence, the best form of treatment need to be selected Seldom cured with medication or surgery in difficult cases epilepsy is generally being controlled clinically. The main methods for epileptic seizure control are: ( 1 ) the use of antiepileptic drugs (AED) ; ( 2 ) surgical removal of the seizure focus; ( 3 ) electrical stimulation If a continuing tendency to have seizures is diagnosed, a regular use of AED will be prescribed usually Other methods including a special diet, brain surgery, or vagus nerve s timulation will be tried i f AED are not effective and successful. Especially, if seizure s are localized in one area and is caused by an underlying correctable brain condition, surgery may be able to alleviate the severity of epilepsy or even stop seizure s T he goal of epilepsy treatment is to prevent seizures, or reduce the frequency and alleviate the severity, while minimizing the side effect, to increase the quality of life of epileptic patients. Patients with refractory epilepsy may feel stressed or li ve in depression with the uncertainty of seizure occurring time in daily lives, since most epileptic seizures come all of a sudden. An automatic seizure prediction method can help anticipate the PAGE 24 24 occurrence of forthcoming seizures and alleviate the uncertai nties, hence improving the quality of life of epileptic patients. PAGE 25 25 CHAPTER 3 PATIENTS CLASSIFICATION BETWEEN GENERALIZED NCSE AND TME USING EEG QUANTIFICATION FOR DIAGNOSIS OF STATUS EPILEPS TICUS Significance Non convulsive status epilepticus (NCSE) is usually defined as an epileptic state lasting 30 minutes or more with some clinically evident change in mental status or behavior (unless comatose) associated with ictal activity on EEG (Brenner, 2002). Nonconvulsive seizures during NCSE do not contain convulsive motor activity and may be repeated discrete seizures or more continuous prolonged seizures. NCSE is particularly difficult to diagnose in patients presenting with stupor or coma as there are often little or no specific clinical symptoms Misdiag nosis or d elay to diagnosis of adult NCSE may engender morbidity and mortality (Bearden et al., 2008). Often NCSE requires an electroencephalogra m (EEG) recording to be examined (Towne et al. 2000 ; Bearden et al., 2008). A n EEG study of 236 patients concluded that 8% of comatose patients without clinical seizure activity had NCSE that would not have been diagnosed without EEG recording (Towne et al. 2000). Nevertheless, e ven when EEG recordings are performed diagnosis of NCSE can still be difficult due to the presence of other disorders with similar EEG waveforms (Brenner 2002). In particular, some non epileptic encephalopathies, such as toxic/metabolic encephalopathy (TME), are indistinguishable by clinical symptoms, and can produce si milar EEG p atterns as NCSE. Most EEG patterns of NCSE are generalized epilepti form waves or focal rhythmic ictal transformations that are not difficult to identify as electrographic seizure activity (Geiger & Harner, 197 8). Unfortunately, EEG waveforms of generalized triphasic like waves or diffuse semi rhythmic delta activity can be seen in patient with either TME or PAGE 26 26 generalized NC SE (Figure 3 1) Treatment with benzodiazepines (BZDs) can reduce electroencephalographic seizure activity and/or improve clini cal state, which may be useful diagnostically. However, the diagnostic utility of BZDs can be limited because BZDs may transiently suppress triphasic waves due strictly to metabolic encephalopathy as well as triphasic like waves seen duri ng NCSE (Fountain & Waldman, 2001) Furthermore, NCSE and TME can coexist, and at present, combining clinical judgment and experience is the only effective method to diagnose NCSE in equivocal cases (Bearden et al., 2008). Rapid diagnosis and treatment of NCSE is desirable as untreated NCSE can cause irreversible damage to the brain (Jirsch & Hirsch, 2007) In the past, quantitative EEG analysis has been investigated as an aid in diagnosis of some psychiatric and neurological disorders (Pardalos et al., 2004; Liu, 2008). One approach to expedite diagnosis would be to introduce a fast EEG classification algorithm. The purpose of this study is to identify EEG characteristics that can distinguish NCSE from non epileptic encephalopathy. This may assist physicians in making a more rapid diagnosis of NCSE or non epileptic encephalopathy. In our study, data mining algorithms are applied to extract complex and intrinsic features from large datasets in multi channel EEG recordings, since the EEG patterns of NCSE can be visually similar to those of non epileptic encephalopathy. Approach Data description: quantitative analysis was performed on EEG recorded from male and female adult patients having been diagnosed with NCSE (generalized or focal) or with a non epileptic encephalopathy, specifically toxic/metabolic encephalopathy (TME). The clinical diagnosis and EEG interpretation were conducted by a board certified PAGE 27 27 electroencephalographer. Eleven human subjects were involved in this study, in which six were diagnosed wit h NCSE, whereas the other five were diagnosed with TME. Length of the recordings for each patient can be seen in Table 3 1 EEG recordings were acquired using the international 10 20 electrode placement. For each subject, there are 21 recording channels (F p1, Fp2, F3, F4, C3, C4, A1, A2, P3, P4, O1, O2, F7, F8, T3, T4, T5, T6, Fz, Cz, and Pz). An effective approach in order to shed light into the underlying non linear structure of time series, in our case EEG recordings, is analysis in the state space (also known as phase space). In this approach EEG recordings are mapped from the time domain onto an n dimensional space of states. In this study, method of delays was used for state space reconstruction (Takens, 1981). If we consider the time series (ro w vector of length m ) we can construct the m by n delay matrix (where is the delay parameter and n is the embedding dimension). Then, every delay vector in A is a point on the n dimensional state space.. In the analysis of dynamic systems the attractor is a very useful visualization tool of the signal because it can reveal structure and patterns that are well hidden in the time domain. In order to quantify state changes and attractor geometrical patterns one can employ several well defined nonlinear dynamic measures. For the purpose of our study we used three such measures, namely: Short term Lyapunov exponents (STLmax): Is a nonlinear state space based metric that quantifies the instability of a dynamic system. S hort term Lyapunov exponents have been used extensively in nonlinear EEG analysis mostly for seizure prediction (Iasemidis et al., 2001, 2004, 2005). Phase of the attractor: Also called angular/frequency phase, estimates the rate of change of the stability for a dynamic system. Complementing the Lyapunov PAGE 28 28 exponents phase is a broadly used nonlinear state space measure that is employed to characterize chaotic systems. Approximate entropy (ApEn): Is another phase space statistic that quantifies the complexity of a dynamic system (Pincus, 1991). In order to analyze the continuous EEG recordings we segmented the whole time series into epochs of a finite time window of 10 seconds at a sample rate of 200 Hz. For each 10 second epoch, values for short term Lyapunov exponents, phase of attractor and approximate entropy were computed. Points were plotted in a three dimensional space (should not be confused with the n dimensional embedding space used for calculating the three measures). Plots for the three nonlinear dyn amical measures show very good separation betw een patients with NCSE and patients with TME. In Fig ure 3 2 we can see an example 3D scatter plot of the three measures (for one channel of EEG recordings from two different patients). X axis, Y axis and Z axi s correspond to approximate entropy (ApEn), phase of attractors (phase/angular frequency), and short term Lyapunov exponents (STL max) respectively. The formation of two clearly distinct and well separated clouds of points led us to suggest an automated al gorithm for real time detectio n of NCSE based on the centroids and the mean radius of the clouds, defined by the 3 dimensional scatter plot points. We define as the centroid of N points ( , (3 1) PAGE 29 29 Figure 3 1 10 second EEG recording of generalized NCSE and TME : A) an equivocal year old male patient with subdural hematoma, who began having focal NCSE which evolved into generalized NCSE ; B ) a 10 second EEG recording of triphasic waves in a 56 year old male TME patient with severe acute liver failure. Visual similarities regarding the triphasic waves of TME and tri profound. PAGE 30 30 Table 3 1. List of EEG recordings of generalized NCSE and TME patient s in this study Patient ID 1 2 3 4 5 6 7 8 9 10 11 Diagnosis NCSE NCSE NCSE NCSE NCSE NCSE TME TME TME TME TME Recording length(min) 27 33.5 53.5 25 29 63.8 56.16 13 21.66 22.16 25.33 T otal ly 14 patients we re recorded and involved in this study 6 with generalized NCSE, 5 with TME, and 3 with focal NCSE. Since focal NCSE is not difficult for clinical diagnosis, only the information of t he first two types of patients is shown in the table. PAGE 31 31 The mean Euclidean radius of the N point cloud is: (3 2) Figure 3 2 Scatter plot of the three non linear measures ( STL max, Phase ApEn ) for two different patients. Eac h circle (NCSE) or star (TME) corresponds to a 10 second EEG epoch Example shows recordings of a single channel. In order for each of the features to contribute equally, all three features are normalized in scal e [0, 1]. Based on these, the real time algorithm can be described in the following well defined steps: Every 10 seconds compute one point on the 3D space spanned by the three nonlinear measures; Update the centroid and the mean radius distance using eq ua tions 3 1 and 3 2 correspondingly; If the centroid coordinates do not change more than a predefined quantity issue diagnosis based on the last centroid coordinates. Results Analysis carried out in the whole datasets (11 patients, 21 EEG channels) shows results similar with the preliminary example in Fig ure 3 2 In Fig ure 3 3 one can see a PAGE 32 32 cumulative scattered plot with the three nonlinear dynamical measures for all the 11 patients plotted together, using 21 channel averages in the 3 D feature space. Although the cloud of points is very dense, one can see that there is clear separation between the two classes. Figure 3 3 Plot of non linear dynamical me asures (STL max, Phase ApEn) for all 11 patients together. E very circle (NCSE) or star (TME) corre sponds to a 10 second EEG epoch sample Plot uses 21 channel average. We further observed that STLmax and phase of attractor play a major role in this separation, whereas approximate entropy contributes very little. This makes it possible to quantify the s described in Eq uation 3 1 and Eq uation 3 2. First, for each patient, we computed the centroids taking into consideration epoch samples generated from all 21 EEG channels. Scatter diagram of the centroids for all 11 patients in 3 D feature space is shown in Fig ure 3 4 ( A ) It is worth mentioning that even if the centroids are projected into two dimensional subspace spanned by STLmax and phase of attractor, there is still clear separation between generalized NCSE and TME patients, as shown in Fig ure 3 4 (B). Besides, the PAGE 33 33 centroids of the six NCSE patients are much more concentrated than those of the five TME patients. This indicates that NCSE can be detected with even fewer (two instead of three) dynamical measures and computational effort. (A) (B) Figure 3 4. Centroids of each patient (circles correspond to NCSE and stars to TME). We can see that two classes are easily separable either by A) using all three nonlinear dynamical measures, or by B) just using STLmax and phase of attractor. Centroids are computed using the full length of the EEG recording available for each patient. Next, for each patient, the mean distance from the already recorded EEG epochs to their centroid in t he feature space is computed, as defined by Eq uation 3 2. In Fig ure 3 5 for each patient, the evolution curve of the mean distance is plotted. After 20 minutes of EEG recording, the mean distances of the NCSE patients all drop below 0.06, converge and sta y there, whereas the mean distances of the TME patients do not converge rapidly and keep above 0.08 after starting EEG recording for a while. This indicates that epoch samples of TME are more scattered than those of NCSE in the feature space. The significa nt difference can be used as an effective criterion to judge whether a patient should be classified into generalized NCSE or TME. Numerical results of the experiments are summarized in T able 3 2. PAGE 34 34 Figure 3 5 Plot of the evoluti on curve of the mean distances, f rom EEG samples to the centroid, for 11 patients. Each curve corresponds to one patient (blue thick: NCSE; red thin: TME). The distance is calculated in the two dimensional space spanned by STL max and phase of attractor. M ean distances of NCSE patients (blue thick) after a time frame of 10 20 minutes all converge to some low value. On the other hand, mean distances of TME patients (red thin) tend to increase or oscillate and does not converge fast. Table 3 2 and Fig ure 3 4 already indicate the STLmax and phase of attractors are sufficient to separate those two classes of patients. We then performed a statistical t test in order to further justify this observation. In Fig ure 3 6 one can see the results of paired t tests in e ach of the 21 EEG channels between every pair of patients using short term Lyapunov exponent, phase of attractor, and approximate entropy, respectively. In the paired t test, the minimum length of the EEG recordings of the two patients is used for the comp arison. The brightness of each cell in Fig ure 3 6 represents how many rejections (of null hypotheses) there are out of the 21 EEG channels: the white color corresponds to 21 rejections, and the black to zero. The null hypothesis of the paired t test is tha t: two matched sample sets from any two patients, in the vectors X and Y, come from distribution with equal means. X Y, the difference of X and Y, are assumed to come from a normal distribution with unknown variance. We expect the null PAGE 35 35 hypothesis will be a ccepted for patients who belong to the same group and rejected for patients that belong to different groups. Because t tests between patient pair (A, B) and (B, A) have the same results, all subplots in Fig ure 3 6 are diagonally symmetric. Table 3 2 Summ ary the numerical values of the centroids and the mean distance to centroid of all 11 subjects of the study. Patient id diagnosis (ApEn) (STLmax) (Phase) 1 NCSE 0.53 0.08 0.17 0.05 2 NCSE 0.46 0.08 0.16 0.09 3 NCSE 0.55 0.12 0.13 0.05 4 NCSE 0.61 0.07 0.14 0.07 5 NCSE 0.60 0.09 0.13 0.05 6 NCSE 0.59 0.04 0.13 0.05 7 TME 0.50 0.58 0.56 0.09 8 TME 0.50 0.34 0.30 0.20 9 TME 0.52 0.23 0.43 0.13 10 TME 0.54 0.63 0.46 0.10 11 TME 0.59 0.57 0.29 0.16 In th is table, Euclidean distance (norm 2) i s used for calculating the mean distance. For the two nonlinear dynamical measures: short term Lyapunov exponent and phase of attractor (Fig ure 3 6 upper left and right), the number of null hypotheses rejections in pairs within the six NCSE patients is very small (colors are dark in the 6 by 6 sub matrix in upper left corner); whereas most of the 21 null hypotheses in any patients pair between NCSE and TME are rejected (colors are light in the 5 by 6 sub matrix in lower left corner). Th is reveals the reason why in Fig ure 3 2 and Fig ure 3 3 cloud of epochs of the NCSE patients are grouped together and are apart from epochs cloud of TME. Furthermore, in pairs within the five TME patients, the number of null hypotheses rejections is larger than that of pairs within the six NCSE patients. This implies the TME patients are more scattered by themselves than the NCSE patients, observed from the PAGE 36 36 two dimensional feature space spanned by short term Lyapunov exponent and phase of attractor. This imp lication is consistent with results in Fig ure 3 5 Figure 3 6 Pairwise comparison using nonlinear dynamic measures of all 6 NCSE (Patient id 1~6) and 5 TME (Patient id 7~11) patients. Upper left is STLmax; upper right is phase of attractor; lower left is approximate entropy. For each pair of patient i and j t tests are performed in 21 EEG channels respectively. The value (brightness) of each cell ( i j ) is the number of rejections of null hypothesis, ranging from 0 (darkest black) to 21 (brightest white). The null hypothesis is expected to be accepted for patients wh o belong to the same group and rejected for patients that belong to different groups. Patient 8 (TME) and Patient 9 (TME) have EEG recor dings of very short length (Table 3 1). As is mentioned before, the minimum length of the EEG recordings of two patients is used for the comparison in the paired t test. However, when recording time is not long enough, epoch samples of some NCSE patients may still be scattered in the three dimensional feature space of nonlinear dynamical measures, and may be close to epoch samples of TME Patient 8 or 9. This gives an explanation why the t test in pairs between NCSE Patient (id 1~6) and TME Patient 8 or 9 does not reject the null PAGE 37 37 hypothesis in some EEG channels. As shown in Fig ure 3 6 in the upper left and upper right subfig ure, some cells in Row 8 or Row 9 (Column 8 or Column 9 accordingly) are slightly dark. For the third nonlinear dynamic measure: approximate entropy (Fig ure 3 6 lower left), the t test indicates that its differentiating effect between NCSE and TME is not optimal Discussion a nd Conclusion Distinguishing generalized equivocal EEG patterns (semi rhythmic delta and triphasic like waves) which can be seen in either NCSE or TME (Toxic/metabolic encephalopathy) is a hard clinical problem. Currently, the only effective method to distinguish generalized NCSE from TME in equivocal cases is combining clinical judgment and experience (Bearden et al., 2008). From the neuroscience perspective it would be interesting to uncover and model the basic mechanism that under lines the evolution of mental status changes of generalized NCSE patients, though at present the properties of epileptic brain are still not well understood. Nonlinear dynamics have been used in seizure detection and prediction in epilepsy over the last de cade. The maximum Lyapunov exponent has been proposed as an indicative metric of the abrupt transient drop in chaoticity in EEG recording before seizure onset Another example, approximate entropy, which quantifies the unpredictability of fluctuations in a time series, has been proposed for EEG epileptic seizure detection (Srinivasan et al., 2007). In this study, three nonlinear dynamic measures (STLmax, phase and ApEn) were applied in order to distinguish between generalized NCSE and TME patients. Findings based on a dataset of 11 subj ects not only suggest that two of those nonlinear PAGE 38 38 dynamics, STLmax and phase, could potentially be used as a bed side assistive tool for generalized NCSE diagnosis in intensi ve care units; but also provide additional evidence (through application other than seizure prediction or detection) that STLmax and phase of attractors can be useful in quantifying neurophysiological states related to epilepsy. Hopefully, this may trigger fruitful discussion among neurologists, clinicians and scientists about the exact role of nonlinear dynamics in the field of epileptic disorders. Although groups of researchers relating epilepsy with dynamic systems (Milton & Jung, 2003) are active, at present there is no commonly accepted theory, and the mechanism of this connection remains unknown. Based on findings in this study, we suggest that epileptic activity is highly associated and can be modeled u sing dynamic system analysis. As mentioned above, STLmax and phase of attractor were found most likel y effective for differentiating NCSE patients from TME. In the feature space spanned by those two nonlinear dynamics, for each patient, the mean distance from the centroid to all EEG samples is a useful metric to measure the degree of concentration of EEG epochs. In the 11 subjects in our study, samples of a TME patient were always more scattered than samples of a NCSE patient. Thus, t he mean distance can be utilized to assist diagnosis of NCSE The paired t tests help quantify the efficacy of the three non linear dynamics (STLmax, phase and ApEn) when differentiating NCSE from TME. The paired t test verifies that STL max and phase of attractor can be the base for an automatic machine learning classifier. Such a classifier would use dynamical support vector machine or artificial neural network, and could be trained to do the classification of EEG recordings. PAGE 39 39 For this purpose a larger training set is needed so that classification error can be minimized. Though the small number of patients (11) used in this pi lot study precludes us from reaching statistically significant conclusions, we are encouraged that the two discriminators we identified were able to accurately distinguish all of our patients. These patients were selected because they had particularly dif ficult equivocal generalized EEG patterns that experienced electroencephalographers found very difficult to classify as NCSE or TME. We anticip ate following th is study with a p erspective study using larger number of subjects to allow testing for statistica lly significant results. This algorithm may be useful during bedside real time EEG recording or during post recording analysis. The three dynamics and derivative metrics proposed in this study were computed using an ordinary laptop computer. The computatio n time is very close to real time computation and the computation can be done simultaneously as the EEG recording is done. This may assist physicians in generating a more rapid accurate diagnosis of NCSE or non epileptic encephalopathy (Zhang et al., 2010 ) PAGE 40 40 CHAPTER 4 MINIMUM PREDICTION E RROR MODELS Minimum error p rediction is crucial in time series analysis. Especially, prediction based on multiple time series, such as EEG and ECoG, plays an important role, in the application of epileptic patients monitoring, seizure control and treatment A prediction is a statement that a particular event will occur in the future, which can refer to the estimation of unknown situations in time series, cross sectional or longitudinal data. Prediction implies at lea st two factors: importance and difficulty (Stevenson, Howard, ed., 1998). Risk and uncertainty are crucial to prediction, thus it is important to minimize the prediction errors. There are no specifically assumed statistical distributions in the prediction models in this article. In the first part of this chapter, two prediction models are discussed: the Least Square Estimation (LSE) model, and the Minimum Mean Square Error (MMSE) model. Both of these two models aim at minimizing the square errors in their r espective senses: (1) sum of squares, and (2) mean square. LSE is known as a linear deterministic model that can be used to obtain approximate solutions of over determined systems. The basic problem is: a limited length time series {Z (i)} is being predict ed from another limited length time series {X (i)}, which can be generalized to infinite length time series case. B y minimizing the sum of error squares in time series prediction, the principle of orthogonality is obtain ed as a necessary condition of reach ing the global minimum error. It is further shown that: the energy of the observed time series {Z (i)}, is the summation of the energy of the LSE estimation time series {Z LSE estimation (i)} and the energy of the LSE error process time series {e (t)}. Furt hermore another way of understanding the LSE model is exhibited: a linear projection onto the column space of the data matrix A PAGE 41 41 where the columns of A are composed of the time delay vectors of the time series {X (i)}. I n addition a stochastic model is c onverted to : MMSE, which is also linear, to minimize the mean square error in prediction We further show that the MMSE model is a special case of ARMA / ARIMA models. LSE Model Basic LSE Model At the beginning of this article, the Least Squares Error (LSE ) model is presented. In a causal system, consider two time series X ( i ) and Z ( i ) ( i = 0, 1, 2, 3, ...), where X ( i ), X ( i 1), X ( i i M + 1) are the unobserved underlying variables and will influence the observed time series Z ( i ) in a linear way, expressed as following: (4 1) w here the 1) are parameters of the LSE model and is the model error. The measurement error is unobservable and is used to count the Let: (4 2) (4 3) The measurement error process in the LSE model is assumed white with zero expectation and the same variance: (4 4) PAGE 42 42 A nd (4 5) w here In the LSE model, given the obser ved data, we select / design the tap weights ( k M 1) to minimize the sum of error squares (viewed as error energy). The objective function is as following: (4 6) w here s and t ( ) are the index limits at which the error minimization occurs. Once the tap weights ( k M 1) are selected, they are fixed as parameters in the interval Not losing generality, suppose index range of the observed data is [1, N ], and set s = M t = N Then, we get the observed input data matrix (Makhoul, 1975; Markel and Gray, 1976): w here the arrays above, from the first column (M) to the last column (N) correspond to the unobserved input variable vectors. Hence, the cost function in Eq uation 4 6 becomes: (4 7) PAGE 43 43 Assume that is differentiable, then the necessary condition for minimizing is: (4 8) Apply Eq uation 4 7 to Equation 4 8 we get: (4 9) Therefore, a necessary condition for to reach its global minimum is: (4 10) Let the tap weights that is optimized to operate the LSE condition have the special values Then, we define as: (4 11) Multiply both sides of Eq uation 4 10 by and then sum the results over the value of 1), then we get: (4 12) Interchange the order of summation, we get: PAGE 44 44 (4 13) Use Eq uation 4 11 to substitute the inside summation item in Eq uation 4 13, we get: (4 14) Eq uation 4 14 is a necessary condition for global minimization of which is called the principle of orthogonality. In Eq uation 4 14, (defined in Eq uation 4 11 ) is considered the least square (LS) estimate of the desired (observed) response Z (i). Hence can be written as then Eq uation 4 14 can be reformulated as: (4 15) Taking the time average of the left hand side of Eq uation 4 15, we find that the cross correlation of two time series and is 0. So, the LSE estimate of the observed response Z (i) represented by and the LSE model error process are orthogonal over time. (4 16) If we define: PAGE 45 45 B y using the principle of orthogonality (Equation 4 15 ), we get: (4 17) Equation 4 17 tells: The energy of the observed time series Z (i) is the summation of the energy of the LSE estimation time series, and the energy of the LSE error process time series. All the three items in Equation 4 17 are nonnegative. Another A pproach of U nderstanding the LSE Model Since in Eq uation 4 11 is also written as in Eq uation 4 16, according to Eq uation 4 11 and Eq uation 4 16, we get: (4 18) Substitute Eq uation 4 18 into Eq uation 4 14, we get: (4 19) Reorganize we get: PAGE 46 46 (4 20) Here we define: (4 21) (4 22) Then Eq uation 4 20 can be rewritten as: (4 23) Let the M by M matrix be : (4 24) Let the vector (M by 1) be: (4 25) Let the vector (M by 1) be: (4 26) Then, E quation 4 23 can be rewritten as: (4 27) Assume is nonsingular, then: PAGE 47 47 (4 28) Eq uation 4 28 is the design of a linear least square estimation (LSE), where is the inverse of the time average cross correlation matrix and is the cross correlation vector between unobserved time series X and observed time series Z When noise is assumed to be white Gaussian distributed, its counterpart filter design is the Wiener Hopf Equation. According to Eq uation 4 22 and Eq uation 4 24, can be decomposed into : (4 29) w here (4 30) Let (4 31) T hen, according to Eq uation 4 30 and 4 31 and the definition of in Eq uation 4 25 and 4 21 (4 32) Eq uation 4 27 can be rewritten as the following using the decomposition expression in Eq uation 4 29 (4 33) Integrating Eq uation 4 32 and Eq uation 4 33 (4 34) Therefore, PAGE 48 48 (4 35) Let (4 36) Substitute Eq uation 4 35 into Eq uation 4 11, we get: (4 37) Noting that is the projection operator onto the linear space spanned by the columns of the data matrix A Figure 4 1 T he LSE prediction can be considered as the observed variable vector Z projected from the ( N M +1) dimensional whole space onto the M dimensional subspace spanned by columns of data matrix A Hence, the LSE estimation vector, which is expressed by or is the orthogonal projection of observable variable vector onto the linear space spanned by the columns of the data matrix A Suppose M < N M + 1 (Stewart, 1973) (Fig ure 4 1). If data matrix A is full rank, then, its column space is an M dimensional sub space of the (N M+1) dimensional full space. Assume that, the data matrix A is already known with no uncertainty. In pract ice, A could be unknown. In this case, a set of training data with limited length will be used to estimate the LSE parameters. Some further properties of LSE are PAGE 49 49 investigated by other literatures (Miller, 1974; Goodwin & Payne, 1977; Hanson, 1995; Hayes, 2 008). In an alternative way of thinking, X ( t ) can be considered as another time series observed. Then, we interpret the problem as: the present states of time series Z ( t ) can be predicted by the present and previous states of another time series X ( t ) using the LSE model. MMSE Model Prediction Next, we introduce the Minimum Mean Square Error (MMSE) prediction, described by the conditional expectation (Hamilton, 1994; Priestley, 1981): (4 38) If we want a linear form p rediction, should have the following form (Whittle, 1983; Priestley, 1981): (4 39) Suppose we identify a particular model for a given time series, we need to estimate the model parameters and wish to compute forecasts from the fitted model ( ) instead of the true model ( ) which we do no t know If we use a quadratic loss function, the best way to compute a forecast is to choose 1) to be the conditional expected value of on the model, and use information available at time N (4 40) Then, the MMSE forecast with MA model of possibly infinite order would be (Box et al., 1994): PAGE 50 50 (4 41) where the future values are replaced by zero. More generally, the MMSE forecast from an Auto Regressive Moving Average (ARMA) model (Shamway, et al., 2000) can be computed by: (4 42) where the future value of Z is replaced by zero and the future values of X is replaced by the conditional expectation (Box et al. 1994). Similarly, if we use the Auto Regressive Integrated Moving Average (ARIMA) model (Shamway, et al., 2000), then the M MSE forecast would be: (4 43) w here the future value of Z is replaced by zero and the future values of X is replaced by the conditional expectation of ARMA and ARIMA models can convert to each other, though generally ARMA models are used for accumulative value forecast, while ARIMA models are used for differential value forecast. For special cases, e.g. ARIMA (0, 1, 1), it is a simple exponential smoothing. Recursive calculation is commonly used for the q step MMSE forecast at time p using ARMA or ARIMA models. PAGE 51 51 Other Prediction Methods To handle non stationarity and seasonality, the prediction process need to select a suitab le model for a given time series. The prediction based on the more general class of ARIMA or Seasonal Auto Regressive Integrated Moving Average (SARIMA) models is called the Box Jenkins prediction (Box et al., 1994). The model involves in an iterative proc edure with : (1) F ormulating, (2) fitting, (3) checking and adjusting. Depending on how the first few observations are treated, there exist different prediction procedures for fitting ARIMA model (e.g. maximum likelihood and conditional least square). But only for short series (no more than several hundred points), the choice of procedures is important and makes a significant difference. It is shown that, for short series even when asymptotically unbiased, parameter estimates are likely to be biased (Ansley and Newbold, 1980). Furthermore, different software package, using different estimation routines, can produce model parameter estimates with non trivial differences (Newbold et al., 1994). Hence it is a wise choice to use software which specifies exactly what estimation procedures is adopted. The way that trend is removed before applying the Box Jenkins approach can be vital for non stationary time series. Especially, the order of differencing can be crucial if differencing is applied. Alternative methods of removing trend, prior to fitting an ARIMA model, may lead to better forecasting ( Makridakis and Hibon, 1997). There are other categories of prediction methods not included in this chapter, like the ensemble prediction, simulation methods or the judgment al methods incorporating intuitive judgments, opinions and subjective probability estimates (e.g. composite prediction; Delphi method; Scenario building). PAGE 52 52 Different models have their respective limitations and advantages in difference scenarios. When selec ting the fit prediction models for epileptic seizure prediction, comprises need to be made to tradeoff the robustness, performance, feasibility and complexity (Zhang et al., 2011) PAGE 53 53 CHAPTER 5 CAUSAL RELATIONS BET WEEN TWO OR MULTIPLE TIME SERIES The causal re lations are useful in the modeling of directed network s of epileptic brain and in the epileptic focus localiza tion for epilepsy surgery In this chapter, a wide sense stationary, purely nondeterministic time series is studied, for an issue closely related to prediction: the causal relations between multiple time series. We start from two variable (time series) and three variable (time series) models. The causal relations in two or multiple time series are discussed in comparison, in statistical senses. W e first disc uss ( 1) the simple causal model and ( 2) the instantaneous causal model, between two time series {Xt} and {Yt } (Granger 1969). We assume that: the time series are wide sense stationary and purely nondeterministic. In the simple causal model, the previous states of {Xt} and {Yt} are used, to predict the present state of {Xt} and {Yt}, respectively. In the instantaneous causal model, the previous states of {Xt}, as well as the previous and present states of {Yt}, are used to predict the present state of {Xt}; and vice versa. Though this is in stochastic sense, the principle of orthogonality is still kept for minimizing prediction error. Then t he three variable models are discussed briefly in frequency formulation (Granger 1969) Finally causal relations between one group of time series and the other group of time series are studied in frequency domain (Geweke 1982) The underlying principle of orthogonality is kept for error minimization. Two Variable (Time Series) Models Let and be two wide sense stationary, purely nondeterministic time series with zero means. The simple causal model (Granger, 1969) is (5 1) PAGE 54 54 (5 2) w here and are two zero mean serially uncorrelated noise process (time series) with for any t and s. Since Since In the simple causal model, the previous states of {X t } and {Y t } are used to predict the present state of {X t } and the present state of {Y t }, respectively. A more general model: instantaneous causal model (Granger, 1969), is (5 3) (5 4) In the instantaneous causal model, the previous states of {X t } as well as the previous and present states of {Y t }, are used to predict the present state of {X t }. Similarly, the previous states of {Y t }, as well as the previous and present states of {X t }, are used to predict the present state of {Y t }. Comparing Equation 5 1 and Equation 5 3 if the instantaneous causality is occurring ( ), then, knowledge of the present state of {Y t } will improve the prediction for the present state of {X t }, and will decrease. Similarly, comparing Equation 5 2 and Equation 5 4 if the instantaneous causality is occurring ( ), then, knowledge of the present state of {X t } will improve the prediction for the present state of {Y t }, and will decrease. PAGE 55 55 To obtain the frequency formulation of the simple causal model, we can rewrite Equation 5 1 and Equation 5 2 by using the time shift (delay) operator U, where UX t = X t 1 (5 5) (5 6) where COEF (U) (COEF = a, b, c, d) are power series in U with coefficient of U 0 being zero. Next, we use Cramer representation of the series to get (Granger, 1969): (5 7) (5 8) Note: Cramer representation is widely used in spectrum estimation for wide sense stationary process; similar to the role of Fourier transform in time frequency analysis in deterministic time series. Because k units of time shift (delay) (represented by U k ) in the time domain, is equivalent to being multiplied by e in the frequency domain, also noting that Cramer representation is a linear transformation, thus we get: (5 9) (5 10) (5 11) PAGE 56 56 (5 12) Integrate Equations 5 7 to 5 12 into Equations 5 5 to 5 6 we get (for every t): (5 13) (5 14) Because Equation ( 5 13 ) and ( 5 14 ) hold for every t, it further implies: (5 15) where [ ]. Therefore, (5 16) Then, we obtain the spectral and cross spectral matrix: (5 17) w here PAGE 57 57 When we study the causality between the two time series { X t } and { Y t }, if { Y t } is NOT causing { X t }, then ; similarly, if { X t } is NOT causing { Y t }, then Therefore, the cross spectrum can be decomposed into and { X t } causing { Y t Y t } causing { X t and can be considered as two arms of the feedback mechanism. In Equation 5 18 and Equation 5 19 two measures of strength of causality are defined X t } to { Y t Y t } to { X t (5 18) (5 19 ) w here PAGE 58 58 In Equation 5 20 and Equation 5 21 the causal phase diagrams are further defined, to study the phase lag against frequency cross spectrum: (5 20) (5 21) Similar to the frequency formulation of the simple causal model (refer to Equations 5 1 and 5 2 ), we consider the frequency formulation of the instantaneous causal model (refer to Equa tions 5 3 and 5 4 ). A spectral and cross spectral matrix can be obtained similar to what is given by Equation 5 17 as follow: (5 22) w here PAGE 59 59 Note: and are the coefficients of instant causality of { Y t } in Equation 5 3 and { X t } in Equation 5 4 respectively Summary: the spectral methods can be applied to analyzing causal relations between two time series { X t } and { Y t }. Furthermore, the spectral methods in the frequency domain are robust in their interpretation, compared with the equation models in the time domain. Th ree Variable (Time Series) Models Extend the two variable model (simple causal model) in Equation 5 5 and 5 6 to three variable model (Granger, 1969), we get (5 23) (5 24) (5 25) Here we keep using the time shift operator U ( ). Since there is no instantaneous causality, the coefficients of are always zero. Applying the frequency formulation (using Cramer representation) to the three variable model, the spectral and cross spectral matrix is obtained: PAGE 60 60 (5 26) w here Similar to the calculation in the two variable model, the power spectrum of , and as well as the cross spectrum: , and can be calculated in the three variable model. Further, the Partial cross spectrum b etween and given is defined as where are derived from Equation ( 5 26 ). When partial cross spectrum is considered, more useful results arise PAGE 61 61 (5 27) w here As shown in Equation 5 27 the partial cross spectrum between and can be decomposed into , and where the first component represents the relations of and through and the second and third components together represent generalizations of the causal cross spectra (between and ) which come from the two variable case. Multiple Variable Models Generalizing the two and three variable causal relation models to multiple variable models brings in depth insights into the causal relation and feedback problem (Geweke, 1982). In this sub s ection, causal relations between one group of time series PAGE 62 62 and another group of time series are studied. The discussion in this sub section is still based on wide sense stationary, purely nondeterministic time series. Let be a zero mean wide sense stationary, purely nondeterministic multiple time series. And suppose that, has been partitioned into two sub vectors: We try to find out the cau sal relations as well as the measurement of linear dependence and feedback between and By assuming the existence of the spectral density matrix at all frequencies we have a corresponding partition of the spectral density matrix (5 28) We assume there is a moving average representation of as follow : (5 29) In which, is serially uncorrelated; for any holds; and ( is the square root of the largest eigenvalue of ; is the square root of the determinant of .) The existence of moving average representation is equivalent to the existence of the spectral density matrix at almost all frequencies (Doob, 1953). Furthermore, a lower bound on the mean squared error of one step ahead minimum mean is given by using : (5 30) PAGE 63 63 If there exists a constant such that for almost all we have (Rozanov, 1967) : (5 31) Which means both and are positive semi definite, and that the spectral density matrix is bounded uniformly away from zero almost everywhere in which is w eaker than the condition that, is an autoregressive moving average (ARMA) process of finite orders with inv ertible AR and MA parts. Then, the relationship in Equation 5 27 can be inverted so that become a function of and (5 32) Predict from the P ast of I tself Considering has been partitioned into two sub vectors if the condition in Equation 5 30 holds, then, and each has an autoregressive (AR) representation, denoted by (5 33) (5 34) w here is one step ahead error when is predicted from its own past ; similarly is one step ahead error when is predicted from its own past PAGE 64 64 Recall Eq uations 4 15, 4 16 and 4 17 (the Principle of Orthogonality) in the LSE model in the previous chapter By generalizing the Principle of Orthogonality to infinite orders, to minimize the one step prediction error, is ort hogonal (uncorrelated) to (the projection of on its own past), and is orthogonal (uncorrelated) to (the projection of on its own past). Thus, is by itself serially uncorrelated, and is by itself serially uncorrelated. But and may be correlated with each other contemporaneously and at various leads and lags. Predict from the P ast of I tself and the P ast of O thers A further step to take is the projection of on its own p ast and the past of and the projection of on its own past and the past of as follow (5 35) (5 36) By generalizing the Principle of Orthogonality in the LSE model to infinite orders, is orthogonal (uncorrelated) to (the projection of on and ), and is orthogonal (uncorrelated) to (the projection of on and ), to minimize the one step prediction error. Thus, is by itself serially uncorrelated, an d is by itself PAGE 65 65 serially uncorrelated. It is worth mentioning that and ma y be correlated with each other contemporaneously and at various leads and lags. The linear prediction model in Eq uation 5 34 and 5 35 is another expression of the linear prediction model in Eq uation 5 31 The two expressions are the same intrinsically, considering Therefore, and we have the partition of as follow (5 37) Predict from the P ast of I tself and the P resent S tate and the P ast of O thers A third step to take is to pre multiply the following matrix on Eq uation 5 34 and Equation 5 35 Then, the first k equations can be considered as a new linear system of projected on and ; similarly, the next l equations can be considered a new linear system of projected on and as follow (5 38) PAGE 66 66 (5 39) w here (5 40) (5 41) Hence and similarly can prove Note that: each of and is uncorrelated (orthogonal) to both and as already discussed above. According to Equations 5 39 and 5 40 each of and is uncorrelated to both and According to Equation 5 39 is also uncorrelated to : Hence, is uncorrelated to all and According to Eq uation 5 35 is uncorrelated to Thus, is uncorrelated to and Therefore, Eq uation 5 37 is a linear (orthogonal) projection of on and Similarly, according to Equation 5 40 is also uncorrelated to : PAGE 67 67 Hence, is uncorrelated to all and According to Eq uation 5 34 is uncorrelated to Thus, is uncorrelated to and Therefore, Eq uation 5 38 is a linear (orthogonal) projection of on and Predict from the P ast of I tself and the T emporally C omplete I nformation of O thers A fina l step is to consider the linear projection of on and and the linear projection of on and Recall the partition of the spectral density matrix in Eq uation 5 27 Let (5 42) Note: a ssume eq uation 5 30 is true, then, the spectral density matrix exits and is bounded uniformly away from zero almost everywhere in ) Use inverse Fourier transform of : (5 43) Let (5 44) PAGE 68 68 Then, from the spectral representation of (Eq uation 5 27 ), it is be proven that is uncorrelated (orthogonal) with all Therefore, if we rewrite Eq uation 5 43 in then, Eq uation 5 43 can be considered as an expression of the linear projection of on From Eq uation 5 27 and 5 43 it can be further shown that: (5 45) Considering the assumption in Eq uation 5 30 : because consists of the first k rows and columns of and Eq uation 5 30 holds, we also have (5 46) Note: f or some constant and for all in above. Therefore, has an autoregressive representation, as follow (5 47) w here is uncorrelated (orthogonal) with Use from Eq uation 5 43 to substitute all in Eq uation 5 46 after re organizing, we get (5 48) PAGE 69 69 By Eq uation 5 46 is a linear function of and since any is uncorrelated (orthogonal) with all therefore, is uncorrelated with all According to Eq uation 5 43 is a linear function of and for any Also, is uncorrelated (orthogonal) with as already discussed above. Therefore, is uncorrelated with Hence, Eq uation 5 47 is the linear projection of on and Similarly (symmetrically), Eq uation 5 48 is the linear projection of on and (5 49) It is worth mentioning that: all the prediction models in this chapter are based on the assumption that the time series analyzed are a zero mean wide sense stationary, and purely nondeterministic. The analysis could be extended to some non stationary and purely nondeterministic time series, and similar autoregressive representations can be obtained (Zhang et al., 2011). PAGE 70 70 CHAPTER 6 EPILEPTIC BRAIN NETWORK MODELI NG AND BRAIN STATE TRANSITION STU DY Network M odeling of B rain s In the last decade, the study of complex networks has dramatically expanded across diverse s cientific fields, including the study and investigation of the brain. There is a tendency that science is concerned more about the structure, behavior, and evolution of the brain, which is an extremely complex system using modern network approach A brain includes a huge number of components : molecules, cells neurons functional regions and systems To understand the brain, not only the knowledge of elementary components is required, but also the knowledge of the ways in which those components interact wi th each other and the corresponding emergent properties. Generally, normal brains display organized and characteristic diverse patterns, due to the highly structured and selective coupling between the components, for example cells neurons or functional brain regions. The connectivity also has many different types, such as molecular interactions coupling between neurons, and the driving force between different brain zones in nervous systems or other biological systems B rain networks span multiple spatia l scales, from micro scale of individual cells, neurons, and synapses, to the macro scale of cognitive systems embodied organisms and brain functional regions It is impossible to fully understand brain functions as an integrated system unless the brain i s approached on all of those multiple scales and knowledge of the network interactions on and across multiple scales of organization is obtained It is worth mentioning that no ne of the scale s in the hierarchy is privileged over any other s Through multipl e scale network interactions, molecules and cells give PAGE 71 71 rise to behavior and cognition. T he q uantitative study of brain network s requires sophisticated mathematical and statistical techniques. The network modeling and analy tics of brain s ha ve opened new experimental and theoretical directions in many areas of neuro science It plays a significant role in neuro anatomy, neurodevelopment, electrophysiology, functional brain imaging, and the neural basis of cognition. Especially, t he analysis of ne twork connectivity topology and architecture illuminates many problems that concern diagnosis and treatment of brain disorders, the brain state transition study, and integrative brain function s : investigat ion of the brain network interactions among groups of a number of nerve cells between brain functional regions or connected in local brain circuits; quantitatively study of the brain network to shape brain anatomy for architectural principles such as localize the focal areas of an epilepsy ; Trace the ext ent and location of structural brain network damage for biomarkers to measure the nature and severity of insults like brain trauma, brain traumatic injury, and/or other cognitive dysfunctions Neocortical Epilepsy and Its Clinical Challenges Neocortical epilepsy is a type of seizure disorder that can be either partial (focal) or generalized, in which the seizures originate in the neocortex, part of the external surface layer of the brain. Seizures caused by neocortical epilepsy do not respond well to medication, and therefore epilepsy surgery is often one of the very few options that patients with neocortical epilepsy have. Unfortunately, surgery for neocortical epilepsy does not have the same success rates as that for other kinds of epilepsy. In patients with other forms of epilepsy, such as mesial temporal lobe epilepsy (MTLE), in which the lesion that causes the seizure is clearly defined, surgery can provide seizure freedom in 70 to 80 percent of cases. On the other hand, patients with neoco rtical PAGE 72 72 epilepsy see lower success rates. This could be because (1) it may be difficult to isolate the region that is causing the seizures with the currently available technologies ; (2) there may be multiple initiation points ; or (3) the initiation points m ay be related to important functional sites in the brain. Because of the complications, patients with neocortical epilepsy often need individually tailored surgeries. A recent study ( Lee et al. 2005) in 89 patients with nonlesional neocortical epilepsy sh owed that, based on the epileptogenic focus location, only 47% of the patients were seizure free (Engel Class I) after the surgery, and additional 7% experienced a significant reduction in seizure frequency (Engel Class II). Therefore, a successful develop ment of a more reliable method that can be used to identify suitable epilepsy surgery candidates with high confidence for patients with neocortical epilepsy would be a major contribution to advance the care for those patients. Existing Solutions in Clinica l Application s High resolution magnetic resonance imaging (MRI) is the most important diagnostic technique for epilepsy surgery in neocortical epilepsy. It offers a high predictive value for surgical outcomes (Cascino et al., 1993). However, MRI is ineffe ctive in 29% of patients with partial epilepsy (Semah et al., 1998), and many patients referred to epilepsy centers for surgery have normal MRI results. Previous studies report that surgical outcome is poor for patients with neocortical epilepsy with norma l MRI, but these conclusions were based on limited numbers of patients (Cascino, 2004). Intracranial EEG monitoring is indispensable for neocortical epilepsy with normal MRI. Depth and subdural EEG monitoring could provide accurate localization of the reg ion of seizure onset. However, sampling error can lead to false or missed localization. False and missed localization during intracranial EEG monitoring is PAGE 73 73 particular common with extratemporal seizure origin (Williamson et al., 1992a; 1992b). Nevertheless, several studies have demonstrated the usefulness of focus localization in neocortical epilepsy with quantitative intracranial EEG analyses. Worrell described the finding, in six patients with neocortical epilepsy, of high frequency (60 100 Hz) epileptifor m oscillations are highly localized in the seizure onset zone (Worrell et al. 2004) They further reported that focal neocortical seizures often occur during periods of increased high frequency activity in the seizure onset zone, indicating that high freq uency activities are involved in neocortical ictogenesis. Ochi also reported high frequency oscillations (HFOs) were observed in 78 out of 79 seizures studied ( Ochi et al. 2007 ) In addition, in four postoperatively seizure free patients, more electrodes recorded higher frequency HFOs inside the resection area than outside before and after clinical seizure onset, while in five patients with residual seizures, electrodes recorded more HFOs that were of higher or equal frequency outside the surgical area tha n inside examples to demonstrate that basic dynamic changes in focal epilepsy of neocortical origin may be useful in localizing the origin of seizures. The findings of these stu dies and many others support the hypothesis that the dynamics of intracranial EEGs are associated with the behavior of epileptogenic focus in neocortical epilepsy. Epileptic Brain Network Modeling for Neocortical Epilepsy Recent observations in humans with MTLE and in the animal models for this condition (Spencer and Spencer 1994; Bertram, 1997) suggest multifocal seizure onset (i.e. seizures that begin focally within different limbic structures with each seizure) as well as synchronized regional seizure onset (i.e. presumed simultaneous seizure initiation). These observatio ns, together with multifocal physiological and anatomical PAGE 74 74 changes in the animal models, have raised the possibility of a widely distributed neural network (e.g., specific cortical and subcortical networks) in the gensis and expression of partial and gener alized onset seizures and anatomically connected, bilaterally represented, set of cortical and subcortical brain structures and regions in which activity in any one part affects activity in all the others ( Spencer, 2002 ) One central observation on which this definition is based is that vulnerability to seizure activity in any one part of the network is influenced by activity everywhere else in the network, and that the network as a whole is responsible for the clinical and electrographic phenomenon that we associate with human seizures. She further distinguishes between the network and the "seizure onset zone". The seizure onset zone for neocortical seizures is defined by the ictal EEG. It is the area where the seizure discharge is first seen on ECoG (subdural electrodes). It may involve one to many electrode sites. In some patients, seizures (based on clinical manifestations) can start with EEG discharges beginning in different "onset zones" and spreading al ong different routes as the seizure progresses. At least three human epilepsy networks have been reported based on the extensive observations on epilepsy patients (Spencer, 2002). The first one is the medial temporal/limbic network which includes the hippo campi, the amygdalae, the entorhinal cortices, lateral temporal neocortices, and extratemporal components of the medial thalamus and the inferior frontal lobes. The other two networks are less commonly identified, even in their component parts: the medial occipital/lateral temporal network and the superior parietal/medial frontal network. These lines of evidence based on clinical observations, intracranial EEG, functional neuroimaging, anatomic observations, PAGE 75 75 and the response of seizures to specific invasive treatments, support the existence and importance of the networks in the genesis of human epilepsy, and further suggest that the network structures are essential to the development of the seizure and thus the existence and maintenance of the epileptic diso rder. Although much of the evidence of epileptic network has been most consistent in MTLE, extratemporal epilepsy also is associated with multifocal, sometimes extensive, reductions in metabolism demonstrated on interictal PET studies, but this has not bee n as reproducible or as well documented as that in the medial temporal / limbic syndrome (Juhasz et al., 2000). Interestingly, despite high sensitivity to extratemporal regions of epileptogenicity, which show hypometabolism on interictal studies, PET scans are not sensitive to initial areas of seizure propagation in extratemporal epilepsy, another demonstration of the nonequivalence of network and propagation areas (Juhasz et al., 2000). The existence of epileptic network can also be suppor ted by the responses to the invasive therapy. If human epilepsy is the expression of specific, abnormally active, intrinsically defined and connected cortical/subcortical/bilateral networks, then one could theoretically alter seizure expression by interven ing in any part of the specific network. Published reports document 60 90% excellent response, meaning cessation of seizures, after any kind or extent of temporal lobe resection in patients identified as having medial temporal lobe epilepsy. Operations in volving anterior temporal lobe, medial structures only, lateral structures only, or more or less extensive lateral temporal resection can cure this disorder. Procedures with no anatomic overlap are similarly successful This cannot be explained unless the multiple areas are all critical in the PAGE 76 76 production of the intractable seizures of this disorder. Then interruption of the network in any one of those areas would be (and apparently is) sufficie nt to alter the seizures. The similarly excellent response with cessation of seizures after temporal lobe resection in well selected patients who have bilateral independent medial temporal lobe origin of seizures is another example of the existence of a network, interference with which at any site alters the expression of the intractable seizures (Hirsch et al., 1991a; 1991b). Based on the above observations, it may be pertinent to consider study of other kinds of phenomena in individual patients, which may define the network in better terms than we have sought in the past because of our single minded attention to defining regions of so called seizure onset. For example, quantitative intracranial EEG analysis, background patterns, sleep effects on interictal and ictal activity, and other types of functional assessments may contribute considerably to our understanding of the role of networks in the expression of the epilepsies. Studying broad regions of brain structures related by the presence of such networks, using quantitative EEG analyses and sophisticated approaches, may detect alterations in the behavior of the network before manifestation clinically or on traditional EEG. EEG Analysis on Networks in Epileptic Brains Compared wi th scalp EEG s intracranial EEGs or ECoGs generally have higher signal to noise ratio and/or signal to artifacts ratio and are of higher recording quality. Intracranial EEGs exhibit the most significant observations that support the network hypotheses (Spe ncer, 2002). Because the entire network participates in the expression of the seizure activity and can be entrained from any of its various parts, initial electrical PAGE 77 77 occurre nce within the network. The initial area of apparent seizure involvement is not might even vary from seizure to seizure in a given patient. This locational variability m ay only one part of the network (Spencer and Spencer, 1994) This may be the main reason why, with the surgery candidates and the resection brain region determined by the current focus localization techniques, the chance of a patient with neocortical epilepsy becoming seizure free after an invasive brain surgery is still less than 50%. A reliable quantitative method based on intracranial EEGs for determining the spatial dis tribution of the focal region(s) and their possible ne twork may assist epileptologist s to more accurately identify surgery candidates that will result in better surgery outcomes for patients with neocortical epilepsy. Monto studied epileptic brain networks by quantifying the long range temporal correlation (LRTC) in subdural human EEG recordings ( Monto et al. 200 7) Their observations on the spread of abnormally large LRTCs suggested that the epileptic focus is associated with significant changes in networ k behavior even in the cortical areas immediately surrounding th e clinically determined focus. By applying synchronization and graph analysis to intracranial EEG recordings, Ponten investigated the hypothesis that functional neuronal networks during tempor al lobe seizures change in configuration before and during seizures ( Ponten et al. 2007 ) The study found that the functional brain networks change from a more random configuration during an interictal recording period to a more ordered configuration durin g seizures, especially during seizure spreading. The authors further suggested that the findings supported the PAGE 78 78 theory that a random network (during interictal periods) even had a stronger tendency to synchronize (Netoff et al., 2004), which could cause sei zures. More recently, Zaveri calculated a magnitude squared coherence (MSC) on background intracranial EEG as a measure of functional connectivity to investigate the network effects within and outside the seizure onset area. The analysis demonstrated an i nverse relationship between the connectivity strength and the distance from the seizure onset area (especially during the frequency band) ( Zaveri et al. 2009 ) Another similar work by Ortega et al. (2008) showed the presence of clusters of increased syn chronization in different locations on the lateral temporal cortex in patients with temporal lobe epilepsy. Although remain unproven, it is possible that increased connectivity ( as can be defined by coherence, synchronization clusters, or nonzero functiona l connectivity ) helps initiate seizures. If increased brain connectivity or alteration in network topology helps initiate seizures, then network nodes and pathways could serve as targets for resective or disconnective surgery implantable devices, and inve stigations of seizure anticipation. To u se EEGs to e xplore the role of the cortex (e.g. frontal cortex) in the epileptic brain network in a quantitative way is a more promising approach compared with t raditional ways to identify areas of the cortex involved in an epilepsy network by look ing for areas that demonstrate interictal spikes. It is generally accepted that the areas that show interictal spikes do not consistently correspond to the seizure onset zone. Nevertheless, it is very likely that areas with interictal spikes and the seizure onset zone are both within the same epileptic brain network. One important issue is to identify the most critical nodes or critical regions within the network that are responsible for the epileptogenic process PAGE 79 79 CHAPTER 7 BRAIN STATE T RANSITION STUDY FOR SEIZURE DETECTION AND PREDICTION, USING MULTIPLE FEATU RE EXTRACTION, PRINC IPAL COMPONENT ANALY SIS AND CLASSIFICATION Purpose and Scope In epileptic seizure predi ction, a question need s to be answered : whether a consistent trend and quantitative measure can be discovered, from inter ictal period to ictal period, based on the EEG recording from epileptic patients. Technically, it is an epileptic brain state transition study, which has the potential to b e applied to real time epileptic seizure detection and seizure prediction. Epileptic patients as subjects get involved for continuous EEG recording in several days. For each subject, several hours of inter ictal recording, and at least three seizures of ic tal recording and the corresponding pre ictal recording ( mostly start from one hour before the onset of seizure), will be included. The whole recording will be divided into few seconds epochs for feature extraction, principal component analysis, and SVM le arning and testing. A quantitative score will be applied to each epoch for classification, depending on (1) the type of applications (e.g. seizure prediction, seizure detection) ; and (2) specifications (e.g. seizure prediction horizon). Procedure Design an d Performance Evaluation Criteria In order to serve as a real time seizure detection and seizure prediction (warning) system or as a part of an epileptic brain state transition study, the experiment on the algorithm will need to demonstrate its capability of classifying the inter ictal, pre ictal, and ictal states of the epileptic brain and differentiating through continuous EEG recording. The usefulness of the algorithm is dependent on the quality of EEG recordings (signal to noise ratio, type and strength of artifacts) and the specification of PAGE 80 80 application. The task of investigating and demonstrating (if useful) this capability with patients of epilepsy will be undertaken in three steps : Step 1, Epileptic Brain State Transition Study Using Multiple Features Extraction on Different Lengths and Overlap of Data Epochs Step 2, Preprocess and Feature Space Dimensionality Reduction Step 3 Classification Algorithm Using Different Types of Support Vector Machines: Training (Learning), Testing and Performance Evaluation Consistent separation of the inter ictal, pre ictal, and ictal states of the epileptic brain is the key of usefulness of this algorithm. For this reason, a considerable sized database (10 15 subjects) need be setup for testing the algorithm. Th e dataset in this study is mainly from BIH MIT database of scalp EEG recordings, most of which are with 23 channels. T he algorithm is patient specific and need to be trained by a segment of continuous EEG recording including inter ictal, pre ictal, and icta l periods. Some definitions are put as in the following: S ENSITIVITY (S E ) Percentage of events correctly identified by the algorithm in relation to the events actually marked by the subject. F ALSE A LARM R ATE (FAR) Frequency of events incorrectly detected by the algorithm per hour. P OSITIVE P REDICTIVITY V ALUE (PPV) T he proportion of true positives on the total of true positives and false positives, which reflects the probability of a positive test. F EATURE E XTRACTION (FE) A special form of dimensionality reduction in classification and pattern recognition. P RINCIPAL C OMPONENT A NALYSIS (PCA) A mathematical procedure that transforms a number of possibly correlated variables into a smaller number of uncorrelated variables. So me criteria and design input s are applied to evaluate the performance of the algorithm such as SE FAR and PPV, as in the following: PAGE 81 81 Consistency of EEG recording: each subject needs to be in a consistent state (e.g. sleeping / awake) during inter ictal, p re ictal and ictal periods for each recording. Ambulatory Artifacts: EEG segments that include very strong artifacts (the amplitude is one order higher than the background EEG) need to be annotat ed in one or more channels. Prediction horizon: better between 3 min and 50 min Sensitivity (Se) : milestone sets Se at >75% False Alarm Rate (FAR): milestone sets FAR better at < 1 / hour also depending on the seizure occurrence frequency and the prediction horizon. S ubjects each in clud e at least three seizures in EEG recording through inter ictal, pre ictal and i ctal periods. For each subject, a randomly selected one third of seizures are included in the training dataset, and the other two thirds of seizures are included in the tes ting dataset. For each seizure, the pre ictal period better starts from one hour before the onset or earlier. The sensitivity and positive predictive value also depend on the level of artifacts and the consistency of state (e.g. sleep / awake; sit / inten sive movement) of the subjects being recorded. There are a great number of types of epileptic seizures. Three major types of epileptic seizures need be included : Tonic chronic epileptic seizures Absence epileptic seizures and Simple or complex partial se izures T he performance of the algorithm needs to be evaluated. Improvement can be made either by parameter tuning / optimization, or by introducing innovative / supplementary methodology. Further, the efficacy of the algorithm can be re assessed after the tuning phase. PAGE 82 82 Brain State Classification of Different Epileptic Stages Conjecture The epileptic seizures occur aft er continuous and gradual brain stat e changes, which are generally difficult to be noted by the behavior of the patients, but can b e recognized and identified by EEG recordings. Method s Scalp EEGs are recorded during daily life ( in ambulatory environment) of epileptic patients. Seven measures (features) (1.STLmax, 2.PM R S, 3.FMX, 4.STM, 5.STX, 6.TEG, 7.STD) are calculated on the epoch samples in each of the EEG channels to evaluate brain stat e in different epileptic stages (inter ictal, pre ictal, and ictal) and help real time detect and/or predict epileptic seizures. For each of the seven measures (features), principal component analys is (PCA) is applied to extract the largest component from feature values of all the channels and to reduce the number of dimensions from the number of total recording channels to one (generally the largest component has more than 90% energy of all the eigenvalues). Then, using only the largest component of all the channels, principal component analysis (PCA) is applied a second time to extract the three largest components (as t hree indices) from the seven measures (features), so as to further reduce the number of dimensions from seven to three, and to visualize the results. Finally, kernel support vector machine (Kernel SVM) is applied to the three (or two) indices to design a n on linear classifier to separate / distinguish pre ictal and inter ictal periods so as to provide indications of different kinds of forthcoming epileptic seizures. PAGE 83 83 Results and Conclusion The algorithm was preliminarily tested on epileptic patients. Each of them had ictal recording ictal recording and recording in the ictal period with the patient awake in both stages. 8% of the total recordings were randomly selected for SVM training and the remaining 92% recordings were used as testing samples in SVM classification. In one of the patients, we had very excellent results in the separation of pre ictal and inter ictal stages, whereas in the other two patients, we had good and acceptable results in the classification of diffe rent epileptic stages. By twice applying PCA respectively on all the recording channels and on the seven measures (features) of the scalp EEG s some effective quantitative indices can be found to help recognize and identify the gradual brain status changes in different epileptic stages. Furthermore, those indices can be used to help design a classifier of pre ictal and inter ictal stages to support the accurate prediction of epileptic seizures. In Figure 7 1 to Figure 7 6, the plot of the two largest PCA co mponents and the three largest PCA components of Patient 48, 51, 52 are shown each from the interictal period to the ictal period, respectively. In Figure 7 7 to Figure 7 12, the plot of the classification effects based on kernel methods are shown, in whic h Figure 7 7 to Figure 7 9 use polynomial kernels and Figure 7 10 to Figure 7 12 use quadratic kernels. PAGE 84 84 Figure 7 1 Plot of the two largest components of patient 48 from interictal to ictal period: Interictal in red color; preictal period are divided into consecutive 5min segments and are plotted in different colors and shapes. Figure 7 2. Plot of the three largest components of patient 48 from interictal to ictal period: Interictal in red color; preictal period are divided into consecutive 5min segments and are plotted in different colors and shapes. PAGE 85 85 Figure 7 3. Plot of the two largest components of patient 51 from interictal to ictal period: Interictal in red color; preictal period are divided into consecutive 5min segments a nd are plotted in different colors and shapes. Figure 7 4. Plot of the three largest components of patient 51 from interictal to ictal period: Interictal in red color; preictal period are divided into consecutive 5min segments and are plotted in differe nt colors and shapes. PAGE 86 86 Figure 7 5. Plot of the t wo largest components of patient 5 2 from interictal to ictal period: Interictal in red color; preictal period are divided into consecutive 5min segments and are plotted in different colors and shapes. Fig ure 7 6. Plot of the three largest components of patient 52 from interictal to ictal period: Interictal in red color; preictal period are divided into consecutive 5min segments and are plotted in different colors and shapes. PAGE 87 87 Figure 7 7. Plot of kernel classification effect between interictal period and preictal period using polynomial kernel, of P atient 48 : i nterictal in magenta color ; preictal in cyan color Figure 7 8. Plot of kernel classification effect between interictal period and preictal per iod, using polynomial kernel, of Patient 51: interictal in magenta color; preictal in cyan color. PAGE 88 88 Figure 7 9. Plot of kernel classification effect between interictal period and preictal period, using polynomial kernel, of Patient 52: interictal in magen ta color; preictal in cyan color. Figure 7 10. Plot of kernel classification effect between interictal period and preictal period, using quadratic kernel, of Patient 48: interictal in magenta color; preictal in cyan color. PAGE 89 89 Figure 7 11. Plot of kernel classification effect between interictal period and preictal period, using quadratic kernel, of Patient 51: interictal in magenta color; preictal in cyan color. Figure 7 12. Plot of kernel classification effect between interictal period and preictal period, using quadratic kernel, of Patient 52: interictal in magenta color; preictal in cyan color. PAGE 90 90 Events D etection U sing S ub B ands D ivision, M ultiple F eature E xtraction, and C lassification in H igh D imensional S pace on A mbulatory EEG In an experimental pr otocol for supplemental materials for NIH fast track application entitled: seizure mitigation through continuous EEG with responsive Vagus nerve stimulation T ypes of Events being detected are : ( 1 ) Eye Close ; ( 2 ) Eye Flutter ; and ( 3 ) Eat and Chewing S ubjects being recorded are : (1) Sub ject 1 with 282 min EEG record for test ; and (2) Sub ject 2 with 243 min EEG record for test The electrode placement of EEG r ecording is: a subset of International 10 20 electrodes, Fp1, Fz, Pz, Fp2, C3, O1, F3, C4, O2, F 4, Cz (Use P4 as reference). And the sample rate is: 240 Hz Event Type : Eye Close Set threshold of classification score at 1.00 For Sub ject 1 the Probability of Detection (PD) is 100% (detect 5 out of 5) and the False Alarm Rate ( FA R ) is 0 per hour For Sub ject 2 the Probability of Detection (PD) is 100% (detect 5 out of 5) and the False Alarms Rate (FA R ) is 0 per hour Figure 7 1 3 Eye Close Detection Results of Subject 1. (Blue: Classification Score of each 2 sec epoch (with 1 sec overlap); Red : reference annotation of onset of eye close; Yellow: reference annotation of offset of eye close) PAGE 91 91 Figure 7 1 4 Eye Close Detection Results of Subject 2 (Blue: Classification Score of each 2 sec epoch (with 1 sec overlap); Red: reference annotation of onset of eye close; Yellow: reference annotation of offset of eye close) Event Type: Eye Flutter Set threshold of classification score at 1.25 For Sub ject 1 the Probability of Detection (PD) is 75% (detect 3 out of 4) and the False Alarm Rate ( FA R ) is 0.638 per hour (detect 3 false positives in 282 min record for test) For Sub ject 2 the Probability of Detection (PD) is 100% (detect 4 out of 4) and the False Alarms Rate ( FA R ) is 0.741 per hour (detect 3 false positives in 243 min record for test) Figure 7 15 Eye Flutter Detection Results of Subject 1. (Blue: Classification Score of each 2 sec epoch (with 1 sec overlap); Red: reference annotation of onset of eye flutter; Yellow: reference annotation of offset of eye flutter) PAGE 92 92 Figure 7 16 Eye Fl utter Detection Results of Subject 2. (Blue: Classification Score of each 2 sec epoch (with 1 sec overlap); Red: reference annotation of onset of eye flutter; Yellow: reference annotation of offset of eye flutter) Event Type: Eat & Chewing Set threshold of classification score at 1.08 For Sub ject 1 the Probability of Detection (PD) is 100% (detect 3 out of 3) and the False Alarm Rate ( FA R ) is 0.213 per hour (detect 1 false positive in 282 min record for test) For Sub ject 2 the Probability of Detection (PD) is 57.14% (detect 4 out of 7) and the False Alarm Rate ( FA R ) is 0.494 per hour (detect 2 false positives in 243 min record for test) Figure 7 17 Eat & Chewing Detection Results of Subject 1. (Blue: Classification Score of each 2 sec epoch (with 1 sec overlap); Red: reference annotation of onset of eat & chewing; Yellow: reference annotation of offset of eat & chewing) PAGE 93 93 Figure 7 18 Eat & Chewing Detection Results of Subject 2. (Blue: Classification Score of each 2 sec epoch (with 1 sec overlap ); Red: reference annotation of onset of eat & chewing; Yellow: reference annotation of offset of eat & chewing) Patient S pecific D etection of Intractable Seizures U sing S ub B ands D ivision, F eature E xtraction, and C lassification in H igh D imensional S pace on S calp EEG R ecorded from P ediatric S ubjects Database and Background Th e CHB MIT database, collected at the Children's Hospital Boston, consists of EEG recordings from pediatric subjects with intractable seizures. Subjects were monitored for up to several days following withdrawal of anti seizure medication in order to characterize their seizures and assess their candidacy for surgical intervention. ( Goldberger et al., 2000) Recordings, grouped into 23 cases, were collected from 22 subjects (5 males, ages 3 22; and 17 females, ages 1.5 19). (Case chb21 was obtained 1.5 years after case chb01, from the same female subject.) Each case contains between 9 and 42 continuous .edf (European Data Format) files from a single subject. In most cases, the .edf files c ontain exactly one hour of digitized EEG signals, although those belonging to case chb10 are two hours long, and PAGE 94 94 those belonging to cases chb04 chb06 chb07 chb09 and chb23 are four hours long; occasionally, files in which seizures are recorded are shorter. All signals were sampled at 256 samples per second with 16 bit resolution. Most files contain 23 EEG signals (24 or 26 in a few cases). The International 10 20 system of EEG electrode positions and nomenclature was used for these recordings. In a few records, other signals are also recorded, such as an ECG signal in the last 36 files belonging to case chb0 4 and a vagal nerve stimulus (VNS) signal in the last 18 files belonging to case chb09 In all, these records include 198 seizures (182 in the original set of 23 cases) Results of C ase 1 as A n E xample Take case 1 as an example There are seven *.edf (Euro pean Data Format) file has one seizure A ll the other 39 files are seizure free. Three randomly selected seizure events are included in the training dataset on this patient Apply SVM supervised learning, and run the SVM to classify and detect seizures on the remaining part of recording which includes all the other seizures As shown in Figure 7 19 to Figure 7 22, 100% sensitivity and 0 / hr false alarm rate if set threshold of c lassification score at 0.5 for this patient. The patient specific seizure detection results of other cases are also very good with high sensitivity approaching to 100% and very low false alarm rate. It indicates that, for the design of patient specific aut omatic seizure detection algorithm high performance is expectable though a machine training session needs to be designed for real world practical use (Shoeb 2010) PAGE 95 95 Figure 7 19 Detection Result of Intractable Epileptic Seizure in chb01_04.edf (Blue: Classification Score of each 2 sec epoch (with 1 sec overlap); Red: reference annotation of seizure onset; Yellow: reference annotation of seizure offset ) Figure 7 20 Detection Result of Intractable Epileptic Seizure in chb01_16.edf. (Blue: Classification Score of each 2 sec epoch (with 1 sec overlap); Red: reference annotation of seizure onset; Yellow: reference annotation of seizure offset) PAGE 96 96 Figure 7 21 Detection Result of Intractable Epileptic Seizure in chb01_18.edf. (Blue: Classification Score of each 2 sec epoch (with 1 sec overlap); Red: reference annotation of seizure onset; Yellow: reference annotation of seizure offset) Figure 7 22 Detection Result of Intractable Epileptic Seizure in chb01_26.edf. (Blue: Classificat ion Score of each 2 sec epoch (with 1 sec overlap); Red: reference annotation of seizure onset; Yellow: reference annotation of seizure offset) PAGE 97 97 CHAPTER 8 R WAVE/BEAT DETECTION FROM SURFACE ECG WIT H MUSCLE ARTIFACTS FOR CARDIAC BASED SEIZURE DETECT I ON Background and Purpose s For some types of partial epilep tic seizures such as temporal lobe epilepsy the phenomenon of an accompaniment of ictal tachycardia exists in a high percentage of seizure events partly because the epileptic discharges may influence portions of the central automatic network during the seizure occurrence Th us t he cardiac alteration s at the transition from the preictal state to the ictal state can be u sed to assist the automatic seizure detection, by identifying the changes i n heart rate and cardiac arrhythmias. Since the long term ambulatory ECG monitoring is technically feasible and common compared with the long term ambulatory EEG, it would be beneficial if a reliable R wave (heart beat) detection algorithm can be developed to monitor the heart rate and cardiac arrhythmias accurately for further partial epileptic seizure detection In the real application in the neuromodulation industry an accurate instantaneous heart rate may need to be detected based on ECG for a mbulatory use Furthermore, the ECG electrode placement may be on surface of the skin and may not us e the typically used lead 2 which has a very favorable signal to noise ratio. As is known, the h eart beat detection algorithm using lead 2 ECG ha s been well developed in the last few decades. However, in applications whe re the recording electrode placement is not using the lead 2 and is on surface of the skin the design of a reliable R wave (heart beat) detection algorithm is much more challeng ing A maj or challenge is that, the skelet al muscle artifacts are so strong sometimes (e.g. in body movements, exercises, or falls) that its amplitude could be one order higher than the desired ECG signal Furthermore, the muscle artifacts generally are not PAGE 98 98 separable from ECG signal in frequency domain or other commonly used transform domains Another challenge is that: when a surface ECG recording system is not using lead 2 electrode placement, instead us ing a patch on the chest for example the recorded ECG signal from heart would be much weaker than that of using lead 2. Hence the signal to noise ratio is low, even without very strong muscle artifacts. One of the projects I have been enrolled in a neuromodulation company is design ing a product aiming to lo g seizure events automatically detected by an algorithm built in the product. The seizure detection algorithm therein is based on the changes of heart rate and heart rate variability, both necessitates heart beat detection from acquired ECG with adequate a ccuracy. A new beat detection algorithm i s designed developed and implemented to overcome those challenges and has been verified and tested that it has satisfactory results. Compared with some open source algorithms and a nother beat detection algorithm (referred to as IMEC) designed and developed in a large group in Belgium, Europe, the proposed beat detection algorithm is superior on standard testing performance and is more robust in extreme cases. The new algorithm meets the performance requirements de fined by various standards established by ANSI/AAMI. I am the major contributor (more than 80% of ideas and work) to the design and development of the proposed algorithm. Design Features of the Proposed Beat D etection Algorithm The propo sed beat detection algorithm has taken the role of the originally used IMEC algorithm for ECG based beat detection and heart rate monitoring. It is a n improved and successful algorithm of R wave detection under strong noises, including PAGE 99 99 muscle artifacts, impedance noise, power line noise, and baseline wandering noise It has many features embedded in the architecture of design and in different modules Those features have help ed improve the performance of the algorithm such as detection accuracy, sensitivity, false positive rate and positive predictivity value (PPV) S ome of th e features are described selectively in this section The flow chart of the algorithm is shown as in Figure 8 1. Multiple features are extracted from the filtered a nd preprocessed ECG input, such as the shape and waveform of R wave, the morphology of QRS complex, the first and second order differential values, the innovative sharpness feature etc. Among those, features are extracted from the time domain, the frequency domain and other transformed domains. All th e features interact with each other, and are infused for a decision which would generate a pool of heart beat candidate s in a given time window Meanwhile, an individualized likelihood score is generat ed for each of the heart beat candidate s In a given time window, t he heart beat candidate s are then put together for further screening and selection in order to output the true heart beats. In the selection process, individualized neighborhood widths are partly determined by the likelihood score of each heart beat candidate as in Figure 8 2. In addition, the individualized neighborhood widths in a given time window are multiplied by a scaling factor based on previous heart rate sta tistics, as in Figure 8 3. For each heart beat candidat e the likelihood score is a measure of the likelihood th at it is a true heart beat ; meanwhile, the likelihood score is inversely related with the value of the individualized neighborhood width. PAGE 100 100 The previous heart rate statistics can influence the selection of true heart beats from heart beat candidate s pool, dynamically. Generally, the higher the previous heart rate, the better chance the heart beat candidate s in the current time window are to be selected. The method of calculating the previous heart rate from previous heart beat data needs to be designed according to the actual application in real world though not unique. Figure 8 4 is an example. In feature extraction process, new features are explored and designed for reducing the false positive rate of the beat detector. An example is the sharpness feature, which can assist distinguishing some types of false heart beats from the true heart beats as in Figure 8 5. Performance and Standard Testing Results of the Proposed Algorithm Th e Test ing is designed to characterize the performance of the new and the proposed beat detection algorithm in processin g not only Human ECGs from MIT BIH D atabase s of Arrhythmia with or without stressed noise, wh ich are widely used as testing ECGs, but also the simulated ECG with major ECG features, such as heart rate, R wave amplitude, QRS duration, and T wave amplitude, vary within their normal dynamic ranges. Since ECG recordings with limited period of time and from limited number of patients can not cover the broad dynamic ranges of major ECG parameters, such as heart rate, R wave amplitude, QRS duration, and T wave amplitude, FDA requires that a beat detection algorithm be tested not only with human ECGs but a lso with simulated ECG with various ECG parameters. PAGE 101 101 Figure 8 1 Flow chart of multiple features extraction in time domain, frequency domain and other transformed domains and feature infusion PAGE 102 102 Figure 8 2 P lot of the selection process that individualized neighborhood widths are partly determined by the likelihood score of each candidate heartbeat, in a given ti me window. A ) B ) C ) D ) are four cases In each case, f ive heart beat candidates are competing with each other on axis of time within a given time window. Different colors represent different c andidates. Each candidate has a likelihood score (the height) and a corresponding individualized neighborhood width (the width). For a candidate the width of the neighborhood is inversely related to likelihood score. If a candidate has the highest likelihood score within its own neighb orhood, then it can be selected (with a yellow ring) PAGE 103 103 Figure 8 3 P lot of dynamical scaling of the neighborhood widths according to prev ious heart rate statistics, in a given time window. A), B), C) are three cases with f ive heart beat candidate s competing on the axis of time. Each candidate has a likelihood score (the height) and an individualized n eighborhood width (the width). T he neighborhood width is partly determined by the likelihood score, and is later on multiplied by a positive scal a r from the previous heart rate data. The value of th e scalar is the same for all candidates in a given time window, and is inversely related to the previous heart rate. If a candidate has the highest likelihood score within its own neighborhood, then it can be selected. Some documents applicable to the design, development and testing of the proposed algorithm are as in the following: ANSI/AAMI EC13, 2002, Cardiac Monitors, Heart Rate. ANSI/AAMI EC38, 2007, Ambulatory electrocardiographs. ANSI/AAMI EC57, 1998 Testing and reporting performance results of cardiac rhythm and ST segment measurement algorithms. EN 60601 2 25, 1996, Medical electrical equipment. Particular requirements for safety. Specification for electrocardiographs. PAGE 104 104 EN 60601 2 27, 1995, Medica l electrical equipment. Particular requirements for safety. Specification for electrocardiographic monitoring equipment. EN 60601 2 49, 2001, Medical Electrical Equipment Part 2 49: Particular Requirements for the Safety of Multifunction Patient Monito ring Equipment. Figure 8 4 Plot of calculating previous heart rate statistics using exponentially decreasing weighting (an example) for a given time window. Figure 8 5 Plot of a new f eature built into the proposed beat detection algorithm that can help distinguish the true heart beats from a class of false ones. C yan colored circles are true heart beats, according to reference annotations; while the red colored circles are false heart beats that a re ruled out by the sharpness feature design in the algorithm. PAGE 105 105 Test : Beat Detection Accuracy with Human ECG According to EC57 3.1.1 General description of available databases, the MIT BIH Arrhythmia ECG database with beat annotation is used to test the acc uracy of beat detector. Performance of beat detector with match window of 5 10 20 50 100 150 millisecond were calculated. The results of the test were presented in Table 8 1 Table 8 1. Results of t he p roposed b eat d etection a ccuracy with h uman ECG from MIT BIH Arrhythmia Database Test Condition Total Beats Total True Positive Total False Negative Total False Positive Sensitivity (%) Positive Predictivity (%) As for Fiducial Point 109494 13057 96437 96110 11.92 11.96 As for Fiducial Point 5 ms 94657 14837 14510 86.45 86.71 As for Fiducial Point 10 ms 104176 5318 4991 95.14 95.43 As for Fiducial Point 20 ms 106378 3116 2789 97.15 97.45 As for Fiducial Point 30 ms 107601 1893 1566 98.27 98.57 As for Fiducial Point 50 ms 108498 996 669 99.09 99.39 As for Fiducial Point 150 ms 109040 454 127 99.59 99.88 Considering the inaccuracy of the reference annotations (the majority is within two samples from the R wave or 5.6 ms) of MIT BIH arrhythmia database, 20 ms match window is selected. The testing results showed that the proposed beat detection algorithm has adequate performance in detecting the beat of the ECG database, with sensitivity of 9 7 15 % and positive predictivity of 9 7 45 % with 20 ms match window It has significant improvement in beat detection accuracy compared to IMEC algorithm. In an applicat ion where the displayed heart rate is the average heart rate of the most recent beats (but not real instantaneous heart rate), the proposed beat detection algorithm can work very well in most cases In the application, where heart rate PAGE 106 106 variability (HRV) ma y need be used, it is better to have both the sensitivity and positive predictivity above 98% with 20 ms match window on MIT BIH arrhythmia database Test : Beat Detection Accuracy with Noise Stressed Human ECG According to EC57 3.1.1 General description of available databases, the MIT BIH noise stressed Arrhythmia ECG data set with beat annotation is used to test the accuracy of the beat detection. Performance of beat detection by comparing the detection with annotated fiducial points of 20 millisecond as match window is calculated. Table 8 2. Results of the proposed beat detection algorithm (20 ms match window) Recording Total TP FN FP Sensitivity (%) PPV value (%) 118e24 2278 2276 2 2 99.91 99.91 118e18 2273 5 10 99.78 99.56 118e12 2255 23 119 98.99 94.99 118e06 2100 178 501 92.19 80.74 118e00 1745 533 927 76.60 65.31 118e_6 1490 788 1127 65.41 56.94 119e24 1987 1958 29 23 98.54 98.84 119e18 1952 35 31 98.24 98.44 119e12 1933 54 124 97.28 93.97 119e06 1881 106 413 94.67 82.00 119e00 1629 358 867 81.98 65.26 119e_6 1336 651 1120 67.24 54.40 Total Beats: 25590 Total TP: 22828 Total FN: 2762 Total FP: 5264 TOTAL Sensitivity (%) PPV value (%) Average 25590 89.24 82.52 Gross 25590 89.21 81.26 The testing results showed that the proposed beat detection algorithm has significant improvement in beat detection performance in noise stressed ECG data set, compared to IMEC algorithm. In the application, the artifacts (including noises) may be stronger a nd SNR may be lower than 6 db. PAGE 107 1 07 Table 8 3 Results of the proposed beat detection algorithm (5 ms match window) Recording Total TP FN FP Sensitivity (%) PPV value (%) 118e24 2278 838 1440 1440 36.79 36.79 118e18 836 1442 1447 36.70 36.62 118e12 833 1445 1541 36.57 35.09 118e06 785 1493 1816 34.46 30.18 118e00 652 1626 2020 28.62 24.40 118e_6 561 1717 2056 24.63 21.44 119e24 1987 1778 209 203 89.48 89.75 119e18 1775 212 208 89.33 89.51 119e12 1768 219 289 88.98 85.95 119e06 1709 278 585 86.01 74.50 119e00 1461 526 1035 73.53 58.53 119e_6 1169 818 1278 58.83 47.60 Total Beats: 25590 Total TP: 14165 Total FN: 11425 Total FP: 13927 TOTAL Sensitivity (%) PPV value (%) Average 25590 56.99 52.53 Gross 25590 55.35 50.42 Table 8 4 Results of IMEC beat detection algorithm (5 ms match window) Recording Total TP FN FP Sensitivity (%) PPV value (%) 118e24 2278 113 2165 2762 4.96 3.93 118e18 124 2154 2635 5.44 4.49 118e12 123 2155 2400 5.4 4.88 118e06 117 2161 2193 5.14 5.06 118e00 118 2160 2159 5.18 5.18 118e_6 103 2175 2717 4.52 3.65 119e24 1987 779 1208 1880 39.2 29.3 119e18 888 1099 1552 44.69 36.39 119e12 889 1098 1270 44.74 41.18 119e06 867 1120 1131 43.63 43.39 119e00 867 1120 1119 43.63 43.66 119e_6 644 1343 2011 32.41 24.26 Total Beats: 25590 Total TP: 5632 Total FN: 19958 Total FP: 23829 TOTAL Sensitivity (%) PPV value (%) Average 25590 23.25 20.45 Gross 25590 22.01 19.12 PAGE 108 108 Furthermore, in the application, the QRS complex in ECG recording may have biphasic morphology, instead of the QRS morphology of normal ECGs. Therefore, for the proposed algorithm, the performance in detecting heart beats from the noise stressed ECG recordings need be further improved. The noise stressed testing results of the proposed beat detection algorithm on the two ECG recordings (#118, #119) are ACCEPTABLE basically. Test : Beat Detection on Rang e and Accuracy of Heart Rate with Simulated ECG According to EC13 S ection 4.2.7 Range and accuracy of heart rate meter, minimum allowable heart rate meter range shall be 30 to 200 BPM with an allowable readout error of no greater than 10 percent of the input rate or 5 BPM whichever is greater. Cardiac monitors labeled for use with neonatal/pediatric patients shall have an extended heart rate range of at least 250 BPM In addition input ECG signals at rates less than the disclosed lower limit of the rate meter range shall not cause the meter to indicate a rate greater than this lower limit. Input signals at rates above the disclosed upper limit of the rate meter range, up to 300 B PM (350 BPM for monitors labeled for use with neonatal/pediatric patients), shall not cause the meter to indicate a rate lower than this upper limit. To test the performance of the beat detection algorithm, a simulated AAMI NSR signal of 1 millivolt (mV) a s the amplitude of R wave and 25% of the R amplitude as the amplitude of T wave is used to test the range and accuracy of heart rate. For the rate of the simulated signal from 30BPM to 180BPM, the simulated signal is of 100 millisecond (ms) as QRS d uration, 180ms as the T wave duration (d T ), and 350ms as the QT interval (d QT ); for the rate from 200BPM to 30 0BPM, the corresponding parameters are 40ms, 60ms and 150ms, respectively. QRS amplitude is PAGE 109 109 defined as a r + a s, as described in EC13. The sample r ate of the simulated ECG is 250 Hz. The testing results showed that the proposed algorithm can detect heart rate correctly when the heart rate ranges at 30, 40, 50, 60, 80, 120, 160, 200, 240, 260, 280, 290 and 300 BPM of the simulated ECG. The error of de tected heart rate by the proposed is very small (<0.1%). The proposed algorithm passed the test On the other hand, except when the heart rate of simulated ECG ranges at 50, 60, 80, 100, 120, 160 and 200 BPM, the IMEC beat detection algorithm exhibits much bigger errors than the application allows in both ends of low and high heart rates. T he IMEC algorithm failed the test Table 8 5. Results of the t est of the p roposed b eat d etection algorithm on r ange and a ccuracy of h eart r ate with s imulated ECG Simulated Rate ( BPM ) 30 40 50 60 80 120 160 200 240 260 280 290 300 Measured Rate ( BPM ) 30 40 50 60 80 120 160 200 240 260 280 290 300 Error (Simulated Measured) ( BPM ) 0 0 0 0 0 0 0 0 0 0 0 0 0 Table 8 6. Results of the t est of IMEC b eat d etection algorithm on r ange and a ccuracy of h eart r ate with s imulated ECG Simulated Rate ( BPM ) 30 40 50 60 80 120 160 200 240 260 280 290 300 Measured Rate ( BPM ) 69 51 50 60 80 120 160 200 120 102 125 131 127 Error (Simulated Measured) ( BPM ) 39 11 0 0 0 0 0 0 120 158 155 159 173 PAGE 110 110 Test : Beat Detection on Range of R Wave Amplitude and QRS Duration with Simulated ECG According to EC13 S ection 4.2.6.1, range of QRS wave amplitude and duration, for a continuous train of simulated ECG pulses, the device shall meet the heart rate range and accuracy requirements of 4.2.7. The minimum range of QRS amplitude (a r + a s ) is 0.5 to 5.0 mV, and the duration of the QRS wave is between 70 and 120 ms (40 and 120 m s for neonatal/pediatric monitors). For monitors set for adult pa tients, the heart rate meter shall not respond to signals having a QRS amplitude of 0.15 mV or less, or a duration of 10 ms or less with an amplitude of 1 mV. Response to either or both of these types of signals is permitted in monitors set for neonatal/pe diatric Table 8 7. Results of the t est of the p roposed b eat d etection algorithm on r ange of R wave a mplitude and r ange of QRS d uration with s imulated ECG QRS Duration (ms) 20 30 40 50 60 80 100 120 140 R wave Amplitude (mV) 0.15 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 0.17 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 0.20 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 0.25 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 0.30 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 0.35 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 0.40 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 0.45 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 0.50 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 0.60 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 0.70 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 0.80 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 0.90 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 1.0 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 2.0 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 3.0 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 4.0 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 5.0 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 80.001 To test the range of R wave amplitude and QRS duration, the simulated AAMI NSR signal of minimum range of R wave amplitude from 0.15 to 5.0 mV, and of the PAGE 111 111 duration of the QRS complex between 10 and 140 ms is used. The heart rate meter should not respond to signals having QRS amplitud e of 0.15 mV or less, or duration of 10 ms or less with amplitude of 1mV. The amplitude of the T wave is 25% of the QRS amplitude, and the rate of the AAMI NSR signal is 80 BPM. The sample rate of the simulated ECG is 250 Hz. The testing results showed tha t both the proposed and IMEC beat detection algorithm s presented consistent accurate measurement of the heart rate of the simulated ECG. Both the proposed and IMEC algorithm s passed the test Table 8 8. Results of the t est of IMEC b eat d etection algorithm on r ange of R wave a mplitude and r ange of QRS d uration with s imulated ECG QRS Duration (ms) 20 30 40 50 60 80 100 120 140 R wave Amplitude (mV) 0.15 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 0.17 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 0.20 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 0.25 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 0.30 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 0.35 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 0.40 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 0.45 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 0.50 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 0.60 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 0.70 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 0.80 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 0.90 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 1.0 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 2.0 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 3.0 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 4.0 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 5.0 80.000 80.000 80.000 80.000 80.001 80.001 80.001 80.001 80.001 Test : Beat Detection on Tall T W ave Rejection Capability with Simulated ECG According to EC13 S ection 4.1.2.1 Disclosure of performance specifications, c) T all T wave rejection capability wave amplitude (a T ) for which heart rate indication is within the error limits specified in 4.2.7. A QRS test signal of 1 millivolt (mV) amplitude and 100 millisecond (ms) duration, with PAGE 112 112 a heart rate of 80 beats per minute ( BPM ), shall be used; the T wave duration (d T ) shall be 180 ms, and the QT interval (d QT ) shall be 350 ms. QRS amplitude is defined as a r + a s in EC13 A 20 second (sec) monitor s tabilization period shall be allowed before testing. The sample rate of the simulated ECG is 250 Hz. The testing results showed that the proposed beat detection algorithm can detect the heart rate correctly, and passed the test On the other hand, except for the small T wave amplitude (50% and below of that of R wave) of the simulated ECG, the IMEC algorithm presented double detections (detection of both R wave and T wave) for tall T wave ECGs. Thus, t he IMEC algorithm has no ad equate tall T wave rejection capability, and failed the test Table 8 9. Results of the t est of the p roposed b eat d etection algorithm on t all T wave r ejection c apability with s imulated ECG Amplitude of T wave (mV) 0.3 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1 .4 1.5 Simulated Rate (BPM) 80 80 80 80 80 80 80 80 80 80 80 80 Measured Rate (BPM) 80 80 80 80 80 80 80 80 80 80 80 80 Error (Simulated Measured) (BPM) 0 0 0 0 0 0 0 0 0 0 0 0 Table 8 10. Results of the t est of IMEC b eat d etection algorithm on t all T wave r ejection c apability with s imulated ECG Amplitude of T wave (mV) 0.3 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 Simulated Rate (BPM) 80 80 80 80 80 80 80 80 80 80 80 80 Measured Rate (BPM) 80 80 167 167 167 167 167 167 167 167 167 167 Error (Simulated Measured) (BPM) 0 0 87 87 87 87 87 87 87 87 87 87 PAGE 113 113 LIST OF REFERENCES Ansley CF Newbold P (1980). Finite sample properties of estimators for autoregressive moving average models. J. Econometrics 13 : 159 183. Bearden S, Eisenschenk S, Uthman B (2008) Diagnosis of n onconvulsive s tatus e pilepticus (NCSE) in a dults with a ltered m ental s tatus: clinic e lectroencephalographic c onsiderations. Am J E lectroneurodiagnostic Technol 48:11 37 Berkovic SF, Scheffer IE (19 99 ) Genetics of the epilepsies. Curr Opin Neurol 12:177 182 Bertram EH (1997) Functional anatomy of spontaneous seizures in a rat model of epilepsy. Epilepsia 38: 95 105 Box GP, Jenkins GM, Rein sel GC (1994). Book: Time s eries a nalysis, f orecasting and c ontrol, 3rd edn. Englewood Cliffs, NJ: Prentice Hall. [Chapter 1 5, 7, 8] Brenner RP (2002) Is i t s tatus? Epilepsia 43 ( Suppl 3 ) :103 113. Cascino GD, Jack CR, Parisi JE, Sharbrough FW, Schreiber CP, Kelly PJ, Trenerry MR (1993) Operative strategy in patie nts with MRI identified dual pathology and temporal lobe epilepsy. Epilepsy Research February, Vol. 14, Issue 2: 175 182. Cascino GD (2004) Surgical treatment for epilepsy. Epilepsy Research July, Vol. 60, Issue 2: 179 186. Chaov alitwongse WA, Prokopyev OA, Pardalos PM (2006) Electroencephalogram (EEG) time series classification: a pplications in epilepsy. Annals of Operations Research 148(1):227 250. Charles LL Richard JH (1995) Solving l east s quares p roblems SIAM publishing Philadelphia, PA [Chapter 2, 3, 4] Chris C (2001). Book: Time s eries f orecasting, Chapman &Hall / CRC Press New York, NY [Chapter 4, 5] Doob J (1953). Book: s tochastic p rocess, John Wiley New York [Chapter 3, 4] Engel JJ (1989). Book: Seizures and e pilepsy, F.A. Davis, Philade lphia PA Engel JJ Pedley TA (1998) Book: Epilepsy: a c omprehensive t extbook Lippincott Raven publishing Philadelphia, P A Fountain NB, Waldman WA (2001) Effects of Benzodiazepines on t riphasic w aves. J Clin Neurophysiol 18(4):345 352. PAGE 114 114 Geiger LR, Harn er RN (1978) EEG patterns at the time of focal seizure onset. AMA Arch Neurol 35:276 286. Geweke J (1982), Measurement of l inear d ependence and f eedback b etween m ultiple t ime s eries. J. American Statistical Association June, Vol. 77, No. 378 : 304 313. Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG, Mietus JE, Moody GB, Peng C K, Stanley HE (2000) PhysioBank, PhysioToolkit, and PhysioNet: c omponents of a n ew r esearch r esource for c omplex p hysiologic s ignals Circulation 101(23): e21 5 e220 Goodwin, GC, Payne RL (1977). Book: Dynamic s ystem i dentification: e xperiment d esign and d ata a nalysis, Academic Press, New York. Granger CW (1969), Investigating c ausal r elations by e conometric m odels and c ross s pectral m ethods. Econometrica Vol. 37, No. 3 Jul. : 424 438 Hamilton, JD (1994). Book: Time s eries a nalysis. Princeton Univ Press Princeton, NJ [1.5, 3.1, 3.4, 4.1, 5.2, 5.4] Hayes, MH (2008). Book: Statistical d igital s ignal p rocessing and m odeling, John Wiley & Sons Inc. New York N Y Heinemann SH; Rettig J; Graack HR, Und Pongs O (1996) Functional characterization of K V channel beta subunits from rat brain, J. Physiol. Lond. 493 (Pt 3) Jun. 15 : 625 33 Hirsch LJ, Spencer SS, Williamson PD, Spencer DD, Mattson RH (1991a) Comparison of bitemporal and unitemporal epilepsy defined by depth electroencephalography. Ann Neurol 30:340 346 Hirsch LJ, Spencer SS, Spencer DD, Williamson PD, Mattson RH (1991b) Temporal lobectomy in patients with bitem poral epilepsy defined by depth electroencephalography. Ann Neurol 30:347 356 Iasemidis LD (1991) O n t he d ynamics o f h uman b rain in t emporal l obe e pilepsy. Doctor of Philosophy Dissertation, University of Michigan. Iasemidis LD, Barreto A, Gilmore RL, Uthman BM, Roper S, Sackellares JC (1993) Spatio temporal evolution of dynamical measures precedes onset of mesial temporal lobe seizures. Epilepsia 35S ( suppl 8 ) : 133 Iasemidis L D, Pardalos PM, Shiau DS, Chaovalitwongse WA, Narayanan K, Prasad A Tsakalis K Carney PR, Sackellares JC (2005) Long term prospective on line real time seizure prediction. J Clin Neurophysiol 116(3):532 544. PAGE 115 115 Iasemidis LD, Principe JC, Czaplewski JM, Gilman RL, Roper SN, Sackella res JC (1997) Spatiotemporal transition to epileptic seizures: a nonlinear dynamical analysis of scalp and intracranial EEG recordings. Book: In Silva FL, Principe JC, and Almeida LB, editors, Spatiotemporal Models in Biological and Artificial Systems 81 88. IOS press p ublishing Amsterdam Netherland. Iasemidis LD, Principe JC, Sackellares JC (2000) Measurement and quantification of spatiotemporal dynamics of human epileptic seizures. Book: wwIn Akay M, editor, Nonlinear biomedical signal processing 2:29 4 318. Wiley IEEE press Publishing New York, NY Iasemidis LD Sackellares JC (1991) The temporal evolution of the largest Lyapunov exponent on the human epileptic cortex. Book: Measuring Chaos in the Human Brain. Singapore: 49 82, World Scientific Publis hing Hackensack, NJ Iasemidis LD, Sackellares JC, Zaveri HP Wil liams WJ (1990) Phase space topography of the electrocorticogram and the Lyapunov exponent in partial seizures Brain Topogr 2: 187 201. Iasemidis LD, Shiau DS, Pardalos PM, Sackellares JC (2001) Phase entrainment and predictability of epileptic seizures. Book: In Pardalos PM and Principe JC, ed itors, Biocomputing. 59 84. Kluwer Academic Publishers Dordrecht Netherland Ia semidis L D, Shiau D S Sackellares JC Pardalos PM, Prasad A (2004) Dynamical resetting of the human brain at epileptic seizures: application of nonlinear dynamics and global optimization tecniques. IEEE Trans Biomed Eng 51(3):493 506. Intern ational Leagu e Against Epilepsy (1993) Guidelines for e pidemiologic s tudies on e pilepsy: c ommission on e pidemiology and p rognosis Epilepsia 34(4) : 592 596. Jirsch J, Hirsch LJ. (2007) Nonconculsive seizures: developing a rational approach to the diagnosis and management in the critical ill population. Clin Neurophysiol 118(8):1660 16 70 Juhasz C, Chugani DC, Muzik O, Watson C, Shah J, Shah A, Chugani HT (2000 ) Relationship between EEG and position emission tomography abnormalities in clinical epilepsy. Journal of Clinical Neurophysiology January, Vol. 17, Issue 1: 29 42. Kaplan DT, Furman MI, Pincus SM, Ryan SM, Lipsitz LA, Goldberger AL (1991) Aging and the complexity of cardiovascular dynamics. Biophys J 59:945 949. Kluger Y, Basri R, Chang JT, Gerstein M (2003) Spectral biclustering of microarray data: coclustering genes and conditions. Genome Res 13(4):703 7 16. Kolmogorov AN (1958) A new metric invariant of transient dynamical systems and automorphisms in lebesgue space. Dokl Akad Nauk SSSR 119:861 864. PAGE 116 116 Lebedev MA, Nicolelis MA (2006) Brain machine interfaces: past, prese nt and future. Trends Neurosci 29(9):536 5 46. Lee SK, Lee SY, Kim KK, Hong KS, Lee DS, Chung CK (2005) Surgical outcome and prognostic factors of cryptogenic neocortical epilepsy. Annals of Neurology Vol. 58, Issue 4: 525 532 Liu CC (2008) Brain d ynamics, s ystem c ontrol a nd o ptimization t echniq ues w ith a pplication in e pilepsy. Doctor of Philosophy Dissertation, University of Florida. Makhoul Proc. IEEE vol. 63, pp. 561 580. Makhoul IEEE Trans. Acoust. Speech Signal Process vol. ASSP 25, pp. 423 428. Makridakis S Hibon M (1997). ARMA models and the Box Jenkins methodology. J. Forecasting 16, 147 163. Markel JD Gray JR (1976). Book: Linear Prediction of Speech, Springer Verlag, New York NY Miller KS (1974) Book: Complex Stochastic Processes: An Introduction to Theory and Application, Addison Wesley, Reading, MA. Milton J, Jung P (2003) Book: Epilepsy as a Dyna mic Disease Springer Publishing New York, NY Monto S, Vanhatalo S, Holmes MD, Palva JM (2007) Epileptogenic neocortical networks are revealed by abnormal temporal dynamics in seizure free subdural EEG. Oxford Journals: Cerebral Cortex Vol. 17, Issue 6: 1386 1393. Nehlig A, Daval J, Debry G (1992) Caffeine and the central nervous system: m echanisms of action, biochemical, metabolic and psychostimulant effects. Brain Res Rev 17: 139 170. Netoff TI, Clewley R, Arno S, Keck T, White JA (2004) Epilepsy in s mall world networks. The journal of neuroscience September, 24(37): 8075 8083. Newbold P, Agiakloglou C, Miller J (1994) Adventures with ARIMA software. Int. J. Forecasting, 10, 573 581. Ochi A, Otsubo H, Donner EJ, Elliott I, Iwata R, Funaki T, Akizuki Y, Akiyama T, Imai K, Rutka JT, Snead OC (2007) Dynamic changes of ictal high frequency oscillations in neocortical epilepsy: using multiple band frequency analysis Epilepsia 48:286 296. PAGE 117 117 Orteg a GJ, Prida LM, Sola RG, Pastor J (2008) Synchronization clusters of interictal activity in the lateral temporal cortex of epileptic patients: intraoperative electrocorticographic analysis. Epilepsia. September, Vol. 49, Issue 2: 269 280. Pardalos PM, Sack ellares JC, Carney PR, Lasemidis LD (2004) Book: Quantitative Neuroscience. Kluwer Academic Publisher s New York, NY Pincus SM ( 1991 ) Approximate entropy as a measure of system complexity. Proceedings of the National Academy of Sciences of the United Stat es of America 1991 : 2297 2301 Pincus SM (1995) Approximate entropy (apen) as a complexity measure. CHAOS 5 (1):110 117. Ponten SC, Bartolomei F, Stam CJ (2007) Small world networks and epilepsy: graph theoretical analysis of intracerebrally recorded mesi al temporal lobe seizure. Clinical Neurophysiology. April, Vol. 118, Issue 4: 918 927. Priestley MB (1981). Spectral Analysis and Time Series, Vols. 1 and 2. Academic Press London [Chapter 1 5, 10] Rafael CG Richard EW (2007) Book: Digital Image Processing, Prentice Hall publishing New Jersey [Chapter 5]. Sanes JN Donoghue JP (2000) Plasticity and Primary Motor Cortex, Annual Review of Neuroscience Vol. 23: 393 415 Schwartzkroin PA, Moshe SL, Noebels JL, Swann JW (1995) Book: Brain Development and Epilepsy Oxford University Press New York NY Semah F, Picot MC, Adam C, Broglin D, Arzimanoglou A, Baz in B (1998) Is the underlying cause of epilepsy a major prognostic factor for recurrence? Neurology 51: 1 256 62 Shoeb A (2009) Application of m achine l earning to e pileptic s eizure o nset d etection and t reatment Doctor of Philosophy Dissertation Massac huset ts Institute of Technology Shumway RH Stoffer DS (2000) Book: Time s eries a nalysis and i ts a pplications, Springer publishing, New York NY [Chapter 2]. Simon H (2002). Book: Adaptive Filter Theory, Prentice Hall Press, Hamilton, Ontario, Canada [Chapter 8] Spreafic R, Arcelli P, Frassoni C, Canetti P, Giaccone G, Rizzuti T, Mastrangelo M, Bentivoglio M (1999) Development of layer I of the human cerebral cortex after midgestation: architectonic findings, immunocytochemical identification of neurons and gl ia, and in situ labelling of apoptotic cells. J Comp Neurol 410: 126 142 PAGE 118 118 Srinivasan V, Eswaran C, Sriraam N (2007) Approximate entropy based epileptic EEG detection using artificial neural networks. IEEE Trans Inf Technol Biomed 11(3): 288 295. Stewart GW (1973) Book: Introduction to Matrix Computations, Academic Press, New York NY [Chapter 3] Spencer SS (2002) Neural networks in human epilepsy: evidence and implications for treatment. Epilepsia March, Vol. 43, Issue 3: 219 227. Spencer SS, Spencer DD ( 1994) Entorhinal Hippocampal interactions in medial Temporal Lobe epilepsy. Epilepsia July, Vol. 35, Issue 4: 721 727. Takens F (1981) Detecting strange attractors in turbulence. Book: Dynamical Systems and Turbulence, Lecture m in Mathematics, Springer V erlag Publishing Heidelberg Germany 898:366 381. Th ulasidas M, Guan C, Wu J (2006) Robust classification of EEG signal for brain computer interface. IEEE Trans Neural Syst Rehabil Eng. 14(1):24 9. Towne AR, Waterhouse EJ, Boggs JG, Garnett LK, Brown AJ, Smith JR, DeLorenzo RJ (2000). Prevalence of Nonconvulsive Status Epilepticus in Comatose Patients Neurology 54(2):340 345. Whittle P (1983). Prediction and Regulation, 2 nd edn., revised. University of Minnesota Press Minneapolis [Chapter 4] Willi amson PD, Thadani VM, Darcey TM (1992a) Occipital lobe epilepsy: clinical characteristics, seizure spread patterns and results of surgery Ann Neurol 31:3 13. Williamson PD, Boon PA, Thadani VM (1992b) Parietal lobe epilepsy: diagnostic considerations and results of surgery. Ann Neurol 31:193 201 Wolf A, Swift JB, Swinney HL, Vastano JA (1985) Determining lyapunov exponents from a time series. Physica D 16:285 317. Worrell GA, Parish L, Cranstoun SD, Jonas R, Baltuch G, Litt B (2004) High frequency osci llation and seizure generation in neocortical epilepsy. Oxford Journals: Brain Vol. 127, Issue 7: 1496 1506. Young LS, Rand D (1982) Book: Dynamical Systems and Turbulence, Lecture Notes in Mathematics. 898:366 381, Springer Verlag Publishing New York Zaveri HP, Pincus SM, Goncharova I, Duckrow RB, Spence DD, Spence SS (2009) Localization related epilepsy exhibits significant connectivity away from the seizure onset area. Neuroreport June, Vol. 20, Issue 9: 891 895. PAGE 119 119 Zhang J, Xanthopoulos P Liu CC Bearden S Uthman BM Pardalos PM (2010) Real time differentiation of nonconvulsive status epilepticus from other encephalopathies using quantitative EEG analysis: a pilot study. Epilepsia February, V ol 51, I s sue 2 : 243 250. Zhang J Xanthopoulos P Chie n JH Tomaino V Pardalos PM (2011) Minimum Prediction Error Models and Causal Relations between Time Series, Wiley Encyclopedia of Operations Research and Management Science, January, Vol 5 : 3271 3285, John Wiley & Sons Inc. Publishing New York NY PAGE 120 120 BIOGRAPHICAL SKETCH Jicong Zhang is currently a Ph.D. candidate in the Department of Industrial and Systems Engineering at University of Florida He defended his doctorate dissertation on March 15 and is scheduled to gr aduate in 2011 H e was born in the City of Wuhan in China and spent most of his earl y life gr o w ing up in that city He earned his B .S. degree majoring in electronics informatics and M .S. degree majoring in signal processing and information systems both from the Department of Electronic Eng ineering at Tsinghua University, the t op one u niversity in China. Then he joined the Ph.D. program in the Department of Industrial and Systems Engineering at University of Florida in fall 2006 semester. H e ha s once worked on shipping problem by ships, empty container allocation problem, and perishable inventory (blood) control problem and h e earned the highest score in the Ph.D. general exam in 2007 From the second year of his Ph.D., he selected the topic, da ta mining and optimization in health care under the supervision of Distinguished Professor Panos M. Pardalos. Mr. Zhang ha s more than ten publications in journal s books and proceeding s up to date and he is the first author in five of them In addition, h e has a lot of presentations and posters in conferences and workshops Especially, as a Ph.D. student in engineering school, Mr. Zhang published a full length paper as the first author in the peer reviewed medical journal Epilepsia From 2008 to 2011, he collaborated with Optima Neuroscience Inc., VA hospitals and Cyberonics Inc. to work on data mining, machine learning pattern recognition and prediction problems i n healthcare and biomedical applications S pecifically he develop ed EEG quantification met hods of patient classification and fast diagnosis of generalized NCSE together with his colleagues He also developed improved methods for seizure detection and prediction in epileptic patients Furthermore he has worked on biomarkers indicative of PAGE 121 121 some neurological disorders such as traumatic brain injury. There are many interesting and significant results from Mr. Zhang s work which give strong indications or conclusions for real world application s in biomedicine and healthcare Mr. Zhang is very happy working with his Ph.D. academic advisor Dr. Panos M. Pardalos in the past four years In 2004, Mr. Zhang worked full time at Intel China Research Center for three months In 2010, h e worked in R&D department of Cyberonics Inc. on two major projects and one minor project for eight months He is a member of IEEE and INFORMS Anyone interested in Mr. Zhang CV publicatio n list or other information can contact him at : gatorjicong@gmail.com 