Sistema de predicción epileptogenica en lazo cerrado basado en matrices sub-durales

Sistema de Predicción Epileptogenica en Lazo Cerrado Basado en Matrices Sub-Durales James Brian Romaine Memoria presentada para optar al grado de Doctor por la Universidad de Sevilla Director,Tutor: Manuel Delgado Restituto Director: Ángel Rodríguez Vázquez Departamento de Electrónica y Electromagnetismo Universidad de Sevilla May 2019 I would like to dedicate this thesis to my family, friends and co-workers ... Declaration I hereby declare that except where specific reference is made to the work of others, the contents of this dissertation are original and have not been submitted in whole or in part for consideration for any other degree or qualification in this, or any other university. This dissertation is my own work and contains nothing which is the outcome of work done in collaboration with others, except as specified in the text and Acknowledgements. James Brian Romaine May 2019 Acknowledgements Firstly, I would like to thank Manuel Delgado Restituto and Ángel Rodríguez Vázquez for their constant collaboration and intuitive guidance throughout my PhD. Secondly, I would like to thank Miguel Angel Lagos Florido and Rafaella Fiorelli Martegani for their design guidence and help with test diagnostics and technical designs. Last but not least, I would also like to thank Verónica Siles Santana and Marco Trevisi for their constant support, friendship and patients, throughout the course of this PhD. This work has been supported by the Spanish Ministry of Science and Innovation under grant TEC2016-80923-P and the FEDER Program. Abstract The human brain is the most complex organ in the human body, which consists of approximately 100 billion neurons. These cells effortlessly communicate over multiple hemispheres to deliver our everyday sensorimotor and cognitive abilities. Although the underlying principles of neuronal communication are not well understood, there is evidence to suggest precise synchronisation and/or de-synchronisation of neuronal clusters could play an important role. Furthermore, new evidence suggests that these patterns of synchronisation could be used as an identifier for the detection of a variety of neurological disorders including, Alzheimers (AD), Schizophrenia (SZ) and Epilepsy (EP), where neural degradation or hyper synchronous networks have been detected. Over the years many different techniques have been proposed for the detection of synchronisation patterns, in the form of spectral analysis, transform approaches and statistical based studies. Nonetheless, most are confined to software based implementations as opposed to hardware realisations due to their complexity. Furthermore, the few hardware implementations which do exist, suffer from a lack of scalability, in terms of brain area coverage, throughput and power consumption. Here we introduce the design and implementation of a hardware efficient algorithm, named Delay Difference Analysis (DDA), for the identification of patient specific synchronisation patterns. The design is remarkably hardware friendly when compared with other algorithms. In fact, we can reduce hardware requirements by as much as 80% and power consumption as much as 90%, when compared with the most common techniques. In terms of absolute sensitivity the DDA produces an average sensitivity of more than 80% for a false positive rate of 0.75 FP/h and indeed up to a maximum of 90% for confidence levels of 95%. Table of contents Appendix B Code listings 129 B.1 MatlabCodeDDA .............................129 B.2 Matlab Code Threshold . . . . . . . . . . . . . . . . . . . . . . . . . . 132 B.3 VHDLCodeDDA .............................133 B.4 RTL-compiler................................148 B.5 MBED-C++ test platform code . . . . . . . . . . . . . . . . . . . . . . 149 Appendix C Pseudo codes 155 C.1 Pseudo code DDA hardware flow . . . . . . . . . . . . . . . . . . . . . 155 C.2 Pseudo code for training and calculation hardware flow . . . . . . . . . 157 Appendix D Control signals for future 16-channel testing platform 159 Appendix E James Romaine, Publication list and curriculum 163 List of figures 1.1 Main lobes of the larger brain . . . . . . . . . . . . . . . . . . . . . . . 7 1.2 Action potential and anatomy of the neuronal cell . . . . . . . . . . . . 10 1.3 Overview of common neuronal oscillatory bands . . . . . . . . . . . . . 12 1.4 Example of phase rotations around a stationary wave . . . . . . . . . . 15 1.5 Graph theory, application to real data . . . . . . . . . . . . . . . . . . 18 2.1 Transition periods and approximated phases of two signals . . . . . . . 34 2.2 Synchronisation indexes obtained through the DDA algorithm . . . . . 35 2.3 Results from DDA and Hilbert transform approaches . . . . . . . . . . 36 2.4 Effects of quantisation on algorithm . . . . . . . . . . . . . . . . . . . . 38 2.5 Effects of sample distribution on algorithm . . . . . . . . . . . . . . . . 39 2.6 Effects of window size on algorithm . . . . . . . . . . . . . . . . . . . . 41 2.7 FR1084Block0161, electrodepositionsM5................. 43 2.8 FR1084Block0161, electrodepositionsGA2................. 44 2.9 SNReffects ................................. 45 2.10 patient #7, SI results from sensitivity test. . . . . . . . . . . . . . . . . 48 2.11 Zoom synchronisation indexes for patient #7 . . . . . . . . . . . . . . . 50 2.12 Average Sensitivity against FPR for DDA + SI and PLV + HT . . . . 50 3.1 2 channel synchronisation processor block diagram . . . . . . . . . . . . 56 3.2 Minama detection logic and application . . . . . . . . . . . . . . . . . . 58 List of figures 3.3 Exponential filter block diagram . . . . . . . . . . . . . . . . . . . . . . 62 3.4 VHDL simulation of parameters and verification . . . . . . . . . . . . . 65 3.5 VHDL simulation of minima detection . . . . . . . . . . . . . . . . . . 65 3.6 VHDL simulation showing the results for phi and detection . . . . . . . 65 3.7 VIVADO schematic of the prototype synchronisation processor . . . . . 66 3.8 Detailed VIVADO schematic of the prototype synchronisation processor 66 3.9 CADENCE layout synchronisation processor . . . . . . . . . . . . . . . 67 3.10 2 channel synchronisation processor correlation and power trade-off’s . 70 3.11 Bar chart showing the FPGA resource allocation for the designed DDA algorithm .................................. 71 3.12 AMS-0.18-micro microscopic view . . . . . . . . . . . . . . . . . . . . . 72 3.13test-setup .................................. 72 3.14spimeasurement............................... 73 3.15mbed-setup ................................. 76 3.16 FR1084Block0161, electrodepositionsM5................. 77 3.17 EEG connectivity maps for patient #7 . . . . . . . . . . . . . . . . . . 81 3.18 Graph theory results for patient #7 . . . . . . . . . . . . . . . . . . . . 82 3.19 Graph theory contour map results for patient #7 . . . . . . . . . . . . 83 4.1 16-Channel ASIC implementation top level block diagram overview . . 87 4.2 Microscopic view of 16-Channel ASIC . . . . . . . . . . . . . . . . . . . 88 4.3 1 to 160 SPI for neural input data . . . . . . . . . . . . . . . . . . . . . 90 4.4 160 to 1 SPI for captured synchronisation values . . . . . . . . . . . . . 90 4.5 Control logic for 16-channel device . . . . . . . . . . . . . . . . . . . . 93 4.6 16-channel synchronisation processor flow . . . . . . . . . . . . . . . . . 94 4.7 Neural input multiplexing scheme and dedicated control logic . . . . . . 96 4.8 Modified exponential filter block diagram . . . . . . . . . . . . . . . . . 98 4.9 Training and calculation module block diagram . . . . . . . . . . . . . 103 List of figures 4.10 Dedicated control logic overview for memory latching . . . . . . . . . . 104 4.11 Array of 15 synchronisation processors synthesised using VHDL code . 106 4.12 Single synchronisation processor VHDL synthesis . . . . . . . . . . . . 106 4.13 Training logic VHDL simulation . . . . . . . . . . . . . . . . . . . . . . 107 4.14 Training logic VHDL simulation2 . . . . . . . . . . . . . . . . . . . . . 107 4.15 Output results from VHDL simulation of 16-channel processor . . . . . 109 A.1 graphtheory basic concepts . . . . . . . . . . . . . . . . . . . . . . . . . 127 List of tables 1.1 Basic overview of brain functions for specific lobes . . . . . . . . . . . . 7 1.2 Synchrony in different neurological disorders. . . . . . . . . . . . . . . . 20 1.3 State of the art comparison table . . . . . . . . . . . . . . . . . . . . . 28 2.1 Records used in the verification of the phase synchronisation processor 42 2.2 DDA and PLV results sensitivity test . . . . . . . . . . . . . . . . . . . 49 3.1 Power breakdown of the synchronisation processor based on the blocks inFig.3.1 .................................. 75 3.2 Normalised performance summary . . . . . . . . . . . . . . . . . . . . . 75 4.1 Table showing which neural data streams are connected to specific processor during each of the 8 epochs of time . . . . . . . . . . . . . . . 92 4.2 Table showing the memory allocation during calculation phase of the patient specific thresholds . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.3 Table showing the signals for processor 1 during VHDL simulation of 16-channeldesign..............................108 4.4 Table showing all design resource allocation for the DDA hardware implementation................................108 Preface This thesis was written to give a broad overview of the work conducted during the doctoral scholarship over a period of approximately 5 years. The main goal during the life of the project was the design and implementation of low powered and lower resource consuming VLSI hardware, for the efficient and reliable detection of synchronisation between neural assemblies. The work was conducted in several key stages: Firstly, through the research of state of the art in the field of neurological disorders and synchronisation, we gathered a thorough understanding of the current limitations and design challenges faced in this field of study. Using this information, we then targeted the design of a new mathematical algorithm, which focused on solving the most commonly faced problems. Furthermore, to test the algorithmic design a MATLAB implementation was introduced which was compared against the state of the art in terms of comparable features and accuracy. At this point, important aspects are tested such as the signal-to-noise ratio, algorithmic sensitivity, sliding window size and sampling rate distribution to name a few. To follow, the algorithm was implemented into a field programmable gate array (FPGA) using a hardware description language (VHDL), which allowed us to verify the algorithm running at the logic gate level. This also provided us with prototype metrics which could then be used to validate the scalability of the device compared to the state of the art. The FPGA implementation of the design was later mapped into CADENCE as a register technology level (RTL) circuit. This design, operating at 1.2V, gave us a good estimation as to the real life functionality of the circuit and algorithmic design. Once the 1.2V design was verified at the RTL level in CADENCE, a 0.5V sub-threshold biased design was implemented using custom layout techniques. Afterwards, the ASIC design was tested in the laboratory using a mixture of test apparatus such as a logic analyser for basic variable verification and micro-controllers for longer in-depth neural signal simulations. Later, a 16-channel mass VLSI implementation was integrated. This design incorporated 15 main synchronisation processors based on the previously integrated 2-channel processor. In order to integrate this design a similar method to that of the previous design was adopted. This included the FPGA implementation and verification before a CADENCE based custom layout and verification. NOTE: Due to a change in the MPW policy of the selected foundry, AMS AG, regarding its 0.18 µ m CMOS process (MPWs were finally cancelled in 2018), the 16- channel VLSI DDA processor was received 13 months after the tape-out. The limited time that was available for the completion of the thesis meant that this design could not be tested in the laboratory. Nevertheless, FPGA simulations have been incorporated into the thesis as a verification of functionality. Organisation This thesis has been separated into five chapters. Chapter 1 describes the most relevant principles needed for the rest of the Thesis. Starting with a brief overview of the anatomy of the brain at the system level ("the larger brain"), including its distinct regions and basic functions, we later go on to describe the anatomy of the brain at the cell level ("the smaller brain") and the electrical signals, such as action potentials (APs) and local field potentials (LFPs), which are involved. Moving forward, we introduce the concept of phase and synchronisation in both stationary signals and neural signals. This includes the introduction of the most common methods of detection and how synchronisation is being treated as one of the most fundamental underlying principles for the encoding and transfer of information between neuronal clusters. Afterwards, we give a small overview of some of the most common neurological disorders which suggests that synchronisation could be used as a tool for the detection of pathological states. Finally, we introduce the most recent advances in hardware implementations for synchronisation detection algorithms. Chapter 2 presents the proposed Delay Difference Analysis (DDA) algorithm for phase synchronisation computation. The chapter discusses the mathematical premises behind the algorithm and illustrates its functionality in real world scenarios. It also shows how parameters such as quantisation, sample distribution or window evaluation length, affect synchrony estimations. Next the algorithm is verified in the Matlab environment using neural recordings from 10 epilepsy patients. These recordings sum in total over 1800 hours of neural data List of tables 3 and contain over 60 annotated seizures. The noise tolerance and sensitivity performance of the DDA algorithm is evaluated and compared against more conventional techniques to give a true impression of the algorithms capabilities. Chapter 3 introduce the design of a VLSI phase synchronisation processor based on the algorithm described in chapter 2. The chapter firstly gives a broad overview of the design including its strengths and weaknesses. Next, the steps followed towards the implementation of the VLSI design are presented. They included the design in an FPGA, followed by the synthesis with an RTL compiler and finally the custom layout in the CADENCE environment. Then, the experimental results from the prototype are introduced. After describing the laboratory test setup, a breakdown of the design’s power consumption is presented. Then, several measurements using neural recording data from the mentioned epilepsy dataset are presented. They include functional connectivity and graph theoretical results. Chapter 4 includes the design and implementation of a multi-channel phase synchronisation processor based on the initial design in chapter 3. This section is laid out similar to the previous section in that firstly, all of the individual digital blocks and digital flow functionality is introduced. This is accompanied by many illustrations of the digital blocks as well an illustration of the overall flow and the implemented ASIC micro-photograph. Next, a FPGA implementation of the VLSI circuit is introduced and system simulations are presented. Finally, Chapter 5 provides some final concluding comments and draw possible future work which will be conducted using the designed algorithm and VLSI implementations. 10 Introduction at the input of the cell. If there are more excitatory signals, the membrane potential increase towards the AP threshold. For the signal to end, the synaptic cleft must then be cleaned of neurotransmitter. Electrical synapse have a direct connection in the form of a gap junction, which allows current to flow extremely fast from one neuron to the next. This has the advantage of allowing neural circuits to synchronise rapidly. One of the benefits of an electrical synapse is the ability to allow current flow in both directions, meaning a bi-directional flow of information. It is believed that these types of synapses are used in fast response actions such as the ’fight or flight mode’ in humans [5]. Axon terminalsExtracellular fluid Flow, Action potential Nucleus Myelin sheath Dendrites Soma Information in Information out -55mv 0mv +30mv -70mv Sodium channels open More sodium channels open Sodium channels close Potassium channels open Potassium channels close Depolarisation repolarisation Hyperpolarisation Action potential Neuronal cell Time Potential Fig. 1.2 Overview of neural connectivity, left, shows an example of an AP including approximate channel opening and closing times, right, shows the basic structure of a single neuron 1.2.2 Neural oscillations The term neuronal oscillation refers to periodic variations in neural activity. These variations can be measured at different scales such as the microscopic scale, which could consist of the measuring of direct synaptic currents which produce spike trains of AP’s, the mesoscopic scale, which consists of the summation of many synaptic currents from a Local network, named a local field potential (LFP), or the macroscopic scale, in the form of electroencephalography (EEG) data, which is the measure of multiple local networks [6]. The history of neural oscillations started with the first ever recordings of individual microscopic recordings. First recorded by Richard Caton in 1875, a physician from 1.2 Connectivity and the smaller brain 11 Liverpool, England, who used electrodes to probe the exposed surface of the brain in animals. Using a reflecting galvanometer, he was successful in monitoring signals in the microamperes range [7]. Nevertheless, it was the work of Berger in 1929, who first reported that he had identified oscillations of approximately 8-12Hz coming from the brain. Berger named this the alpha band since it was the first oscillation detected from such a source. This was early on in his pioneering work on the study of EEG. During one of his experiments Berger noticed that while the patients eyes were closed the prominent oscillations of the alpha band were visible. In the power spectrum the alpha waves produced a lone prominent peak of at approximately 12 Hz. However, when the patient opened their eyes the oscillations seemed to diminish in amplitude and become faster and more erratic, in a way the signals seemed to desynchronise and the prominent peak disappeared. This is known as the ’Berger effect’. Noting this, it was apparent that the alpha oscillations were attenuated by other more complex frequencies due to the visual stimuli originating from the occipital lobe. Further observations of Hans Bergers works were carried out by W. Grey Walter in 1934, who devised an experiment entitled ’Thought and Brain: A Cambridge experiment’, in Which one set of results was quoted as: ’When the eyes of a subject were open the lines were irregular, but when the subject’s eyes were closed they showed as regular sinusoidal type signals occurring approximately every 10 seconds. Interestingly, When the subject closed their eyes and was asked to perform a simple arithmetic problem, the regular patterns disappeared, similar to that of when the patient’s eyes were open. The lines continued to appear irregular until the simple problem was solved at which time the regular patterns reappeared’. The fundamental observation by these pioneers revealed that the electrical fields recorded, must have been the combination of thousands of neurons firing in synchrony as stated by [8]. Since then, further advances in technology have led to the identification of multiple active frequency bands in which oscillations tend to occur. Namely, • Infraslow ( < 0 . 2 Hz ) and δ = 0 . 2 − 3 . 5 Hz . These wave patterns are usually seen during sleep of newborn infants or adults with serious brain diseases. •θ = 4 − 7 . 5 Hz . These waves have been identified to mainly oscillate in the parietal and temporal regions of children’s brains, and in healthy adults. Nonetheless, these slow varying oscillations are absent or indistinguishable under normal 12 Introduction conditions. Nevertheless, these oscillations do appear in healthy adults who show signs of emotional stress or frustration. •α, µ = 8 − 13 Hz . These oscillations tend to occur during the period when a person is awake but in a very relaxed state, mainly with the eyes closed. The introduction of complex thought was proven to also diminish alpha waves causing alpha blocking due to complex thought. •β = 14 − 30 Hz . The lower spectrum of the beta band usually falls into the same category as the alpha band, in that it can be disturbed by small amounts of complex activity or tension. The upper spectrum of the beta band is only affected by extreme thought processes, such as complex arithmetic. The beta band is usually considered as the awake and alert band, as we constantly process visual and auditory senses, as well as perform cognitive tasks throughout the day. Other oscillatory bands include, γ = 30 − 90 Hz and the upper limit of High Frequency Oscillations HFO > 90 Hz . Recent research is, now starting to identify the gamma band as a communication layer between individual neural circuits [ 9 ]. Fig. 1.3, shows the decomposition of an original neural signal measured using EEG, into some common respective oscillatory bands. Fig. 1.3 Example of the most common neuronal oscillatory bands The interconnection of neurons forms several states of connectivity, three main types are as follows [10] : 1.3 Phase and Synchronisation 13 1. structural connectivity, which refers to the anatomical connectivity via the structure of major tracts (i.e, white matter) connecting local and remote cortical regions of the brain. 2. functional connectivity, which is a statistical measure of the amount of connectivity in one region vs another and their dependence upon each other. 3. effective connectivity, the amount of connectivity in one cortical region which directly affects another cortical region by either activation or depression 1.3 Phase and Synchronisation 1.3.1 Basic concept Phase is defined as the angular change in a signal over time. This can be represented as the anticlockwise rotation around the unit circle. Each adjustment is measured in terms of degrees, ranging from 0 ◦ to 360 ◦ , which indicates a full rotation around the unit circle. The rotation can also be measured in the SI unit radians, where 360 ◦ = 2 π , radians. A sinusoidal wave can be expressed as: x(t) = A·sin(ωt +θ)(1.1) where, A is the amplitude, θ (this has no relation to the theta oscillatory band described above) is the initial phase of the signal and ω = 2 πf denotes the rate of change of the vector rotation around the unit circle. The rate of change of phase increases or decrease as the frequency of the signal. If the frequency remains constant, the rate of change in phase will also remain constant. These types of signals are known as stationary signals. Non-stationary signals, however, refers to a signal in which frequency changes over time, such as a synthetic sinusoidal chirp in which the rate of change in phase is constantly increasing. Another non-stationary source is that of EEG signals as LFP’s tend to incorporate a wide range of frequencies. 14 Introduction Fig. 1.4, shows an example of the change in phase angle from a typical sinusoidal waveform and a cosine waveform. In this case, the first peak in the sine wave represents 0 radians and the first minimum in the cosine wave represents, π radians. As the two signals evolve, by the time the first minimum has appeared in the sine wave the phase angle has increased to πradians and the cosine has evolved to 2πradians. Since the sine wave and cosine wave are both stationary at 10Hz, the phase of both signals increases at a fixed rate. In comparative situations, it is desirable to take the phase difference between the two waveforms. This can give rise to information on their frequency offsets with respect to one another and can be calculated as: θ=ϕsignal1−ϕsignal2(1.2) Using specific markers in each signal such as minima or maxima, if the first derivative d(ϕsignal1−ϕsignal2) dt , is constant it can be said that the signals have a constant phase difference. Since the two signals in Fig. 1.4 are of a constant 10Hz, the phase difference between the two is constant. In this case, the constant phase change is π radians [11]. 1.3.2 Synchronisation and phase locking In general terms, synchronisation is defined as any two-coupled (phase locked) oscillating systems, which have a relationship in both phase and frequency, but at the same time are independent of amplitudes. In non-disturbed and noise free systems synchronisation is usually described as the phase locking between instantaneous phase angles such that n·ϕs1−m·ϕs2 = const , where ϕs1 and ϕs2 are the instantaneous phase angles of signal 1 and signal 2, respectively. m and n, are integer based weights which need adjusting based on the source of the oscillating system. In the example in Fig. 1.4, although the two signals are in anti-phase there does exist a relationship between the two, in terms of phase locking. In the case of non stationary signals, the phase locking condition should be considered as: |n·ϕs1−m·ϕs2−δ|< const . Where δ , is some average phase shift. This condition indicates that although the rate of change is not completely constant, it should reside (fluctuate) around a common constant [12]. 1.3 Phase and Synchronisation 15 θ θ θ θ θ θ 0,0º π/2,90º π,180º 3π/2,270º Phase difference 2π,360º π/2,90º π,180º 3π/2,270º 2π,360º π/2,90º π,180º 3π/2,180º 2π,360º π/2,90º π,180º 180º-0º 270º-90º 360º-180º 180º-0º90º-270º180º-360º Fig. 1.4 Phase: example of how phase rotates anti-clockwise around the unit circle in a stationary sine wave 16 Introduction In many cases of phase locking, the values m and n are usually set to a fixed 1:1 ratio. This is especially true if the signals originate from the same source such as the brain. Nevertheless, this phase locking attribute also holds true for signals of different frequencies, such as a stationary sine wave of 10Hz and a stationary sine wave of 20Hz. Although the frequencies are not the same there does exist a relationship in phase known as cross-frequency phase locking [13, 14]. 1.3.3 Synchronisation in the brain The human brain seems to effortlessly communicate large quantities of information over multiple regions of the brain. However, one of the questions still unanswered is how neural communication between assembles throughout the larger brain actually works. That‘s to say, what decides what information is accepted or declined at the input of neural assemblies. For many years it was assumed that general information processing was carried out by changes in firing patterns of so-called ’smart neurons’. which were dedicated to certain tasks and functions. Nevertheless, the work derived from Hans Berger and others has led to the question as to whether synchronisation of neural assemblies is the key to information processing and learning. The fact that neural oscillatory patterns exist, indicates that to some degree, neural synchrony at the given frequency must be intrinsic of some type of neural assembly collaboration. Over recent years, new evidence is emerging that precise synchronisation and desynchronisation of neural assemblies may, in fact, be the key to the encoding of information transfer within the brain. This can be seen in various reviews and experimental results from [15–17] Two main hypotheses have been suggested for the efficient transfer of information between neural circuits and brain regions. Namely, coherence, in the form of phase synchronisation and gating by inhibition. Phase synchronisation could be used as an identifier between neural assemblies, similar to that of two radios set to the same frequency channel. If the presynaptic response from one neural circuit is oscillating at the correct excitability rate of the postsynaptic neurons, then the information is accepted. However, if the presynaptic response is oscillating at a different rate (different frequency), then the postsynaptic inputs of the second neural circuit then the information is blocked or diminished, as with a radios set at different frequencies. 1.4 Synchronisation in pathological brain states 17 In general terms, it is the phase of an oscillatory neuron population which determines their excitability, hence their spike timings. It is, therefore, the phase relationship between neuron populations which determines routing. This has been proposed to occur in the gamma band [18]. The flow of information is controlled by an increase of alpha oscillations which inhibit the neural population which is not firing in synchrony [19]. 1.4 Synchronisation in pathological brain states The topic of the synchrony between neural circuits for proper functional neuronal communication is a trending topic in neuroscience today. When comparing brain states we like to consider the differences between normal and pathological brains, in order to identify the connectivity differences which may exist between them. Therefore, the data of pathological brain state patients is usually compared with that of a control subject (a patient with no medical history of neurological disorders). Nevertheless, in the case of neurological disorders such as epilepsy, the patients inter-ictal epileptic period can be associated with a non-pathological brain state, hence used in further analysis. In order to distinguish between these two brain states, we use a common mathematical tool known as graph theory. This allows us to gather insight into the structural, functional and effective connectivity of the brain. It allows us to visualise the overall connectivity as a function of time. For an in-depth look into the most common mathematical calculations and an overview of the steps required to calculate complex graph networks. Please see appendix A. It is widely accepted that the brains neural communication is characteristic of a small world network where neural clusters are well connected and transfer information over small distances. Therefore comparing the small-worldness of a non-pathological brain state with that of pathological states can identify the changes which occur in terms of overall functional connectivity. Fig. 1.5, shows an overview of the steps required for identifying small-worldness in epilepsy when taken from EEG sources. 18 Introduction , dd CL Synchronisation values Randomised Synchronisation values Binary matrix Threshold Threshold Graph Graph Binary matrix , rand rand CL , dd rand rand CL CL Fig. 1.5 Example of how graph theory methods can be applied to real data. In this case, we compare an epileptic patients ictal period to that of their non ictal period. We then calculate M and A using an binary/undirected graph and calculate metrics as per needed 1. Firstly, all synchronisation values for every possible combination of electrode positions are computed, creating a O = NXM matrix. Where N corresponds to the electrode combination and M corresponds to each synchronisation value for the given electrodes over time. 2. Next, two matrices G 1 = NXN For the non ictal period and G 2 = NXN for the ictal period are constructed. This is done by creating a single mean representative value for each electrode combination from Osuch that: Gx(i,j)=1 M M X z=0 Oi,j (1.3) 3. The next step is to apply thresholds to each matrix G, for the purpose of creating two binary based adjacent matrices Ax . The threshold which is applied is subject to the experimental data and the related subject at hand. In many cases, various thresholds are tested. This allows for the identification of the strongest connections and builds 3d topologies which give rise to information on C, L, k. 1.4 Synchronisation in pathological brain states 19         ϕ11 ϕ12 ϕ13 . . . ϕ1n ϕ21 ϕ22 ϕ23 . . . ϕ2n . . .. . .. . ..... . . ϕd1ϕd2ϕd3. . . ϕdn         =        0 0 1 . . . 0 1 0 1 . . . 1 . . .. . .. . ..... . . 0 0 1 . . . 1         4. Finally, we can calculate λ , γ and η , by calculating the mean clustering coefficients and the mean shortest path lengths.to distinguish the small-world characteristics and small-worldness. 1.4.1 Synchrony in neurological disorders Although each neurological disorder has distinct causes and effects, synchrony has been found to manifest itself as an identifier, to distinguish between the normal and neurologically damaged brain. Table 1.2 [ 20 ], shows an overview of the most common neurological disorders as described above and how long and short-range communication is affected by increases or decreases in synchronisation. Other numerical methods such as graph theoretical analysis (please see appendix A) have also been shown to be a useful tool in the identification of both functional and structural connectivity. In Alzheimer Disease (AD) synchronisation has been observed to increase in the θ and δ bands whilst, conversely, reducing in both the α and β bands at resting state. Using a synchronisation likelihood (SL) measure, it was determined that these bands decreased in synchronisation over both long and short range recording sites [ 21 – 23 ]. Whilst global efficiency measures are higher in AD patients compared to the normal brain, the brain’s ability to synchronise decreased significantly over multiple EEG bands in AD patients when their eyes were closed according to [ 24 ]. Synchronisation in AD patients has also been analysed in patients during memory working tasks. These include EEG and functional Magnetic Resonance imaging (fMRI) analysis of functional connectivity. Once again it was shown that a decrease in synchronisation was present when compared to control patients specifically in the αand βbands [25]. The way our brain synchronises does not just affect cognitive functions. In fact, there is evidence which suggests that synchronisation in the β band is related to pre-motor movements (i.e preparation for movement). Nevertheless, once the actual movement takes place the β activity is replaced by γ activity [ 26 , 27 ]. Moreover, there is evidence to suggest that an increase in β activity could be responsible for tremors in Parkinson Disease (PD) [ 28 ]. Where in patients who do not suffer from PD a short 26 Introduction certain bin l, the index is given by, Sse = 1 + 1 ln(L) L X l=1 Pse(l)·ln[Pse(l)] (1.24) where L is the number of bins. Note that Sse is comprised in the interval [0 , 1], where Sse = 0 corresponds to a uniform distribution (no synchronisation) and Sse = 1 corresponds to a distribution localised in one point. The estimate of this index heavily depends on L[44]. Another synchronisation index is based on the conditional probability Pcp ( k, l )to find the phase of a first signal in a bin k when the phase of the second signal falls in the l -th bin, l = 1 , . . . , L . Denoting as Ml the number of hits of the second signal in the lbin, a synchronisation index can be defined as, Scp =1 L L X l=1  1 Ml Ml X k=1 exp[jPcp(k, l)] (1.25) It has been found that this index is particularly suitable for revealing weak interactions between signals [44–47]. 1.6 VLSI implementations of PS estimators As mentioned in Sec.1.5.1, one popular approach for the calculation of the instantaneous phase of a signal relies on the construction of the analytic representation sa ( t ) = s ( t ) + j˜s ( t ), where ˜s ( t )is the Hilbert Transform of s ( t )[see 1.4]. Using time-frequency domain techniques, the Hilbert Transform can be derived by [ 48 ] (i) obtaining the Fourier transform of s ( t ), (ii) nulling the negative frequency components, (iii) calculating the inverse Fourier transform and (iv) taking the imaginary part of the result. Given the non-stationary nature of neural signals, the Fourier transform should be applied over short intervals. This suggests the use of the Short Time Fourier Transform (STFT) algorithm instead of the classical Fast Fourier Transform (FFT) [36]. 1.4 and 1.19 A simpler method for the generation of the analytic signal representation sa ( t )uses Hilbert Transformers based on digital filtering. These transformers can be implemented in complex domain or formed by the combination of two sub-components; a Hilbert filter with ideal transfer function −jsgn ( f )for the generation of the imaginary part of sa ( t ), and a delay equal to that of the Hilbert filter for the real part. In both cases, 1.6 VLSI implementations of PS estimators 27 the total group delay of the signal components should be kept constant to not cause phase errors. The vector analyzer presented in [ 49 ] follows the second implementation approach and employs 16-tap Hilbert finite impulse response (FIR) filters for obtaining ˜s ( t ), and 16-tap All-Pass FIR filters for implementing the delay. In a more recent implementation [ 50 ], the in-phase and quadrature components are directly provided by the mixed-signal acquisition front-end, and four 64-tap FIR filters are used for band selection. Another transform-based approach for deriving the instantaneous phase of a signal uses Wavelets (see 1.18). In this case, the on-chip implementation of the method requires Multiply-and-ACcumulate (MAC) processors to correlate the input time series with the wavelet templates. As an example, [ 51 ] proposes a charge-recycling mixedsignal MAC processor for epilepsy detection in which time series of 1024 samples are serially loaded and correlated in parallel with all the Morlet wavelet templates (coded in 32 frequency bins with 4-bit coefficient resolution) stored in the on-chip memory. The trigonometric functions involved in 1.7 and 1.17 as well as the computation of the phase locking value in 1.19 can be readily implemented by CORDIC (COordinate Rotation DIgital Computer) units. This is actually done in [ 52 ] where a digital signal processing (DSP) block formed by three 16-bit CORDIC cores and two 32-tap moving average FIR filters is used for computing the phase locking value ( PLV ) from the analytical signal sa(t). Numerical methods for computing PS, such as the one based on probabilities of recurrence, pr ( τ ), described in Sec. 1.5.2, typically use simpler logic than transformbased techniques at the price of increased latency and more memory resources. Thus, for instance, the evaluation of pr ( τ )in 1.22 requires the collection and combination of phase space trajectories, demanding for a large memory allocation. Nevertheless, it is worth observing that, depending on τ , there are coincident terms in the sum 1.22. This opens doors for reducing the computational complexity of the algorithm by avoiding recalculations. This is actually the approach followed in [ 53 ], where a differential scheduling for the computation of the correlation integral of a signal (similar to the calculation of recurrences) allows for a substantial reduction in complexity. Following the calculation of pr ( τ ), a MAC unit should be used for the calculation of the cross-correlation in 1.23 to serve as a synchronization index. As shown in 1.24 and 1.25, numerical methods can also be used for the computation of synchronization indexes. Although the computational complexity is relaxed as compared to PLV (in particular for the index based on Shannon entropy), long sliding 28 Introduction Table 1.3 Performance summary of proposals related to on-chip PS calculation Process Supply Core Area Input rate Clk Freq. Power cons. Computation capability [49] 0.13µm 1.2 V 1.86 mm27.2kHz (8-bit) 10MHz 400 µW 32 input-pairs multiplexed phase synchronization processor [50] 0.13µm 1.2 V 1.27 mm2NA NA 260 µW 32 input-pairs multiplexed phase synchronization processor [51] 0.35µm 1.65 V 16 mm2NA 50Hz NAa128×1024 binary MAC operator [48] 0.13µm 1.2 V 1.18 mm240kHz (16-bit) 10.24MHz 3.64.4 µW 16-channel multiplexed 16-bit Hilbert Transform block [52] 0.13µm 0.85 V 0.178 mm21.7kHz (8-bit) 2.5MHz 102 µW 10-bit PLV Processor [53] 90nm 1.0 V 0.143 mm2256Hz (9-bit) 4096Hz 2.34 µW 16-channel Correlation Integral Processor [55] 90nm 1.0 V 1.02 mm2256Hz (16-bit) 522.24 kHz 57.3 µW 1-channel EMD processor extracting 5 IMFs and residue a Power consumption mainly due to the 4 × 128 8-bit Σ∆ algorithmic ADCs. At 15 kHz parallel sample rate, the bank of ADCs dissipates 6.3 mW from a 3.3 V supply [54]. windows are needed to populate the histograms and obtain meaningful distribution of probabilities. Hence, circuit simplicity is counterbalanced with an increase in latency and memory depth. Table 1.3 illustrates the performance of the different integrated circuits presented in this section. The two first rows are for complete PS processors while the rest corresponds to building elements which can be eventually used depending on the particular PS implementation strategy. The table clearly shows that the algorithms used for PS calculation, although easily programmable in computers, tend to occupy large area and consume significant power when dealing with silicon integration. The problem is even more acute in the context of functional brain connectivity, where multiple phase synchronisation indexes from different sites are needed. For instance, if functional connectivity has to be measured between 16 recording sites, as it is typically done in a standard 10-20 EEG electrode system, 120 different combinations of neural signal pairs should be simultaneously addressed. This means application specific integrated circuits (ASIC’s) for functional connectivity quickly become far to big and/or power hungry, even if sophisticated multi-thread sharing DSP techniques are employed. 1.7 Thesis Objectives This thesis aims to address three major aspects: • Development of an algorithm for the detection and quantification of phase synchronisation between pairs of neural signals which significantly alleviate the constraints of power consumption and area occupation as compared to state-of-art methods. The algorithm should exhibit little or no degradation 1.7 Thesis Objectives 29 on performance metrics such as seizure detection sensitivity, noise tolerance or functional connectivity computation. • Implementation of a low-power, low area silicon prototype of the proposed algorithm to confirm its suitability for VLSI, and performance validation with human recordings in a realistic application scenario, namely, the detection of seizures in patients with different epileptic pathologies. • Demonstration of the scalability of design through the implementation of a multichannel processor for surface or intracranial EEG signals, and confirmation that the proposal represents a low-cost, low-complexity, programmable and portable solution for computing functional connectivity in wearable/implantable systems. Chapter 2 Delay Difference Analysis for Phase Synchronisation Computation 32 Delay Difference Analysis for Phase Synchronisation Computation In this chapter, an algorithm for phase synchronisation computation, which aims to alleviate the area/power consumption burden of silicon integration of previous approaches such as in Sec.1.6 is proposed. The synchronisation algorithm, called Delay Difference Analysis (DDA), is based on the piece-wise linear approximation in [ 12 ], however, it has been modified by us to achieve a more hardware friendly implementation. Furthermore, a new synchronisation index is introduced which allows for a more reliable identification of synchronisation peaks. The two components of the algorithm, were then subjected to tests, including: simulation of nearly two thousand hours of real neural data, noise analysis and sensitivity analysis. Moreover, we include the effects of sampling rate, window length and sample distribution. All of the following was carried out by the author of this thesis. 2.1 Proposed algorithm Similar to the transform-based methods, this proposal is also based on the estimation of instantaneous phases in band-limited signals, however, instead of using complex filtering or wavelet convolutions, phases are extracted through the identification of distinct marker events in wave forms. Without loss of generality, in this work, it is assumed that such events correspond to the local minima of the signal at time instants tn . The time interval between two consecutive minima correspond to one complete cycle and, therefore, the phase increment during this period is exactly 2 π . Hence, the signal phase can be approximated for any arbitrary time instant tn< t < tn+1 as, φs(t)≈2πt−tn tn+1 −tn + 2πn (2.1) where the term 2 πn unwraps the phase angle and Ts n = tn+1 −tn measures the duration of the transition period between the minima, as illustrated in Fig.2.1. Note that this approximation notably simplifies phase estimation as the problem basically reduces to a measure of time intervals. For this reason, the proposed algorithm is denoted as Delay Difference Analysis. Based on the approximation in (2.1), an index for quantifying the synchronisation level between two signals has been proposed in [ 56 ]. It relies on creating histograms on how many times the events in one signal are preceded by events in the other. The approach has a low computational cost, however, similarly to the numerical techniques 2.1 Proposed algorithm 33 described in Sec. 1.5.2, it requires long time series to build meaningful statistics and, in some sense, it wastes the time information between events. The proposed approach for measuring synchronisation is similar to the PLV technique in (1.19) but avoids the use of trigonometric functions. DDA analyses the time difference between the transition periods of both signals, not their relative phase △φ , and, remarkably, it operates on-the-fly with no need to store long signal sequences. The procedure simply consists in pairing transition periods from both signals as soon as they are detected. If a pair is formed between the i -th transition period of signal s1 ( t )and the j -th transition period of s2 ( t ), the absolute difference △Tn = |T1 i−T2 j| is calculated, where T1 i and T2 j denote, respectively, the durations of said transition periods. If two or more transition periods in one signal are detected before a transition period is found in the other, only the last transition period of the former signal is accounted for in the computation of △Tn . Fig.2.1 shows an example of the procedure. This, in turn, increases the phase error distribution. Let us denote by K the number of absolute differences △Tn computed in the k -th sliding window of length N , and let us assume with no loss of generality that N = 2 m . A synchronisation index between signals s1(t)and s2(t)can be defined as, S(r) dda(k)=1−2r Nmin K−1 X n=0 |△Tn| − Tos,N 2r!(2.2) where r∈ [0 , m ]is a user-defined selectivity control parameter and Tos is a positive integer offset for adjusting the synchronisation index range. Fig.2.2(a) plots S(r) dda for Tos = 0 and different selectivity configurations assuming that s1 ( t )is a sine wave at 20Hz, and s2 ( t )is a chirp whose frequency linearly sweeps from 10 to 30Hz. Both signals had a sampling rate of 1024S/s and were quantised at 10-b of resolution. Every calculated index was obtained assuming a non-overlapping sliding window of length N = 1024. Clearly, the higher the selectivity parameter, the narrower the band for which the synchronisation index obtains non-null values. No matter the r parameter, the plots in Fig.2.2 are not symmetrical around the tone frequency of s1 ( t ); they decay linearly for chirp frequencies below said tone but fall away more smoothly for higher frequencies. This is because the number of absolute differences K is limited by the smaller frequency, fmin , of both signals as K = N·fmin/fs . In Fig.2.2(a), it is also observed that the peak value of the synchronisation index decreases with the selectivity parameter. This can be compensated, without drastically increasing the sampling frequency or the number of samples per observation window, by means of the offset 34 Delay Difference Analysis for Phase Synchronisation Computation 1 1 T 1 2 T1 3 T 1 4 T 1 5 T AB CD E 2 1 T 2 2 T 2 3 T 2 4 T 0 2 0 2  1s 2s  1 ()st 2 ()st Fig. 2.1 Transition periods and approximated phases of two signals. Starting at point A, the DDA algorithm generates the following absolute differences between transition periods: |T1 1−T2 1|,|T1 2−T2 2|,|T1 3−T2 3|and |T1 5−T2 4|. parameter Tos . This is illustrated in Fig.2.2(b) where Tos = 6 making the synchronised indexes for all selectivity parameters to mostly cover the range between 0, indicating no coupling, and an approximate 1 for perfect coupling. Note that too large of an offset may eventually make the synchronisation index larger than 1, hence, a limiter, which can be easily implemented in hardware by an overflow detector, should be used together with (2.2). 2.1 Proposed algorithm 35 (b) (a) 0 0.2 0.4 0.6 0.8 1 10 12.5 15 17.5 20 22.5 25 27.5 30 Frequency, Hz DDA Synchronization index r=0 r=1 r=2 r=3 r=4 0 0.2 0.4 0.6 0.8 1 10 12.5 15 17.5 20 22.5 25 27.5 30 Frequency, Hz DDA Synchronization index r=0 r=1 r=2 r=3 r=4 Fig. 2.2 Synchronisation indexes obtained through the DDA algorithm for different selectivity parameters, using a sine wave and a swept-frequency cosine as input signals. (a) Tos = 0; (b) Tos = 6. 42 Delay Difference Analysis for Phase Synchronisation Computation Table 2.1 Records used in the verification of the phase synchronisation processor Id Patient Dominant pattern Electrodes aSync. Pairs Sample Rate Seizures Hours #1 male 22y rhythmic sharp waves 1/114/0/0 SCL7 - {SCL8, SCR9, TPL1} 1024 Hz 26 113.20 #2 female 29y rhythmic beta waves 0/121/0/0 HAR1 - {HPR12, HPR1, HPR2} 1024 Hz 9 183.10 #3 female 16y low amplitude fast activity 21/0/12/48 IHB2 - {IHA1, IHA2, GB6} 256 Hz 9 229.30 #4 male 18y repetitive spiking 21/20/46/32 BLB1 - {BLB3, TRB1, FP1} 1024 Hz 13 245.20 #5 male 35y rhythmic beta waves 21/10/46/0 TBLB1 - {TBLB3, TBRB1, FP1} 256 Hz 26 180.00 #6 female 53y rhythmic beta waves 19/25/16/0 HL2 - {HL1, HL4, TBB1} 256-512 Hz 6 164.40 #7 male 5y low amplitude fast activity 21/20/50/0 TBA1 - {TBA2, HRA4, HRA5} 1024 Hz 22 170.60 #8 male 15y rhythmic beta waves 51/72/52/58 GG3 - {TL2, GG4, GH2} 256 Hz 19 199.8 #9 female 32y rhythmic beta waves 36/34/35/37 HL4 - {HL2, HL3, HL5} 1024 Hz 9 162.6 #10 female 11y rhythmic alpha waves 60/57/59/4 HR12 - {HR9, HR11, TBA4} 1024 Hz 14 155.0 a(S/D/M/G) for surface, depth, strip and grid electrodes. contains recordings of multiple patients including 225 scalp recordings, 50 intracranial recordings and over 100 annotated sets using both intracranial and surface electrodes which were organised into a standard 10-20 format. In this work, simulation experiments including the DDA and the PLV with 10 recording sets included in the European Epilepsy database have been conducted. In total, 1803.2 hours of recordings and 61 annotated seizures have been analysed. Simulations were carried out in Matlab using mathematical models of the DDA and the PLV processors. In the DDA case, the model closely follows the block diagram of Fig. 3.1 for N = 1024, r = 1, M = 10 and Q = 2 (see chapter 3, subsection 3.2.1 for definition of M and Q ) and the Matlab code can be seen in Appendix B.1 while, in the PLV case, the model uses the architecture presented in [ 50 ]. In both cases, an input signal resolution of 10-b has been assumed. Table 2.1 shows the records used for the verification of the proposed phase synchronisation calculation approach. For each patient the table shows the electrode configuration (both EEG and invasive), sampling rate, number of seizures, recording duration and dominant seizure pattern. The EEG data were acquired using a Neurofile NT digital video EEG system, while the invasive measurements were obtained with depth, strip or grid electrodes from Ad-Tech. For each patient, three different phase synchronisation indexes, involving four intracranial electrodes, were calculated. In all cases, the three contact pairs have one electrode in common, located in the proximity of the epileptic focus, while the other three electrodes, not necessarily from the same array, were used as reference. These electrodes are identified in Table 2.1 as #F-{#R1, #R2, #R3}, where #F denotes the electrode at the seizure focus and #Rx denote the references. Prior to the calculation of the synchronisation indexes, records were 2.3 Verification 43 Fig. 2.7 (a) Simulated DDA algorithm index and simulated PLV value from recordings at positions M5 and M8 of a 1 × 8EcoG strip implanted in patient FR_1084. Signals were pre-filtered in the β -band. The input rate of the chip was adjusted to the sample rate of the recordings, i.e., 1024S/s. (b) Correlation between both plots. band-pass filtered in the frequency bands where higher activities were observed. This gave the selection of β -band for all patients but patient #4, for which the α -band was most insightful. In the case of patient #10, the β - and θ bands offered similar results. 2.3.1 DDA-PLV correlation For illustration purposes, 2 examples of the correlation results between the DDA SI index and the PLV can be seen in Figs. 2.7 and 2.8, respectively, the synchronisation results obtained from two recording blocks measured in patient FR_1084 (48 year old female). In both blocks (6h long) three low amplitude fast activity (lafa) seizures were annotated (red dotted lines). For each figure, the plots on the left show the simulated synchronisation index (labeled as "SI") and the calculated PLV value (labeled as "SIM PLV"), while the plots on the right represent the correlation between both measures calculated at the time instants in which the experimental synchronisation indexes were recorded. The DDA algorithm detects all major changes in synchrony and indeed follows a similar trend as the PLV algorithm. This is verified by the high correlation results of approximately 90% on average. Figs. 2.7 and 2.8 also show that seizures manifest with a remarkable increase in the synchronisation measures. This feature is in fact on the basis of the detection 44 Delay Difference Analysis for Phase Synchronisation Computation Fig. 2.8 (a) Similar as for Fig. 2.7 using the electrodes GA2 and GA5 of a 8 × 8 subdural grid implanted in patient FR_1084. (b) Correlation between both plots. mechanism presented in [ 52 ] and [ 50 ] in which seizures are detected by applying thresholds on the PLV indicators calculated between pairs of channels. 2.3.2 SNR The simplicity of the DDA circuit comes at the expense of higher susceptibility to noise on the recorded signals as compared to the PLV approach. This is illustrated in Fig. 2.9 which shows the correlation between the ideal PLV values obtained with two noise-free 60 minutes long neural signals, and the synchronisation indexes, derived from the DDA approach, when the same signals are contaminated with white Gaussian noise. The correlation is plotted against the spot Signal-to-Noise Ratio ( SNR ), assuming 1Hz bandwidth, at the maximum frequency of the β band, i.e. 30Hz. PS estimations are derived on that frequency band as well. Spot SNRs in the range between 10dB and 70dB at 5dB steps have been considered and, for each SNR , a Montecarlo analysis (100 instances) has been carried out. The error bars show both the mean value and the ± 3 σ of the synchronisation indexes. It is worth observing that for spot SNRs above 30dB the average correlation between the DDA approach and the noise-free PLV technique is larger than 90%. Also note that for large SNRs , i.e. when the noise contamination is negligible, the correlation is in the order of 95%. This "intrinsic" deviation can be attributed to the structural differences between the DDA and the PLV algorithms. Fig. 2.9 also shows the correlation between the noise-free PLV values and those obtained with a 16-tap Hilbert transformer with passband ripple of 2.3 Verification 45 SNR(dB) 0 10 20 30 40 50 60 70 80 Correlation to ideal PLV 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 SI [M=10, Q=2] PLV Fig. 2.9 Correlation of the DDA synchronisation index and the PLV processor in [ 52 ] with an ideal noise-free PLV calculation when the input signals are contaminated by Gaussian noise. 0.01dB when input signals are corrupted by noise, following the implementation in [ 52 ]. Similar as for the DDA case, performance worsens for SNRs below 20dB. To assess the impact of the noise behaviour shown in Fig. 2.9, it must be emphasised that surface and intracranial EEG signals exhibit 1 /fx -like power spectra at low frequencies [ 58 , 59 ]. This makes low frequency bands more tolerant to the aggregated noise contributed by biological tissue, electrodes and recording front-ends [ 60 , 61 ]. Indeed, it has been shown that even with modest amplifier noise specifications, without suppressing flicker contributions, SNRs larger than 30dB can be obtained at low frequencies [ 62 ]. This suggests that PS calculations can be relevant as a neurological biomarker in the β and lower frequency bands although the reliability decreases for the highγ band and higher frequencies. This observation is indeed confirmed in Section 2.3.1, where it is also demonstrated that DDA deviations compared to PLV can be tolerated without sensibly degrading the quality of the results while the DDA approach 46 Delay Difference Analysis for Phase Synchronisation Computation offers light and modular hardware implementations, easily scalable to large neural recording arrays. 2.3.3 Sensitivity In order to estimate the sensitivity (number of seizures which are correctly detected divided by the total number of real positive cases) of the DDA and the PLV approaches in terms of the False Positive Rate, FPR (number of wrong positive detections divided by the total number of actual negative events), different thresholds per contact pair and patient have been applied over the corresponding synchronisation indexes. Thresholds for the DDA and the PLV algorithms do not necessarily coincide even if applied on the same electrode pair. Thresholds are obtained by scanning the synchronisation results and identifying the index value leading to a specific FPR performance. To identify each of the 5 thresholds for each individual contact pair and patient the following procedure is enforced. Firstly, by identifying the total length of the resulting SI and PLV values (the length is in hours and does not include the pre, post and ictal periods). Secondly, by multiplying the total hours by the desired FPR we can extract the exact amount of false positives allowed during the resulting SI and PLV values. Since we are looking for peaks (false positives), we then search the SI and PLV values for the top x peaks in descending order. From there, using the last peak from the x peaks selected we set a threshold. The Matlab code for the thresholding can be seen in Appendix B.2. In our experiment, a detection is produced whenever one out of three synchronisation indexes (or PLV values) rises above their corresponding thresholds. Hence, decisions are taken based on the outcomes from three contact pairs, not just one. Additionally, the following considerations have been adopted [ 63 ]. An observation window of 15 min before the annotated seizure has been defined for the detection of true positives. Hence, if a crossing is detected any time within this window, which can be associated to a preictal state, a true seizure detection is accounted for. Further, in order to exclude effects from postictal states, recording periods within 30min after the onset of a seizure were discarded for threshold calculations. If two successive seizures are separated by less than 45min, the preictal window before the onset of the following seizure is preserved in the analysis. Fig. 2.11 illustrates a true positive detection identified in patient #7. In this case, the synchronisation indexes calculated with DDA from the three contact pairs pass their respective thresholds during the preictal period. 2.3 Verification 47 Fig. 2.10, shows an example of the full SI index results obtained from patient #7. The 5 horizontal red lines indicate the different thresholds for each FPR and the vertical red boxes indicate the pre/post and ictal periods. All though the results can be somewhat noisy, it is apparent that there do exist very Strong prominent peaks in synchronisation. Fig. 2.12 shows the sensitivity variation with respect to FPR , averaged over all the 10 patients in Table 2.1, according to the one out of three detection rule described before. Five different FPR values, in units of number of events per hour, have been considered. Both DDA and PLV detectors were analysed. Error bars indicate the 95% confidence interval limits per method and FPR value. The average detection threshold variation along the FPR range (from 0.15 to 0.75 FP/hr) was 8.8% for DDA and 14.5% for PLV . Table 2.2, shows all of the true positives, false negatives and sensitivity results which where used for the sensitivity analysis. As can be seen in Fig. 2.12, there are little differences between the two approaches (only 4.6% sensitivity variation averaged over all FPR ), with DDA slightly outperforming PLV detection, excepting at 0.45 FP/hr for which both methods obtain similar results. As the constraint is relaxed in terms of FPR , it is observed that the sensitivity rises and the confidence levels get tighter. In fact, at 0.75 FP/hr, DDA achieves an 84% sensitivity whilst its upper confidence level reaches as much as 92%. These results are comparable to those in [50]. Although results can be hardly generalised, it was observed in the analysis that certain seizure patterns tend to offer better sensitivities than others. In particular, those patients showing low-amplitude, fast activity or rhythmic β waves, obtained some 5% better sensitivity in β band. This was observed both for DDA and PLV PS detection. 48 Delay Difference Analysis for Phase Synchronisation Computation Fig. 2.10 patient #7, SI results from sensitivity test. The horizontal red lines represent the threshold levels and the vertical red boxes represent the pre/post and ictal periods. 2.3 Verification 49 Table 2.2 DDA and PLV results showing the number of true positives (TP), false negatives (FN) and sensitivity for each of the 10 patients shown above. The FPR rate ranges from 0.15 through to 0.75 Patient Metric DDA PLV FPR 0.15 0.3 0.45 0.6 0.75 0.15 0.3 0.45 0.6 0.75 TP 4 5 5 6 6 4 4 5 5 5 #1 FN 3 2 2 1 1 3 3 2 2 2 Sensitivity 0.57 0.71 0.71 0.86 0.86 0.57 0.57 0.71 0.71 0.71 TP 7 8 9 9 10 9 10 10 10 10 #2 FN 3 2 1 1 0 1 0 0 0 0 Sensitivity 0.70 0.80 0.90 0.90 1.00 0.90 1.00 1.00 1.00 1.00 TP 2 3 3 3 3 3 3 3 3 3 #3 FN 1 0 0 0 0 0 0 0 0 0 Sensitivity 0.67 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 TP 2 2 2 2 2 2 2 2 2 2 #4 FN 1 1 1 1 1 1 1 1 1 1 Sensitivity 0.67 0.67 0.67 0.67 0.67 0.67 0.67 0.67 0.67 0.67 TP 1 1 1 1 1 1 1 1 1 1 #5 FN 0 0 0 0 0 0 0 0 0 0 Sensitivity 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 TP 4 4 4 4 4 3 3 3 3 3 #6 FN 0 0 0 0 0 1 1 1 1 1 Sensitivity 1.00 1.00 1.00 1.00 1.00 0.75 0.75 0.75 0.75 0.75 TP 1 6 6 7 10 0 3 7 8 10 #7 FN 11 6 6 5 2 12 9 5 4 2 Sensitivity 0.08 0.50 0.50 0.58 0.83 0.00 0.25 0.58 0.67 0.83 TP 3 3 3 3 3 2 3 3 3 3 #8 FN 4 4 4 4 4 5 4 4 4 4 Sensitivity 0.43 0.43 0.43 0.43 0.43 0.29 0.43 0.43 0.43 0.43 TP 3 3 3 3 4 1 3 3 3 4 #9 FN 1 1 1 1 0 3 1 1 1 0 Sensitivity 0.75 0.75 0.75 0.75 1.00 0.25 0.75 0.75 0.75 1.00 TP 5 5 6 7 7 1 6 7 7 7 #10 FN 5 5 4 3 3 9 4 3 3 3 Sensitivity 0.50 0.50 0.60 0.70 0.70 0.10 0.60 0.70 0.70 0.70 50 Delay Difference Analysis for Phase Synchronisation Computation Threshold 0.15 FP/h PREICTAL POSTICTAL Seizure detection Seizure detection Seizure detection PREICTAL POSTICTAL ICTAL PREICTAL POSTICTAL ICTAL Threshold 0.15 FP/h Threshold 0.15 FP/h DDA (SI) DDA (SI) DDA (SI) 23 23.1 23.2 23.3 23.4 23.5 23.6 23.7 23.8 23.9 0 0.2 0.4 0.6 0.8 1 Hours 23 23.1 23.2 23.3 23.4 23.5 23.6 23.7 23.8 23.9 0 0.2 0.4 0.6 0.8 1 Hours 23 23.1 23.2 23.3 23.4 23.5 23.6 23.7 23.8 23.9 0 0.2 0.4 0.6 0.8 1 Hours (c) ICTAL (b) (a) Fig. 2.11 Zoom synchronisation indexes for patient #7 during a clinically annotated seizure and pre/post ictal padding at a false positive rate of 0.15 per hour. (a) Pair TBA1-HRA4, (b) pair TBA1-HRA5, (c) pair TBA1-TBA2. 0.15 0.30 0.45 0.6 0.75 0.4 0.5 0.6 0.7 0.8 0.9 1 False Positives per Hour Average Sensitivity DDA PLV Fig. 2.12 Average Sensitivity against FPR for DDA + SI and PLV + HT including 95% confidence analysis obtained via 1803.2 hours of recordings and 61 annotated seizures. 2.4 Conclusions 51 2.4 Conclusions In this chapter, a hardware friendly algorithm for the detection of phase synchronisation has been introduced. The algorithm achieves a very high absolute sensitivity of approximately 80% for an FPR of 0.75 FP/h. Indeed, it produces an absolute sensitivity of close to 90% at 0.75 FP/h at the upper bound of the 95% confidence levels. The effect of noise has also been assesed which shows for SNR values greater than 30dB the algorithm performs almost perfectly with over 90% correlation to an ideal PLV. The correlation between other proven methods of phase detection has been assessed giving on average a 90% correlation to the PLV. 58 Synchronisation Processor OX<3> OX<2> OX<1> OX<0> OX<3> OX<2> OX<1> OX<0> OX<2> OX<3> OX<1> OX<0> OX<1> OX<3> OX<2> OX<0> OX<0> OX<3> OX<2> OX<1> OX<9> OX<8> OX<7> OX<6> OX<9> OX<8> OX<7> OX<6> OX<8> OX<9> OX<7> OX<6> OX<7> OX<9> OX<8> OX<6> OX<6> OX<9> OX<8> OX<7> OX<5> OX<4> Minimum detection logic 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 1 0 0 01 1 1 1 1 1 1 Minimum detection Example Fig. 3.2 Minima detection logic (left) and example of bubble detector(right). The example shows an example of a real EEG data stream pre filtered into the beta band. Here we can note that the first minimum in the signal is detected as O = 0000011111. As for the second minimum it is also detected due to the relaxed outlier condition such that O= 0100011111 3.2 Design 59 In the case that a minimum is detected, M1 activates a 10-b array of reading logic gates which are implemented as standard 2 input AND gates and the current value in the up counter Co (9 : 0) is passed into a 10-b latch and stored. simultaneously, since a new minimum has been detected the counter is reset and it once again begins its count towards the next minimum in the signal. 3.2.3 Calculation The architecture of the calculation layer starts with a signed 10-b subtracter with one-bit overflow. This is used to calculate the difference △Tn = |T1 i−T2 j| , between the two signals detected minima. The subtracter runs continuously subtracting values every cycle. However, since the algorithm simplifies the calculations to a minima detection basis only, rather than a sample by sample basis, a small piece of extra logic is needed. This logic comes in the form of two 1 bit latches and a logic AND gate. Once T1 n has been detected in S 1the logic value of M 1gets latched. This indicates that S 1, is ready for the calculation stage, Similarly, M2, is latched if T2 nis detected in S2. If M 1 = M 2 = 1, then a signal is sent to read logic block situated at the output of the previous latching stage. This allows the values from both of the counter modules to be subtracted. Simultaneously the small piece of logic resets the latches holding M 1 and M 2such that the output is now of logic state 0. In terms of the read logic if they are inactive then the output of the read logic is a constant logic 0, for all 10-b. Therefore, although the subtracter is in continuous operation, when no relevant data is ready for calculation it will simply subtract 0 − 0. Since non-switching is necessary for this specific case, the dynamic power consumption of the subtracter during non-crucial data is almost null and only consumes static power. The initial subtraction can create two different outputs at Sub (9 : 0), firstly, if L 1(9 : 0) > L 2(9 : 0), then the result will be a positive unsigned 10-b result. However, in the reverse case that L (9 : 0) < L 2(9 : 0), a 2’s complement negative result is produced. In this case, the MSB (the eleventh bit of the subtracter(S)), is set to 1 indicating the result is negative. Since the algorithm desires |△Tn| , the 2’s complement result needs converting back to standard unsigned notation. To achieve this the design incorporates an array of 10 bit XOR gates which perform the logic conversion to 1’s complement and a full adder adds 1 to Sub(9 : 0) producing Abs(9 : 0). 60 Synchronisation Processor The next stage requires the accumulation of errors for the value of Sdda ( k ), This is achieved using a 10-b add and accumulate block, where the window size K , consists of a variable up counter with a maximum size of 2048 samples. When the window counter is less than its maximum value. Adderout (9 : 0), consistently calculates Abs (9 : 0) + feed (9 : 0), However, once again, since Abs (9 : 0), is only non zero when relevant data is available the resulting sum is 0 + feed (9 : 0), which requires no switching. Once the counter reaches its maximum value. A pulse is generated which resets the feedback register in the accumulator creating a zero-sum, at which point read logic at the output of the accumulator captures Sdda ( k ) = ϕ ( i )and is passed to the post-processing stage for smoothing. This happens every K samples and avoids irrelevant data being passed into the next stage. 3.2.4 Filter Due to fast fluctuations of ϕ ( i )at the processors output, the detection of high and low values of ϕ ( i ), becomes troublesome. To counteract this problem an unconventional exponential filter was adopted, which is based on an autoregressive IIR filter. Mainly used for time series, the filter acts as a low pass filter and provides a substantial amount smoothing with very little cost in terms of hardware. the basic equation for the exponential filter is given by: y(i) = α·ϕ(i) + (1 −α)·y(i−1) (3.1) Where, y ( i ), is the new filtered value, ϕ ( i ), is the new incoming value (in this case new synchronisation value Sdda ( k )), y ( i− 1), is the previously filtered value, and α , is a weighting for the function, also known as a smoothing factor. At this point, it is clear to see that this function produces an approximation of what the next value in the sequence should be, based solely on the current synchronisation value and the previously filtered value. This is why only currently relevant data can be passed into the filter, otherwise, its predictions will be erratic. The smoothing factor, can be set such that 0 ≤α≤ 1, where a α = 1 provides no smoothing, hence y(i) = ϕ(i), and α= 0, provides maximum smoothing. 3.2 Design 61 At first glance, the terms α·ϕ ( i )and (1 −α ) ·y ( i− 1), may seem complicated leading to floating point multiplication. However, re-writing the formula we get the equation: y(i) = y(i−1) −y(i−1) ·α+ϕ(i)·α(3.2) By choosing, α , such that it is always a power of 2, we can replace α , with 2 pos such that: y(i) = y(i−1) −y(i−1) ·2pos +ϕ(i)·2pos (3.3) Where pos , is an integer value. Since a bit shift right, results in a division, a drastic hardware reduction can be made. This approach limits our range of smoothing, nevertheless, for the generalised purpose of smoothing ϕ ( i ), it is sufficient. As an example, if pos is set to 1, then that would correspond to: y(i) = y(i−1) −y(i−1) 2+ϕ(i) 2(3.4) It is also important to note that since the filter relies heavily on its previous value, if y(i), is initially set to one, a certain time frame must pass before stabilisation. Fig.3.3, shows the basic architecture of this block. Firstly, a 10-b shift right register, which is externally adjustable to provide the smoothing factor and a 10-b full subtracter are used to calculate y ( i− 1) −y ( i− 1) · 2 pos . Secondly, a 10-b shift right register is used to calculate ϕ ( i ) ·α , which is added to y ( i− 1) −y ( i− 1) ·α , via a full adder to produce the final value of y(i). 62 Synchronisation Processor φ (i) - + x xShift right register ( 1)yi ()yi ( 1) 2pos yi ( ) 2pos i   ()i  Shift right register -3dB cutoff Fig. 3.3 Block diagram of the implemented exponential filter with variable smoothing capabilities (left). Filters magnitude response (right) 3.2.5 Detect The detection block was designed for the detection of large changes in synchronisation. Namely, this would prove useful for such neurological disorders as epilepsy. In which large changes in synchronisation occur during different periods before, during and after the event. it is important, to stress the fact that since Sdda ( k ), is now effectively equal to PK−1 n=0 △Tn , then a large drop in synchronisation is equivalent to better synchrony, whilst the opposite is true increases in Sdda ( k ). Therefore if ϕ ( i )< Thh , then an alarm is triggered indicating there has been a significant increase in synchrony, contrarily if ϕ ( i ) > Thl , an active high alarm is triggered indicating a significant reduction in synchrony. This is performed by 2 10 bit comparators. Nonetheless, due to the inherent noisy output, 2 small counters of 4 bits are also placed in the detect block. These counters place the constraint that Thl and Thh , must be active for a set number of consecutive outputs. This diminishes the possibility of false positive detections. This block relies heavily on the previous filtering stage such that insufficient smoothing can increase the number of false detections. Thl , Thh and the set number of consecutive outputs are all injected into the device from an outside source and are not set dynamically within. This gives more control for testing. 3.3 Methods 63 3.3 Methods The hardware implementation of the test was built and verified at several keys stages via various platforms. 3.3.1 Design flow FPGA implementation and verification The digital design was firstly realised in the VIVADO design environment using a VHDL. After successful simulation of the design within the environment, it was later synthesised and downloaded onto a XILINX ARTIX 7 AC701xc7a200tfbg676-2 evaluation board. For ease of testing, a hardware Co-Simulation of the DDA algorithm was ran through the MATLAB DSP SIMULINK environment along with a MATLAB based calculation of the HT and PLV index. Fig.3.7, shows an overview of the top-level RTL compiled schematic from the VHDL code. This consists of 4 main blocks. Firstly, the synchronisation processor as described above (a view of the synchronisation processor blocks can be seen in Fig. 3.8). Secondly, the SPI interface, which used to convert serial input streams to parallel data for the processor. The SPI has 3 main inputs ’datain1’, ’datain2’ and ’loadparameters’. The first two inputs convert neural signals 1 and 2, from a 1-b binary serial stream into a 10-b parallel unsigned binary sample. Since the processor requires 2x10-b unsigned parallel data at the input (1x10-b for each signal), the conversion of ’datain1’ and ’datain2’ are made via a 2x10 shift registers which operates 10 times faster than the core clock. Due to the 10 times increase in speed for the datain SPI interface the ’datain1’ and ’datain2’ samples from the SPI are latched at the output until 1 clock cycle of the main processor is complete. The final input ’loadparameters’, allows for all of the static variables of ’K’, ’Thh’, ’Thl’, ’2 pos ’ and others, to be uploaded via a single binary stream. This allows us to reduce the total number of PADS used in the ASIC implementation, hence, reducing overall area consumption of the final design. The third block, ’Check parameters’, was used as a test strategy in order to verify that all of the uploaded variables were correct. This, block simply re-downloads the data sent onto the chip for verification. 64 Synchronisation Processor Lastly, the PSI interface is the download link for the data generated by the synchronisation processor. Similar to the SPI interface this block is used to reduce area consumption of the device. Nevertheless, this block works in the opposite fashion converting parallel data streams into a single serial data stream. Assuming a clock frequency of approximatley 1024Hz operating frequency for the main processor this output PSI can run at the standard processor frequency. Fig. 3.4, shows an example of how ’loadparameters’ and the data stream ’datain1’ are uploaded in a serial fashion and ultimately re-downloaded and checked via ’Parameter check’ and ’checkdata1’. Fig. 3.5, shows the minima detected by the system for signals S1 and S2. Here we can see that there exists a delay between the actual minimum and the detection pulse (the delay is indicated by the blue box). This delay exists because we need to gather 10 samples of history before an accurate decision can be made. Nevertheless, this delay in no way disturbs the systems ability to calculate phase as the exact same delay is present in both signals. Furthermore, we can see that no minima were detected during the section of signal encompassed in the red boxes. This is due to the outlier logic described previously which correctly dismisses these very small dips for minima. Fig. 3.6, shows an example of the output ϕ ( i )using a filtering of α = 1 and two detection thresholds Thh = 170 and Thl = 590. During the simulation 1 drop in synchrony was detected (i.e when ϕ ( i )> Thl ’dropdetected’ = 1) and one increase in synchrony (i.e when ϕ ( i )< Thh ’risedetected’ = 1) (both of which are highlighted by blue boxes). There do exist false positive at the beginning of the output due to the time needed for the filter to build reliable data. RTL compiler verification In the second phase of the design, the VHDL code was implemented using the RTL compiler from CADENCE (An example of the RTL compiler script used to build the RTL schematic can be seen in Appendix B.4. This compiler maps all of the logic gates from the synthesised FPGA design to a specific technology in this case 0.18 µm . The RTL compiler delivers a complete schematic of the design at the register level. Once obtained the schematic is then uploaded to the CADENCE environment and simulated for error detection. If the design produces errors, certain attributes and optimisations can be enforced in the RTL compiler and re-simulated as per needed. 3.3 Methods 65 Fig. 3.4 VHDL simulation showing the uploading of parameters and data stream S1 along with the download verification ’Parameters check’ Fig. 3.5 VHDL simulation showing the minima detection of a small section of EEG recording for S1 and S2 Fig. 3.6 VHDL simulation showing the results for phi for real EEG data and the detection 66 Synchronisation Processor Fig. 3.7 Synthesised block RTL schematic of the prototype processor designed in VHDL. Where, the labled blocks correspond to figure 3.1. Fig. 3.8 Synthesised block RTL schematic of the prototype processor designed in VHDL. Where, the labled blocks correspond to figure 3.1. 3.3 Methods 67 Minima detection Time stamp Filter/calculation/detect Fig. 3.9 Layout view of the 2-channel 0.5V synchronisation processor from CADENCE CADENCE layout Fig 3.9, shows the circuit layout from the CADENCE design environment. Due to a sub-threshold operating voltage the design was implemented via a manual layout approach. The design was fully tested using extracted layouts including parasitic components. 74 Synchronisation Processor other test points of the ASIC through a pod connector. For two tones at 25 and 30Hz (100k samples each), deviations from the ideal synchronisation index were noticeable for supply voltages below 0.35V. At this point, the index starts to drop and, for 0.2V supply, the system breaks down completely with no observed activity. The reset signal generated in the ASIC after and-ing the latched outputs of the two minima logic blocks (see Fig. 3.1) was also monitored along this transition. This signal goes high every time a transition period is computed. For supply voltages above 0.35V, the train generated by the ASIC is quite regular (only minor variations due to the quantisation of the inputs) in good agreement with the expected response, however, pulses are eventually missed below 0.35V and, at 0.2V supply, they are no longer observable. Electrical simulations show that this behaviour is compatible with data transfer failures from the SPI input ports of the ASIC. This is illustrated in Fig. 3.14 which plots: (i) the experimentally observed degradation of the synchronisation index as the supply voltage decreases according to Smeas/Sid , where Smeas and Sid denote, respectively, the measured and ideal values of the index and, (ii) the simulated SPI register position, normalised to the vector resolution, at which input data fails to be shifted. In both cases, 1 indicates no error while 0 means complete breakdown. It is worth observing the good agreement between both curves. Similar behaviour was obtained for other tones in the βband. At a signal input rate of 1024S/s and 0.5V supply voltage, the power consumption of the chip core is 15nW. Table 3.1 shows the power budget of the chip. Event detection and time stamp calculations by the finite state machine are the two most power demanding tasks of the processor due to the higher switching activity. Table 3.2 shows the normalised performance of the ASIC compared to other designs in Table 1.3. Variables have been normalised to a single channel and, in the case of core area occupation, an additional normalisation has been applied by arbitrarily assuming all chips are fabricated in CMOS 65nm. Accordingly, the original area per channel has been scaled by (65 nm/Wmin ) 2 , where Wmin is the feature size of the fabrication technology. Note that the proposed implementation exhibits the smallest normalised area and it is among the designs with lower energy consumption. 3.5.1 Experimental Setup with Neural Recordings The processor has been also verified and validated using neural recording data available in the European Epilepsy database [ 57 ]. This database contains recordings of over 3.5 Experimental results 75 Table 3.1 Power breakdown of the synchronisation processor based on the blocks in Fig.3.1 Block Average Power (nW) Event detection 4.01 Finite state machine 4.52 Counter 1.46 Indexing 2.64 Filtering 2.86 Table 3.2 Normalised performance summary Norm. area Power Energy (mm2/ch) (µW/ch) (pJ/ch/bit) [49] 0.0073 6.25 108 [50] 0.0050 4.06 NA [51] 0.0043 49.22 NA [48] 0.0184 0.23 0.36 [52] 0.0223 51 3.75k [53] 0.0047 0.15 63.5 [55] 0.5320 57.3 14k This work 0.0033 0.0075 0.73 275 patients including 225 scalp recordings and 50 intracranial recordings. As very long recording sessions, often lasting for more than 6h, are included in the database, a specific experimental set-up has been devised. As shown in Fig. 3.15, it uses an mbed LPC1768 MCU based on a 32-bit ARM ® Cortex ™ -M3 core from NXP Semiconductors. The test platform has been expanded with an application board which includes an RJ45 Ethernet connector and a USB-A Host/Device port (in addition to the built-in USB drag’n’drop FLASH programmer). The MCU defines settings and transmits driving stimuli to the DDA chip through the available GPIO ports and other I/O interfaces. Neural input signals are stored as CSV files in a 32 GB Fat-32 formatted stick memory plugged into the USB-A port of the test platform. These files are sequentially read by the MCU and serially transmitted to the SPI ports of the ASIC. Long testing sessions have been run uninterruptedly using this setup with no need to split the database recordings into smaller blocks (this would clear the internal states of the 76 Synchronisation Processor MBED LPC17688 MICROCONTROLLER CUSTOM LEVEL SHIFTING BOARD ASIC 32GB FAT 32 FORMAT USB DRIVE 5V USB 3.0 POWER SUPPLY AND SERIAL DATA ACQUISITION CABLE 0.2V – 1.2V DATA INPUT REGULATORS 5V SUPPLY DATA OUTPUT FROM ASIC DATA INPUT FROM MBED MBED APPLICATION EXTENSION BOARD Fig. 3.15 Test setup showing the ASIC mounted on a dedicated PCB, mbed LPC1768 micro controller and dedicated level shifters. ASIC). The synchronisation index generated by the chip is transferred back to the MCU and steered through the USB port to a host computer where data are further analysed with MATLAB®. The C++ code to send data from the MBED device to the ASIC and retrieve data from the ASIC can be seen in appendix B.5. Neural recording signals from the database have to be pre-processed prior to being used as stimuli in the described platform. This involves the following steps: (i) bandpassfiltering for selecting the frequency band of interest (fourth-order Butterworth filters have been used), (ii) resampling for matching the frequency rates of the recordings and the ASIC (only if needed), (iii) adding noise for a given input-referred SNR value, (iv) digitising the samples to the resolution of the ASIC (10-b), and (v) casting data in CSV format files for storing in the stick memory. As mentioned, files as large as 10GB have been handled without any observed problem. 3.5 Experimental results 77 0 10 20 30 40 50 60 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Time (mins) SI and PLV values ASIC (SI) SIM (PLV) Seizure onset Fig. 3.16 Experimental synchronisation index and simulated PLV value from recordings at positions M5 and M8 of a 1 × 8EcoG strip implanted in patient 1. Signals were pre-filtered in the β -band. The input rate of the chip was adjusted to the sample rate of the recordings, i.e., 1024S/s. The Sdda curve (labeled "ASIC SI") was obtained at a sample rate of 1S/s. 3.5.2 Epilepsy Detection Results Different blocks of intracranial data from the referred database were used for assessing the epilepsy detection capabilities of the DDA chip. These blocks were selected for their high confirmed seizure rate. No other consideration was taken into account. For each block, signals were first pre-filtered in MATLAB ® into conventional neuro-physiological bands. After a preliminary computer analysis of the records, the frequency band with the largest synchronisation index activity was finally selected. For each experiment, the set of valid values were plotted alongside with the theoretical PLV value obtained through the Hilbert Transform method. For illustration purposes, Fig. 3.16 shows the synchronisation results obtained from two channels of a recording block measured in a 48 years old, female patient. In this block, one low amplitude fast activity (lafa) seizure was annotated (dotted line). Before driving the ASIC, Gaussian noise were added to the signals for a spot SNR of 30dB at 30Hz. The plot shows the experimental synchronisation index (labeled as "ASIC SI") and the calculated PLV value (labeled as "SIM PLV"). The DDA algorithm detects all major changes in synchrony and, indeed, follows a similar trend as the PLV algorithm. 78 Synchronisation Processor This is verified by the high correlation measure of approximately 91% (in line with Fig. 2.9) . Fig. 3.16 also shows that the seizure manifests with a remarkable increase in the synchronisation measures. This feature is in fact on the basis of the detection mechanism presented in [ 52 ] and [ 50 ] in which seizures are detected by applying thresholds on the PLV indicators calculated between pairs of channels. In some cases, as in Fig. 3.16, seizures are preceded by a sudden drop in synchronisation [67], however, we have not verified this feature in all the analysed blocks. 3.5.3 Functional Connectivity Results The proposed DDA prototype has been also used for estimating functional brain connectivity. As in the previous section, chip results were later compared with computer simulations of the Hilbert Transform approach with PLV indexing. Recordings from patient #7 have been considered. This patient suffers from simple and complex partial seizures originated from temporal lobe. Among the recordings, a block including one rhythmic β wave episode has been analyzed. This has been the only aspect taken into account for the selection. To assess the functional connectivity strength between neural assemblies, a mapping approach has been followed in which synchronisation indexes for every possible electrode combination at a given observation window were arranged into a symmetric matrix. A 1h long recording block including inter-ictal, pre-ictal, ictal and post-ictal periods was examined. Signals were band-pass filtered in the βband. Both EEG and ECoG recordings are available from patient #7, however, only the results from EEG analysis are herein presented. The EEG recordings consisted of 16 EEG captures selected from a standard 10:20 system. The resulting 120 pairwise combinations were arranged in a 16 × 16 matrix. For easy visualisation, index values were colour-coded between yellow (index close to 1) and dark blue (no coupling). Figure 3.17 shows the connectivity maps obtained through the Hilbert Transform (first column) and the DDA chip (second column) for observation windows in the four mentioned time segments. Both sets are clearly correlated as verified by the absolute differences between the two methods (third column). This is particularly noticeable during the ictal stage (first row), where the difference across all neural combinations is almost null. This shows that for high values of synchronisation, DDA retains the detectability of high synchronisation changes when compared to other more 3.5 Experimental results 79 complex algorithms. For smaller synchrony levels, as occurs in the pre-ictal, post-ictal and inter-ictal periods, the variations between both methods become a little more pronounced. That said, the biggest difference can be seen in the post-ictal period (third row) between electrodes T3-F8 at approximately 10%. Similar conclusions were drawn from the analysis of ECoG maps on the same recording block (not shown). Good correlation between both approaches was also observed. Nevertheless, it was noticed an overall increase in synchrony for all the observation windows. This is consistent with the closer proximity of ECoG electrodes compared to surface electrodes. In any case, the ictal-period was still clearly identifiable. Next, we analysed the graph theoretical metrics as described in Appendix A. Using non-directional, un-weighted binary matrices were the thresholds were set for each individual partition. To maintain consistency each threshold was set as 25% less than the maximum swing between the maximum and minimum respective values, such that: Thresholdpartition =max −max−min 4. In order to calculate the small world metric we introduce a comparison between the inter-ictal and ictal periods. Notably, we use the graph metrics from the inter-ictal period as a random reconnection of the ictal period such that λ = Cictal Cinter−ictal and γ=Lictal Linter−ictal . Conversely, for the inter-ictal period λ=Cinter−ictal Cictal and γ=Linter−ictal Lictal . Fig.3.18, shows the various average graph theoretical metrics over the ictal and inter-ictal partitions. As expected the average nodal degree and clustering increases substantially during the ictal period. Conversely, the average path length drops over the ictal period. This indicates an overall more heavily synchronised and connected network, as can be seen in the ictal connectivity graph. Fig. 3.18, also shows that the ictal partition tends towards a higher small world connectivity than the inter-ictal period. Nevertheless, this highly connected ictal period as compared to the inter-ictal period indicates a mass synchrony between multiple nodes. In a normal brain network, however, these characteristics can be considered undesirable as neural clusters should operate with a certain degree of independence and not oscillate at the same frequencies. Further analysis, was made for patient #7 in terms of the variance in nodal degrees and clustering coefficients. Using, ’topoplots’ in the form of colour graded contour plots projected onto a 2D representation of the brain. Fig. 3.19, show the variance of these metrics for patient #7. Over the ictal period, there exists a high nodal degree towards the right temporal lobe which would indicate higher connectivity from other nodes to this area. This is indeed consistent with the 80 Synchronisation Processor clinical prognosis for this patients type of epilepsy. Nevertheless, The clustering during the ictal period tends to shift to the left hemisphere, indicating that although the right temporal lobe is heavily connected during this period the connections tend to be longer range and less clustered than the left hemisphere, Nevertheless they remain high at approximately 0.6. Over the next three periods, noticeable reductions in nodal degree in the right temporal lobe can be seen, indicating that this area actually has a lower connectivity than other regions. In fact, a reduction in both the inter-ictal and pre-ictal periods accounts for an approximate 23% reduction from in the epileptic event time. Interestingly, the highest connectivity during the non-ictal period shifts towards the P4 parietal electrode which is consistent with the correlation colour maps seen in Fig. 3.17. Another interesting, characteristic for this patient is the reduction in clustering during the pre-ictal period, indicating a sudden drop in clustering around the epileptic centre before the ictal event is triggered. Once again the gathered DDA data shown in the plots was compared with that of the PLV and HT method for a comparative visual aid. In all cases, the two algorithms show a strong correlation with all shifts in clustering and nodal degrees trending towards the same pattern. In conclusion, the ictal period proves itself to evolve from a non-hyper synchronous network to a hypersynchronous network. Where the epileptic event originates from one point of origin and incorporates new clusters into the firing process. This is indicated by the short path lengths, meaning that long-range synchrony does not necessarily take place, instead, the synchrony evolves from one cluster to another. On the contrary, during the inter-ictal period, long rage synchrony is more apparent and the average nodal degree is low. 3.5 Experimental results 81 PLV ICTAL fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 DDA ICTAL fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 PLV - DDA ICTAL fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 PLV PRE ICTAL fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 DDA PRE ICTAL fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 PLV - DDA PRE ICTAL fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 PLV POST ICTAL fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 DDA POST ICTAL fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 PLV - DDA POST ICTAL fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 PLV INTER ICTAL fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 DDA INTER ICTAL fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 PLV - DDA INTER ICTAL fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 fp1 fp2 f3 f4 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Fig. 3.17 EEG connectivity maps for patient #7. Columns 1 through 2 are colour maps representing the mean values for both the PLV and DDA for all possible combinations of EEG signals. Column 3, shows the absolute difference between columns 1 and 2. EEG data were organised into a standard 10-20 format. 82 Synchronisation Processor PLV-ICTAL DDA-ICTAL PLV-INTER-ICTAL DDA-INTER-ICTAL 0 2 4 6 8 10 12 Average degree ICTAL/INTER-ICTAL PLV-ICTAL DDA-ICTAL PLV-INTER-ICTAL DDA-INTER-ICTAL 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2Small worldness ICTAL/INTER-ICTAL PLV-ICTAL DDA-ICTAL PLV-INTER-ICTAL DDA-INTER-ICTAL 0 0.5 1 1.5 2 2.5 Average path length ICTAL/INTER-ICTAL PLV-ICTAL DDA-ICTAL PLV-INTER-ICTAL DDA-INTER-ICTAL 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Average clustering ICTAL/INTER-ICTAL ICTAL PLV/DDA fp1 fp2 f3 c2 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 fp1 fp2 f3 c2 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 INTER PLV/DDA fp1 fp2 f3 c2 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 fp1 fp2 f3 c2 c3 c4 p3 p4 o1 o2 f7 f8 t3 t4 t5 t6 Fig. 3.18 Graph theoretical results for patient#7. Row 1, shows the average metrics for both the PLV and DDA approach over both the ictal and inter-ictal periods. Row 2, shows the functional connectivity representation of the brain using a standard 10-20 EEG electrode set-up for both the ictal and inter-ictal partitions. The blue solid line represents the connectivity of the PLV value and the dashed red line represents the connectivity of the DDA approach 3.5 Experimental results 83 Fig. 3.19 Graph theoretical coloured contour maps for patient 442012 during each partition period ictal,inter-ictal,pre-ictal and post-ictal. Where each row shows both the variance of nodal degree and clustering for both the PLV and DDA approaches. 90 Multi-channel synchronisation processor LSBMSBLSBMSB L1L2L3L4L5L6L7L8L9L10 LSBMSB L151L152L153L154L155L156L157L158L159L160 L11L12L13L14L15L16L17L18L19L20 B1B2B3B4B5B6B7B8B9B10B11B12B13B14B15B16B17B18B19B20B151B152B153B154B155B156B157B158B159B160 10 bit unsigned neural data stream 1 sample 110 bit unsigned neural data stream 2 sample 1 10 bit unsigned neural data stream 16 sample 1 10 bit unsigned neural data streams 15 to 3 sample 1 10 bit unsigned neural data streams 15 to 3 sample 1 10 bit unsigned neural data streams 15 to 3 sample 1 Datain Clk2 Enable Pulse 160 Pulse 1 LATCH OLD DATA MASTER CLK NEW DATA Enable Clk2 Clk Fig. 4.3 1 to 160 SPI for the conversion into 16-10bit unsigned digitised neural signals LSBMSBLSBMSB LSBMSB Data in 10 bit correlation values processor 15 Data in 10 bit correlation values processor 15 to 3 Master clk Data in 10 bit correlation values processor 2 Data in 10 bit correlation values processor 1 Dataout Fig. 4.4 160 to 1 PISO converter for the conversion of 16-10bit unsigned correlation results into a single digital stream 4.2 Design 91 4.2.3 Epoch selection Using 16 main electrodes, a combinational total of 120 neural combinations need to be assessed. Normally, this would lead to a large array of 120 synchronisation processors, where each processor is capable of handling two data streams. To combat this, we use a custom epoch based multiplexing of neural data to reduce the array of processors from 120 down to just 15. The ’epoch’, on which the method is based, can be considered as a small time segment window, where 1 epoch lasts for a total window size of x samples. The actual time of the window is dependent on the sampling frequency of the master clock, however, for argument’s sake, let’s assume for now that it is set to a standard sampling rate of 1024S. During an epoch, each processor will receive 1024 samples from any two given neural data streams. Therefore, during 1 epoch the 15 processors can calculate the correlation between 15 neural data combinations. In order to reach the total of 120 combinations, 8 epochs are required. The epoch system, therefore, reduces the number of processors from 120 to 15, at the cost of losing small amounts of data between epoch cycles. In terms of samples, this means that each neural data combination will get updated every 8192 samples or in terms of time every 8 seconds. The tradeoff between data loss and realistic implantation schemes is questionable. Nevertheless, it is clear that as a diagnostic tool for pathological brain states such as Alzheimers, Parkinsons or schizophrenia, the need for constant data monitoring is unnecessary. Mainly due to the constant and stable change or deterioration in the neural network. That’s to say, the instability in the neural network does not change rapidly or fluctuate from a better to a worse state or visa versa. In epilepsy, however, abrupt rapid changes in the neural network are the most common attribute. A paper by et. al Sigmund Jenssen [ 69 ], reports that on average a Tonic seizures last for a median time of 18.5 seconds, simple partial seizures last for 28 seconds, primary generalised tonic-clonic seizures lasted approximately 66 seconds and complex partial seizures lasted for as long as 78 seconds. From these results, it is clear that an 8-second refresh rate for neural synchronisation values is sufficient. In the fastest seizures of 18.5 seconds there would be a minimum of two attempts at detection of the epileptic event and in the slowest seizures, there would be almost 10 attempts at identification. 92 Multi-channel synchronisation processor Processor Epoch 1 Epoch 2 Epoch 3 Epoch 4 Epoch 5 Epoch 6 Epoch 7 Epoch 8 1 1,2 1,3 1,4 1,5 1,6 1,7 1,8 1,9 2 2,3 2,4 2,5 2,6 2,7 2,8 2,9 2,10 3 3,4 3,5 3,6 3,7 3,8 3,9 3,10 3,11 4 4,5 4,6 4,7 4,8 4,9 4,10 4,11 4,12 5 5,6 5,7 5,8 5,9 5,10 5,11 5,12 5,13 6 6,7 6,8 6,9 6,10 6,11 6,12 6,13 6,14 7 7,8 7,9 7,10 7,11 7,12 7,13 7,14 7,15 8 8,9 8,10 8,11 8,12 8,13 8,14 8,15 8,16 9 9,10 9,11 9,12 9,13 9,14 9,15 9,16 10,1 10 10,11 10,12 10,13 10,14 10,15 10,16 11,1 11,2 11 11,12 11,13 11,14 11,15 11,16 12,1 12,2 12,3 12 12,13 12,14 12,15 12,16 13,1 13,2 13,3 13,4 13 13,14 13,15 13,16 14,1 14,2 14,3 14,4 14,5 14 14,15 14,16 15,1 15,2 15,3 15,4 15,5 15,6 15 15,16 16,1 16,2 16,3 16,4 16,5 16,6 16,7 Table 4.1 Table showing which neural data streams are connected to specific processor during each of the 8 epochs of time Table 4.1, shows a grid of the neural signals being computed for each processor during any given epoch. As an example processor 1 (P1), calculates the synchrony between neural data channels, [1-2,1-3,1-4,1-5,1-6,1-7,1-8,1-9] at epochs 1 through 8 respectively. Whilst processor 15 (P15), calculates the synchrony between [15-16,16- 1,16-2,16-3,16-4,16-5,16-6,16-7]. 4.2.4 Window, epoch and memory addressing logic Structurally speaking, the control logic for the entire system is simple and consists of basic logic gates and counters. Firstly, as can be seen in Fig. 4.5 (Window counter), the window logic consists of a simple 10-b up counter with self-reset. This counts from 0 to 1024, when the maximum of 1024 is reached a single clock cycle positive pulse ’K’ is created. From observations in the previous design discussed in chapter 3, a fixed window of 1024 was chosen. This allows for a reduction in hardware and power consumption by omitting the variable shifting aspect. The epoch selection of the neural signal inputs and addressing of memory elements is driven by a 3 bit up counter, built from D-type flip-flops as seen in Fig. 4.5 (Epoch selection and memory addressing logic). Where the output of each flip-flop is hard-wired to 8 distinct 3 input AND gates in combinations which create an up sequence ranging from 1 through 8 (a1 through a8). The counter has an enable line connected to ’K’, 4.2 Design 93 hence the counter only counts up once ’K’ is introduced. Contrarily, when ’K’ is low, the counter is paused in its current state meaning that the count is kept constant until the next ’K’. This combination allows us to produce logic high signals of exactly 1024 samples long. Using the hard-wired AND gates at the output has the advantage of hardware reduction over ring counters as a reduction of five D-type flip-flops can be introduced. The long pulses of a1 through a8 are used as a reference for the current epoch as well as address memory locations, which will be described later. Fig. 4.5 shows a basic timing diagram, where the 1024 sample pulses of a1 through a8 rotate continuously. In the zoomed view of a1, we can see how a single cycled pulse of ’K’ activates the epoch a1. D E Q D E Q D E Q c1 c2 c3 Epoch selection and memory addressing logic D Q D Q D Q D Q D Q D Q D Q D Q D Q D Q K Clk Window Counter 1024 Samples a1 a2 a3 a4 a5 a6 a7 a8 NOT(c1) NOT(c2) NOT(c3) c1 NOT(c2) NOT(c3) NOT(c1) c2 NOT(c3) c1 c2 NOT(c3) c1 NOT(c2) c3 NOT(c1) c2 c3 NOT(c1) c2 c3 c1 c2 c3 a1 a2 a3 a4 a5 a6 a7 a8 K Clk a1 1024 Clk cycles Zoom Epoch control selectors Fig. 4.5 Dedicated control logic for the window length, epoch selection and memory addressing. 4.2.5 Flow From Fig. 4.6, we can see an overview of the hardware implemented for processor 1 (P1) (all 15 processors are identical). Firstly, during any given epoch, 2 of the 9 possible channels of parallel neural data are multiplexed into the synchronisation processor as per, table 4.1. During epoch 1 signals S1 and S2 are evaluated, the phase error between these two signals is then continuously monitored and updated until the window has finished (K = 1) producing a phase value ϕepoch1 . Simultaneously, the Epoch has increased to 2, hence new data streams S1 and S3 are feed into the 94 Multi-channel synchronisation processor synchronisation processor for assessment. Working as a pipelined system ϕepoch2 , is now being processed whilst ϕepoch1 continues to be post-processed. In general ϕepoch(n) , is always being processed whilst ϕepoch(n−1), is always being post-processed. ϕepoch1 , is then passed into a dedicated custom low pass filter and stored in a custom dedicated memory for future retrieval. After filtering, ϕepoch1 , is then passed to a dedicated training and calculation module which uses the values of synchrony to produce patient-specific epileptic detection and graph theoretical thresholds for online epileptic detection and offline graph analysis respectively. The dedicated thresholds are then used as a comparison in a detection block which compares new synchrony values to determine whether a rise or fall in a patient’s specific synchrony values has occurred. Custom Tristate multiplexing Phase synchronisation processor 10 bit custom exponential filter Memory bank 1 8x10 Memory bank 2 8x10 Memory bank 3 8x10 Processor (P) Training and threshold calculations Training block Threshold calculations block Detection block 10 bit upcounter 3 bit counter with AND logic out Bypass 0.5v Train Clk Logic unit Select K Bypass Select memory location S1 S2 N1 N2 Clk 0.5v Clk Epoch 1:8 Slection logic Fig. 4.6 Overview of the 16-channel processors control and flow 4.2.6 Neural multiplexing To accommodate all 120 neural data combinations over 8 epochs of time, each processor must manage the phase calculations of 9 independent neural combinations. In order to minimise the large area occupation enforced by standard multiplexing schemes, our design incorporates a non-conventional hard-wiring via tristate buffers. 4.2 Design 95 Fig. 4.7 shows the basic hard-wiring technique. Where select processor input channels are connected via a multi-dimensional 8x10 array of tristate buffers. Each of the 8 rows of 10 tristate buffers connects to an individual neural data stream of 10-b, which all have a dedicated enable line connected to a specific epoch. This allows us to connect 8, 10-b neural data streams to a single bus. When a current epoch is active (a1 through a8, see section 4.2.4), the output of that given tristate array is set to low impedance, Hence, the data can pass onto the main line and feed into the synchronisation processor unimpeded. Contrarily, when an epoch is not active the output is set to a state of high impedance (Z), meaning that the neural data is blocked from passing onto the main line. From Table. 4.1, for synchronisation processors 1 through 8, the first channel is always hard connected to neural data streams 1 through 8, respectively. Whilst channel 2, is connected to each processor’s dedicated tristate array. As an example in Fig. 4.7, channel 1 (Ch1) of processor 1, is hard-wired to the neural data stream 1 (D1), whilst channel 2 (Ch2), is connected to the tristate array, the first row in the array is connected to neural data stream 2 (D2), the second row to neural data stream 3 (D3) and so on until the 8th row, which is connected to neural data stream 9 (D9). Synchronisation processors 9 through 15 incorporate an extra 2-row tristate array at channel 1. As per table 4.5, eventually the input at channel 1 will have to change, however, this is not a uniform process so a small piece of distinct logic is added to control the switching of these tristate rows. As an example during epoch 8 of processor 9, the Ch1 input will need to switch from D9 to D10. To achieve this a simple inverter is added to the enable line of the first row, such that if epoch 8 is not active then neural data D9 can flow into Ch1, else D10 is allowed to pass unimpeded. The additional logic for each of the processors 10 through 15 can be seen in Fig.4.7. The timing diagram shows the currently active rows for processors 1 and 9 during epochs 1 (highlighted in green) and 8 (highlighted in blue). 96 Multi-channel synchronisation processor a8 ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB 10 bits 10 bits a1 a2 a8 D2 D3 D9 10 bits Ch2 φ1 D1 Ch1Ch1 φ9 Ch2 ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB D9 D10 Example processor 1 connection Example processor 9 connection ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB ENB a1 a2 a8 D10 D3 D1 a1 a2 a3 a4 a5 a6 a7 a8 a1 a2 a3 a4 a5 a6 a1 a2 a3 a4 a5 a1 a2 a3 a4 a1 a2 a3 a1 a2 TRISTATE ARRAY2 TRISTATE ARRAY1 D10 D11 TRISTATE ARRAY2 TRISTATE ARRAY1 D11 D12 TRISTATE ARRAY2 TRISTATE ARRAY1 D12 D13 TRISTATE ARRAY2 TRISTATE ARRAY1 D13 D14 TRISTATE ARRAY2 TRISTATE ARRAY1 D14 D15 TRISTATE ARRAY2 TRISTATE ARRAY1 D15 D16 a1 Channel 1 input logic processor 10 Channel 1 input logic processor 11 Channel 1 input logic processor 12 Channel 1 input logic processor 13 Channel 1 input logic processor 14 Channel 1 input logic processor 15 Fig. 4.7 Example of the neural input multiplexing scheme and dedicated control logic. The top figures show how tristate buffer arrays are used to control 8 separate neural streams to one processor. The bottom figures shows the extra control logic needed in order to control the first input channel 4.2 Design 97 4.2.7 Post processing The epoch based method chosen means that temporary variables and data need to be stored for future access. The exponential filtering is an example of this, where the original low pass filtered ϕ(i), from Fig. 3.3 of section 3.2.4 in chapter 3: Now needs to store the previous value of y ( i )for each epoch, until a full rotation of the signals has been complete. The equation can then be expressed as: y(i)(P)=y(i−1)(P) epochx −y(i−1)(P) epochx ·2pos +ϕ(i)(P) epochx ·2pos Were y ( i− 1) epochx , represents the previously stored filtered value for that current epoch and P, represents the processor. From herein, we will omit the term epochx and (P) for clarity of reading, hence we assume that all further references are indicative of all 15 processors and epochs Since each processor manages multiple neural data streams, a dedicated 8x10 memory, which is illustrated in Fig. 4.8 (a), is used to store the previous values of y ( i− 1), until the epoch selection module makes a full round and returns back to the chosen epoch. The memory itself consists of a 10x8 array of D-type latches all interconnected at the output via tristate buffers (T1 through T8) onto a single 10-b bus. This allows for both the read and write operations to be carried out simultaneously. The latching of the data to the correct address requires a single cycle positive pulse in this case ’a1Lf’ through ’a8Lf’. Which are activated by ’anding’ the current epoch with the single cycled window pulse of ’K’ as illustrated in Fig. 4.8 (b). Once the new data is latched to the correct address, the output of the memory is read constantly throughout the life of each individual epoch via epoch selection signals a1 through a8. Fig. 4.8 (c), shows an overview of the hardware requirements for such a filter. Which is similar to the hardware described in from Fig. 3.3 of section 3.2.4 in chapter 3. Nevertheless, by assuming a constant value of pos , in this case, pos = 2 a reduction of two 10-b shift registers can be omitted by hard wiring the shift into the dedicated hardware. Firstly, a 10-b subtracter calculates y ( i− 1) −y ( i− 1) · 2 pos (collected from memory), where, y ( i− 1) · 2 pos is hard-wired to the subtracter by taking the most significant bits from 10 down to 2 and grounding out the first 2 most significant bits and concatenating bits 10 to 2 in position 8 down to 1 98 Multi-channel synchronisation processor The result from the subtracter, is then added to ϕ ( i ) · 2 pos (once again hard-wired), using a 10-b adder to create the final result y ( i ), which is simultaneously passed on to the next stage and stored into memory as y(i−1). D E QD E QD E QD E Q D E QD E QD E QD E Q D E QD E QD E QD E Q 1 2 3 10 1 2 8 10x8 Memory T1 T2 T8 Tristate arrays 10 bit 10 bit 10 bit Common bus 10 bit a1Lf a2Lf a8Lf - ( 1)epochx yi ( 1) 10:2 epochx yi   + () epochx yi a1 a2 a8 Exponential filter 10x8 filter memory x10 a1 a2 a3 a4 a5 a6 a7 a8K a1Lf a2Lf a3Lf a4Lf a5Lf a6Lf a7Lf a8Lf fIlter memory logic ( ) 10:2 epochx i   ab c Fig. 4.8 Adapted exponential filter from processor 1, using a hard-wired smoothing factor and memory array 4.2.8 Training and Thresholds overview It is well known that there are a wide variety of different forms of epilepsy which can produce different varying quantities of synchrony between different regions. Additionally, each and every patient will produce varying amounts of synchrony due to the uniqueness of every brain. These variations make it difficult to set reliable static thresholds for epileptic detection as seen in Fig. 3.6 of section 3.2.5 in chapter 3, where the thresholds were set based on pre simulation. Therefore, in this design, a patient-specific epileptic threshold calculator was implemented. To begin, an initial training phase is set, such that 2 separate arrays are built, 1 containing the maximum filtered synchronisation values and the other the minimum values, from each processor at each epoch (i.e the maximum and minimum filtered values for each neural data combination). In an ideal situation, the training period would be performed during a patient’s epileptic episodes by clinical specialists. The arrays, give an accurate starting point for 4.2 Design 99 the patient’s pathological synchronisation state. Once training is complete, the data accumulated is then used to set a patient specific low and high threshold voltage such that: Thh =max(y(i)) −∆(4.1) Thl =min(y(i)) + ∆ (4.2) Where, p , is the processor in question and ∆is a margin calculated as the voltage swing between the maximum and minimum values stored in the array divided by a constant such that:: ∆ = "max(y(i)) −min(y(i))#·1 2n(4.3) A threshold low (Thl) is set to produce alerts when the patient synchrony levels are increasing towards the maximum recorded value in the array, warning the patient to a possible upcoming episode. The threshold high (Thh) is set as a diagnostic tool to track a possible upcoming epileptic seizure. The time for an accurate array to be built depends greatly on the patient and the amount of successfully induced episodes recorded. Take for example a patient who, did not reach a positive seizure state during the training period, this patients Thh and Thl thresholds will not give an accurate gauge of his/her hypersensitive network detection. This dynamic approach to thresholds gives rise to several advantages over fixed threshold based architectures. Firstly, since the thresholds are based from a patient’s specific source, the thresholds are dynamically adjusted to account for a synchronisation offset which may occur from patient to patient. Secondly, In the fixed threshold hold based algorithmic designs a normalization between 0 and 1 is usually applied. This makes setting the thresholds easier, however, the thresholds are still susceptible to large error margins due to offsets of patient-specific synchrony. The dynamic thresholds also provide an all-important feature towards the graph theoretical analysis. The complexity of threshold setting for the binary adjacent matrices (see appendix A for definition of graph theory), usually leads researches to handle multiple data sets at varying thresholds. 106 Multi-channel synchronisation processor Fig. 4.11 Array of 15 synchronisation processors synthesised using VHDL code SYNC_BLOCK_1gen[0].sync8 SYNCBLOCK CALC TRAININ a1 a2 a3 a4 a5 a6 a7 a8 ce clk flagin rst rst2 alpha[3:0] data_in1[9:0] data_in2[9:0] offset[2:0] sumout[9:0] Calcs Calculations clk flag1 flag3 flag5 ready1 ready2 rst Inst_Calcflags Calcflags clk flag1 flag3 flag5 rstorage countin1[9:0] countin2[9:0] countsubed[9:0] offset[2:0] Inst_Indexing Indexing clk flag flagin rst cal[9:0] sum[9:0] Inst_Smooth_Detect Smooth_Detect clk flag flagwrite rst alpha[3:0] sumin[9:0] suminfrommem[9:0] sumout[9:0] Mainone Main clk rst valley counto[9:0] data_in[9:0] Maintwo Main clk rst valley counto[9:0] data_in[9:0] expmem1 exp_mem a1 a2 a3 a4 a5 a6 a7 a8 clk flagwrite rst memin[9:0] memout[9:0] train1 TRAIN CALC K TRAININ a1 a2 a3 a4 a5 a6 a7 a8 clk rst high[9:0] low[9:0] SYNC_BLOCK_1gen[1].sync8 CALC TRAININ alpha[3:0] SYNCBLOCK rst rst2 Fig. 4.12 Single synchronisation processor including all of the blocks from Fig. 3.8 in chapter 3 and an additional train and calculation block 4.3 FPGA Implementation 107 Fig. 4.13 Simulation of the basic control logic and threshold memory storage Fig. 4.14 Simulation of the basic calculation stage 108 Multi-channel synchronisation processor Signal (P1) Static or chirp frequency range 1 chirp 26-51Hz 2 static 37Hz 3 chirp 27-52Hz 4 static 37Hz 5 chirp 28-53Hz 6 static 37Hz 7 chirp 29-54Hz 8 static 37Hz 9 chirp 30-55Hz Table 4.3 Table showing the signals for processor 1 during VHDL simulation of 16- channel design Resource Utilisation LUT’s as logic 4035(3%) Register as Flip Flop 4418(2%) Input/output 25(9%) Clocking 2(6%) Table 4.4 Table showing all design resource allocation for the DDA hardware implementation. frequencies converge at an earlier stage as the processors get higher due to the odd chirp signals increasing in frequency the higher the signal number. As another example the first result for epoch 1 is 292 which corresponds to the synchronisation results for S1(26Hz) and S2(37Hz), since their frequencies are distant from one another we get a large value at the output. However, during epoch 2 signals 1(26Hz) and 3(27Hz) are compared hence producing a small value of 43. Table 4.4, shows the resources used in the FPGA implementation of the 16-channel synchronisation processor. 4.3 FPGA Implementation 109 Fig. 4.15 Output results from VHDL simulation of 16-channel processor 110 Multi-channel synchronisation processor 4.4 Conclusions In this chapter a 16-channel VLSI phase synchronisation processor was implemented The design only occupies 3.64 mm2 . Due to a change in the MPW policy of the selected foundry, AMS AG, regarding its 0.18 µ m CMOS process, the processor was received 13 months after the tape-out and finally it was not tested. Nevertheless, the design has been fully verified in an FPGA environment. Chapter 5 Conclusions 112 Conclusions 5.1 Conclusions In this thesis a new algorithm and synchronisation index for the detection of phase synchronisation between neural signals has been introduced. The algorithm shows a high correlation of approximately 90% when compared to more complex transformbased methods with only a slight decrease on noise performance. In addition, an exhaustive sensitivity analysis with human neural recordings shows that the algorithm produces average values of over 80% for FPR values of 0.6 and 0.75. Two VLSI phase synchronisation processors have been introduced. The first circuit, implemented with a custom 0.5V subthreshold digital library in a 0.18 µ m standard CMOS process, allows for the detection of phase synchronisation between two independent neural signals. The circuit occupies just 0.05 mm 2 , including SPI AND PSO interfaces. The design consumes on average 15nW of power at 2KHz operating frequency which has been proven to be much lower than current state-of- the-art. In fact, when normalised to a 65nW technology, the design consumes just 0.0075µW per channel. That is a remarkable 20 times less power. The functionality of the prototype has been extensively tested with neural recordings from an epilepsy database and its correct operation demonstrated both in seizure detection tasks and functional connectivity computation. The second multi-channel processor, using the same library and technology node as above, is capable of calculating phase synchrony between 16 independent neural signals. The circuit occupies 3.64 mm 2 in total and incorporates novel patient-specific training and on-the-fly seizure detection mechanisms. The design has been implemented in FPGA and shows very small resource allocations of just above 4000 LUTs and 4400 registers. 5.1 Conclusions 113 En esta tesis se ha introducido un nuevo algoritmo e índice de sincronización para la detección de la sincronización de fase entre señales neuronales. Los resultados del algoritmo propuesto mostraron una alta correlación de aproximadamente el 90% en comparación con el método PLV basado en la transformación más complejo. También se demostró que cuando se realizó un análisis de SNR de ruido puntual en el algoritmo, se produjeron valores de correlación con una señal ideal libre de ruido superior a 90% para valores de SNR de 25 a 70 dB. Esto es solo un poco más bajo que el del método más complejo de HT + PLV. Además, un fuerte análisis de sensibilidad absoluta mostró que el algoritmo produjo valores promedio de más del 80% para valores de FPR de 0,6 y 0,75. Increíblemente, se demostró que el algoritmo DDA superó el método más complejo de HT + PLV en todos los valores de FPR. Además, se han introducido dos procesadores de sincronización de fase implementados por VLSI. El primero de ellos permitió detectar la sincronización de fase entre 2 señales neuronales independientes y tenía una frecuencia de operación principal de 2 KHz. El VLSI se implementó utilizando una biblioteca digital por debajo del umbral de 0.5V y una tecnología estándar de 0.18 µm. La implementación de FPGA del diseño mostró que el diseño usaba hasta un 90% menos de recursos de FPGA en comparación con otras 2 implementaciones avanzadas. El VLSI ocupaba solo 0.05 mm 2 que incluía interfaces SPI Y PSO. Además, cuando se normalizó en torno a una tecnología de 65 nm, el diseño ocupó solo 0,0033 mm2 por canal, que es 1,5 veces más pequeño que el diseño más pequeño del estado de la técnica. El diseño consumió un promedio de 15 nW de potencia a 2 KHz de frecuencia de funcionamiento, que se ha demostrado que es mucho más bajo que el estado actual de la técnica; de hecho, cuando se normaliza a una tecnología de 65 nW, el diseño consume solo 0,0075 µ W por canal. Eso es un notable 20 veces menos poder. Además, a través de una prueba de laboratorio, se demostró que el diseño puede operar tan solo a 0.3 voltios antes de que la corrupción de datos comience a infringir el circuito. Los resultados tomados de la implementación VLSI a través de la configuración de prueba de laboratorio mostraron que el ASIC estaba funcionando correctamente, lo que se indica por un aumento repentino en la sincronización justo después del inicio de la crisis. Esto fue respaldado por la simulación del método HT + PLV que muestra un patrón muy similar. 114 Conclusions Además, el análisis estadístico mostró altos aumentos en la sincronía alrededor del evento epiléptico en términos de conectividad funcional y análisis teórico de gráficos. El segundo procesador multicanal fue capaz de calcular la sincronía de fase entre 16 señales neuronales independientes. Una vez más, se implementó utilizando una biblioteca digital por debajo del umbral de 0.5V y una tecnología estándar de 0.18 µ m. El VLSI ocupó 3,64 mm 2 en total e incorpora el entrenamiento específico del paciente y el umbral de cálculo, así como la detección de la mosca. Los umbrales específicos del paciente son de última generación y no se han integrado hasta este punto. El diseño se implementó en FPGA y mostró asignaciones de recursos muy pequeñas de poco más de 4000 LUT y 4400 registros. References [1] Cold Spring Harbor Laboratory. 3D Brain. DNA Learning Center, AXS Studio, Vivid Apps, 2009. URL http://www.brainfacts.org/3d-brain#intro=truehttps: //itunes.apple.com/nl/app/3d-brain/id331399332?mt=8. [2] Wulfram Gerstner. Neuronal dynamics: From single neurons to networks and models of cognition, 2014. URL http://neuronaldynamics.epfl.ch/online/Ch1.S1. html. [3] Khan Academy. Neuron action potentials: The creation of a brain signal (article) | Khan Academy. pages 1–13, 2017. URL https: //www.khanacademy.org/test-prep/mcat/organ-systems/neuron-membran e-potentials/a/neuron-action-potentials-the-creation-of-a-brain-signal. [4] 12.4 The Action Potential – Anatomy and Physiology. URL https://opentextbc.c a/anatomyandphysiology/chapter/12-4-the-action-potential/. [5] Chemical and Electrical Synapses Two Kinds of Synapses. URL http://cbm.msoe .edu/markMyweb/ddtyResources/documents/synapseTypes.pdf. [6] Tomas Ros, Bernard J Baars, Ruth A Lanius, and Patrik Vuilleumier. Tuning pathological brain oscillations with neurofeedback: a systems neuroscience framework. Frontiers in human neuroscience, 8:1008, 2014. ISSN 1662-5161. doi: 10. 3389/fnhum.2014.01008. URL http://www.ncbi.nlm.nih.gov/pubmed/25566028h ttp://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC4270171. [7] R Caton. The Electric Currents of the Brain. Brit Med J, 2(1):278–, 1875. ISSN 09598138. doi: 10.1136/bmj.2.765.257. URL https://www.tandfonline.com/doi/ abs/10.1080/00029238.1970.11080764. [8] THOUGHT AND BRAIN: A CAMBRIDGE EXPERIMENT » 5 Oct 1934 » The Spectator Archive. URL http://archive.spectator.co.uk/article/5th-october-1934/ 10/thought-and-brain-a-cambridge-experiment. [9] Fernando Lopes da Silva. EEG and MEG: Relevance to neuroscience, 12 2013. ISSN 08966273. URL http://www.ncbi.nlm.nih.gov/pubmed/24314724http://link inghub.elsevier.com/retrieve/pii/S0896627313009203. [10] Mikail Rubinov and Olaf Sporns. Complex network measures of brain connectivity: Uses and interpretations. NeuroImage, 52(3):1059–1069, 9 2010. ISSN 1053-8119. 122 References 1–1, 2014. ISBN 978-1-4799-4982-3. doi: 10.1109/ICGCCEE.2014.6922342. URL https://pdfs.semanticscholar.org/23f6/05ebba13b1ac014bbf3e6d737200ae3f 5c15.pdfhttp://ieeexplore.ieee.org/document/6922342/. [67] Florian Mormann, Thomas Kreuz, Ralph G Andrzejak, Peter David, Klaus Lehnertz, and Christian E Elger. Epileptic seizures are preceded by a decrease in synchronization. Epilepsy Research, 53(3):173–185, 2003. ISSN 09201211. doi: 10.1016/S0920-1211(03)00002-0. URL https://pdfs.semanticscholar.org/6506/ 728a434aa128ac9961c0104dd5914c2848f6.pdf. [68] Sparsh Mittal. A Survey of Techniques for Approximate Computing. ACM Computing Surveys, 48(4):1–33, 3 2016. ISSN 03600300. doi: 10.1145/2893356. URL http://dl.acm.org/citation.cfm?doid=2891449.2893356. [69] Sigmund Jenssen, Edward J. Gracely, and Michael R. Sperling. How long do most seizures last? A systematic comparison of seizures recorded in the epilepsy monitoring unit. Epilepsia, 47(9):1499–1503, 9 2006. ISSN 00139580. doi: 10. 1111/j.1528-1167.2006.00622.x. URL http://doi.wiley.com/10.1111/j.1528-1167. 2006.00622.x. [70] Leonhard Euler. Solutio problematis ad geometriam situs pertinentis. Commentarii Academiae Scientiarum Imperialis Petropolitanae, 8:128–140, 1736. [71] Ed Bullmore and Olaf Sporns. Complex brain networks: graph theoretical analysis of structural and functional systems. Nature Publishing Group, 10(3):186–198, 3 2009. ISSN 1471-0048. doi: 10.1038/nrn2575. URL http://www.nature.com/artic les/nrn2575http://www.ncbi.nlm.nih.gov/pubmed/19190637. [72] Hunter C Gits. Relating Connectivity and Graph Analysis to Cognitive Function in Alzheimer’s Disease. pages 45–65, 2016. ISSN 2470-9727. doi: http://dx.doi.o rg/10.3998/mjm.13761231.0001.111. [73] Bernadette C. M. van Wijk, Cornelis J. Stam, and Andreas Daffertshofer. Comparing brain networks of different size and connectivity density using graph theory. PLoS ONE, 5(10):e13701, 10 2010. ISSN 19326203. doi: 10.1371/journal.pone .0013701. URL http://dx.plos.org/10.1371/journal.pone.0013701. [74] Duncan J. Watts and Steven H. Strogatz. Collectivedynamics of ’small-world’ networks. Nature, 393(6684):440–442, 1998. ISSN 00280836. doi: 10.1038/30918. URL http://www.nature.com/nature/journal/v393/n6684/abs/393440a0.html %5Cnhttp://www.nature.com/nature/journal/v393/n6684/pdf/393440a0.pdf% 5Cnhttp://www.nature.com/doifinder/10.1038/30918. Chapter A Graph Theory concepts Graph theory is common mathematical approach which provides interesting metrics for the determination of networks. The applications for graph theory have a strong hold in many disciplines including computer science, biology and sociology to name a few. A historical mathematical problem named Seven bridges of Konigsberg written by Leonhard Euler in 1736, is considered to be the first introduction to graphs and topologies, where Euler proved using graphs that it was impossible to traverse the city of Königsberg exactly once and return to the starting point using seven different interconnected bridges [70]. In neuroscience graph theory is used as a tool to further identify underlying characteristics of the brains structural, functional and effective connectivity. Figure A.1, shows an overview of some of the underlying concepts behind graph theory. Firstly, every object in the network is assigned as a node (identified as a grey dot). Where each node in the network are connected via edges, which are the abstract connections between the nodes. As an example, in a typical EEG/MEG scalp recording set-up the electrodes would be considered the nodes and all the interconnections between the electrodes would be considered the edges of the network. In the 10-20 system this would typically consist of 16 nodes and 120 edges (excluding isolated nodes). In general, graph networks are complex networks which should follow the following logical steps[71, 72]: 1. Define all nodes in the network: As stated above, the definition of the nodes could be in the form of EEG/MEG electrodes or as anatomically defined regions via the interpretation of MRI or DTI data. 124 Graph Theory concepts 2. The second step, is to identify a continuous measure of association between nodes. This is usually done between all combinations of nodes. This continuous measure could be time based correlation, spectral analysis or DTI analysis. 3. The next step is to derive an ’association matrix’, M , which consists of a square NxN, map of all the continuous measures from the previous step. 4. The next step is to derive an ’adjacency matrix’, A , from the association matrix. These can be constructed in several ways depending on the application intended. Nevertheless, the most common A , is a binary matrix, in which a static or dynamically varying threshold is applied the M . In this case anything above the threshold is set to 1, conversley, anything below is set to 0, this produces a symmetrical Amatrix filled with zeros and ones. 5. The final step in the process is to conduct analysis on A , using graph theoretical based metrics to build information maps on the network. Other types of networks do exist beyond the non-weighted non-directional network, such as weighted/directed graphs, binary/directed graphs, however functional connectivity is usually set as binary/undirected graphs. In general there are three main desirable metrics which can prove useful to the identification of structural and functional connectivity in brain networks[73]: 1. Degree The degree k, is a metric which identifies the number of edges connected to a given node such that: ki= N X i Ai,j (A.1) Where A , is the adjacency matrix and i and j are the row and coulomb in the matrix respectively. In the case of binary threshold matrices, this gives a clear indication of how heavily each node in the network is connected with other nodes. More often than not the averaging of the edges per node is calculated such that: ⟨k⟩=1 N N X i ki(A.2) 125 This gives an overall value for the brains connectivity patterns (i,e) the total number of nodes which passed the threshold. This metric is especially useful for time varying analysis and dynamic thresholds as time lapse maps can be built. 2. Path length The average path length L, indicates the minimum number of edges between any two nodes within a network, using global efficiency Eto account for isolated nodes. The path length indicates how well information or other entities in the network can be transferred from one region to another. Given as: E=1 N(N−1) N X i,j,i=j 1 di,j (A.3) L=1 E(A.4) Where, d, is the distance between two nodes. 3. Clustering coefficient The clustering coefficient C, measures the occurrence of clusters in a given network by showing the probability that a nodes neighbours are also connected. This measure is an index of local structure and is used as a resilience factor, i.e if a node is lost (doesn’t pass threshold), what’s the probability that its nearest neighbours are still connected. C=1 N N X i 2ni ki(ki−1) (A.5) Where ni, represents the number of edges between the neighbours of node i. 4. Edge density The edge density m, identifies the fraction of existing edges in a network out of the total number of edges in the network: 126 Graph Theory concepts m=1 N(N−1) N X i,j Ai,j (A.6) figure A.1, shows an example of how graph theory can be applied to a network. In this case the network (shown in the upper left hand corner), comprises of 5 nodes, A,B,C,D and E. To calculate the average nodal degree of the network the mean of the total number of edges per node is taken. In this case, A = 3, as it is connected to 3 other nodes via edges, B = 4, etc... To grade the overall network the average is taken such that: avg(ki) = 1 5·(3+4+2+2+1), giving a mean total of 2.4. The average minimum distance between any two nodes is such that: E = 1 5·4· (3 . 5+ 4 + 3 + 3 + 2 . 5), giving a mean total of 0 . 8, this is known as the global efficiency. The Path length is calculated as L=1 0.8= 1.25. As an example in the case of the distance between nodes A and E, there exists a shortest route of 2 edges, therefore, 1 di,j =1 2= 0.5. For the clustering coefficient node A = 2·2 3·(3−1) = 1 3 , the denominator signifies the amount of edges connected to the node A, whereas the numerator signifies the number of connections made between other nodes connected to these edges. In A the numerator ni= 2, since C is connected to B and D is also connected to B. In the case that isolated nodes are located such as node E, the clustering coefficient is set to 0, because the denominator ki(ki−1), is set to 0, and is hence undefined. From the metrics above it is possible to construct special networks topologies which derive an association to brain connectivity and can be seen in figure A.1. These are,ordered, small-world and random networks. The ordered network is a such that each node is connected to its K nearest neighbours. In general the ordered network has a high C and high L (neighbouring nodes are well connected and the path length is long). The random network is such that all edges are randomly connected to nodes, in this case C tends to be very small L very small (neighbouring nodes are not well connected and the path length is short). 127 The final network, small world, is a form of the ordered network, where some nodes are disconnected and reconnected based on probability, Small world tends to have a high C but a very low L. The brains connectivity tends to show characteristics similar to that of a small world network, where neuronal clusters are well connected via short connections [74]. The small-worldness of a given network H, which comprises of N nodes can be obtained by normalising its C and L with that of a random network such that: λ = CH Crandom and γ = LH Lrandom . If λ > 1and γ≈ 1, the network is said to be characteristic of a small world network. η = λ γ , defines the small worldness but must strictly be greater than 1 [73]. It is however important to stress that the probability factor allows for many networks between p= 0 and p= 1, By increasing the rewiring probability. Node Edge Degree Clustering Path length Ordered Small World Random P = 0P = 1 Nodes Edges 0.66 0.33 1 1 0 A B C DE 0.59 Mean A A B B C C D D E E 3 4 2 2 1 Total 2.4 A B C D E 0.66 1 1 1 0.5 A B C DE 3.5 Mean 1 0.33 1 1 1 4 1 1 1 0.5 0.5 3 1 1 0.5 1 0.5 3 0.5 1 0.5 0.5 0 2.5 A B C D E FROM TO 0.8 1.25 , o CL , o CL , o CL Average Degree Average path length Average clustering Fig. A.1 Graph theory: Overview of some of the basic concepts in graph theoretical analysis. A node is represented as a grey dot and an edge is represented as the connection between two nodes. The degree defines the number of connections connected to a node, in this case the number of connections connected to the blue dot. The path length, is the shortest distance between two nodes and clustering refers to the probability that neighbouring nodes are connected to each other. In the top example, which consists of nodes A through E, the degree of node A is 3, the minimum path length between nodes A and E is 2. The clustering coefficient of A is 2 3 . The three network topologies show an ordered network, which has a high C and long L, the random network has a small C and a short L, the small world network has a high C and shot L, which are properties similar to that of the human brain. In the ordered case the number of nodes are 8 N=8 and k=4 Chapter B Code listings B.1 Matlab Code DDA 1clc; 2clear all; 3%Inital setup of signals 4samples = 1024; 5bins = 2048; 6t = 0:1:samples*bins; 7fs = 2048; 8fc = 20; 9%Signal 1 10 S1 = 10*sin(2*pi*fc/fs*t); 11 f1 = 10; 12 f2 = 30; 13 %Signal 2 14 S2 = 10*chirp(t,f1/fs,t(end),f2/fs); 15 r = 4; 16 % Initiate variables 17 j = 0; p = 0; l = 0; n = 0; k = 0; storagea1 = 0; storageb1 = 0; 18 count1 = 0; count2 = 0; v = 1; Deltas(1) = 0; Deltas(2) = 0; 19 wincount = 0; RS = 0; window = samples; 20 21 22 for i = 2:length(s1) 23 windowcount = windowcount + 1; 24 25 %Signal 1 event detection based on variable M and Q 130 Code listings 26 %Toggeling comparator 27 if s1(i) < s1(i−1) 28 M(end) = 0; 29 elseif s1(i) > s1(i−1) 30 M(end) = 1; 31 end 32 %Shift register minima 33 M = circshift(M,[0 1]); 34 35 %Logic outlier decision 36 if sum(M(1:(length(M)/2))) == length(M)/2 && ... sum(M((length(M)/2)+1:end)) == 0; 37 j=j+1; 38 elseif nnz(¬M(1:length(M)/2)) < Q && ... nnz(M((length(M)/2)+1:end)) < Q 39 j=j+1; 40 end 41 42 %Same set−up for signal 2 43 if s2(i) < s2(i−1) 44 M2(end) = 0; 45 elseif s2(i) > s2(i−1) 46 M2(end) = 1; 47 end 48 M2 = circshift(M2,[0 1]); 49 50 51 if sum(M2(1:(length(M2)/2))) == length(M2)/2 && ... sum(M2((length(M2)/2)+1:end)) == 0; 52 p=p+1; 53 elseif nnz(¬M2(1:length(M2)/2)) < Q && ... nnz(M2((length(M2)/2)+1:end)) < Q 54 p=p+1; 55 end 56 57 58 59 % signal1 minimum found start count 60 if j == 1 61 count1 = count1 + 1; 62 % signal1 second minimum found store count and reset 63 elseif j == 2 64 count1 = count1 + 1; 65 storagea1 = count1; B.1 Matlab Code DDA 131 66 count1 = 0; 67 j = 1; 68 n = 2; 69 end 70 71 % signal2 minimum found start count 72 if p == 1 73 count2 = count2 + 1; 74 % signal2 second minimum found store count and reset 75 elseif p == 2 76 count2 = count2 + 1; 77 storageb1 = count2; 78 count2 = 0; 79 p = 1; 80 l = 2; 81 end 82 83 % Minimum trasition period found in both signals abs(diff) 84 if n==2&&l==2 85 Deltas(i−3) = abs(((storagea1) −(storageb1))); 86 storagea1 = 0; 87 storageb1 = 0; 88 n = 1; 89 l = 1; 90 v=v+1; 91 else 92 Deltas(i−3) = 0; 93 end 94 95 % accumulate error for N samples (window) 96 if wincount < window 97 RS = Deltas(i−3) + RS; 98 else 99 %If window reached store phi(i) 100 k=k+1; 101 RSfinal(k,:) = [k*window RS]; 102 wincount = 0; 103 RS = 0; 104 end 105 end 106 107 108 %Synchronisation index 109 SI = RSfinal(2,:);