ALGORITHMS TO ASSESS STRUCTURAL AND FUNCTIONAL REMODELING SURROUNDING NEURAL PROSTHESIS By Akash Saxena A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Electrical Engineering - M a ster of Science 2020 ABSTRACT ALGORITHMS TO ASSESS STRUCTURAL AND FUNCTIONAL REMODELING SURROUNDING NEURAL PROSTHESIS By Akash Saxena Recorded signals from implanted electrodes in the brain can be used to restore motor function for patients suffering from neurological injuries or neurodegenerative diseases, either through decoding the signal to control exterior assistive devices or as a trigger for downstream neurostim ulation. However, an ongoing challenge in the field is the limited and unpredictable ability to record from neural prosthesis for longer periods of time due to the biological response from the brain, the technical issues related to the signal processing al gorithms used, and the stability of the microelectrode itself. Our lab has recently uncovered multiple new observations of structural and functional remodeling of neurons surrounding implanted electrodes, which may contribute to the instability in recordin g quality over time. Some of these observations indicate losses, or changes, in syna ptic connectivity of individual local neurons with the surrounding network. In this thesis, I have introduced signal processing tools which will help us to understand and c haracterize the structural and functional remodeling happening around the neural pro sthesis, through the incorporation of multi - unit activity, coherence analysis, spike - triggered average of the local field potential, and spike triggered covariance of the l ocal field potential. A better understanding of these changes can provide us with ne w insights into how the network is remodeling itself over time and the major factors which influence these processes, which will help us to enable the design of better and more biocompatible neural prostheses. iii ACKNOWLEDGEMENTS I would like to thank Dr. Eri n Purcell, for always encouraging me, providing her excellent guidance towards bridging the gap between electrical engineering and neuroscience and making it a fun experien ce, this research would not have been possible without her. I would also like to tha nk the members of the lab for always supporting me and taking care of me both inside and outside the lab. iv TABLE OF CONTENTS LIST OF TABLES iv LIST OF FIGURES v CHAPTER 1 I NTRODUCTION 1 1. Overview of neuronal structure and function 1 1.1. Signal generation of a neuron 1 1.2. Transmission of neuronal signals 4 2. Extracellular s ignal d etection by i mplanted e lectrodes 6 2.1. Neuro nal geometry and architecture 8 2.2. Temporal scaling properties 10 2.3. Role of volume conduction 10 3. Implanted electrodes in brain: Progress and challenges 11 3.1. Silicon based micro electrodes 12 3.2. Challenges to effective integration 14 3.2.1. Biological response 1 5 3.2.1.1. Synaptic remodeling 1 6 3.2.1.2. Modulation of network activit y 1 6 3.2.1.3. Remodeling of subtype - specific markers: excitatory (VGLUT1) and inhibitory (VGAT) synaptic markers 1 7 3.2.1.4. Effect o n ion channel expression 1 7 3.2.2. Technical issues 18 4. Summary and organization of the T hesis 19 CHAPTER 2 MULTIUNIT ACTIVI TY, LOCAL FIELD POTENTIAL, AND THE SPIKE TRIGGERED AVERAGE 2 0 1. What is Multiunit Activity? 2 1 1.1. Relationship bet ween MUA and LFP 2 3 1.2. Results 2 4 1.3. Interpretation of results 2 5 2. What is Spike Triggered Average? 2 5 2.1. Steps of computation in the code used for spike triggered average 2 7 2.2. Results 3 0 2.3. Interpreting STA 3 3 2.4. Applications for future use 3 3 2.5. Disadvantage of using Spike triggered average 3 4 3. STA - LFP and LFP - MUA coherence 3 5 CHAPTER 3 SPIKE TRIGGERED COVARIANCE AND CORRELATION 3 6 1. What is Spike triggered Covariance - Correlation? 3 7 1.1. Covariance and Correlation Matrix 38 1.2. Spike triggered Covariance and Correlation 4 1 v 2. Results 4 2 2.1. Discussion and interpretation 4 5 3. Futu re use of STC as a tool 51 CHAPTER 4 CONCLUSION 5 2 APPENDICES 5 5 APPENDIX A: Coherence analysis 5 6 APPENDIX B: Local field potential - spike triggered average 59 APPENDIX C: Spike triggered covariance and correlation 67 REFE RENCES 7 4 vi LIST OF TABLES Table 1. Major contributors influencing LFP [12] 6 Table 2. Examples of m icro electrode arrays and their applications 1 4 Table 3. Covariance Matrix for cases in Fig. 2 8 4 7 vii LIST OF FIGURES Figure 1 . The F igure shows an action potential recorded from a pyramidal neuron in t he CA1 region of a rat hippocampus, illustrating commonly measured parameters [10] 2 Figure 2 . (A) B asic subunit of Shaker B. Four of these subunits are assembled into a functional channel [11]. (B) a schematic view o f the channel from the outside showing the assembled four subunits or domains (or subunits, in the case of potassium channels) [11] 4 Figure 3 . (A) E xample of a neuron with the structure for a closed field (B) Neuron structure that would generate a n open field (bottom left to right) (C) O rganization of neurons that would generate an open field (D) Organization of neurons that would generate a closed field [13] 9 Figur e 4 . (A) 8 * 8 Utah multielectrode array (B) Single electrode site in th e Utah multielectrode array (C) Labelled diagram of single shank Michigan style probe (D) Single shank Michigan style probe [26] 13 Figure 5 . Steps for computation of MUA [56] 22 Figure 6 . Coherence plot between MUA and LF P 2 4 Figure 7 . LFP in frequency domain 2 4 Figure 8 . The waveform represents Local Field Potential and the vertical lines below it shows the spikes recorded simultaneously [5 9 ] 2 6 Figure 9 . The black waveform is the LFP - STA waveform [5 9 ] 2 7 Figure 1 0 . Steps of computation of LFP - STA 29 Figure 1 1 . Shows the overlap of snippets used to construct the STA in F igure 1 2 , the spike occurs at 0 msec 3 0 Figure 1 2 . LFP - STA for the sn ippets in F igure 1 1 3 1 Figure 13 . LFP - STA for a wind ow size of 8 msec (4 msec on each side) 3 1 Figure 14 . Shows the LFP - STA for a window size of 12 msec (6 msec on each side) 3 2 Figure 1 5 . Shows the LFP - STA for a window size of 20 msec (10 msec on e ach side) 3 2 Figure 16 . Covariance matrix for a microelectrode with 16 channels 40 viii Figure 17 . Correlation matrix for the covariance matrix in Figure 16 4 0 Figure 18 . The LFP Spike triggered a verage for the mentioned animal 4 2 Figure 19 . Overlap of al l the 194 LFP snippets used for construction of STA 4 3 Figure 20 . LFP - STA Covariance Matrix 4 3 Figure 2 1 . A contoured version of the LFP - STA covariance matrix in Fig ure 20 4 4 Figure 2 2 . LFP - STA Correlation Matrix 4 4 Figure 2 3 . Contoured versi on of the LFP - STA Correlation matrix in Fig ure 2 2 4 5 Figure 2 4 . Different scenarios that can exist when covariance is calculated (A) close to zero phase shift, (B) A random phase shift x, (C) Close to 180 degrees phase shift 4 6 Figure 2 5 . (A) LFP snippe ts from 66 to 70 with LFP - STA (B). LFP snippets from 101 t o 105 with LFP - STA 48 Figure 26 . (A) LFP snippets from 42 to 48 with LFP - STA (B). LFP snippets from 78 to 8 2 with LFP - STA 49 1 CHAPTER 1 INTRODUCTION This thesis pertains t o the development of new signal processing algorithms to unders tand the tissue response to implanted electrode arrays in the brain. In particular, these tools are incorporated into the data analysis to unmask th e effects of synaptic remodeling [1], changes in local dendritic arbors, and alterations in intrinsic excitability surrounding devices on the recorded extracellular potential [2] [3] . Many of the techniques described herein pertain to understanding the relationship between the recorded spikes of indiv idual neurons (unit activity) and the activity of the broader synaptic network (the local field potential). A brief overview of neuronal structure and function is provided to provide contex t for the subsequent chapters, as well as a primer on implanted ele ctrode arrays . 1. Overview of neuronal structure and function The brain is made up of individual nerve cells or neurons which are the basic functional cellular unit of the brain. T he structure of a typical neuron includes projections which serve as the locati on of signal input and output (dendrites and axons, respectively). The fibrous structures known as dendrites, together known as the dendritic tree, are covered with even smaller structures known as spines. A single axon projects from the soma, and as it a pproaches its target (which can be other neurons, muscles or gland cells), it branches into a number of smaller axons known as terminals or knobs. These points of functional contacts are kn own as synapses, whereby an electrical signal in a neuron can be tr ansmitted to a downstream target. 1.1. Signal generation of a neuron electrophysiology. The signal is largely shaped by the intracellular influx of sodium ions followed 2 by the e fflux of potassium ions ( Figure 1 ). In 1952, Hodgkin and Huxley wrote a series of five papers [4][5][6][7][8] that described the seminal experiments they conducted to determine the laws that govern the movement of ions in a nerve cell during an action pote ntial. The first paper examined the function of the neuron membr ane under normal conditions and outlined the basic experimental method pervasive in each of their subsequent studies. The second paper examined the effects of changes in sodium concentration o n the action potential as well as the resolution of the ionic cu rrent into sodium and potassium currents. The third paper examined the effect of sudden potential changes on the action potential (including the effect of sudden potential changes on the ionic conductance). The fourth paper outlined how the inactivation pr ocess reduces sodium permeability. The final paper put together all the information from the previous papers and turned them into a mathematical model. Critical to understanding the under pinnings of the action potential (AP) was the development of the voltage membrane potential (V m) at a desired level. The difference between the measured variable and the set poi tem that decreases Figure 1. The Figure shows an action potential recorded from a pyramidal neuron in the CA1 region of a rat hippocampus, illustrating commonly measured parameters [1 0 ] 3 the magnitude of the error. An amplifier measures the potential difference between an intracellular electrode (V - in) and an extracellular electrode (V - out). The outp ut of this amplifier (V m) is the controlled variable and is compared wi th the command potential (V - command), or set point, by an amplifier called the control amplifier. If V m is not equal to V - command, an error signal causes a current to flow through an axial wire that is connected to the output of the control amplifier and inserted longitudinally through the axon. The current then flows out through the membrane to a grounded electrode to complete the circuit. The current passing through the axial wire ra pidly and continuously causes V m to remain equal to V - command. One way to measure the transmembrane current (I m) is simply to measure the current flowing out of the control amplifier. The transmembrane currents are carried by ion channels embedded in th e neuronal cell membrane, which are large macromolecular proteins that o ften consist of several peptide subunits. These channel proteins extend across the lipid bilayer and are in contact with the aqueous environment on both sides of the membrane. The chan nel protein forms a water - filled pore that allows ions to cross the memb rane by electro - diffusion. Channels can open and close spontaneously and in response to various stimuli. Voltage - gated ion channels have an open probability (the fraction of time the c hannel is open) that depends on V m. A typical Na+ voltage gated ion cha nnel consists - subunit contains all the functional characteristics of Na + channels, including the pore - forming region and - subunit contains four homologous domains (desig - helical membrane - spanning segments (S1 to S6). Segment S4 has a positively charged amino acid at every third residue and is the voltage sensor [9] . A sequence of residues connecting S5 to S6 on the extracellular side of the channel, called the P region or P loop, dips partway into the membrane to form part of the channel pore. Figure 2 (A), shows a similar arrangement but in the case of 4 shaker B. K+ channels are similarly structured, with the exception that four i ndividual subunits (rather than domains) aggregate to form the tetrameric channel. Figure 2 . (A) B asic subunit of Shaker B. Four of these subunits are assembled into a functional channel [ 11 ] . (B) a schematic view of the channel from the outside sho wing the assembled four subunits or domains (or subunits, in the case of potassium channels) [ 11] The four voltage - sensing domains (VSDs) from each S4 are a rranged symmetrically around the central pore [11] as shown in Figure 2(B) . At a negative V m, S4 is electrostatically attracted to the inside of the membrane, and in that position it holds the channel pore in the closed configuration. Upon depolarization, S4 moves in the outward direction and may also rotate. That movement is c oupled to the opening of the channel. The movement of charges on S4 generates a measurable current, known as the gating current. 1.2. Transmission of neuronal signals Once the action potential is generated, how does the spread of electrical activity occur, and across what time course? The three passive membrane properties that are important to understand the spread of electrical activity are the membrane resistance, the membrane capacitance, and the internal resistance of long thin processes or cells. By modelli ng the membrane as an el ectrical circuit, the electrotonic potential can be analyzed. Many ion channels behave like conductors and is open. Most permeant ion s are distributed asymme trically across the plasma membrane. This (A) ( B ) 5 results in a chemical driving force that tends to push the ion through the open channel (E). This chemical force functions as a battery. Membranes usually contain several different types of ion channels that are ea ch present in large numbers. In electrical terms, single channels in the membrane represent conductors arranged in parallel, i.e. the total conductance would be equal to the multiplication of, No (number of open ion channels), and t he conductance of a sing le channel, connected in parallel with the elements representing the ion channels (the membrane capacitance ). Membrane ionic curre nts cause V m to change gradually, because of membrane capacitance. In response to a pulse of constant current, V m follows an exponential time course to a new level. - m, an d is equal to the membra ne capacitance multiplied by the membrane resistance. The passive spread of small - amplitude, subthreshold signals along the surface of an axon is affected by membrane and axoplasmic resistances. The amplitude of these signals decrea ses as an exponential fu nction of distance along the axon, which can be modeled by cable theory. The neuron makes efficient use action potential at locations of highly concentrated volt age - gated sodium channel expression amongst intervening periods of passive propagation along the length of the axon. If the signal is successfully propagated to the synaptic terminal, transmission to a downstream neuron may occur th rough release of a neuro transmitter across the synaptic cleft. Neurotransmitters may have excitatory or inhibitory effects on the downstream neuron, depending on the nature of the receptor channels expressed in the post - synaptic membrane. In a chemical syn apse , t he axon terminal is the top knob - like structure and the spine of the receiving neuron underlies it. There is a presence of a large number of vesicles clustered in the 6 presynaptic axon terminal. These vesicles contain chemical synaptic transmitters f or a particular synapse. In between the pre - and post - synaptic terminals , there is a space known as synaptic cleft. When a synapse is active, the vesicles in the presynaptic terminal fuse with the presynaptic membrane and release their content of transmitt er into the synaptic spa ce. The transmitter molecules diffuse across the narrow synaptic space and attach to specific chemical receptors on the surface of the post - synaptic membrane, which activates the post - synaptic target cell. 2. Extracellular signal dete ction by implanted elect rodes Electric current contributions from all active cellular conductances within a volume of brain tissue superimpose at a given location in the extracellular medium and generate a potential, Ve (a scalar measured in Volts), with r espect to a reference po tential. The difference in Ve between two locations gives rise to an electric field (a vector whose amplitude is measured in Volts per distance) that is defined as the negative spatial gradient of Ve. The low frequency band of Ve, w hen recorded by a microe lectrode in the brain, is referred to as the local field potential (LFP). The LFP is of the electrode site [1 2 ]. Other examples of cont ributors to the measured potential include sodium action potentials, calcium - mediated spikes, intrinsic currents and resonance, spike afterhyperpolarizations, and neuroglial interactions (Table 1). Contributors to Ve Description Fast action Potentials (Na +) The strongest current contribute to the LFP that much because the fields generated are of very short duration. Table 1 . Major contributors influencing LFP [12] 7 Table 1 (cont d ) Calcium spikes (Ca2+) Long lasting (10 - 100 ms, amplitude varies from 10 - 50 mV) spikes can actively propagate within the cell and due to this their contribution to the LFP can be substa ntial Intrinsic currents and r esonance Some neurons possess resonant properties i.e. when the intracellular depolarization is large enough, self - sustained oscillations are generated. For example: currents flowing through hyperpolarization de - inactivated c yclic nucleotide - gated channels (I - h) and low - threshold hyperpolarization - induced transient Ca2+ currents, which often lead to burst firing (I - t) contribute to intrinsic resonance and self - sustained oscillations. Since resonant properties depend on both fr equenc y and voltage, its impact varies with LFP in a complex manner. Spike after polarizations and down states Elevation of the intracellular concentration of a certain ion may trigger influx of other ions through activation of ligand - gated channels, and this w ill, in turn, contribute to LFP (the amplitude and duration of such burst - induced afterhyperpolarizations can be as large as synaptic events) Gap Junctions and Neuroglia interactions Direct electrical communication between neurons through gap juncti ons (a also can affect neuronal excitability and contribute indirectly to the extracellular field. The signal that is detected from the electrode placed inside the brain is the summa tion of activities from all the sources listed in the table above and many other sources. In order to obtain 8 the LFP or the units/spikes, signal processing techniques are ap plied: a first step is to filter the signal into relevant frequency bands. For spik ing activity, a high pass filter is applied at the desired range, which will attenuate all the frequency content less than that value. This high pass filtered signal can be further processed to obtain the spikes. Often, a threshold is set as several times the standard deviation above or below the mean of the signal to initially extract the spikes as threshold - crossing events. Principal component analysis can further distingui sh individual spikes based on their similarity in waveform shape. Similarly, if we are interested in obtaining the LFP, which are slow travelling waves, a low pass filter is applied which attenuates all the high frequency content above a chosen value. Anal yzing the LFP requires a broader understanding of the factors that influence it: 2.1 N euronal geometry and architecture Neuronal architecture and geometry influence the LFP significantly. For example, neurons that bution to the extracellular potential, whereas the 2 ]. This, in turn, depends on the distance between the current sources and the sinks (i.e., open conductances along the cell membrane), which is decided by th e geometrical features of a neuron. Figure 3 (A) shows an example of a neuron where the source and the sink are close together giving rise to a closed field, whereas Figure 3 (B) shows a neuron where the source and the sink are at a larger distance (giving rise to an open field). The measured potential depends on how these neurons are arra nged spatially in different regions, which can lead to synchronous superposition (e.g., in Figure 3 (C) neurons lie parallel to each other, contributing significantly to LFP). A different arrangement might lead to the waves cancelling each other out (e.g., in Figure 3 (D) the neurons individually generate open fields, but the way they are arrange d gives rise to a closed field). In other words, a small LFP 9 amplitude could occur even when the magnitude of the population activity might be substantial, if the arr angement of neurons occurs in a manner similar to Figure 3 (D). This is an important consid eration when assessing LFP amplitude as a reflection of signal quality, and it also is relevant for understanding field contributions in the context of dendritic loss and remodeling surrounding implanted electrodes. Figure 3 . (A) E xample of a neuron with the structure for a closed field (B) Neuron structure that would generate a n open field (bottom left to right) (C) O rganization of neurons that would generate a n open field (D) Organization of neurons that would generate a closed field [13] (A) ( B ) ( C ) ( D ) 10 2.2 Temporal scaling properties In a ddition to cytoarchitecture, a second critical factor in determining the magnitude of the extracellular current is the temporally synchronous fluctuations of the membrane potential in large neuronal aggregates. A quantitative feature of the LFP is that the magnitude of LFP power (that is, the square of the Fourier amplitude) is inversely related to temporal frequency f ; that is, there is 1/ f ^n scaling (the exact value of n depends on various factors). There are several reasons that this might happen. On e of them might be due to the low - pass frequency filtering property of dendrites [1 4 ], typically owing to a serial capacitance which depends on the distance between the soma and the location of the input, and on the membrane time constant. As the electroto nic length and input resistance of neurons can be effectively altered by synaptically induced excitatory and inhibitory conductance changes [1 5 ][1 6 ], the frequency filtering performance of neurons depends not only on the geometric characteristics of the ne uron s but also on their physiological state. Another reason is related to the capacitive nature of the extracellular medium itself [1 7 ],[1 8 ]. Network mechanisms also contribute to the 1/ f ^n feature of the power spectrum [19] . In a brief time window, only a lim ited number of neurons can be recruited in a given volume, whereas in longer time windows, the activity of many more neurons can contribute to the LFP, therefore generating larger amplitude LFP at slower frequencies. 2.3 Role of volume conduction The elect ric field specifies the forces acting upon a charged particle. The field is defined at be transmitted through volume (for example, through brain tissue); a phe nome non known as volume conduction. In a volume, the Ve induced by a current dipole depends on the magnitude 11 and location of the current source, and on the conductivity of the extracellular medium. Conductivity in the medium depends on the degree of isotro py a nd homogeneity of the medium and is therefore a function of a number of factors, including the geometry of the extracellular space. The relationship between Ve and the current source density (CSD) J (measured in amperes per meter square) at a particula r po of electromagnetism. In their simplified form (that is, when the magnetic contributions can be neglected), these equations dictate that in S m 1) is the waveform and functionality of the spatiotemporal Ve deflections. Measurements of the extracellular medium in the re levant frequency range (<10 kHz) have not yet fully resolved thi s issue, with some experiments concluding that the extracellular medium is anisotropic and homogeneous[ 20 ],[ 21 ], and others suggesting that it is strongly anisotropic, inhomogeneous [2 2 ],[2 3 ] and may even possess capacitive features. 3 Implanted electrodes i n the brain: Progress and challenges A microelectrode is an extremely fine wire or microfabricated electrode that en ables the measurement of the electrical potential following implantation int o nervous tissue. In recent years, numerous fabrication strategi es encompassing a wide range of materials and architectures have permeated the neural engineering field [2 4 ]. Historically, silicon microtechnology has been a mainstay approach to meet the dem ands of neural applications, allowing multiple functions to be i ntegrated on a single implantable microsystem. Besides reproducibility and low variability of silicon - based microelectrodes, combination of various functionalities like standard electrophysiol ogy, integrated signal processing, local drug delivery, neuroche mical detection and optogenetic stimulation is also possible using these microsystems. This ability makes silicon 12 microelectrodes good candidates to provide high - resolution recording and stimu lation in the electric, fluidic, chemical or optical domain in m ore complex neurophysiological experiments in the future [ 2 5 ]. 3.1 Silicon based microelectrodes Some of the most widely used silicon based microelectrodes are the Utah multielectrode array and Michigan style probes (Figure 4 ) . The Utah array is a comm ercially available intracortical electrode array consisting of up to 100 silicon needle - shaped electrodes, which are produced via microscale fabrication techniques such as thermom igration, a combination of mechanical and chemical microm achining, metal deposition, and encapsulation with a polymer made of imide monomers The Utah array enable s users to acquire neural activity from different cortical areas. However, t hey are limited in their ability to target deep neural structure s in the axial direction. This is where Michigan probes come into play where the electrode sites are laid across the shank and can be used to record from multiple layers at the target site. Mi chigan probes have been successfully used in several neuroscienc e applications, but they also suffer from some disadvantages related to limited thickness (>15 um) and length ( < 8 mm). Nonetheless, these devices remain important tools in neurophysiology stud ies. Examples of the application of Michigan probes and Utah arrays are provided in Table 2. 13 Figure 4 . (A) 8 * 8 Utah multielectrode array (B) Single electrode site in the Utah multielectrode array (C) Labelled diagram of single shank Michigan style probe (D) Single shank Michigan style probe [2 6 ] ( A ) ( B ) ( C ) ( D ) 14 Multielectrode - Array (MEA) E xample A pplications Utah multielectrode array Neuroprosthetic arm control via cortical motor activity in rhesus monkeys (Macaca mulatta) using Utah arra ys [2 7 ]. Electrically stimulate the spinal - cord according to decoded neural signals from the Utah array in the motor cortex of monkeys [2 8 ] . Induced tactile sensations in the hand via targeted neural stimulation with Utah arrays implanted in the somato sensory cortex, as illustrated in [2 9 ] . Michigan style probe I nvestigated neural activity in the hindlimbs of rats with spinal cord injuries using a 16 - channel Michigan probe and explored treatment methods based on neuro - stimulation [ 30 ]. I nvestigated t he use of a neural interfacing system to recover neural function after brain injuries i n small animals (i.e. rats) d emonstrating that the missing brain functions can be restored using a brain machine brain interface [ 31 ] . Table 2 Examples of silicon m icro electrode arrays and their applications . 3.2 Challenges to effective integration The probes we discussed above are useful in a wide array of applications and have been very useful, but an ongoing problem with these electrodes is that t hey often do not provide reliable recordings when implanted for long periods of time. Signal amplitud es vary on a daily basis 15 [3 2 ], which can lead to difficulties analyzing, decoding, and interpreting the signals recorded. There are a wide array of reason s why this may happen, but two major contributors are: (1) the biological response to the implantation of the device (neuronal and non - neuronal response), and (2) the technical issues that arise with implantation of the electrode (signal processing algorit hms, electrical and mechanical failures). 3.2.1 Biological response In terms of the biological response, wor k by Biran et. al, used neuronal loss and gliosis as parameters to investigate device integration as a function of distance using immunohistochemistry (3 3 ). Significant neuronal loss was observed at 4 weeks within the first 100um compared to stab control, which did not fully resolve until ~500um, and a significant loss of neurofilament extended out beyond 200um. Non - neuronal cells in the brain include spec ialized supportive cells (such as astrocytes, microglia, oligodendrocytes and NG2 - glia) and neurovascu lar cells. Microglia are the resident macrophage cells of the brain, which initiate the foreign body response by the release of inflammatory factors. They rapidly react to inflammation in the central nervous system, leading to encapsulation of the device a t the site. Astrocytes are activated by the microglial signaling and likewise contribute to the device encapsulation. The barrier nature of gliosis has tr aditionally been assessed through in vivo measurements of the impedance of the tissue/electrode interf ace and modelled using static circuit elements. However, the electrode/tissue interface, especially in the presence of reactive gliosis, cannot be fully d efined by these traditional methods [3 4 ]. More specifically, the barrier nature of gliosis is reflecte d in models of the effective volume of tissue activated, where greater gliosis reduces the number of neurons stimulated. The stimulation paradigm affects the impact of the glial barrier: in constant - voltage stimulation, 16 voltage is controlled, and the actua l current delivered to tissue varies as the tissue response evolves (increased impedance due to gliosis reduces the stimulation delivered). Meanwhile, the impact of glial encapsulation on the quality of signals recorded in vivo remains ill - defined. A few l ines of indirect evidence support the idea that glial encapsulation acts as a barrier to signal detection by implanted electrodes. For instance, astroglio sis, as identified by increased glial fibrillary acidic protein (GFAP) immunoreactivity, was associated with reduced recording quality of Utah - style arrays implanted in the rat cortex in a study that investigated the relationship between histology and recording quality. Additional possible biological impacts to the local tissue are summar ized below: 3.2.1.1 Synaptic remodeling Synaptic remodeling is largely affected by its neurochemical environment, which, in turn, is affected by the device implan tation. Due to the electrode implantation, cellular membranes are punctured, and astrocytes and microg lia are activated. This results in increase of neurotransmitters in the extracellular environment [3 4 ]. This increase in neurotransmitter leads to the cre ation of a gradient which may attract and reinforce gliosis [3 5 ][3 6 ]. There is also the evidence of in creased markers of inhibitory synaptic transmission surrounding devices [1]. All of these changes in the surrounding environment of the implant have the p otential to lessen the effect of stimulation and add noise to the recorded signal. Nonetheless, these adaptations may reduce the spread of excitotoxicity, neuronal loss and also the likelihood of excessive synchrony. 3.2.1.2 Modulation of network activity Interco nnected astroglia are able to orchestrate synchrony through the integration of signaling within neuron al circuits and across functional regions of the brain [3 7 ]. The 17 stimulated actions of a single astrocyte could dictate functional consequences on an enti re network of neurons [3 7 ]. Given evidence for astrocyte coordination of neuronal networks, reactive gliosis likely impacts not only the generation and transmission of action potentials between single neurons, but also the broader population of activity de tected and stimulated by electrodes implanted in the brain. 3.2.1.3 Remodeling of subtype - specific m arkers: excitatory (VGLUT1) and inhibitory (VGAT) synaptic markers Neurons, astrocytes, and microglia each have multiple unique subtypes which may be affected by in jury in different ways. Neural circuitry in the brain is complex, where individual cells may receive thousands of connections from other cells, and neurons contain a large diversity of form and function. In recent work from our lab [1], a progressive shift from VGLUT1 to VGAT predominance at the device interface was observed over time. The result excitation to a later elevation in inhibition. 3.2.1.4 Effect on ion channel expression Ion channel expression and f unction is highly dynamic and modulated by many factors [2], including changes to the surrou nding environment caused by injury [3 8 ] and inflammation [3 9 ]. Channel modulation can impact not only the signal generation capabilities of single neurons, but also their frequencies, patterns, and waveform characteristics that underlie information encodin g. Additional work from our lab [ 2 ] investigated four voltage - gated ion channels (Kv1.1, Kv4.3, Kv7.2, and Nav1.6) based on their roles in regulating action potenti al generation, firing patterns, and synaptic efficacy. To understand whether shifts in the e xpression of selected ion channels 18 occurs at the interface of implanted single - shank microelectrode arrays, quantitative immunohistochemistry was performed over 6 w eeks (with time points at 1 day, 1 week and 6 weeks). These observations were analyzed based on their expression, both spatially and temporally. Based on the spatial expression, the results support a shift toward a decrease in sodium channel expression and an increase in potassium channel expression over the chronic 6 - week time course. When the t emporal expression was depression at 6 weeks, coupled with the sustained elev ation in all Kv channels at both time points, indicates a temporal shift from hyper - to hypo - excitability within the recordable radius of the device. 3.2.2 Technical issues While biological factors are important determinants of signal quality, non - biological con tributors are also important to consider. Technical issues include poor signal detection due to problems with noise, signal processing techniques, and electrical and mechanical failures. Signal processing affects the signal based on the spike sorting algor ithms used, the features used for classification of spikes, how that data is pr ocessed after spike sorting, what kind of model is being used and how reliable it is. For example, as discussed by Michelson et. al. [ 40 ], the accuracy of a machine learning dec oder (MLD) depends upon the features that are used for decoding after spike sor ting (e.g., PCA - based single unit activity, PCA - based multi - unit activity, or a constant threshold). Additional observations underscore how different signal processing technique s affect the data and the results [40]. Electrical and mechanical failure likew ise can drastically affect the recording. For example, a sudden increase in impedance caused by lead breakage or electrode site delamination, 19 insulation failures which might res ult in a drastic decrease in impedance, [41] [42] tissue - device mechanical misma tch, device variability due to subtle differences during the microfabrication process, and motion artifacts. 4 Summary and o rganization of the T hesis In th is introductory chapter , we have reviewed the source of the signal we obtain for analysis and what feat ures might be affecting the signal, which can inform our signal processing techniques. The second chapter introduces two tools for assessing the connectivity between neuronal sp iking and the network activity , namely local field potential - spike triggered av erage (LFP - STA) and s pectral c oherence analysis for m ultiunit activity (MUA) and the LFP . We talk about computation of these two methods and how the results can be interpreted. The last part of the chapter discusses potential future applications of these tools. The third chapter introduces another tool , which works with LFP - STA , known as spike triggered covariance (STC) . STC looks further into the STA and the non - linear features responsible for STA. In this chapter , we begin by describing how simple correlation and covariance can be interpreted, and then move towards how covariance and correlation matrix can be calculated for LFP - ST A. We then propose the possible interpretations t he results might point towards and then discuss additional potential applications this tool can be used for in the future. The concluding chapter briefly summarizes the work, and how the tools introduced in each chapter integrate to reveal the changes in t he tissue response around the electrode, an d how these tools c ould be improved further and used in conjunction with other analysis techniques to reveal even more features. 20 C HAPTER 2 MULTIUNIT ACTIVITY, LOCAL FIELD POTENTIAL AND THE SPIKE TRIGGERED AV ERAGE The extracellular data recorded from the electrode can be processed and used to investigate how neurons behave as an assembly (LFP) or as single neurons (unit activity). Like the LFP, which represents the behavior of a neuronal population, multiunit activity (MUA) is yet another way to look at a neuronal population, but in a different way. Whereas the LFP is often viewed as representing the activity of the local neuronal assembly within the first ~300 m icrons of the electrode (and, predominantly, the synaptic conductances [1 2 ]), MUA is used to represent the superimposed activity of the individual neurons in a relatively smaller region around the electrode (~up to 100 microns). The computation performed b y a neuron can be formulated as a combination of dimensional reduction in stimulus space and the nonlinearity inherent in a spiking output. White noise stimulus or spike triggered average (STA) is one of the tools used to understand which dimension in the stimulus space the neuron is more sensitive to. L ikewise, for the spike triggered average of the local field potential (LFP - STA), it can reflect the coupling between individual neurons and the local network [ 4 3][ 4 4][ 4 5] [46] . LFP - STA represents the average LFP waveform surrounding individual spiking event s. The earliest attempts at neural characterization, including classic electrophysiological experiments [ 4 7 ], presented the neural system with simple, highly stereotyped stimuli with one or at most a few free parameters; this makes it relatively straightfo rward to map the input - output relationship, or tuning curve, in terms of these parameters. While this approach has been invaluable and is still often used, it has the shortcoming of requiring a stron g assumption about what the neuron is ding on this concept, white noise was used as an alternative stimulus, which 21 4 8 ]. The stimulus feature best correlated with spiking is then the spike - triggered aver age (STA); that is, the average stimulus history precedin g a spike [ 4 9 ] [ 50 ]. Spike triggered average is used for variety of applications. Some use it to study the stimulus right before the spike in order to understand the characteristics of what the wavef orm that makes the neuron spike is like with goals of cre ating a model mimicking the behaviour [ 5 1 ], some use it to study the relationship between spiking activity and the extracellular potential to understand the effects of membrane potential on extracell ular data [ 5 2 ], and some studies employ it as a tool to s tudy connectivity between different regions based on the increase or decrease in the value of STA for different regions [ 5 3 ]. In this chapter, we explore what the relationship between MUA and LFP tel ls us about connectivity within the network and how LFP - S TA can also be used to reveal the connectivity within a network. While both of these relationships tell us similar things, they do so in unique ways and reveal different characteristics about network connectivity. As such, they are useful tools to compleme nt our assessment of the impact of the tissue response on the quality of signals collected from implanted electrode arrays. 1. What is Multiunit Activity? Multiunit activity is used to represent the com bined activity of neurons in the vicinity of the electrod es in the order of (~100 um). Unlike unit activity or LFP, the frequency content present in respectively. It co ntains contains frequency content spanning the interim fr equency range (i.e., 300 - 6000 Hz). Likewise, MUA does not require the additional signal processing steps required to discriminate individual unit waveforms through sorting algorithms (often, through the use of thresholding, principal component analysis, or wavelets) [ 5 4 ][ 5 5 ]. The flow chart below shows how MUA is calculated in the program used in the analysis [ 56 ]: 22 Figure 5 . Steps for computation of MUA [5 6 ] As seen from Figure 5 , the frequency of MUA is between (300 - 6000 Hz). The (1/f) noise or pink nois e is a signal with a power spectral density (energy or power per frequency interval) inversely proportional to the frequency of the signal. This makes the spiking activity (> 500 Hz) less prone to noise because it carries only high frequency content, and f or this reason MUA is less affected by the noise. Additionally, since MUA does not require a thresholding step, its detection is less sensitive to underlying noise than unit detec tion . MUA also has been shown to have more independent features when looking at a particular brain area, because while the LFP (<300 Hz) at multiple electrodes tends to be highly correlated, comparatively little correlation exists in the MUA for the same r ecordings [ 5 6 ]. In the same study, the MUA predicted movement (reach and gras p) for a brain computer interface (B C I) with greater accuracy in comparison to the LFP and multiple spikes (MSPs). Apart from MUA being less affected by noise and having more feat ures for characterization, an additional advantage in the context of BCI is t hat it can be used to replace the for prediction require several computational ly intense steps, which in turn requires more processing power and more compl ex circuitry. These factors might also lead to an increase in the 23 size of the device being used for decoding or limit wireless transmission; computation of MUA is an alternative, simplified way of assessing spiking activity. 1.1. Relationship between MUA and LF P coherent electrical activity can communicate more effectively with each other [ 57 ] . In other words, neurons which fire in phase with one another will have greater functiona l connectivity and information transfer. In the case of the LFP, synaptic currents dominate; i.e., it represents the excitability of the broader neural network. The L to the system, and the local neuronal spiking calculating the coherence of the MUA and LFP may reflect the communication strength between the local population of in dividual neurons within 100 microns of the device and the broader neural network. As such, the coherence may serve as an indicator of the level of neuron - network coupling. However, previous studies which calculated low coherence values between MUA and LFP suggest that this interpretation must be assessed carefully, as the relationship is likely to be complex and may require a more granular analysis within specific frequency bands [ 5 8 ]. We can perform the coherence analysis using equation 1 where deno tes the Fourier transform of (a), and denotes the power in the frequency spectr um of (a) . The numerator is the product of fourier transforms of MUA and LFP and the denominator is the product of power in the frequency spectrum of LFP and MUA: Equation 1 24 1.2. Results We perform coherence analysis using E quation 1 on recordings from the primary motor cortex in the right hemisphere of a rat, recorded by a 16 channel Michigan probe. The recordings used are part of a larger data set used in a pre viously published paper [ 2 ]. Figure 6 . Coherence plot between MUA and LFP Figure 7 . LFP in frequency domain 25 1.3. Interpretation of results From Figure 6 , we see that maximum coherence exists in the gamma band (20 Hz - 90 Hz), and the maximum coherence exists at 82 Hz. This tells us that at that frequency there was the least phase difference between MUA and LFP, i.e maximum transfer of information occurs at that frequency. When coherence analysis is applied to the LFP and MUA, the frequencies over which the LF P and MUA have the least phase difference will reflect the frequencies at which the highest information transfer can occur ( the frequency at which they lock with each other ) . We can then inspect the LFP power spectrum and determine the relative signal stre ngth at that frequency (Figure 7 ) , which might inform us about: (1) the strength of connectivity between MUA and LFP in that network, and (2) the amount of information transfer between LFP and MUA. As such, this calculation may contextualize parallel work in our lab which indicates decoupling of neurons from the broader network via image analysis of dendritic arbors and spines. 2. What is Spike Triggered Ave rage? An alternative approach to assess neuron - network coupling is through investigation of the spike tr iggered average of the LFP (LFP - STA). If we record spikes from a neuron (either intracellularly or via sorted units recorded extracellularly) and also r ecord the local field potential near that neuron simultaneously, we can look at the relationship between these two recordings. In the LFP - STA, snippets of the LFP recorded within an affixed time window is centered at the timestamp of each detected spike, a nd the average of these snippets yields a final waveform known as the spike triggered average of the LFP . 26 A summary of the general computational process of the LFP - STA follows: 1. Once we have the LFP and spike data, the first step is to align them in the time domain and index the times at which we observe the spikes [ 5 9 ]. Figure 8. The waveform represent s Local Field Potential and the vertical lines below it show the spikes recorded simultaneously [ 5 9 ] 2. After labelling these points on the extracellular waveform, we need to select the duration of the time window, which depends on various factors including the firing rate and the nature of the frequency response observed in the frequency domain. 3. Once we decide the size of this fixed time wind ow, we need to determine what information we want out of the waveform. Are we interested in information just before th e spike or information when the spike is happening? Based on the answer, in the first case, we take snippets from the extracellular wavefo rm just before the time points recorded in the first step, or we place the window centered around the time point in th e second case. Based on the results, we can see how the two responses are related in terms of phase (synchrony) and amplitude (strength of the response). In turn, this can yield new information on the relationship between individual neuronal spikes and the broader network activity [ 60 ] [61] . 27 Figure 9. The black waveform is the LFP - STA waveform [ 5 9 ] 2.1. Steps of computation in the code used fo r spike triggered average A brief summary of the code generated to calculate the LFP - STA is provided below. The MATLAB script is printed in Appendix A of this thesis. This code builds upon previous scripts used for LFP and spike sorting analysis [ 6 2 ][2]. The workflow is depicted in Figure 1 0 . 1) The code needs two files to start running. The first file is the raw extracellul ar time series data and the second file has the information about the sorted units of the raw data obtained by processing the raw data through another program (TDT analysis). O nce we have this data we can input them into the code. 2) The sorts file has a fold er which contains the unit waveforms for every channel. This file is useful to see which channel has units and which just recorded noise as units. The TDT analysis code does th is by high pass filtering the data at 500 Hz and creating a threshold. 28 The final step is applying principal component analysis (PCA) to the threshold crossings (2 - 3 msec snippets which are removed and stored), which clusters waveforms with similar shapes t ogether for each channel, hence giving us the units waveforms. As desired, a fin al step involves visual inspection of the data to confirm legitimate unit waveforms. In the code written for LFP - STA in this chapter, we are processing all detected units irres pective of the waveform shape. In the future, we may analyze waveform shape manu ally or automate this process. 3) The second part of the code takes the LFP data files of the channels raw data and applies a notch filter at 60 HZ to eliminate any contaminating hum noise. This information is stored in a separate matrix with respect to its channel until further use. 4) The third part of the code is used to index the time points at which these spikes occurred to mark those points in the LFP data and isolate a snippet of the desired duration second from that data centered around the time point. W e do this by storing the time points at which spikes occurred in a different matrix for each of the channels. Next, we use these time points with respect to their channel and m ark these time points on the LFP data acquired in second step. Once these time p oints are marked, a window is applied centered around that time point and that data is isolated. This is done for every spike in the channel, so in the end we will have a matri x of number of spikes , times the number of windows . 5) In the last part of the cod e we obtain the STA waveform. Each of the snippet s obtained in the previous step are added together, and then divided by the number of spikes (hence, averaging the mean LFP wav eform). Therefore, we obtain one STA waveform for each of the channels that were included in the processing after Step 2 for the recording time period of interest. 29 Data with Principal component analysis data Looking at the features of the data to see if it is a unit or not Information about channels with units stored Raw electrop hysiology data Applying Filters for separating spiking and LFP 300 Hz Low Pass Filter S piking data from TDT analysis Noisy LFP data 60 Hz notch filter Hum noise Removal Processed LFP Data Extracting the timepoints at which spikes occurred in a matrix LFP with marked time points Window Applied (duration selected) Snippets of desired duration at each marked timepoints A dding all the available snippets and dividing by the number of spikes Spike Triggered Average Figure 10 . Steps of computation of LFP - STA 30 2.2. Results Here , we look at the STA of recordings from the primary motor cortex in the right hemisphere of a rat, recor ded by a 16 channel Michigan probe. The recordings used are part of a larger data set used in a previously published paper [2]. For the recording used in the following examples, the duration was 300 seconds with 194 spikes detected on the selected channel within that time. I have analyzed the STA for different window lengths (4, 8, 1 6 , and 20 msec). The snippet overlap is only shown for a window of 4 msec ( Figure 1 1 ). Widening the window was performed to visualize LFP modulation surrounding the spike occurr ence over a broader time scale (thus capturing lower frequencies). Figure 11 . Shows the overlap of snippets used to construct the STA in Figure 1 2 , the spike occurs at 0 msec 31 Figure 12 . LFP - STA for the snippets in Figure 1 1 Figure 13 . LF P - STA for a window size of 8 msec (4 msec on each side) 32 Figure 14 . Shows the LFP - STA for a window size of 12 msec (6 msec on each side) Figure 15 . Shows the LFP - STA for a window size of 20 msec (10 msec on each side) 33 2.3. I nterpreting STA The spike trigger ed average of the LFP yields a visual representation of the relationship between the activity of individual neurons and the local field potential, i.e., single neuron activity with the local network. A large amplitude could reflect a high degree of couplin g between the individual neurons and the local network activ ity. In the context of investigating the tissue response to implanted electrodes, this could reflect changes in local synaptic connectivity due to the presence of the device. Likewise, the STA als o illustrates the phase shift between the spiking and LFP: d oes the spike occur on a relative up - or down - phase of the LFP? In summary, the STA can be used a tool to assess functional connectivity between the single unit activity and the LFP [ 6 3 ]. In turn, this may be an additional tool to assess structural remodel ing in neurons surrounding implanted electrodes. Importantly, our recent observations suggest that neurons local to the device lose synaptic connectivity with the network over time. 2.4. Applications f or future use Additionally, there are multiple opportunities for new applications in future use: 1) Its use can be further extended if multiple microelectrodes or electrodes with multiple electrode sites are used; using the spiking activity at one electrode s ite as an index for LFP at another electrode site could deli ver a map of functional connectivity between two regions. Also, when looking at the waveforms with one site as a spiking reference, we can look at phase differences between regions. 2) We can look a t these waveforms across different time points to show how c onnectivity between regions changes due to the electrode implant. We might be able to assess the STA as a reflection of the morphological changes that are happening in the vicinity of the electrod es. 34 3) It can be used for assessing the effects of genetic knoc kdown delivered locally to the electrode, revealing even more specific relations and effects due to the neural implant. For instance, we may use the STA to assess the effects of knockdown of synap tic transporters or ion channels surrounding devices [ 1 ] [2 ]. 4) STA and current source density (CSD) analysis can be used together, which can give us a better picture into the composition of LFP in terms of contributions from different neural populations [ 6 4 ] . 5) Another opportunity is to look at the spike - triggered aver age of the time frequency spectrum of the LFP (STTFA), as used in [ 6 5 ]. Here, instead of averaging over segments of the LFP centered on each spike, the two - dimensional time frequency energy plot o f the LFP is averaged around each spike. 2.5. Disadvantage of usi ng Spike triggered average The STA only provides us with the linear relationship between Spike - LFP activity and reveals nothing about the non - linear relationships between them. What this means is that when we look at the STA, it shows the average changes i n the LFP at every spike for a given time window. As such, we lose information about the variation in the LFP snippets and the LFP e reveal underlying structure in the data set. For instance, assessing the variation in the LFP snippets might reveal new observations that could reflect the underlying dynamics in neuron - network coupling, which could, in turn, be interpreted relative to s tructural and functional remodeling of neurons surrounding d evices. This is a primary motivator for exploring the Spike - Triggered Covariance of the LFP (STC - LFP), in Chapter 3. 35 3. STA - LFP and LFP - MUA coherence In summary, the LFP - STA is one tool to represen t the relationship between neuronal spiking and network activity, where higher amplitudes may reflect greater coupling between individual neurons and the broader network. Likewise, the coherence of LFP - MUA may be interpreted similarly as representing the s ynchronization between LFP and MUA, which can be further assessed within specific frequency bands. As a complementary analysis to the STA, LFP - MUA coherence also provides us with added information on the frequ ency dependence of the relationship between spi king and network activity: does the communication between LFP and MUA increase in a certain band of frequencies? The combination of these tools may reveal new reflections of the tissue response in recorded ext racellular electrophysiology data. 36 Ch apter 3 SPIKE TRIGGERED COVARIANCE AND CORRELATION In the previous chapter, we investigated the spike triggered local field potential average (LFP - STA) and how it can be used as a tool to understand and analyze the influence of neuronal spiking on network activity in different scenarios. However, the LFP - S TA only reveals the linear characteristics of the neural response, while neurons exhibit other non - linear behaviors that are not revealed by the STA. In order to get a better understanding of these non - lin earities, how they influence the neuronal response, and what this might reveal about the functional connectivity, we explored the use of spike triggered covariance of the LFP (LFP - STC). Spike triggered covariance has been used to reveal highly complex mod alities exhibited in the space of sensory stimuli. For instance, Schwartz et al. used STC to reveal and characterize a neuron model with gain control [ 6 6 ]. In a separate paper using non - centered spike triggered covariance, it was shown that neurotrophin - 3 acts as a developmental regulator of receptive fiel d properties in retinal ganglion cells [ 6 7 ]. STC also can be used to extract a low - dimensional subspace of the full stimulus space that is primarily responsible for generation of the neural response [ 6 8 ], including both excitatory and suppressive component s in monkey V1. Spike triggered covariance takes STA a step further and is used to reveal the hidden and most relevant features of the stimuli that might be shaping the neuronal response, where that inform ation may be used to characterize different non - lin ear behaviors by capturing those non - linearities in filters to recreate that neuronal response. In this chapter, we use spike triggered covariance and correlation in a novel application to explore what it may reveal about coupling between the single neuro n (spiking activity) and local network activity (local field potential). 37 1. What is Spike triggered Covariance - Correlation? Before describing Spike triggered Covariance or Correlation, we provide a brief rev iew on how the simple covariance and correlation is calculated. For two variables X and Y, both with N number of elements, the following formula is used to find the covarian ce between these two variables (where and are the mean of the X and Y res pectively): Equation 2 So, every element in the array is subtracted with the mean for both the variables. The numerator in the equation is also known as the sum of cross products, where the elemen ts of both arrays are multiplied elementwise (dot product) and then summed together . The last step is to divide this number by the degrees of freedom which in this case in N - 1, giving us the covariance. The resulting covariance value, if positive, tells us that both the variables are increasing or decreasing together. If negative, it mea ns that one of the variables is increasing and one of them is decreasing. Now that we understand the covariance, we can look at how correlation is calculated. The equation t o be used is as follows: Equation 3 In E quation 3, the covariance between X and Y is divided by the product of their standard deviations denoted by and . By doing this, the covariance between X and Y is standardized between 1 and - 1 and is a dimensionless quantity. The correlat ion matrix, in addition to revealing positive or negative linear relationships, also reveals the strength of the relationship 38 between these variables. A value closer to (+ 1) reveals a stronger positive correlation whereas a value closer to ( - 1) reveals a s trong negative correlation. 1.1. Covariance and Correlation Matrix Covariance and correlation between two variables results in a single value which describes the linear relatio nship and strength of the relationship between them. But what if we want to look at t he covariance between the recordings from various channels of a microelectrode? This is to look at covariance between M channels with N time points, i.e. a N x M matrix. We use the following formula to calculate the covariance matrix, Equation 4 In the above equation, l is used to index the channel and is the mean of that channel. In the numerator, every element in the channel is subtracted by its mean and multiplied by the transpose of the same. These values are then summed and divided by N - 1. This operation gives rise to a square matrix of M x M , where each value in the matrix corresponds to the covariance between t hose two channels. For example, (1,2) will reveal the covariance between channel 1 and channel 2. Also, this covariance matrix is a symmetric matrix; i.e., the values above and below the matrix diagonal are exactly the same. And the diagonal of this matrix is nothing but the covariance with itself which is the variance, so the values at (1,1), STA as the mean. viding the covariance values in the matrix by the product of its corresponding standard deviations. The equation to be used is as follows: 39 Equation 5 This m atrix, like the covariance matrix, is a square and symmetric matrix, where the diagonal values are all ones because we are dividing the variance by variance (product of standard deviations of the channels we are looking at the correlation of). This matrix, in addition to the direction of linear relationship, will also give the strength of the relationship between these channels. While the downstream p urpose of exploring the use of covariance and correlation matrices is to apply it on LFP - STA to reveal hidde n effects of the tissue response on the extracellular signal, looking at the simpler example of the relationship between channels is provided as a s tarting point to illustrate what a covariance and correlation matrix can reveal about the data. As such, I h ave created the covariance and correlation matrix for data from a 16 channel, Michigan - style single - shank probe one day after insertion into the rig ht hemisphere of the primary motor cortex of a rat. Each recording consists of 917248 samples collected at a sampling rate of 48828 Hz [2]. Below are 2 images which show a color plot of a covariance and correlation matrix across the 16 channels: 40 Figure 1 6 . Covariance matrix for a microelectrode with 16 channels Figure 17 . Correlation matrix for the covarianc e matrix in Figure 16 The result is a 16 x 16 symmetric matrix, where each element in the covariance matrix is a covariance between the corresponding channel. As the color bar shows, the highest values are represented by yellow and lowest by blue. Looking at the covariance matr ix, we see that there is an increase and decrease in covariance across a group of channels at the same time. This might be showing the oscillatory behavior of the LFP, i.e. the up and down cycles of the LFP related to excitation and i nhibition in the netwo rk. Inspection of the correlation matrix confirms that the diagonal is all yellow (i.e., +1), as expected. The remaining values between 41 (+1) and ( - 1) indicate how strong the correlation is between pairs of channels and not just their linear relationship. We see that the correlation values are high for several channels, which makes sense because the LFP recorded on nearby electrodes is expected to be highly correlated. 1.2. Spike triggered Covariance and Correlation Building on how covarian ce and correlation matrices can be used together as tools, we can apply these methods to the LFP - STA. The spike triggered covariance matrix is a covariance matrix constructed by using the LFP snippets ( ) that were averaged to obtain the LFP - STA as elements of a matrix, where the mean is the STA itself and N is the number of spikes (or the number of snippets). Using E quation 4 for obtaining the covariance matrix, where every element is a LFP - snippet and the mean is the STA, is shown below: Equation 6 Having calculated the LFP - STC, we can create the correlation matrix by using E quation 4. Here, the covaria nce matrix is replaced with the matrix obtained in E quation 4, and the product of the standard deviation will be the standard deviations of the corresponding LFP - snippets used: Equation 7 42 2. Results As proof - of - principle, I have calculated the LFP - STC and correlation matrix for the data from a single channel of a 16 channel Michigan - style single - shank probe one day after insertion into the rig ht hemisphere of the primary motor cortex of a rat. The total recording time was ~300 seconds within which 195 spikes were recorded. The data used was down sampled, and the new sampling rate is 3051.8 Hz for the LFP. Figure 18 . The LFP Spike triggered a v erage for the mentioned animal 43 Figure 19 . Overlap of all the 194 LFP snippets used for construction of STA Figure 20 . LFP - STA Covariance Matrix 44 Figure 21 . A contoured version of the LFP - STA covariance matrix in Fig ure 24 Figure 22 . LFP - STA Corre lation Matrix 45 Figure 23 . Contoured version of the LFP - STA Correlation matrix in Fig. 2 6 2.1. Discussion and i nterpretation In the covariance matrix, every LFP - snippet is subtracted by the STA(mean), thus reflecting how different the LFP snippet is from the ST A. In other words, the covariance matrix illustrates the similarity between the mean subtracted waveforms to each other (how similar they are to each other in terms of their difference from the STA). There can be three kinds of waveforms after the STA is s ubtracted from the snippet. In order to illustrate key concepts, Fig ure 24 shows artificially generated , phase - shifted s ine waves to represent the three scenarios which may result in different covariance values. First, the snippet may have close to a zero phase shift ( Figure 24 (A) ) and the mean subtraction results in a waveform which has a value close to zero across time. W ith these kinds of waveforms, the covariance will be higher with the snippets which have a similar relation with the STA (i.e. ~ 0 phase 46 shift). The second scenario is when there is a random phase shift with the STA ( Figure 2 4 (B) ). In this case, the covari ance will be high between the snippets which have similar phase shift between them (i.e., ~ x phase shift). The third type is when the s nippet is almost completely out of phase with the STA ( Figure 2 4 (C) ), which also results in a waveform which has values waveform, but it will still show a high covariance, albeit with an opposing effect on the formation of STA waveform. Figure 24 . Different scenarios that can exist when covariance is c alculated (A) close to zero phase shift,(B) A random phase shift x, (C) Close to 180 degrees phase shift 47 1 2 3 1 0.0801 - 0.155 0.15 9 2 - 0.155 2.726 3.986 3 0.159 3.986 7.923 Table 3 . Covariance Matrix for cases in Fig. 24 Table 3 shows the covariance matrix for the cases ( A ),( B ) and ( C ), where (1,1) is the variance of ( A ), (2,2) is the variance of ( B ) and ( 3,3) is the variance of ( C ). And we see that the value at ( A , A ) has the least covariance which makes sense because it has close to zero phase shift, it increases at (2,2) and the variance is highest at (3,3) which also makes sense because the phase shift is close to 180 degrees. But when we look at the value at (1,3) it shows posit ive covariance, even though their contribution to the STA is the opposite. In practice, when the covariance matrix of real data is calculated, the cases will not be as extreme. The diagonal of the covariance matrix becomes useful for interpretation (which, of variance mean for the STA. The higher the variance, the more different the sn ippet is from the STA waveform. For example in Fig. 6, between (100 - 105, 100 - 10 5), we can see a yellow portion telling us that those LFP - snippets have a very different shape from the STA, similarly between (65 - 70,65 - 70), we can see green pixels again showi ng that there is a significant variance when compared to the STA. At all the ot her places in the diagonal we can see that the pixels are blue telling us that the variance is low and they must be positively affecting the shape of the STA. Fig ure 25 plots th e snippets at those spike indexes, which indeed show high variance in comparison with the STA. 48 Figure 25 . (A) LFP snippets from 66 to 70 with LFP - STA (B). LFP snippets from 101 to 105 with LFP - STA Alternatively, in Fig ure 26 , we lo ok at snippets that h ave low variance with the STA. From the covariance matrix, we can identify regions of low variance between (40 - 50,40 - 50) and between (78 - 82,78 - 82). In Fig ure 26 , we can see that at 0 msec (the spike time) all of the LFP snippets follow the trajectory of th e STA, hence further reinforcing the interpretation of high or low variance values can be interpreted in terms of the similarity of individual LFP snippets with the STA. (A) ( B ) 49 Figure 26 . (A) LFP snippets from 42 to 48 with LFP - STA (B). LFP snippets from 78 t o 82 with LFP - STA Now that we know which values have a positive and negative variance with respect to the STA, a correlation matrix will standardize the covariance matrix, revealing the strength of linear relationship between pairs of snippets. In Fig. 21 , the color plot in dicates which values which had a positive relation with the STA. For example, in the range of 78 - 82, green and dark green colors indicate that those regions have relatively high correlation; these snippets represent a correlation that pos itively affects th e STA. By analyzing the covariance and correlation matrix, we may reveal underlying structure in the data beyond the simple LFP spike - triggered average. In turn, these tools may be useful for identifying otherwise hidden manifestations of structural or fun ctional remodeling in local neurons surrounding devices in recorded electrophysiological data. Possible interpretations may include: (A) ( B ) 50 Sources of Spatial Variation . One interpretation can be that the waveforms representing high variance are spikes from a dif ferent neuron but they have less connectivity with the LFP that i s being recorded by the electrode , or it is farther away from the LFP than the one that is influencing the LFP more dominantly. Once we separate the values of the diagonal like this, and loo k at the values that show high correlation in those rows and columns in the correlation matrix, they can be interpreted as snippets which might be generated from the same neuron. Sources of Temporal Variation . The second interpretation c an be that recordi ngs are from the same group of neurons, but across time the neuronal group is modulating its decision, such that the major variation in the LFP is similar to the STA. Because STA is the average of all the LFP snippets a larger number of spikes gave rise t o a shape like the STA. In Chapter 1 we discussed how the LFP is affected by different contributors, the fast action potentials (Na+), the long lasting potentials (Ca2+) and direct electrical communication between gap junctions. So, when the spike occurs w e can see a sharp dip in the STA, but we also can see that the STA is affected in range of 2 msec, hinting towards the contribution due to the fast acting potentials. The synchrony might also be affected by the communication through gap j unctions, which ca n increase neuronal synchrony indirectly, affecting the LFP towards the mean behavior captured by the STA. 51 3. Future use of STC as a tool In the future, we can use STC to investigate the interpretation of structural and functional network remodeling surrounding implanted electrode arrays. Furthermore, we can use it to assess levels of variation in the data which may not be immediately eviden t using traditional metrics (spike counts, LFP amplitude, etc.). Likewise, we can apply these techn iques in the context of gene knockdown to reveal the effects on the STC matrix, across different time points and implant regions. 52 Chapter 4 CONCLUSION The goal of the thesis was to create signal processing tools to investigate the structural and functional remodeling effects caused by implantation of a neural prosthesis. We begin this by understanding the signal we will be analyzing, as described in the first chapter. In this chapter , we describe how a neuron generates an action potential, ho w the neuron integrates the various input s it is receiving, how the signal travels from the soma to the axon terminals, and how these things affe ct the shape of the action potential generated. Once we understand how a single neuron generates a signal, we c onsider the activity associated with a broader interconnected population of neurons (the local field potential). We discuss how the shape of the LFP might be affected by different contributors, the geometry of the neuron, the architecture of the arrangemen t of the neurons, and on the nature of the recording array. We then discuss how the brain reacts to the implant and how different features around the implant are affected, and how this may influence the recorded data. The last part of the chapter discusses the challenges not created by the tissue response, but rather the technical challenges at play. After understanding the major features influenci ng the recorded data, we move into the second chapter where we look at data from the LFP and MUA, and how the r elationship between them might reflect the connectivity within the network. Here, we discuss the concept of communication through coherence, and how there might exist a frequency where LFP and MUA lock with each other providing maximum communication. The s econd part of the chapter looks at another signal processing tool known as the spike triggered average (STA), which is conventionally used to stu dy the relationship between spiking and the corresponding stimuli. Here, we use it to look at the synchrony bet ween the spiking and the local field potential. In the last part of the chapter, we 53 discuss how the results from these tools can be interpreted t o reveal the connectivity within the network. In this chapter, we introduce two tools to reveal changes in tiss ue response by observing the electrical activity: MUA - LFP coherence analysis and LFP - STA. In the third chapter, we discuss how the tools discusse d in the previous chapter reveal only the linear relationship between the signals, which is useful but hides th e non - linear features affecting the response. Therefore, in the third chapter we look at a tool which reveals the hidden non - linearities present known as spike triggered covariance and correlation (STC). In the first part of the chapter, we discuss what co variance and correlation mean, building toward calculation of the STC matrix. Next, we discuss how to interpret this matrix and how the features revealed by STC might be responsible in forming the final response we observe in the LFP - STA, and what this mig ht mean in terms of connectivity within the network. The tools created in the previous chapters are expected to provide us with new insights abou t what been us ed on large data sets to reveal strong conclusions as - of - yet, this thesis serves as an introduction of these tools and a proof - of - principle for t heir use to understand the functional and structural connectivity. We expect that these tools will be further i mproved and have more features incorporated in them moving forward. Nonetheless, while we have only used them on a small dataset, it is already e vident that a lot of features may be revealed just through its straightforward/initial application. In the futu re, we can use these tools on different neuronal populations, or the same tools layered over other analyses. The LFP - STA can be used to look at t he tissue response across different time points, and it also can be used along with current source density anal ysis to get a clearer picture of how the LFP is being constructed. Likewise, it also can be used for spike - triggered average of the time frequenc y 54 spectrum of the LFP (STFFA). All of these analyses, as discussed, can be looked at in parallel with the MUA - L FP coherence. We can also look at how the communication between different neuronal populations is affected within the network. STC analysis can b e applied to reveal the non - linearities, which in the case of MUA - LFP coherence might be very informative, and might reveal information about non - linearities in the cortical network in future work. 55 APPENDI CES 56 APPENDIX A: Coherence A nalysis 1. %% reading the required data files: Raw LFP and Sorts 2. file1 = uigetfile('*.mat'); %LFP raw: g ets directory 3. %% Raw data 4. T=load(file1); % Raw Wave 5. b=T.data.Wave.data{1,1}; 6. [u1,u2]=size(b);% get the length of arrays in LFP data 7. B=zeros(u2,16); % empty matrix 8. for k=1:16 9. B(:,k)=T.data.Wave.data{1,(k)};% store LFP data in B 10. end 11. %% Select the data file 12. d ata=B(:,3)% change value of 3 from 1 to 16 depending on which channel you want to see 13. %% Filter MUA 14. S1=filter(MUA,data); % MUA = band pass filter form 300hz - 6000hz 15. SD1=std(S1); 16. if abs(S1(:,:))6e - 06)) 57. % c(j)=j; 58. % c1(j)=i; 59. % subplot(4,4,j) 60. % plot(p); 61. % ti tle(['Channel ',num2str(j)]); 62. % xlabel('Time(msec)') 63. % ylabel('volts') 64. % hold on; 65. % end 66. % end 67. % hold off; 68. % sgtitle('Sorted meanwaveforms') 62 69. % end 70. % 71. %% Obtain LFP of the channels with units in a matrix 72. 73. % czero = c(c~=0); % channel numbers with units 74. czero=16; 75. T=load(file1); % raw LFP 76. % cz=length(czero) 77. cz=16; 78. b=T.data.LFPV.data{1,czero(1)}; 79. [u1,u2]=size(b);% get the length of arrays in LFP data 80. B=zeros(u2,cz); % empty matrix 81. for k=1: cz 82. B(:,k)=T.data.LFPV.data{1 ,czero(k)};% store LFP data in B 83. end 84. 85. %% apply fourier on channels with units, plots for time and fourier 86. 87. fourier1=fft(B); 88. [rawpower,axis1]=viewfft(fourier1,3051.8);% view fft is a separate function 89. figure; 90. for k1=1:cz 91. s ubplot(4,4,k1) 63 92. plot(B(:,k1)) 93. title(['Channel ',num2str(czero(k1))]); 94. xlabel('Time(msec)') 95. ylabel('volts') 96. end 97. sgtitle('Raw LFP') 98. figure; 99. for k1=1:cz 100. subplot(4,4,k1) 101. plot(axis1,rawpower(:,k1)) 102. title(['Channel ',num2str(czero(k1 ))]); 103. xlabe l('Frequency(Hz)') 104. ylabel('V') 105. end 106. sgtitle('Fourier Transform') 107. %% Create a time array for the LFP 108. 109. srate=3051.8% sampling rate 110. S=zeros(u2,cz); % empty matrix 111. for k2=1:length(czero) 112. % S(:,k2)=filter(mid60,B(:,k2));% notch filter 113. S(:,k2) =B(:,k2);% no n otch filter 114. end 64 115. len=(length(b))/3051.8;% total number of samples/sampling rate= time in seconds 116. K=1/srate;% distance between each sample= 1/3051.8 117. time=1:K:(len+1);% first array for time 118. Time=time(1:(length(time) - 1));% final array for time 119. % % Extract spike information 120. 121. % c1zero = c1(c1~=0); % MWF waveform to be used 122. P = A.sorts.sniptime{czero(1),c1zero(1)};% spike data 123. g=ones(1,length(P));% time axis for spike array 124. q=P(40:end - 40); %excluding the first few spike so th at the number of sample s is >6000 125. u=length(q); % number of spikes 126. 127. %% create window and calculate STA 128. 129. avg=zeros(1,12001); 130. Finalavg=zeros(1,12001); 131. S1=S'; 132. S11=S1(1,:); 133. 134. e11=zeros(u,12001); 135. 136. for w=1:u 137. f=q(w); 65 138. [~,loc]=min(abs(Time - f)); 139. e=(loc - 6000):(loc+60 00); 140. e11(w,:)=S11(e); 141. avg=avg+S11(e); 142. end 143. 144. Finalavg(1,:)=avg/u; 145. 146. %% Time axis for STA: spike occurs at 2.966 msec 147. lensta=(length(Finalavg (1,:)))/3051.8;% total number of samples/sampling rate= time in seconds 148. K1=1/srate;% distance between each sample= 1/3051.8 149. timesta=1:K1:(lensta+1);% first array for time 150. Timesta=time(1:(length(timesta) - 1));% final array for time 151. Timestanorm=Timesta - 2.96 6; 152. %% plot the STA :spike occurs at 2.966 msec 153. 154. figure; 155. plot(Timestanorm,Finalavg); 156. title('Spike Trigg ered Average'); 157. xlabel('Time(msec)') 158. ylabel('V') 159. %% STA snippets 66 160. figure; 161. for d1=1:u 162. plot(Timestanorm,e11(d1,:)); 163. title('Snippet Overlap'); 164. xlabel('Time(msec)') 165. ylabel('V') 166. hold on 167. end 168. hold off; 67 APPENDIX C: Spike trigg e red covariance and correlation 1. %% reading the required data files: Raw LFP and Sorts 2. %% reading the required data files: Raw LFP and Sorts 3. 4. file1 = uigetfile('*.mat'); %LFP raw: gets directory 5. file2 = uigetfile('*.mat'); %Sorts file: gets directory 6. 7. %% Ob tain LFP of the channels with units in a matrix 8. 9. T=load(file1); % Raw Wave 10. b=T.data.LFPV.data{1,1}; 11. [u1,u2]=size(b);% get the length of array s in LFP data 12. B=zeros(u2,16); % empty matrix 13. for k=1:16 14. B(:,k)=T.data.LFPV.data{1,(k)};% store LFP data in B 15. end 16. %% LFP channel to use 17. S5=B(:,3);% can cange value of 3 to 1 to 16 depending on the channel you want to look at 18. %% Time array for MUA 19. srate=3051.8% sampling rate 20. len=(length(S5))/3051.8;% total number of samples/sampling rate= time in seconds 21. K=1/srate;% dist ance between each sample= 1/500 22. time=0:K: (len);% first array for time 68 23. Time=time(1:(length(time) - 1));% final array for time 24. %% Extract Spike information 25. A=load(file2); %load sorts in A 26. P = A.sorts.sniptime{3,1};% spike data 27. g=ones(1,length(P));% time axis f or spike array 28. q=P(40:end - 40); %excludin g the first few spike so that the number of samples is > 6000 29. u=length(q); 30. %% plot spike and LFP together 31. gg=ones(1,length(g)); 32. m1=max(S5); 33. figure; 34. stem(q,gg(40:end - 40)); 35. hold on; 36. plot(Time,100*S5); 37. hold off; 38. %% Sta begins 39. avg=zeros(1,60001); 40. Finalavg3=zer os(1,60001); 41. S=S5'; 42. S11=S(1,:); 43. e13=zeros(u,60001); 44. 45. %% time window for sta 69 46. 47. for w=1:u 48. f=q(w); 49. [~,loc]=min(abs(Time - f)); 50. e=(loc - 30000):(loc+30000); 51. e13(w,:)=S11(e); 52. avg=avg+S11(e); 53. end 54. 55. Finalavg3=avg/u; 56. 57. %% %% Time axis for STA: spike occurs at 2.966 msec 58. lensta=(length(Finalavg3))/3051.8;% total number of samples/sampling rate= time in seconds 59. K1=1/srate;% distance between each sample= 1/500 60. timesta=0:K1:(lensta);% first array for time 61. Timesta=time(1:(len gth(timesta) - 1));% final array for time 62. Timestanorm3=Timesta - (1.965*5); 63. 64. %% STA plot 65. figure; 66. plot(Timestanorm3,Finalavg3); 67. title('Spike Triggered Average'); 68. xlabel('Time(msec)') 70 69. ylabel('V') 70. 71. %% Input the spike indexes you want to co mapre 72. % figure; 73. % subplot(2,1,1) 74. % plot(Timestanorm3',e13(66,:)); 75. % hold on; 76. % % plot(e13(66,:)); 77. % plot(Timestanorm3',e13(67,:)); 78. % plot(Timestanorm3',e13(68,:)); 79. % plot(Timestanorm3',e13(69,:)); 80. % plot(Timestanor m3',e13(70,:)); 81. % plot(Timestanorm3',Fina lavg3,'color','k','LineWidth',1.2); 82. % title('a'); 83. % xlabel('time(msec)') 84. % ylabel('V') 85. % legend('66','67','68','69','70','STA'); 86. % hold off; 87. % subplot(2,1,2) 88. % plot(Timestanorm3',e13(101,:)); 89. % hold on; 90. % plot(Times tanorm3',e13(102,:)); 91. % plot(Timestanorm3 ',e13(103,:)); 71 92. % plot(Timestanorm3',e13(104,:)); 93. % plot(Timestanorm3',e13(105,:)); 94. % plot(Timestanorm3',Finalavg3,'color','k','LineWidth',1.2); 95. % title('b') 96. % xlabel('time(msec)') 97. % ylabel('V') 98. % legend('78','79','80','81','82','STA'); 99. % hold off; 100. % 101. %% sn ippet Sta 102. figure; 103. for d1=1:u 104. plot(Timestanorm3,e13(d1,:)); 105. title('Snippet Overlap'); 106. xlabel('Time(msec)') 107. ylabel('V') 108. hold on 109. end 110. hold off; 111. 112. %% Covariance Begins 113. 114. % subtracting mean/STA from matrix 72 115. rawcov=e13'; % elements cross channels 116. meanrawcov=bsxfun(@minus,rawcov,Finalavg3'); % data - mean: step 1 117. 118. %% step 2: multiplying matrices 119. C=length(Finalavg3) - 1; 120. stacovmat=(meanrawcov'*meanrawcov)/C; 121. figure; 122. imagesc(stacovmat); 123. title('C ovariance Matrix') 124. xlabel('LFP Snippet number'); 125. ylabel('LFP Snippet number'); 126. figure; 127. contourf(stacovmat); 128. title('Contour Covariance Matrix') 129. xlabel('LFP Snippet number'); 130. ylabel('LFP Snippet number'); 131. %% variance 132. % roo=diag(stacovmat); 133. % figure; 134. % plot(q,roo); 135. 136. %% Step 3: correlation matrix 137. newcorrdata=bsxfun(@rd ivide,meanrawcov,(std(meanrawcov,[],1))); 73 138. corr1=(newcorrdata'*newcorrdata)/C; 139. figure; 140. imagesc(corr1); 141. title('Correlation Matrix') 142. xlabel('LFP Snippet number'); 143. ylabel('LFP Snippet number'); 144. figure; 145. contourf(corr1); 146. title('Contour Correlation Matrix') 147. xlabe l('LFP Snippet number'); 148. ylabel('LFP Snippet number'); 74 REFERENCES 75 REFERENCES 1. Functional remodeling of subtype - specific markers surrounding implanted neuroprostheses, JW Salatino, BM Winter, MH Drazin, EK Purcell. Journal of n europhysiology, 118(1):194 - 202 , 2017. 2. Alterations in ion channel expression surrounding implanted microelectrode arrays in the B rain , JW Salatino, AP Kale, EK Purcell. bioRxiv, 518811 , 2019. 3. Genetic Modulation at the Neural Microelectrode Interface: Metho ds and Applications , Winter BM, Daniels SR, Salatino JW, Purcell EK . Micromachines (Basel) , 9(10):476 , 2018. 4. Measurement of current - voltage relations in the membrane of the giant axon of Loligo , A. L. Hodgkin, A. F. Huxley and B. L Katz. Journal of Physi ology, 116(4), 424 - 448, 1952. 5. Currents carried by sodium and potassium ions through the memb rane of the giant axon of Loligo , A. L. Hodgkin and A. F. Huxley. Journal of Physiology. 116(4), 449 - 472, I952. 6. T he components of membrane conductance in the giant axon of L oligo , A . L . H odgkin and A . F . H uxley , J ournal of Physiol ogy . 11 6 (4) , 473 - 496 , 1952. 7. The dual effect of membrane potential on sodium conductance in the giant axon of Loligo , A. L. Hodgkin and A. F. Huxley. Journal of Physiology, 116(4): 497 506 , 1952. 8. A quantitative description of membrane current and its application to con duction and excitation in nerve , A . L . H odgkin and A . F . H uxley . J ournal of Physiol ogy, 11 7 (4) , 500 - 544 , 1952 . 9. Cellular physiology and neurophysiology, Blaustein, Kao, and Ma tteson . 2 nd Edition, 2012. 10. The action potential in mammalian central neurons, Bruce P. Bean . Nature Reviews Neuroscience, 8(6):451 - 465 , 2007 . 11. The Voltage Sensor in Voltage - Dependent Ion Channels , Francisco Bezanilla . P hysiological Reviews Vol. 80 ( 2 ): 55 5 - 592 , 2000 . 12. The origin of extracellular fields and currents -- EEG, ECoG, LFP and spikes. Buzsáki G, Anastassiou CA, Koch C. Nature Reviews Neuroscience, 13(6):407 - 420, 2012. 13. https://neupsykey.com/neurophysiologic - basis - of - the - electroencephalogram - 2/ 14. Nucle us basalis and thalamic control of neocortical activity in the freely moving rat , Buzsáki, G. et al. J ournal of Neuroscience, 8 : 4007 4026 , 1988. 76 15. Sharp wave - associated high - f requency oscillation (200 Hz) in the intact hippocampus: network and intracellula r mechanisms , A Ylinen, A Bragin, Z Nadasdy, G Jando, I Szabo, A Sik and G Buzsaki . J ournal of Neurosci ence, 15 : 30 46 , 1995. 16. Submillisecond synchronization of fast electrical oscillations in neocortex , Barth, D. S. J ournal of Neurosci ence . 23 : 2502 2510 , 2003. 17. A study of nerve physiology , Lorente de Nó, R. Studies from the Rockefeller Institute for Medical Research Part I, 131 , 1947. 18. Electr oencephalography: Basic Principles, Clinical Applications, And Related Fields 5th edn , Niedermayer, Donald L. Scho mer & F.H. Lopes da Silva , 2012. 19. The brain in fractal time: 1/f - like power spectrum scaling of the human electroencephalogram , Pritchard W. S . Int ernational J ournal of Neurosci ence, 66 : 119 129 , 1992. 20. Estimation of population firing rates and current sour ce densities from laminar electrode recordings , Pettersen, K. H., Hagen, E. & Einevoll, G. T. J ournal of Comput ational Neurosci ence, 24 : 291 313 , 2008. 21. Synaptic mechanisms of synchronized gamma oscillations in inhibitory interneuron networks , Bartos, M., Vida, I. & Jonas . P. Nature Rev iews Neurosci ence . 8 : 45 56 , 2007. 22. Synaptic background activity influences spatiotemporal integration in s ingle pyramidal cells , Bernander, O. Douglas, R. J., Martin, K. A. C. & Koch, C. Proc. Natl Acad. Sci. USA , 88 : 11569 11573 , 1991. 23. Theory of current source - density analysis and determination of conductivity tensor for anuran cerebellum, Nicholson, C. & Freeman, J. A. Journal of Neurophysiology, 38 : 356 368, 1975. 24. Regenerative Electrode Interfaces for Neural Prostheses, C ort H. Thompson, Marissa J. Zoratti, Nicholas B. Langhals, and Erin K. Purcell. Tissue Engineering Part B: Reviews , 22 : 125 - 135 , 2016. 25. Fully integrated silicon probes for high - density recording of neural activity, Jun, J., Steinmetz, N., Siegle, J. et al. Nature 551 : 232 236, 2017. 26. Recent advances in silicon - based neural microelectrodes and microsystems: a review, Zoltan Fekete. Sensors and Actuators B: Chemical , 215 : 300 - 315 , 2015 . 27. Restoring the sense of touch with a prosthetic hand through a brain int erface, Tabot GA, Dammann JF, Berg JA, Tenore FV, Boback JL, Vo - gelstein RJ, Bensmaia SJ. Proc Natl Acad Sci U S A , 110:18279 - 18284.43, 2013. 28. Behavioral assessment of sensitivity to intracortical microstimulation of primate somatosensory cortex , Kim S, Cal lier T, Tabot GA, Gaunt RA, Tenore FV, Bensmaia SJ . Proc . 77 Natl Acad Sci U S A , 112:15202 - 15207 , 2015 . 29. Neural populat ion dynamics in human motor cortex during movements in people with ALS , Pandarinath C, Gilja V, Blabe CH, Nuyujukian P, Sarma AA, S orice BL, Eskandar EN, Hochberg LR, Henderson JM, Shenoy KV , 4:e07436 , 2015. 30. Output properties of the cortical hindlimb motor area in spinal cord - injured rats , Frost SB, Dunham CL, Barbay S, Krizsan - Agbas D, Winter MK, Guggenmos DJ, Nudo RJ . J ournal of Neurotrauma , 32:1666 - 1673 , 2015. 31. Restoration of function after brain damage using a neural prosthesis , Guggenmos DJ, Azin M, Barbay S, Mahnken JD, Dunham C, Mohseni P, Nudo RJ . Proc Natl Acad Sci U S A , 110:21177 - 21182 , 2013 . 32. Neural Prosthetic Systems: Current Problems and Future Directions , Cindy A. Chestek, John P. Cunningham, Vikash Gilja, Paul Nuyujukian, Stephen I. Ryu, Krishna V. Shenoy, IEEE EMBS, 3369 - 3375 , 2009 . 33. Neuronal cell loss accompanies the brain tissue r esponse to chronically implanted silicon microelectrode arrays, Biran, R ., Martin, D. C. & Tresco, P. A. Experimental Neurolo g y , 195 ( 1 ): 115 - 126 , 2005. 34. Glial responses to implanted electrodes in the brain JW Salatino, KA Ludwig, TDY Kozai, EK Purcell . Nat ure biomedical engineering, 1(11):862 - 877 2017 . 35. Microthalamotomy effect during deep brain stimulation: Potential involvement of adenosine and glutamate efflux , Y. Shon, Y. M., Agnesi, F. & Lee, K. H. Annual International Conference of the IEEE Engineering in Medicine and Biology Society 3294 3297, 2009 . 36. Glial responses to implanted electrodes in the brain , JW Salatino, KA Ludwig, TDY Kozai, EK Purcell . Nature biomedical engineering, 1: 862 877 , 2017 . 37. Astrocyte roles in traumatic brain injury, Burda , J. E. , Bernstein, A. M. & Sofroniew, M. V. Experimental Neurology, 275 : 305 315, 2016. 38. Emerging role for astroglial networks in information processing: from synapse to behavior, Pannasch, U. & Rouach, N. Trends in Neuroscience, 36 : 405 17, 2013. 39. The up - regula tion of voltage - gated sodium channel nav1.6 expression following fluid percussion traumatic brain injury in rats Mao, Q. et al. Neurosurgery , 66 : 1134 1139 , 2010. 40. Perspective Multi - scale, multi - modal analysis uncovers complex relationship at the brain tiss ue - implant neural interface: new emphasis on the biological interface, Nicholas J Michelson1, Alberto L Vazquez, James R Eles, Joseph W Salatino, Erin K Purcell, J ordan J Williams , X Tracy Cui and Takashi D Y Kozai. Journal of Neural Engineering , 15 (3): 78 033001 , 2018. 41. The relationship between glial cell mechanosensitivity and foreign body reactions in the central nervous system, Pouria Moshayedi , Gilbert Ng, Jessi ca C.F. Kwok, Giles S.H. Yeo, Clare E. Bryant, James W. Fawcett, Kristian Franze, Jochen Guck . Biomaterials Volume 35 (13): 3919 - 3925 , 2014 . 42. Mechanical failure modes of chronically implanted planar silicon - based neural probes for laminar recording , Takashi D.Y. Kozai, Kasey Catt a, Xia Li, Zhannetta V. Gugel, Valur T. Olafsson . Biomaterials, 37: 25 - 39, 2015. 43. Modelling and analysis of local field potentials for studying the function of cortical circuits, Einevoll, Gaute T. et al. Nature Reviews Neuroscienc e 14.11: 770 - 785 , 2013. 44. The gamma cycle , s in neurosciences 30. \ ( 7 ) : 309 - 316 , 2007 45. Neuronal oscillations in cortical networks , Buzsáki, G., & Draguhn, A. Science, 304(5679), 1926 - 1929 , 2004. 46. Modulation of Oscillatory Neuronal Synchronization by Selective Visual Attention, Pascal Fries1, John H. Reynolds, Alan E. Rorie1, Robert Desimone. Science , 291 ( 5508 ): 1560 - 1563, 2001. 47. The basis of sensation , Adrian, E. D , W W Norton & Co. 1928. 48. Nonlinear a nalysis and synthesis of receptive - field responses in the catfish retina. I. Horizontal cell leads to ganglion cell chain , Z Marmarelis , and K I Naka. Journal of Neurophysiology 36 ( 4 ): 605 - 618 , 1973. 49. Triggered Correlation, E. De Boer and P. Kuyper . IEEE Transactions on Biomedical Engineering,15 (3): 169 - 179 , 1968. 50. , H L Bryant ,J P Segundo . The Journal of physiology , 260(2):279 - 314, 1976 . 51. The Spike - Triggered Average of the Integrate - and - Fire Cell Driven by Gaussian White Noise, Liam Paninski. Neural Computation , 18 ( 11 ): 2592 - 2616 , 2006. 52. Dependence of the spike - triggered average voltage on membrane response properties, Laurent Badel, Wulfram Gerstner, Magnus J.E. Richardson. Neurocomputing , 69 : 1062 1065, 2006. 53. Network Rhythms Influence the Rela tionship between Spike - Triggered Local Field Potential and Functional Connectivity, Supratim Ray and John H. R. Maunsell . The Journal of Neuroscience, 31(35): 12674 - 12682, 2011 . 79 54. Unsupervised Spike Detection and Sorting with Wavelets and Superparamagnetic Clustering , R. Quian Quiroga, Z. Nadasdy, and Y. Ben - Shaul . Neural Computation , 16 ( 8 ): 1661 - 1687 , 2004. 55. A review of methods for spike sorting: the dete ction and classification of neural action potentials , Michael S. Lewicki . Medicine, Computer Science Netw ork , 9(4):53 - 78, 1998. 56. Predicting Movement from Multiunit Activity, Eran Stark and Moshe Abeles . Journal of Neuroscience, 27 (31) : 8387 - 8394 , 2007. 57. Rh ythms For Cognition: Communication Through Coherence, Pascal Fries. Neuron, Volume 88 ( 1 ): 220 235, 2015 . 58. Comparisons of the Dynamics of Local Field Potential and Multiunit Activity Signals in Macaque Visual Cortex. Samuel P. Burns, Dajun Xing, and Robert M. Shapley . Journal of Neuroscience, 30(41): 13739 - 13749, 2010. 59. Neural decoding, Naureen Ghani February 2018 . 60. Spike - triggered average electrical stimuli as input filters for bionic vision a perspective, D L Rathbun, N Ghorbani1, H Shabani1, E Zrenner, and Z Hosseinzadeh,Liu, J.K., Schreyer, H.M., Onken, A. et al. Journal of Neural Engineerin g, 15 ( 6 ) : 06300 2 , 2018. 61. Spike - triggered neural characterization Odelia Schwartz; Jonathan W. Pillow; Nicole C. Rust; Eero P. Simoncelli . Journal of Vision, 6 ( 13 ): 484 - 507 , 2016. 62. Flavopiridol reduces the impedance of neural prostheses in vivo without affec ting recording quality, Erin K Purcell, David E Thompson, Kip A Ludwig, Daryl R Kipke . Journal of Neuroscience Methods, 183(2):149 - 157 , 2009. 63. Inference of neuronal functional circuitry with spike - triggered non - negative matrix factorization, Liu, J.K., Schr eyer, H.M., Onk en, A. et al. Nature Communications, 8 ( 149 ) , 2017. 64. Laminar population analysis: estimating firing rates and evoked synaptic activity from multielectrode recordings in rat barrel cortex, Einevoll, G. T. et al. J. Neurophysiol. 97 : 2174 2190 , 2007. 65. Effect of stimulus intensity on the spike local field potential relationship in the secondary somatosensory cortex , Ray, Supratim, et al. Journal of Neuroscience 28 ( 29 ) : 7334 - 7343 , 2008 . 66. Characterizing Neural Gain Control using Spike - triggered Covariance, Odeli a Schwartz, E.J. Chichilnisky, Eero P. Simoncelli. Advances in Neural Information Processing Systems , 14 : 269 - 276, 2001. 80 67. Spike - Triggered Covariance Analysis Reveals Phenomenological Diversity of Contrast Adaptation in the Retina, Jian K. Liu, Tim Gollisch . PLoS Computational Biology, 11(7) : e100442 , 2015. 68. Spike - triggered characterization of excitatory and suppressive sti mulus dimensions in monkey V1, Nicole C. Rusta , Odelia Schwartz, J. Anthony, Movshona, Eero Simoncellia, Neurocomputing , 58 ( 60 ): 793 - 7 99 , 2004.