University of Southern California Dissertations and Theses

Nonlinear modeling of causal interrelationships in neuronal ensembles: an application to the rat hippocampus

(USC Thesis Other)

Nonlinear modeling of causal interrelationships in neuronal ensembles: an application to the rat hippocampus

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content NONLINEAR MODELING OF CAUSAL INTERRELATIONSHIPS IN NEURONAL ENSEMBLES: AN APPLICATION TO THE RAT HIPPOCAMPUS by Theodoros Zanos A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (BIOMEDICAL ENGINEERING) May 2009 Copyright 2009 Theodoros Zanos ii DEDICATION to my brother Stavros, my mentor and companion in science and life to my father Panagiotis and mother Vasiliki my ever-present support and guiding light to my advisor, Vasilis my role model and inspiration in science iii ACKNOWLEDGMENTS Throughout my Doctoral Degree undertaking, I was fortunate enough to interact with brilliant scientists and wonderful people. Their sage advice was and will be my guiding light to my academic journey. I would like to offer my deepest gratitude to my brother, Stavros, for his support in difficult times, his endless patience and enlightening advice that shaped a huge part of who I am, as a scientist but also as a person. I am also extremely grateful to my Ph.D. advisor, Vasilis Marmarelis, who provided generous support throughout my years at USC. His guidance and contribution to my work can only be matched by the inspiration he instilled in me to work in the field of Systems Identification. His pioneering work in the field, coupled with his experience in mentoring students for more than 25 years, provided the perfect environment for a unique learning experience and professional growth. I must also acknowledge the constructive guidance and support of my co-advisor, Dr Theodore Berger and all the post-docs and students at the Berger Lab. I was fortunate to interact and be assisted by Dr Spiros Courellis, with whom I collaborated extensively throughout my first three years of my Ph.D. I would also like to thank Dr Stefan Schaal for his thoughtful suggestions as the outside member of my Ph.D. committee and Dr David D’Argenio, Dr Michael Khoo, Dr Norberto Grzywacz, Dr Bartlett Mel for creating an excellent learning environment at the Biomedical Engineering Department at USC. I would also like to thank the students that collaborated and shared the same working space with me at the Biomedical Simulations Resource, Michalis Markakis, Gopal Erinjippurath, Tushar and Roheed Shaney. I also received valuable help by the BMSR Administrative Coordinator, Marcos Briano, as well iv as the Graduate Student Advisors at the BME Department, Adam Wyatt and Mischalgrace Diasanta. Last, but not least, my academic journey would not have been possible without the irreplaceable and never-ending support and love of my parents, Panagiotis and Vasiliki. v TABLE OF CONTENTS DEDICATION .................................................................................................................... ii ACKNOWLEDGMENTS ................................................................................................. iii TABLE OF CONTENTS .................................................................................................... v LIST OF TABLES ............................................................................................................ vii LIST OF FIGURES ......................................................................................................... viii ABSTRACT ..................................................................................................................... xiii CHAPTER ONE: INTRODUCTION ................................................................................. 1 CHAPTER TWO: NEURONAL SYSTEMS MODELING BACKGROUND.................. 5 Neuronal Modeling Approaches ..................................................................................... 5 Functional Connectivity Studies in Multi-Unit Neural Recordings ............................. 10 Models of CA3 and CA1 of the Hippocampus ............................................................. 13 CHAPTER THREE: NONLINEAR NONPARAMETRIC MULTI INPUT - MULTI OUTPUT MODELS ............................................................................................ 15 Modeling Framework ................................................................................................... 16 Multi-input / Single-Output Modules ....................................................................... 16 Kernel Estimation with the Laguerre Expansion Method ......................................... 18 Selection of Threshold and Evaluation of Model Performance Using ROC Curves ....................................................................................................................... 20 Model Statistics ......................................................................................................... 21 Selection of Significant Inputs .................................................................................. 22 Selection of Model Order and Number of Laguerre Functions ................................ 25 Principal Dynamic Modes......................................................................................... 26 Kernel Interpretation ................................................................................................. 26 Introducing Multi-input/Multi-output Models with an Example .................................. 27 Single-Input/Single-Output Case .............................................................................. 28 Two-Input/Single-Output Case ................................................................................. 29 Three-Input/Single-Output Case ............................................................................... 32 Multi-Input/Multi-Output Simulation Example ............................................................ 33 vi Noise-Free Case ........................................................................................................ 35 Noisy Case ................................................................................................................ 37 Spike Jitter Case ........................................................................................................ 41 Deleted Spikes Case .................................................................................................. 42 Misassigned Spikes Case .......................................................................................... 43 Selection of Model order and Number of Laguerre Functions ................................. 43 CHAPTER FOUR: APPLICATION OF MIMO MODELING TO THE HIPPOCAMPUS............................................................................................................... 45 The Cortical Neuroprosthetic Platform ......................................................................... 45 Data Acquisition ........................................................................................................... 47 Experimental setup and behavioral training ............................................................. 47 Multineuron recording techniques ............................................................................ 49 MIMO Modeling Results .......................................................................................... 50 Functional Connectivity Results ................................................................................... 61 Principal Dynamic Modes Results ................................................................................ 73 CHAPTER FIVE: CONCLUSIONS ................................................................................ 79 ALPHABETIZED BIBLIOGRAPHY .............................................................................. 86 APPENDICES .................................................................................................................. 94 Appendix A: Kernel Computation using Laguerre Basis Functions ............................ 94 Appendix B: Principal Dynamic Modes ....................................................................... 98 Appendix C: ROC Curve / Mann – Whitney Statistic ................................................ 102 vii LIST OF TABLES Table 2.1: Results of MWS-based test for input significance in the simulated noise-free case ................................................................................................................... 35 Table 2.2: Results of MWS-based test for input significance in the simulated 25% added noise case ................................................................................................................ 38 Table 2.3: Results of MWS-based test for input significance in the simulated 50% added noise case ................................................................................................................ 39 Table 2.4: Results of MWS-based test for input significance in the simulated jitter case .................................................................................................................................... 41 Table 2.5: Results of MWS-based test for input significance in the simulated deleted spikes case ............................................................................................................ 43 Table 2.6: Results of MWS-based test for input significance in the misassigned spikes case ......................................................................................................................... 43 Table 2.7: Results of MWS-based test for model order selection .................................... 44 Table 2.8: Results of MWS-based test for number of Laguerre functions selection ........ 44 Table 3.1: Input Selection for the four behavioral tasks ................................................... 52 Table 3.2 : Mann-Whitney Statistics for final MIMO model for LNM, RNM and RS tasks ............................................................................................................................. 56 Table 3.3: Averaged Mann-Whitney Statistics for MISO models estimated from randomly generated independent Poisson spike trains for different MFRs and number of inputs ............................................................................................................... 57 viii LIST OF FIGURES Figure 2.1: (A) Schematic representation of the decomposition of the Multi- input/multi-output (MIMO) model into an array of multi-input/single-output (MISO) modules. (B) Schematic detail of the structure of the MISO module comprising the cascade of the Nonlinear Volterra Transformation (NVT) and the Threshold Trigger (TT) operators. .................................................................................... 16 Figure 2.2: Single-input/single-output model (right hippocampus): (A) spatial location of the input (CA3 neuron) and the output (CA1 neuron) on the hippocampal cross-section, (B) first, second, and third order kernel, (C) ROC curve, (D) Actual response vs. model prediction (segment). ............................................ 29 Figure 2.3: Two-input / single-output model. (A) Computed first-, second-, and third-order self- and cross-kernels. (B) Model prediction with and without cross- kernels. .............................................................................................................................. 31 Figure 2.4: Comparison of the model prediction between a single-input model (A), a two-input model (B), and a three input model (C); the corresponding recorded response (D) and the spatial location of the input and output neurons on the rat hippocampus (E). ................................................................................................... 32 Figure 2.5: The kernels of the second-order simulated MISO system with four inputs and one output: the first column shows the first-order kernels, the second column shows the second-order self-kernels, and the third column shows the cross-kernel between the 1 st and the 4th input (see text). ................................................. 34 Figure 2.6: The estimated kernels of this second-order model from the simulated data for the noise free case, along with ROC curves and theta values of in-sample and out-of-sample prediction (see text). The x axes in the first-order kernels and the x and y axes in the second-order kernels range from 0 to 1000msec, while the ROC curves have in both axes values from 0 to 1. ........................................................... 37 Figure 2.7: (A) The estimated model for the case of 25% spurious spikes in the inputs and output. (B) The estimated model in the case of 50% spurious spikes in the inputs and output. The computed ROC curves are shown in the right compartment of the model for the “in-sample” and “out-of-sample” cases. .................... 40 ix Figure 2.8 The estimated model in the case of jitter in the spike location in the input-output data, along with ROC curves and theta values of in-sample and out- of-sample prediction (see text). The x axes in the first-order kernels and the x and y axes in the second-order kernels range from 0 to 1000msec, while the ROC curves have in both axes values from 0 to 1. .................................................................... 42 Figure 3.1. Diagrammatic representation of the rat brain, showing the relative location of the hippocampal formation on the left side of the brain and a representation of a transverse slice of the hippocampus with the relative locations of the dentate gyrus (DG), the CA3 regions and the CA1 region. .................................... 47 Figure 3.2. Photograph of a live rat performing a lever-press during the experiment......................................................................................................................... 47 Figure 3.3. Illustration of the behavioral tasks the live rat is performing in each trial of the experiment. (A) Right Sample, (B) Delay Phase, (C) Left Non Match. ......... 48 Figure 3.4. Conceptual representation of the multi-electrode array recording from the CA3 and CA1 hippocampal regions of the behaving rat ............................................ 49 Figure 3.5: Illustrative Case of a complete Multiple-Input Multiple-Output (MIMO) model for the Left Sample Event. (A) ROC curve and theta estimate of each output, (B) first-order kernels, (C) second-order self-kernels and cross- kernels, (D) colormap of the meshes used for the second order kernels, (E) schematic of the topology of the input-output neurons in the CA3 and the CA1 respectively, and the causal connections considered in the MIMO model, (F) the model prediction and the actual recorded activity of each output neurons during the Left Sample task ). The x axes in the first-order kernels and the x and y axes in the second-order kernels range from 0 to 500msec, while the ROC curves have in both axes, values from 0 to 1. ....................................................................................... 55 Figure 3.6. Multiple Input Multiple Output model, Left Non-Match Case. (A) ROC Curve of each output, (B) First order kernels, (C) Second Order self and cross kernels, (D) Conceptual representation of the input neurons from CA3, the output neurons from CA1 and the connections considered in this model, (E) The predicted and the actual spatiotemporal activity ............................................................... 58 Figure 3.7. Multiple Input Multiple Output model, Left Sample Case. (A) ROC Curve of each output, (B) First order kernels, (C) Second Order self and cross kernels, (D) Conceptual representation of the input neurons from CA3, the output x neurons from CA1 and the connections considered in this model, (E) The predicted and the actual spatiotemporal activity ............................................................... 59 Figure 3.8. Multiple Input Multiple Output model, Right Sample Case. (A) ROC Curve of each output, (B) First order kernels, (C) Second Order self and cross kernels, (D) Conceptual representation of the input neurons from CA3, the output neurons from CA1 and the connections considered in this model, (E) The predicted and the actual spatiotemporal activity ............................................................... 60 Figure. 3.9. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for each of the behavioral tasks in a single animal case (1052). ............................................................... 63 Figure. 3.10. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for each of the recording days for all behavioral events in a single animal case (1052) .......................... 64 Figure. 3.11. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1052), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. ....... 65 Figure. 3.12. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1029), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. ....... 66 Figure. 3.13. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1036), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. ....... 67 Figure. 3.14. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1040), along with the xi Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. ....... 68 Figure. 3.15. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1042), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. ....... 69 Figure. 3.16. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1051), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. ....... 70 Figure. 3.17. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1053), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. ....... 71 Figure. 3.18. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1054), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. ....... 72 Figure 3.19: The PDMs of the MIMO model for the Left Sample behavioral event. Each triplet of PDMs corresponds to an output in the model. Included in the plots are the singular values for each PDM. ............................................................... 75 Figure 3.20: The PDMs of the MIMO model for the Left Non Match behavioral event. Each triplet of PDMs corresponds to an output in the model. Included in the plots are the singular values for each PDM. ............................................................... 76 xii Figure 3.21: The PDMs of the MIMO model for the Right Sample behavioral event. Each triplet of PDMs corresponds to an output in the model. Included in the plots are the singular values for each PDM. ............................................................... 77 Figure 3.22: The PDMs of the MIMO model for the Right Non Match behavioral event. Each triplet of PDMs corresponds to an output in the model. Included in the plots are the singular values for each PDM. ............................................................... 78 Figure A.1: Graphs of associated Laguerre functions: (A) zero (L=0), second (L=2), and fourth (L=4) order, for fixed α=0.94, and (B) second order (L=2) for three values of α. Note that the order of the Laguerre function is equal to the number of its zero crossings (A) and that as α approaches 1 the “temporal stretch” of a Laguerre function increases (B). ................................................................................ 96 Figure B.1: A single input / single output N-th order Volterra model (A) and its equivalent model based on principal dynamic modes (B). ............................................... 98 Figure B.2: PDM model approximating a second order Volterra model. ....................... 101 xiii ABSTRACT The increasing availability of multiunit recordings gives new urgency to the need for effective analysis of "multidimensional" time-series data that are derived from the recorded activity of neuronal ensembles in the form of multiple sequences of action potentials--treated mathematically as point-processes and computationally as spike-trains. Whether in conditions of spontaneous activity or under conditions of external stimulation, the objective is the identification and quantification of possible causal links among the neurons generating the observed binary signals. Nonparametric, data driven models with predictive capabilities are excellent candidates for these purposes. When modeling input- output relations in multi-input neuronal systems, it is important to select the subset of inputs that are functionally and causally related to the output. Inputs that do not convey information about the actual transformation not only increase the computational burden but also affect the generalization of the model. Moreover, a reliable functional connectivity measure can provide patterns of information flow that can be linked to physiological and anatomical properties of the system. This Doctoral Thesis presents a multiple-input/multiple-output (MIMO) modeling methodology that can be used to quantify the neuronal dynamics of causal functional relationships in neuronal ensembles using spike-train data recorded from individual neurons. Part of this methodological framework is a novel functional connectivity algorithm. Results from simulated systems and for actual applications using multiunit data recordings from the hippocampus of behaving rats are presented. 1 CHAPTER ONE: INTRODUCTION Computational methods have been used extensively in neuroscience to enhance our understanding of information processing by the nervous system – i.e. the dynamic transformation of neuronal signals as they flow through neuronal ensembles. Quantitative neurophysiological studies have produced computational models that seek to describe the functional relationships between observed neuronal variables of the system of interest. Since brain states are inherently labile and characterized by a high degree of complexity, issues such as neuronal coupling, neuronal coding, neuronal feedback, functional connectivity (convergence/divergence) and functional integration continue to be studied extensively because they are not yet adequately understood. The development of quantitative functional models is a growing trend in computational neuroscience [Stark 1968, Marmarelis & Marmarelis 1978, Koch & Segev 1989, Poznanski 1999, Dayan & Abbot 2001, Marmarelis 2004, Berger & Glanzman 2005, Wu et al. 2006]. Extensive literature discusses various types of models that summarize large amounts of experimental data, using linear or nonlinear methods, depending on the specific characteristics of the system and the goals of the study. Generally, the development of reliable computational models for neuronal ensembles has been hindered by the intrinsic complexity arising from the fact that neuronal function is nonlinear, dynamic, highly interconnected and subject to stochastic variations. Modeling neuronal ensemble interactions in the form of point-process systems adds more considerations to the already complex modeling methodology, as presented in the pioneering early [Berger et al. 1988, 2 Sclabassi et al. 1988] and more recent studies [Marmarelis & Berger 2005]. This task is further complicated by the vast amount of time-series data that is required to represent the activity of a neuronal ensemble and by the many possible interactions among the multiple neuronal units. The recent availability of multi-unit data recordings (through multiple electrodes or multi-electrode arrays) presents a unique opportunity to tackle this important problem if effective methodologies can be developed for their proper analysis. This provides the motivation for the work presented herein. The proper analysis of multi-unit recordings requires efficient methodologies for obtaining multi-input/multi-output (MIMO) models in a nonlinear dynamic context [Song et al. 2007, Zanos et al. 2008a]. This presents us with the daunting complexity of incorporating all existing nonlinear interactions among the various neuronal units. While parametric models (that utilize differential equations to define specific mechanisms or compartments) may provide important insight into the biophysical and physiological aspects of neuronal systems, they require proper validation and are faced with rapidly growing complexity in scaling up to many inputs and outputs. On the other hand, nonparametric models (that utilize Volterra-type functional expansions) provide a general approach that is data-based and does not require a priori model postulates, while it is able to scale up to multiple inputs and outputs rather gracefully. Using broadband experimental data, the nonparametric approach yields data-true models that quantify the nonlinear dynamic manner in which multiple “input” neuronal units interact in order to generate the activity of a causally connected output unit. Although issues of computational complexity still arise for large numbers of inputs and outputs, the 3 nonparametric approach allows the methodical search for the important input interactions in a manner that yields complete, robust and accurate models of manageable complexity. In this proposal, we present an extension of the nonparametric modeling approach to the case of multiple inputs and outputs that requires only moderate increase in representational complexity and computational effort [Zanos et al. 2008a]. The presented methodology estimates a Volterra model with multiple inputs and a single output (MISO) [Zanos et al. 2006] and appends to this output a threshold-trigger operator that generates the spikes at each recorded output. The validation of this MISO model is based on its predictive performance for arbitrary inputs, which is assessed with the use of Receiver Operating Characteristic (ROC) curves commonly used in detection systems. The statistical significance of this predictive capability is assessed with the use of the Mann- Whitney statistic. The resulting MIMO model for multi-unit recordings is composed of multiple MISO model modules, one for each output. We also propose a novel method of selection of a distinct subset of inputs that are causally linked to each output (i.e. the relevant inputs that significantly affect it) based on the prediction of the respective MISO models and it’s Mann-Whitney based statistical evaluation [Zanos et al. 2008b]. Finally an alternative method of quantification of the dynamics of the modeled MIMO systems is presented through the Principal Dynamic Modes, a supplemental modeling approach to the Volterra model that reduces the complexity of the resulting model, increases model tractability and the potential to address model interpretation issues. Computer simulated examples are used to illustrate the efficacy of the proposed methodology and it’s robustness in the presence of noise. We also apply our modeling 4 approach to multi-unit recordings from the CA3 and CA1 areas of the hippocampus of live behaving rats, in connection with the prospective development of hippocampal neuroprostheses that may emulate the nonlinear dynamics of neuronal transmission between the hippocampal regions of dentate gyrus and CA3 or CA1 [Gholmieh et al. 2002], or between the CA3 and the CA1 regions that are believed to be involved in the processing of information leading to long-term memory formation [Gholmieh et al. 2007, Gruart et al. 2006]. The rats perform a delayed non-match to sample task and the recordings are spontaneous activity of the neurons on these specific sites. The value of incorporating more than one input to our models is demonstrated, along with the value of cross- interactions between inputs. Multi-input, multi-output models that capture the nonlinear dynamics for different behavioral events are also demonstrated. The application of the significant input and output selection algorithm to the hippocampal data results in functional connectivity diagrams and overall statistics of these connections provide a better understanding of the dynamic nature of the modeled system. Finally, principal dynamic modes derived from the computed models reveal the model’s primitives, that provide a significantly reduced version of the model and provide an alternative description of inherent basic dynamics of the system. 5 CHAPTER TWO: NEURONAL SYSTEMS MODELING BACKGROUND Neuronal Modeling Approaches The quest to understand how neural signals are processed involves modeling the nervous system in different structural scales, from the molecular and cellular level, to circuits, networks and systems. As with every physiological system modeling procedure, the development of quantitative functional representations (models) combines experimental data with knowledge about the anatomy and the underlying physiology of the system analyzed. The evaluation and comparison of such models can be addressed by how successfully the model mimics the capabilities and the behavior of the real system: the prediction paradigm. The number of different methodologies used to address similar questions is large and one must select a specific approach after careful planning and strong consideration of the constraints present in the experimental data, the paradigm, the system itself, and the application of the resulting model. Nevertheless, the proposed models must always maintain a fine balance among the various functional resolution levels of understanding, providing enough detail to address the lower level modeling needs, and clear and interpretable results at the higher level. System identification approaches have been used extensively in studying the properties of individual neurons and neural circuits [Marmarelis 2004, Dayan and Abbott 2001, Friston 2001, Wu et al. 2006]. The main goal has been to quantify the functional relationship between the various observed variables of the studied system, such as neural 6 recordings, EMG, Eye-movement related saccades, and so on. The resulting models can be categorized as parametric or non-parametric, linear or non-linear, static or dynamic, time-invariant or time-varying depending on the properties of the modeled system. Models are also categorized as discrete-time or continuous-time depending on whether sampled input and output data are used or not. All of them have been used extensively in computational neuroscience studies. Although linear methods continue to be used in modeling neuronal systems, the bulk of the applications currently utilize nonlinear methods, because nonlinearities are inherent and essential in most neuronal systems [Carandini et al. 2005, Ito et al]. Nonlinear modeling has been applied successfully to several neuronal systems to date, such as the retina [Citron et al 1981, Citron et al. 1988, Marmarelis & Naka 1973, Marmarelis & Naka 1974], the auditory cortex [McAlpine 2005, Young & Calhoun 2005], the visual cortex [Carandini et al. 1997, Pack et al. 2006, Rapela et al. 2006, Simoncelli & Heeger 1998], the somatosensory system [Lao et al. 2001], the motor cortex [Paninski et al. 2004a, Truccolo et al. 2005], the cerebellum [Yamazaki & Tanaka 2007], and the hippocampus [Berger et al. 1988, Berger et al. 2001, Sclabassi et al. 1988]. At the level of single neurons, much work has been done using multi-compartmental models based on biophysical principles. Such models attempt to capture known attributes of neuronal structure and function that have been revealed by electrophysiological investigations. These models typically contain a fair number of unknown parameters that must be adjusted in each particular case on the basis of experimental measurements. The use of these models has shed light on several aspects of neuronal function, such as 7 synaptic transmission and somato-dendritic integration [Dayan & Abbot 2005, London & Hausser 2005]. Although such parametric models are extremely useful (when accurate), they present considerable challenges in a practical context with regard to the task of model structure specification and the requisite complexity to approximate the actual system (trade-off between model fidelity and tractability). The problem of model complexity becomes especially acute when we need to scale up the model to accommodate multiple inputs and outputs. Limiting the model parameters to a tractable number can lead to oversimplification of the actual system dynamics, while increasing the number of parameters can raise the model complexity to an impractical level. The issue of model complexity and fidelity is fundamental for the particular subject of this paper: the modeling of neuronal systems with multiple inputs and outputs. The use of multi-compartmental models not only makes the problem very complex computationally, but raises also the risk of misrepresentation, since knowledge about the exact way neurons are interconnected and communicate with each other in a nonlinear context is still limited. Thus, inaccurate representation of this interconnectivity (in structural and functional terms) can lead to invalid modeling results and erroneous conclusions. As an alternative to parametric modeling, the nonparametric approach of Volterra-type modeling employs a general model form and avoids the treacherous task of model specification – especially in the case of nonlinear and highly interconnected systems. For this reason, it is considered a general methodology for nonlinear modeling of physiological systems that is based on a rigorous mathematical framework and provides a quantitative description of the nonlinear dynamics of the input-output relation 8 in the form of a functional series that contains unknown kernel functions estimated from the data. This nonparametric model has predictive capabilities for arbitrary inputs [Marmarelis 2004]. The first extension of Volterra-type nonparametric modeling to the case of two inputs was made in the study of motion detection in the fly visual system and of ganglion cells in the catfish retina [Marmarelis & Naka 1973, Marmarelis & Naka 1974]. This approach was extended to spatio-temporal visual inputs and ganglion cell responses in the frog retina [Citron et al. 1981, Citron et al. 1988] or cortical cells [Emerson et al. 1987]. More recent applications to multiple or spatio-temporal inputs have been reported for several neuronal systems [VanKleef et al. 2005, Pack et al. 2006]. The aforementioned cases concern systems with continuous inputs and continuous outputs. However, the nonparametric modeling methodology has also been applied to neuronal systems with point-process inputs [Marmarelis & Berger 2005], two inputs [Dimoka et al. 2008], multiple inputs and a single output [Zanos et al. 2006] and a MIMO system that is modeled using output-triggered feedback [Song et al. 2007]. The latter study presents a similar approach to our proposal, through nonlinear Volterra modeling of MIMO systems, using a probabilistic approach to estimate the Volterra kernels, the noise term of the model, as well as an output-triggered feedback component. This approach presents the efficiency of Bayesian methods in modeling MIMO neuronal systems and showcases the additive value of components such as the noise term and the feedback when modeling such systems. Note that the input-output data required for the effective estimation of the unknown quantities in nonparametric models (Volterra kernels) must 9 have broad spectral content that endows them with predictive capability for arbitrary inputs and makes their estimation robust in the presence of noise. The nonparametric approach has its own limitations. For instance, direct physiological interpretation of the nonparametric model is difficult, because this model is not developed on the basis of biophysical principles or existing knowledge about the system internal structure, but constitutes instead a data-based representation of the input- output transformation/mapping. Another limitation is that this general approach may lead to cumbersome models in the case of highly nonlinear systems – a problem that has been addressed to some extent by recent work [Marmarelis 2004]. Finally, a potential problem with the proposed approach is the utilization of least-squares methods for the estimation of the model parameters (Laguerre coefficients of the Volterra kernel expansions) in the context of point-process outputs, where the binary nature of the output may introduce biases in the kernel estimates under certain circumstances (elaborated in the Discussion section). Nonetheless, computer simulations (where ground truth is available) indicate that these estimation biases are limited in cases with characteristics similar to our application (as demonstrated by the illustrative example presented in the Results section). The obtained models in applications to real systems (where ground truth is not available) can be evaluated on the basis of their predictive performance and can be accepted as reasonable model approximations if such performance is deemed satisfactory. Alternative modeling methodologies for point-processes have been presented in the past [Brillinger 1975, Brillinger et al. 1976, Cox & Lewis 1966, Snyder 1975, Maas et al. 2002, Ogata & Akaike 1982] with Brillinger’s formulation of the maximum-likelihood estimation 10 problem in the binomial context (for Gaussian variations of the postulated stochastic threshold) deserving special note because it yielded excellent results (in the form of linear input-output models) for real point-process data from three Aplysia neurons [Brillinger 1988]. In recent studies, probabilistic formulations of neuronal encoding have led to methodologies utilizing maximum-likelihood techniques [Paninski et al. 2004b, Song et al. 2007, Truccolo et al. 2005] and network likelihood models [Okatan et al. 2005] that may offer advantages in certain applications. Comparison of modeling approaches leads to the conclusion that a carefully modified nonparametric model that could accommodate more inputs but at the same time keep its size in tractable levels is a reasonable choice for multi-input multi-output neural system modeling. The modification of the already well-known Volterra modeling approach as well as specific techniques that will make the model more compact and robust will be presented in later chapters. Functional Connectivity Studies in Multi-Unit Neural Recordings Models of multi-input/multi-output nonlinear systems can grow to intractable complexity as the number of inputs and the degree of nonlinearity increases. In the most general case, each model output is influenced by all the inputs of the model. Practically though, this is not the case since each model output is significantly influenced by a distinct subset of inputs. Reducing the inputs considered for the model estimation has a profound impact on reducing the model complexity. The question translates to identifying 11 functionally connected neural units and tools of assessing functional connectivity between recordings can be used for that purpose. Functional connectivity can be defined as the set of spatiotemporal correlations between neurophysiological events and the influence that one neural system exerts over another either directly or indirectly [Friston et al. 1993]. Much work has been done to characterize functional connections anatomically [Schubert et al. 2001], through functional neuroimaging [Koch et al. 2002], and multiunit electrodes [Snider et al. 1998, Makarov et al. 2005]. At the level of multiunit electrode recordings, functional connectivity can result from stimulus-locked transients (evoked by sensory, behavioral or any other type of event and occurring around the time of the event), a common afferent input (that is, signal input into the neural system after external stimulation), or stimulus- induced coupling of neuronal firing (phasic coupling of neural assemblies mediated by synaptic connections), in a specific brain area, or among different brain areas. Linear and nonlinear measures have been used to assess functional connectivity between two discrete time series measured simultaneously. The most commonly used measures are the linear cross correlation function and the coherence function, which is the normalized value of the cross-spectrum of the two time series (the Fourier transform of the cross-correlation function) by the self-spectrum of each time series. Another measure that does not escape the linear domain but has the advantage of being directed (determining which one of the two time series is useful in forecasting the other) is Granger causality. Granger causality, in the form of partial directed coherence, uses a series of statistical tests (usually F-tests) on lagged values of the input to examine 12 whether these values provide statistically significant information on future values of the output. It has been applied in neural spiking data [Sameshima and Baccala 1999] and continuous local field potential recordings [Bernasconi and Konig 1999]. The directed transfer function is a similar measure that has been shown to be included in the framework of the Granger causality [Kaminski et al. 2001]. Linear measures seem to be a quick, reliable and computationally efficient in assessing functional connectivity, but they fail to capture any nonlinearities that may be present. Nonlinear measures of dependency include mutual information [Schreiber 2000], which quantifies statistical dependencies between the two time series without any assumption of their respective densities and captures linear and nonlinear dependencies, but is only sensitive to static nonlinear dependencies and requires large amounts of data. Nonlinear interdependences [Quiroga et al. 2002], phase synchronization from Hilbert transform [Mormann et al. 2000] and phase synchronization from the wavelet transform [Lachaux et al. 1999] can also be used to assess dependencies and capture both linear and nonlinear aspects of causality. Their main weaknesses are the relatively significant computational burden they add to the modeling procedure and their need for sufficient (generally large) amounts of data. The use of functional connectivity measures for ascertaining significant inputs in multi-input systems has been suggested in the context of brain-machine interfaces [Sanchez et al. 2004] and EMG recordings [Westwick et al. 2006]. As noted by both studies, large models pose not only computational burden but also affect the models’ 13 generalization, thus making the model reduction through input selection an important part of the modeling approach. Models of CA3 and CA1 of the Hippocampus The hippocampus is a brain site directly related with learning and memory processes and is conserved as a structure across species. It plays an important role in the encoding and retrieval of information used to perform various behavioral tasks. It also believed by many researchers to encode spatial information, information used for special navigation [Touretzky and Redish 1996]. The characteristic anatomy of the hippocampus, especially in the areas of DG, CA3 and CA1 suggest the presence of several subsystems that function as a closed-feedback loop system. Input from the neocortex is passed through the entorhinal cortex to the dentate gyrus (DG) through the perforant path; granule cells of the DG project to the CA3 area which in turn projects to the CA1 and back to the entorhinal cortex through the subiculum. This organization makes the hippocampus an intriguing neuronal system from the system identification standpoint and it’s importance in memory function provides an additional motivation for treating damaged cognitive function in that area [Berger et al. 2005]. The dense collaterals between pyramidal cells in the dorsal hippocampus and especially the CA3 to CA1 projections has prompted investigators to associate the area to autoassociative memory [Treves and Rolls 1992], to be the route that cortically processed information follows to form long term memory [Gruart et al. 2006] but also to encoding of spatial and nonspatial nature, related to specific behavioral tasks [Hampson et al. 14 1999]. Parametric modeling methods used to describe the hippocampal neural dynamics [Sargsyan et al. 2004, Levy 1996, Menschik et al. 1999] lead to complex, large-scale representations that often are not evaluated in terms of their predictability. The data used for the modeling studies of the dorsal hippocampus are predominantly single-unit recordings from CA3 and CA1 pyramidal cells, either in vitro, from acute slices of the hippocampus [Dimoka et al. 2003] or in-vivo, from live behaving rats [Hampson et al. 2002]. The activity recorded can be either a result of stimulation of specific sites of the hippocampus with random pulse train stimuli [Gholmieh et al. 2002], or it can be spontaneous activity, as a result of a specific working memory task. Use of nonparametric modeling in order to quantify the underlying dynamics of the various hippocampal subsystems has been a key step towards a neuroprosthetic platform that would provide a reasonable solution to restore lost functionality in case of malfunction in that area [Berger et al. 2005]. Using the Volterra modeling approach predictive models of the dorsal hippocampus were introduced, in the context of a single- input, single output system [Courellis et al. 2004], two-input, single-output systems [Dimoka et al. 2003] and multiple-input, single output systems [Zanos et al. 2006]. Also equivalent modular forms of the Volterra model termed as principal dynamic modes [Marmarelis 2004] were also used, as a way to reduce the complexity of the resulting models [Courellis et al. 2006]. 15 CHAPTER THREE: NONLINEAR NONPARAMETRIC MULTI INPUT - MULTI OUTPUT MODELS Modeling the transformation of spatiotemporal binary neural patterns by individual spatially distributed neural units is a challenge that can be readily and effectively addressed with nonparametric modeling. Nonparametric models do not need to make any explicit assumptions about the biological mechanisms responsible for the transformation neither do they require any a-priori knowledge about the internals of the system they are going to represent. The Volterra modeling approach has been at the forefront of nonparametric modeling. It is based on a rigorous mathematical framework, it provides quantitative description of the nonlinear dynamics of the system it models in the form of kernels, it has predictive capabilities, it is scalable to arbitrary nonlinear order, it is extensible to multiple inputs, and provides means for quantitatively evaluating how closely the model represents the modeled system. The proposed modeling methodology combines the Volterra modeling approach adapted for point-process input and output datasets [Marmarelis and Berger 2005] with a spike generation mechanism and extends it to multiple inputs and multiple outputs. The presented multi-input / multi-output framework focuses on the computation of the Volterra kernels, the determination of an appropriate decision mechanism for firing a spike, and the evaluation of the prediction accuracy of the model. The estimation of Volterra kernels of a single-input/single-output model using Poisson point process input 16 and point process output datasets has been largely investigated theoretically in [Marmarelis and Berger 2005]. Modeling Framework The multi-input/multi-output nonlinear dynamic models we introduce for modeling neural transformations performed on spike sequences are formed as arrays of multi-input/single-output modules (Figure 2.1A). A multi-input/multi-output model with Q-inputs and P-outputs is formed by an array of P Q-input/single-output modules. ` Multi-input / Single-Output Modules Each MISO module represents the neural transformation of the Q point-process inputs into the respective point-process output. This is done in two stages (see Fig. 1B): (1) a “Nonlinear Volterra Transformation” (NVT) of the Q input spike sequences ( ) q x n , q=1,2,3,..,Q into an intermediate variable ( ) p u n that can be interpreted as the transmembrane potential at the axon hillock of the respective p-th output neuron Figure 2.1: (A) Schematic representation of the decomposition of the Multi- input/multi-output (MIMO) model into an array of multi-input/single-output (MISO) modules. (B) Schematic detail of the structure of the MISO module comprising the cascade of the Nonlinear Volterra Transformation (NVT) and the Threshold Trigger (TT) operators. (A) (B) NVT TT 17 (p=1,2,3,..,P); and (2) a threshold-trigger (TT) operating on the intermediate variable u p (n) to generate the spikes (action potentials) of the respective point-process output y p (n) when the intermediate variable exceeds the specified threshold value T p . The NVT stage of each MISO module contains a set of Volterra kernels (characteristic of each input-output mapping) that operate on the Q inputs to generate the intermediate variable ( ) p u n according to the following mathematical expression (2.1) of the continuous multi-input Volterra model: ( ) ( ) ( ) 1 2 1 2 1 2 1 2 1 0 1 1 0 2 1 2 1 2 ( ) , ( ) ( ) (2.1) p p q p q q Q M p y y x q q m y x x q q q q m m u n k k m x n m k m m x n m x n m − = = = + − + − − ∑∑ ∑∑∑∑ where k 0 , k 1 and k 2 denote the zeroth, first and second order kernels respectively. Higher order kernels may generally exist but are not included in these expressions. At the second stage, the threshold-trigger operator for the respective MISO module ( ), p p p TT u n T     is applied on the intermediate variable ( ) p u n and generates a spike when the threshold value T p is exceeded. ( ) ( ) ( ) ( ) , 0,1,2,..., 1 1 ( ) , , (2.2) 0 p p p p p p p p p y n TT u n T p P if u n T T TT u n T otherwise   = = −     > ∈ℜ =     18 where ( ) p y n is the p th -output of the model at time n, TT p is the p th Threshold Trigger function with θ p being the corresponding threshold value, ( ) p u n is the output of the kernel stage at time n, Q is the number of inputs, and P is the number of outputs. Denoting the Q-input / Single-output model for the p th -output as p M the Q- input / P-output model M can be written as { } , 0 ,1 ,2 ,..., 1 p p P = = − M M with p M defined as { } 1 2 1 2 3 0, 1 , 2, 3 , , , , , p p q p q q p q q q p y p y x y x x y x x x k k k k T = M . Kernel Estimation with the Laguerre Expansion Method Given a set of spike (point-process) input and output datasets, the Volterra kernels in equations (2.1) are estimated by expanding the kernels on the associated Laguerre basis functions (Appendix-A) and by performing least squares optimization to compute the Laguerre basis function coefficients. The computed Laguerre coefficients are used to estimate the Volterra kernels. The savings in representational and computational complexity achieved by the Laguerre expansion of the kernels are considerable. For instance, for a second-order model of a single-input/single output system with kernel memory of M lags, the unknown discrete kernel values that must be estimated are: M for k 1 and M(M+1)/2 for k 2 , after taking into account kernel symmetries – for a total of (M+1)(M+2)/2 unknown discrete kernel values, where the zero-order kernel value is also included [Marmarelis 2004]. A typical value of 50 lags would require the estimation of a total of 1,326 unknown discrete kernel values. Obviously, the computational burden of estimating so many unknown 19 kernel values is immense and would require extremely long input-output data sets to avoid high estimation variance. On the other hand, the use of Laguerre expansions for the estimation of the kernels reduces significantly the number of unknown parameters (expansion coefficients) and makes the computational and representational task tractable. If L Laguerre basis functions are used for this purpose, then the total number of unknown expansion coefficients that must be estimated is: (L+1)(L+2)/2. Experience has shown that a maximum of nine Laguerre basis functions (L=9) is required for adequate representation of the kernels of physiological systems, which corresponds to a maximum total of 55 expansion coefficients. Note that L=5 or even L=3 has been found to be adequate in many actual applications, corresponding to a total of 15 or 6 expansion coefficients respectively. Considerable advantages in terms of estimation accuracy, input- output data requirements and model complexity result from this significant reduction in the number of estimated parameters [Marmarelis 1993, Zanos et al. 2006]. Obviously, the number of free parameters increases rapidly in the MIMO case because of the presence of possible nonlinear interaction terms among the various inputs as they affect the outputs. For instance, a second-order MIMO model with Q inputs and P outputs has a total number of free parameters equal to: PQ(QL 2 +3L+2)/2, which grows mainly as the square of the product QL (dominant dependence) and is only proportional to P. Thus, for a MIMO model with 10 inputs and 10 outputs (P = Q = 10), the total number of free parameters, even for L=3, is 5,050. This provides the motivation to screen the various inputs and their interactions regarding the significance of their contributions 20 to the output in order to reduce the total number of free parameters in the model (see the sub-section Selection of Significant Inputs below). Selection of Threshold and Evaluation of Model Performance Using ROC Curves The predictive capability of the obtained models is the basis for the quantitative evaluation of their performance. When the output data are discretized continuous signals, the most commonly used measure of prediction error is the mean-square error. However, in the case of point-process outputs (as in this study), the mean-square error is not appropriate due to the binary nature of the output, and the balance between counts of correct and incorrect predictions of output spikes should be used instead. A spike predicted by the model is termed a “true positive” (TP), if it coincides with an actual output spike, or it is termed a “false positive” (FP) otherwise. Since the model prediction depends on the threshold value T p of the TT operator (which is not estimated through the Laguerre expansion method used for kernel estimation), various values of the threshold are applied upon the continuous NTV model prediction ( ) p u n and the two performance metrics of True Positives Fraction (TPF) and False Positives Fraction (FPF) are computed for each value of the threshold according to the relations: (number of TP) (2.3) (number of actual output spikes) (number of FP) (number of non-spike events) TPF FPF = = 21 When we plot the FPF value in the abscissa and the TPF value in the ordinate for each threshold value, we obtain the Receiver Operating Characteristic (ROC) curve that was developed in the 1950’s as a tool of assessing the performance of detection systems in the presence of noise. The ROC curve is a graphical representation of the competitive relation between sensitivity (TPF) and specificity (1-FPF) of a binary detector/classifier, as its detection threshold changes. The shape of the ROC curve depicts the detection performance and the Area Under the Curve (AUC) is a quantitative measure of this performance (i.e. the model performance is better over all threshold values when the corresponding AUC value is closer to 1, while the worst performance corresponds to an AUC value of 0.5 when the ROC curve coincides with the diagonal). The threshold value that brings the ROC curve closest to the (0, 1) corner-point is often selected as the “optimum” threshold of the model. 2 2 : ( ) min{(1 ) } (2.4) p p p popt p T T T T D T TP FP = = − + Model Statistics In order to establish the statistical significance of each estimated MISO model (as the quantitative representation of a causal link between the inputs and the output) in the presence of noise or other sources of stochastic variability in the data, we utilize the Mann-Whitney two-sample statistic [Hoeffding 1948] that relates to the AUC value of the ROC curve and can be used to test statistically whether a specific model is better than another as a binary predictor. 22 The Mann-Whitney statistic (MWS) represents the probability θ that a randomly selected sample X i from the intermediate variable, u p (n), that corresponds to zero output will be less than (or equal to) a randomly selected sample Y i from the values of u p (n) that correspond to spikes in the output. Essentially the MWS represents how well these two random variables, X i and Y i , are separated. It has been shown that the area under the ROC curve (calculated using the trapezoidal rule) is equivalent to the MWS. The unbiased estimate of the MWS is the average of the samples Ψ: 1 ( , ) (2.5) 0 i i Y X X Y Y X ≤  Ψ =  >  formed by all possible pair combinations of the two sets of samples X i and Y i in the data record. The MWS is a U-statistic and, according to the theory developed by Hoeffding for U-statistics [Hoeffding 1948], it follows an asymptotically normal distribution with unbiased estimates of its mean and variance calculated from the data. Therefore, we can use a t-test to perform statistical testing of significant differences in performance between two models (e.g. with and without certain model component), using the respective MWS. Specifically, the statistical significance of an estimated model can be tested against the null hypothesis of a random predictor (i.e. no causal input-output relationship). Selection of Significant Inputs The complexity of the MIMO Volterra model depends on the number of inputs that are causally linked to each output, since the number of required kernels rises 23 dramatically with increasing number of inputs. This is especially true in the presence of high-order nonlinearities (higher order Volterra kernels) that give rise to numerous nonlinear interaction terms (cross-kernels) and, consequently, to the number of free parameters that must be estimated (i.e. the Laguerre expansion coefficients). As indicated earlier, for a second-order MIMO Volterra model, the total number of free parameters is: PQ(QL 2 +3L+2)/2, where Q is the number of inputs and P the number of outputs. Therefore, it is important to select only the necessary inputs for each MISO module to achieve the minimum model complexity. The significant inputs are those causally linked to each specific output and they are selected on the basis of their contribution to the prediction of the output using the respective MWS. The algorithm that selects the causally-linked inputs to each output, builds successive models with increasing number of inputs and examines whether the inclusion of a specific input (or set of inputs) improves the predictive accuracy of the model in a statistically significant sense. The comparison of the predictive accuracy of two successive models is performed by applying the t-test (with p < 0.01) on the two MWS values that correspond to the compared models. A complete combinatorial search of all possible pair-wise comparisons can be followed in order to select the causally-linked inputs to each output. In order to keep the number of required pair-wise comparisons manageable in the case of a high number of possible inputs, we propose the following procedure that is more efficient: (1) The MWS of each single-input/single-output model is compared with the MWS of a “random” predictor in order to assess whether the input has statistically 24 significant impact on improving the prediction of the respective output (relative to the “random” predictor that serves as the null hypothesis). The random predictor is a single- input/single-output model with the same input as the tested model but with a statistically independent Poisson output having the same mean firing rate as the output of the tested model. The distribution of theta estimates are calculated for the random predictor under these conditions through Monte Carlo runs and a “cut-off” value is established at the 95% significance level that serves as the threshold for the theta estimate obtained from the actual data. The input is accepted as “statistically significant” at the 95% confidence level when the computed theta estimate is higher than the respective cut-off value. After the significant inputs have been selected for each output, we have an initial set of MISO modules comprising the initial MIMO model. (2) We form all possible pairs of the remaining inputs (that were not selected in step 1) with the selected inputs and examine whether the addition of each pair to each initial MISO model has a statistically significant impact on the predictive ability of the model using the MWS. Theta estimates and their variances are calculated for both cases of the model including a pair of inputs and a model without them and used for a t-test with a null hypothesis of equal mean values at a 99% confidence level. This step is done in order to explore any possible modulatory interactions between the initially selected inputs and those that were not deemed significant on their own. The resulting MISO model contains all inputs that are causally-linked with the specific output either directly (step 1) or through modulatory second-order interactions (step 2). One may continue further to examine higher order interactions, but the process 25 becomes rather burdensome and has only limited value (if one posits that any possible higher order interactions are likely to be among inputs that also exhibit second-order interactions already included in the model). This procedure is repeated for each output of the system. The final MIMO model is comprised of the final set of MISO modules thus constructed. Note that the order of nonlinearity of each final MISO module may be determined after the selection of the relevant inputs through a procedure that examines the statistical significance of improvements in model performance (using the MWS) as we extend the order of nonlinearity. Selection of Model Order and Number of Laguerre Functions In order to select the model order (degree of nonlinearity) of the system and the number of Laguerre functions, we employ the MWS in a similar way to the input selection algorithm. For the model order selection, we start by comparing the first-order with the second-order model using the respective MWS estimates for the two models and a t-test. If the second-order model is selected, then the comparison continues between the second-order and the third-order model, until the predictive performance of the higher order model ceases to improve significantly (based on the t-test of the respective MWS). The same procedure is followed in order to select the proper number of Laguerre basis functions. The two procedures should be performed jointly by first increasing the number of basis functions and subsequently increasing the model order, until no significant improvement is observed in the respective MWS. 26 Principal Dynamic Modes Even when Laguerre functions are employed, the Volterra modeling approach may not always provide the most economical representation of our system, when it comes to the implementation of the model in hardware or software. We can determine the minimum set of linear filters that adequately approximate the kernels by employing Principal Dynamic Modes (PDMs) [Marmarelis 2004]. The resulting model is a parsimonious and modular approximation of the discrete Volterra model. Background on PDM based models and the derivation of the Principal Dynamic Modes is described in Appendix A. The advantages of moving to a PDM representation are obvious not only for tractability reasons but also for the potential to address model interpretation issues. Kernel Interpretation Kernel interpretation varies depending on the type of system they model. Kernels modeling nonlinear neural dynamics can be thought of quantities that capture quantitatively facilitatory and depressive effects on the output of the interaction among spikes at the input within a memory window (kernel memory). Positive kernel values indicate facilitation while negative kernel values indicate depression. The magnitude of the kernel value reflects the strength of the facilitatory or depressive effect. The zero- order kernel may capture the output (dc) offset when the input and output are continuous signals, but in the point-processes cases, the zero-order kernel represents the probability of spontaneous firing, that is spikes in the output that are not preceded by any spike in the input during the system’s active memory. The first-order kernel captures the effect on the 27 current output value of past input spikes within the kernel memory window. The second- order kernel captures the effect on the output of pairs of past input spikes within the kernel memory window. The third-order kernel captures the effect on the output triplets of past input spikes within the kernel memory window, and so on. When more than one input is involved in the model, we distinguish second-order and higher-order kernels to self-kernels and cross-kernels. Self-kernels capture the effect on the output of interactions within each input spike sequence (within path interactions). Cross-kernels capture the effect on the output of interactions among input spike sequences (between path interactions). Much like the different order self-kernels, not all cross-kernels need to be present in a model. Unlike self-kernels that are characterized by a number of symmetries which reduce the number of coefficients to be computed, cross- kernels have generally no symmetries and, thus, such savings are not possible. Introducing Multi-input/Multi-output Models with an Example In an effort to demonstrate the extensibility and the scalability of the proposed modeling approach, we will progressively build the multi-input/multi-output model, starting from the single-input/single-output model. We will employ real data recorded from single neuron units in the CA3 hippocampal region (input) and the CA1 hippocampal region (output) of behaving rats performing a Delayed-Non Match-to- Sample Memory task. The resulting models reflect the function mapping between the CA3 and the CA1 hippocampal region. 28 Single-Input/Single-Output Case In the case that we investigate the functional mapping between one CA3 neuron (input) and one CA1 neuron (output), equation (2.1) reduces to: ( ) ( ) ( ) ( ) 1 2 1 2 3 1 0 1 2 1 2 1 2 0 3 1 2 3 1 2 3 ( ) , ( ) ( ) , , ( ) ( ) ( ) (2.2) M m m m m m m u n k k m x n m k m m x n m x n m k m m m x n m x n m x n m − = = + − + − − + + − − − ∑ ∑∑ ∑∑∑ An example of a single-input / single-output model using experimental data from the behaving rat experiments is shown in Figure 2.2. Figure 2.2(A) shows the location of the input and output neurons for the single-input/single-output case coming from the right hippocampus. Figure 2.2(B) shows the computed first, second, and third order kernel, figure 2.2(C) shows the ROC curve, and figure 2.2(D) shows a segment of the actual response vs. the model prediction. The threshold that gave us the largest true positive fraction and the smallest false positive fraction is at 0.018. 29 Two-Input/Single-Output Case In the case that we investigate the functional mapping between two CA3 neurons (2-inputs) and one CA1 neuron (output), equation (2.1) changes to: (A) (B) Figure 2.2: Single-input/single-output model (right hippocampus): (A) spatial location of the input (CA3 neuron) and the output (CA1 neuron) on the hippocampal cross-section, (B) first, second, and third order kernel, (C) ROC curve, (D) Actual response vs. model prediction (segment). (C) Input Output (D) 30 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , , , , , 1 ( ) , ( ), (2.3) 0 , , x x x x x y n TT u n if u n TT u n otherwi e x s θ θ θ θ =     > ∈ℜ   =     + + + + + + + ∑ ∑ ∑ ∑ 1 1 1 1 2 1 1 2 2 2 1 1 2 M -1 M -1 1x 2 m =0 2x 1 2 2 1 2 1x 1 m =0 2x x 2 m m 1 2 1 1 1 2 m m 3x 1 2 3 1 1 1 0 k m (n - m k m (n - m) k m m u n = k + + x (n - m )x (n - m ) k m m ) k m m m x (n - m )x (n - m (n - m )x ) ( ) ( ) ( ) ( ) , , , , , , , , , , , , , x x x x x + + + + + + + ∑ ∑ ∑ ∑ 2 2 2 1 2 3 1 2 1 2 1 1 2 1 2 3 1 2 2 1 2 3 2x x 1 2 1 1 2 2 m m 3x x x 1 2 3 1 1 2 1 3x x x 1 1 2 2 3 m m m 3x x 1 2 3 1 1 2 2 2 2 1 2 2 2 3 m 3 m m 3 m m m m 2 k m m (n - m )x (n - m ) k m m m (n - m )x (n (n - m k m m (n - m )x (n - m )x (n - m ) - m )x (n - m ) k m m m (n - m )x (n - m )x (n + - )x (n - m m ) ) , , ∑ 1 2 3 m m An example of computed first order kernels, and second and third order self- and cross-kernels is shown in Figure 2.3(A). Note that self-kernels are symmetric while cross kernels are not. As a result, implementational savings associated with the symmetry of self-kernels are not present in cross-kernels. Thus, cross-kernels introduce an additional level of complexity to the model with sizable implementational impact. Figure 2.3(B) shows the model response with the inclusion and the exclusion of cross-kernels and Self Kernels Cross Kernels 31 compares it to the actual response. The inclusion of cross-kernels clearly improves model prediction. (A) (B) Figure 2.3: Two-input / single-output model. (A) Computed first-, second-, and third-order self- and cross-kernels. (B) Model prediction with and without cross- kernels. 32 Three-Input/Single-Output Case We will use the three-input case to illustrate how the inclusion of inputs increases the prediction accuracy of the model. We select three CA3 neurons (input) and one CA1 neuron (output) and we create three models including each time, one, two, and three inputs. For simplicity, interactions among the inputs in the two-input and the three-input model have not been included. Figure 2.4 shows the result of each model prediction and compares them to the recorded data from CA1. Clearly, the inclusion of three inputs improves the prediction accuracy of the model. (A) (B) (C) (D) Input Output (E) Figure 2.4: Comparison of the model prediction between a single-input model (A), a two-input model (B), and a three input model (C); the corresponding recorded response (D) and the spatial location of the input and output neurons on the rat hippocampus (E). 33 Multi-Input/Multi-Output Simulation Example In order to illustrate the efficacy of the proposed modeling methodology and its robustness in the presence of noise (spurious spikes) and jitter (small shift of spike location) in the input-output data, we consider a second-order MISO system with four independent inputs and a single output (since each output is treated separately as a MISO module in the proposed MIMO approach – see Figure 2.1). The functional characteristics of this system are defined by the first-order and second-order Volterra kernels shown in Fig. 2.5 that describe the quantitative manner in which each input affects the output. The form of these kernels resembles the ones that have been experimentally observed. Specifically, the 1 st and 2 nd input have first-order excitatory and second-order inhibitory impact on the output, respectively, and the 4 th input has the reverse effect; while the 3 rd input has no impact on the output since all kernels that link it to the output are zero. A cross-kernel that describes the second-order dynamic interactions between the 1 st and the 4 th input as they affect the output is also present (these are the only existing nonlinear interactions among different inputs in this system). The system is simulated for independent Poisson process inputs with different values of the Poisson parameter λ that defines the mean firing rate of each input. 34 Figure 2.5: The kernels of the second-order simulated MISO system with four inputs and one output: the first column shows the first-order kernels, the second column shows the second-order self-kernels, and the third column shows the cross- kernel between the 1st and the 4th input (see text). Since the main application of our modeling methodology involves spike transformations in the hippocampus, we select the mean firing rates of the four inputs to correspond to the average level of activity of four types of hippocampal neurons observed in our experiments (viz. 3 spikes/sec, 10 spikes/sec, 5 spikes/sec and 1 spikes/sec for the 1 st , 2 nd , 3 rd and 4 th input respectively). The memory of this simulated system is set at 1 sec that corresponds to the observed average memory extent of hippocampal neurons. We select a sampling interval (binwidth) of 10 ms which is slightly smaller than the minimum inter-spike interval observed in the data and larger than the refractory period; thus, each kernel dimension has 100 sampled values (bins). To estimate the kernels and T h r e s h o l d x 1 x 2 x 3 x 4 y k 1 k 2 k cross k 1-4 X X 35 apply the proposed input-selection method, we use 6000 samples of input-output data that correspond to 1 min of recorded neuronal activity. The threshold of the simulated system is chosen to be 0.22 in order to produce an output mean firing rate of about 12 spikes/sec, which is consistent with our experimental observations. Validation of the proposed approach is based on the form of the obtained kernel estimates, the input selection and the predictive capability of the model for these input-output data (“in-sample” training set) and an independent set of input-output data (“out-of-sample” testing set). Note that through 500 Monte Carlo simulations of random predictors using the same firing rates and length of data records, we derive a 95% confidence cutoff value to be used as the threshold of significance for each input. Noise-Free Case First, we test the proposed methodology in the noise-free case (i.e. the input- output data have no spurious spikes). The results of our input selection algorithm are shown in Table 1 for the testing data sets. The calculated MWS theta estimates for each of the single-input/single-output models are shown in the same Table 2.1, along with the 95% confidence cutoff value from the “random predictor”. The resulting decisions of the algorithm are all correct regarding the statistical significance of each specific input in terms of its causal effects on the output. Input Estimated θ 95% cutoff value Input Significant x 1 0.833 0.565 Yes x 2 0.752 0.561 Yes x 3 0.545 0.563 No x 4 0.723 0.573 Yes Table 2.1: Results of MWS-based test for input significance in the simulated noise- free case 36 After the input selection is completed, we determine the proper model order and number L of Laguerre basis functions using a similar MWS-based procedure, as described in Methods. Then, we estimate the Volterra kernels for the selected second- order model using the training data set. The estimated kernels for this simulated example are shown in Fig. 2.6 and exhibit close resemblance to the true kernels. Note that all other cross-kernels (except for the existing k 1,4 ) do not contribute to a statistically significant improvement in output prediction when included in the model. Using the estimated kernels, we predict the output of the system for various threshold values of the TT operator and form the ROC curves that are shown in the right compartment of Fig. 2.6 for the training (in-sample) and testing (out-of-sample) data sets. These ROC curves and the corresponding theta values demonstrate the quality of model performance achieved with the proposed methodology. 37 Figure 2.6: The estimated kernels of this second-order model from the simulated data for the noise free case, along with ROC curves and theta values of in-sample and out-of-sample prediction (see text). The x axes in the first-order kernels and the x and y axes in the second-order kernels range from 0 to 1000msec, while the ROC curves have in both axes values from 0 to 1. Noisy Case We now test the proposed methodology in the presence of spurious spikes (noise) in the input-output data. The spurious spikes are independent Poisson processes that are inserted into the noise-free data (all inputs and output). The resulting “noisy” data are analyzed with the previously described algorithms for input selection and kernel estimation. We examined two levels of noise (numbers of inserted spurious spikes): one equal to one-fourth of the total spikes of each input/output process (for example a spike sequence with 100 true spikes will be contaminated with 25 spurious spikes, leading to a total of 125 spikes) and the other equal to half of the total spikes of each input/output process. For the first case (25% added spurious spikes), the results of the input selection 38 algorithm are shown in Table 2.2 and they are not affected adversely. The same is true for the obtained kernel estimates (not shown in the interest of space). Input Estimated θ 95% cutoff value Input Significant x 1 0.829 0.558 Yes x 2 0.758 0.564 Yes x 3 0.549 0.566 No x 4 0.705 0.568 Yes Table 2.2: Results of MWS-based test for input significance in the simulated 25% added noise case After the input selection, we estimate the Volterra kernels using the noisy input- output data and the results are shown in Fig. 2.7A. It is evident that the kernel estimates are very good, in spite of the large number of spurious spikes, demonstrating the remarkable robustness of this approach. Subsequently, we compute the model prediction of the output using the estimated kernels for both training and testing data sets and evaluate its performance using the ROC curves shown in the right compartment of Fig. 2.7A for the “in-sample” and “out-of-sample” cases. For the second case (50% added spurious spikes), we have the results of input selection shown in Table 2.3. After the input selection, we estimate the Volterra kernels using the noisy input-output data and the results are shown in Fig. 2.7B. It is evident that the kernel estimates remain satisfactory, in spite of the large number of spurious spikes, demonstrating the remarkable robustness of this approach even in the presence of a very large number of spurious spikes. Subsequently, we compute the model prediction of the output using the estimated kernels for both training and testing data sets and evaluate its 39 performance using the ROC curves shown in the right compartment of Fig. 2.7B for the “in-sample” and “out-of-sample” cases. Input Estimated θ 95% cutoff value Input Significant x 1 0.684 0.556 Yes x 2 0.611 0.553 Yes x 3 0.528 0.56 No x 4 0.632 0.554 Yes Table 2.3: Results of MWS-based test for input significance in the simulated 50% added noise case 40 Figure 2.7: (A) The estimated model for the case of 25% spurious spikes in the inputs and output. (B) The estimated model in the case of 50% spurious spikes in the inputs and output. The computed ROC curves are shown in the right compartment of the model for the “in-sample” and “out-of-sample” cases. (B) 41 Spike Jitter Case We now test the proposed methodology for the case of random jitter in the location of the input-output spikes that is simulated by introducing stochastic latencies to each one of the spikes, using a normally distributed latency with a mean value of zero and a standard deviation of 2 lags of jitter. The results of input selection are shown in the Table 2.4 and the estimated kernels with the corresponding ROC curves are shown in Figure 2.8. These results demonstrate the remarkable robustness of this approach in the presence of spike jitter. Input Estimated θ 95% cutoff value Input Significant x 1 0.8 0.574 Yes x 2 0.64 0.575 Yes x 3 0.562 0.575 No x 4 0.69 0.574 Yes Table 2.4: Results of MWS-based test for input significance in the simulated jitter case 42 Figure 2.8 The estimated model in the case of jitter in the spike location in the input- output data, along with ROC curves and theta values of in-sample and out-of- sample prediction (see text). The x axes in the first-order kernels and the x and y axes in the second-order kernels range from 0 to 1000msec, while the ROC curves have in both axes values from 0 to 1. Deleted Spikes Case We also test the proposed methodology for the case where some spikes are mistakenly deleted in the input/ output data. In the test performed here, we delete randomly 30% of the total spikes in all inputs and output. The results of input selection are shown in Table 2.5 and they are not affected adversely. The same is true for the obtained kernel estimates (not shown in the interest of space). 43 Input Estimated θ 95% cutoff value Input Significant x 1 0.744 0.579 Yes x 2 0.64 0.578 Yes x 3 0.534 0.585 No x 4 0.66 0.590 Yes Table 2.5: Results of MWS-based test for input significance in the simulated deleted spikes case Misassigned Spikes Case We finally test the proposed methodology for the case where there are misassigned spikes between the four inputs of the model. In the test performed here, we misassign randomly 5% of the total spikes in each input. The results of input selection are shown in Table 2.6 and they are not affected adversely. The same is true for the obtained kernel estimates (not shown in the interest of space). Input Estimated θ 95% cutoff value Input Significant x 1 0.785 0.564 Yes x 2 0.7 0.567 Yes x 3 0.551 0.569 No x 4 0.713 0.57 Yes Table 2.6: Results of MWS-based test for input significance in the misassigned spikes case Selection of Model order and Number of Laguerre Functions As an illustrative example of the selection procedure for the model order and the number of Laguerre functions that we described earlier, we use the case of 25% spurious spikes case to test the proposed model-order selection algorithm by increasing the model order successively from first to second and then to third for the correct L=3. For each model order, we calculate the MWS theta estimates and run a t-test for the significance of 44 the effect of each model order increase on the output prediction. The results are shown in Table 2.7 and yield the correct answer of a second-order model. Model Order Estimated θ (var{θ}) p-value of t-test Increase of Model Order Significant 1st 0.868 (0.001) - - 2nd 0.961 (0.004) 0 Yes 3rd 0.82 (0.002) 1 No Table 2.7: Results of MWS-based test for model order selection We also test the selection algorithm for the number of Laguerre functions by increasing successively the number of basis functions from two to three and then four (for a second-order model), using the same MWS-based approach as with the model-order selection. The results are shown in Table 2.8 and yield the correct answer of L=3. In practice, these two searches should be performed together by incrementing L for each model order until a “No” outcome occurs and then incrementing the model order until another “No” outcome occurs. Number of Laguerre Functions Estimated θ (var{θ}) p-value of t-test Increase of Laguerre Functions Significant 2 0.92 (0.007) - - 3 0.97 (0.003) 0 Yes 4 0.973 (0.003) 0.015 No Table 2.8: Results of MWS-based test for number of Laguerre functions selection 45 CHAPTER FOUR: APPLICATION OF MIMO MODELING TO THE HIPPOCAMPUS The Cortical Neuroprosthetic Platform A major frontier in the biomedical engineering and generally in the neuroscience research is the development of reliable cortical neural prosthetic device. Advances in the development of neural prostheses have been remarkable over the past years, with applications in the auditory cortex (cochlear implant – [Wilson 2000]), the visual cortex (visual prostheses, [Weiland and Humayun 2006]) and the motor cortex [Schwartz et al. 2006] and brain-machine interfaces [Lebedev and Nicolelis 2006]. The rapid growth in these fields enabled devices that are based on these platforms to move on to either manufacturing stage – cochlear implants succeeded clinically and commercially – or clinical trials on human patients – both visual prostheses and brain machine interfaces are tested on human subjects. Approaches for a prosthetic system that would replace higher thought processes like memory and learning are significantly fewer than the aforementioned applications, since the nature of this kind of prosthesis is totally different that any other prostheses approach. The neuroprosthetic platform we consider here is a device that replaces cortical transmission between specific cortical brain regions ([Berger et al. 2001, Berger and Glanzman 2005], the hippocampus and specifically the area between the DG and CA1 [Berger et al. 2005]. Synaptic circuitry that includes the CA3 and the CA1 areas of the hippocampus is believed to be the route that cortically processed information follows to form long term memory [Gruart et al. 2006]. In case of malfunction, due to aging or 46 damage, this cortical neuroprosthetic device provides a reasonable solution to restore the lost functionality. Implementing such neuroprosthetic devices requires quantitative representation of the transformations performed by the malfunctioning area. Thus the most important component of the neuroprosthetic platform is the biomimetic model that will be used to represent the system. Nonparametric, data driven, scalable, compact models with predictive capabilities are excellent candidates to represent hippocampal transformations as they do not require explicit knowledge of the underlying mechanisms, they can be easily adapted to the required size that will provide sufficient predictive accuracy while they remain compact and computationally tractable. A proof-of-concept of the restoration of lost cognitive function through a biomimetic device has been already presented [Berger et al. 2005], using acute hippocampal slice preparations to acquire data and create a nonparametric model that quantified the CA3-CA1 functional mapping [Gholmieh et al. 2002]. It was a single input / single output model, considering only temporal nonlinear dynamics, based on the Volterra modeling approach adapted for point-process input and output signals. In this study, we use data from live, behaving rats recorded contemporaneously from several spatially distinct sites at CA3 and CA1, while the rats are performing a memory task. Consequently, the functional relationship between the CA3 and the CA1 hippocampal region has a spatial and a temporal dimension with nonlinear characteristics, while the data used for the model are from spontaneous activity related to the behavioral events. 47 Data Acquisition Experimental setup and behavioral training Figure 3.1. Diagrammatic representation of the rat brain, showing the relative location of the hippocampal formation on the left side of the brain and a representation of a transverse slice of the hippocampus with the relative locations of the dentate gyrus (DG), the CA3 regions and the CA1 region. Male Long-Evans rats (n = 11, 3-11 months of age) were trained to criterion on a two lever spatial delayed-nonmatch-to-sample (DNMS) task with randomly occurring variable delay intervals of 1-40 sec (Fig. 3.2). Figure 3.2. Photograph of a live rat performing a lever-press during the experiment. 48 The animal performs the task by pressing (SR) the lever presented (SP) in the sample phase (left or right) (Fig. 3.3A). This constitutes the sample/position information on the trial. The lever is then retracted, and the delay phase initiated, during 5 which the animal must nose poke into a lighted device on the opposite wall for the duration of the delay (Fig. 3.3B). Following the termination of the delay, both levers are extended, and the animal must press the lever opposite the sample/position information (nonmatch response, NR) (Fig. 3.3C). The longer the delay between trials, the less likely the animal will be correct on that trial. Correct responses are rewarded with a drop of water delivered to the trough between both levers, and initiation of a 10 sec intertrial interval (ITI). Incorrect responses produce a darkened chamber for 5 sec without reward, followed by a 5 sec ITI. The next trial is always initiated after at least a 5 sec ITI. Daily sessions consisted of 100- 150 consecutive trials, 5 days per week. Only data from correct trials were utilized in this study. Figure 3.3. Illustration of the behavioral tasks the live rat is performing in each trial of the experiment. (A) Right Sample, (B) Delay Phase, (C) Left Non Match. (A) (B) (C) 49 Multineuron recording techniques Arrays of 16 microwires (40 μm) were surgically implanted using both stereotaxic coordinates and spontaneous cell firing activity to position electrode tips in both the CA3 and CA1 cell layers (Fig. 3.4). All arrays had a fixed geometry with 200 μm between pairs (septotemporal axis) and 800 μm between rows (medial-lateral axis) supplied by NB Labs (Denison TX). Electrode tip length was precisely trimmed to follow the longitudinal curvature of the hippocampus. Data was collected from 11 animals and a total of 198 pyramidal cells (average 16- 20 simultaneously recorded per animal) over at least 7 DNMS sessions with stable extracellular action potential waveforms and consistent event specific firing patterns. Neurons that exhibited low firing rate (<0.2spikes/sec) were excluded as unresponsive and neurons with high firing rates (>6spikes/sec) were excluded, as they were considered to be inter-neurons. Figure 3.4. Conceptual representation of the multi-electrode array recording from the CA3 and CA1 hippocampal regions of the behaving rat 50 MIMO Modeling Results Several instances of each behavioral task were considered across a number of different trials, and all trials were recorded from a specific animal during one session (day) of experiments. Each instance of a behavioral task included the neuronal activity within ±1.5 seconds around the time of pressing the lever that is associated with the task. Recorded cells from the CA1 area were considered as outputs of the MIMO model, while the CA3 cells contributing to the activity of each output CA1 cell were considered as inputs to for the respective MISO module. The multi-unit input-output data recorded during each one of the behavioral tasks (left nonmatch, right nonmatch, left sample, right sample) were analyzed with the proposed methodology and second-order Volterra models were obtained for each MISO module, involving first-order and second-order kernel estimates. The input selection algorithm was applied to each MISO module (corresponding to each of the recorded outputs) and the results are summarized for each of the four behavioral tasks in Table 3.1 below. Note that the MWS estimates for all input channels that were not considered significant with regard to a specific output are omitted in the interest of space. Also note that the theta estimate and its variance for the random predictor are the same for a given output channel (based on the firing rate of that specific output channel) regardless of the input channel. Input Selection for Left Sample Task Output Channel Input Channel Estimated θ 95% cutoff value 39 9 0.698 0.580 39 17 0.827 0.585 39 29 0.804 0.579 41 10 0.643 0.629 41 22 0.666 0.635 41 25 0.639 0.628 51 41 29 0.637 0.633 49 3 0.885 0.773 53 25 0.869 0.701 53 29 0.878 0.714 58 3 0.801 0.721 Input Selection for Left Non-Match Task Output Channel Input Channel Estimated θ 95% cutoff value 38 31 0.883 0.777 39 17 0.827 0.662 39 29 0.804 0.664 41 1 0.640 0.583 41 3 0.560 0.556 41 4 0.557 0.554 41 9 0.766 0.598 42 1 0.994 0.629 42 2 0.720 0.632 42 6 0.624 0.621 Input Selection for Right Sample Task Output Channel Input Channel Estimated θ 95% cutoff value 66 121 0.641 0.639 66 122 0.760 0.656 79 121 0.612 0.609 79 122 0.954 0.616 79 125 0.740 0.613 81 114 0.669 0.621 81 125 0.630 0.619 87 101 0.597 0.568 87 117 0.636 0.561 87 126 0.645 0.569 91 117 0.627 0.563 91 118 0.657 0.563 91 122 0.604 0.563 91 125 0.615 0.563 Input Selection for Right Non-Match Task Output Channel Input Channel Estimated θ 95% cutoff value 76 121 0.957 0.788 76 122 0.956 0.781 81 122 0.901 0.561 81 123 0.588 0.569 81 125 0.604 0.568 81 127 0.562 0.560 83 122 0.900 0.611 83 123 0.633 0.598 85 123 0.822 0.685 85 125 0.695 0.684 85 126 0.765 0.688 52 85 127 0.696 0.688 88 121 0.999 0.833 88 123 0.998 0.849 Table 3.1: Input Selection for the four behavioral tasks The resulting selection of CA3 input cells for each CA1 output cell and the causal connections between them for each behavioral task are shown in Fig. 3.5 - 3.8. The results from the analysis of the data collected during the Left Sample task are presented in detail (for illustrative purposes) in Fig. 3.5 that shows the ROC curves (with the corresponding theta value) for each output of the model (A), the estimated first-order kernels (B), the estimated second-order self-kernels and cross-kernels (C), the colormap used to depict the second-order kernels (assigning the red color to the maximum value and the blue color to the minimum value of each kernel (D), a schematic of the topology of the identified input-output causal connections (E), and the MIMO model prediction and the actual recorded activity for each output neuron (F). Note that the Volterra kernels of the obtained second-order MISO modules for all these cases were estimated using three Laguerre basis functions (L=3) and the nonlinear interactions between different inputs were represented in each MISO model by the estimated cross-kernels (when statistically significant). The model predictions shown in Fig. 3.5(F) for all output CA1 units demonstrate the good predictive capability of the obtained MIMO model. We must note that each MIMO model is animal-specific and task-specific. A closer inspection of Fig. 3.5 can provide important information about the functional characteristics of the modeled system and the underlying physiology that are quantified by the Volterra kernels. The vast majority of the first-order kernels and many 53 of the second-order kernels for all outputs and inputs of this system depict a facilitatory characteristic (positive values), an observation that can be explained by the fact that we considered only pyramidal CA3 and CA1 cells (excluding any interneurons), thus modeling predominantly the monosynaptic connections between CA3 and CA1 cells which are mostly excitatory. The “effective memory” of most input pathways of this system is about 500 ms, since all the values of most first-order and second-order kernels have become negligible before or around a lag of 500 ms. A few kernels retain significant (but not large) values beyond the lag of 500 ms but not beyond a lag of 1000 ms. For illustrative purposes, let us examine closer the kernels of output #1 of the Left Sample case (first line of kernels in Fig. 3.5) that correspond to three significant inputs. The first-order kernels of all three inputs have a similar shape, with positive values in the early lags and relaxing to zero around a lag of 150ms. The second-order self-kernels for two (first and third) of the three inputs have a slight depressive effect when the two lags are shorter than 150ms, while the second-order self-kernel of the second input presents a small facilitatory region when both lags are shorter than 100ms. The second-order cross- kernel between the 1 st and 2 nd input is highly asymmetric, revealing the presence of a modulatory effect that the 2 nd input has on the way the 1 st input influences the output -- but not significantly the other way around. Specifically, this modulatory effect is triphasic: initially depressive (for lags less than 30ms), then facilitatory (for lags between 30 and 120ms) and subsequently depressive again (for lags between 120 and 500ms). The other two cross-kernels are more symmetrical and while the one between the 2 nd and 3 rd input is clearly facilitatory for short lags (less than 100ms), the cross-kernel between the 54 1 st and 2 nd inputs is biphasic (initially depressive up to 80ms lags and excitatory for longer lags). For the second input, the shapes of first-order kernels exhibit a variety of functional characteristics for the four different inputs. The effect of the 3 rd input is excitatory, as well as the 4 th (with maximum excitation around 50ms), while the 1 st and 2 nd inputs exhibit a biphasic behavior (initially excitatory and subsequently inhibitory) with the 1 st input dynamics being slower than the 2 nd . Also noticeable is the fact that the first-order and second-order kernels for the third and fifth output (which have the same significant inputs as shown by the schematic in Fig. 3.5(E)) are very similar to each other. 55 Figure 3.5: Illustrative Case of a complete Multiple-Input Multiple-Output (MIMO) model for the Left Sample Event. (A) ROC curve and theta estimate of each output, (B) first-order kernels, (C) second-order self-kernels and cross-kernels, (D) colormap of the meshes used for the second order kernels, (E) schematic of the topology of the input-output neurons in the CA3 and the CA1 respectively, and the causal connections considered in the MIMO model, (F) the model prediction and the actual recorded activity of each output neurons during the Left Sample task ). The x axes in the first-order kernels and the x and y axes in the second-order kernels range from 0 to 500msec and the ROC curves have in both axes values from 0 to 1. θ=0.95 6 θ=0.98 6 θ=0.98 8 θ=0.98 7 θ=0.98 8 (D ) (E) (A) (B) (C) (F) (min) (max) 56 An illustration of the MIMO model predictions and the actual recorded activity of each output neuron for the other three behavioral tasks are shown in Fig. 3.6 -3.8. It is evident that the observed bursts of activity in some output neurons are usually predicted correctly, as well as the intermittent activity. Table 3.2 reports the significant input channels for each output channel and the respective theta estimates for the four behavioral tasks. These results demonstrate the predictive capability of the obtained models and the efficacy of the proposed approach. Left Non Match Task Output Channel Input Channels θ 38 31 0.799 39 17, 29 0.974 41 1, 3, 4, 9 0.913 42 1, 2, 6 0.971 Right Non Match Task Output Channel Input Channels θ 66 121, 122 0.984 79 121, 122, 125 0.998 81 114, 125 0.807 87 101, 117, 126 0.845 91 117, 118, 122, 125 0.807 Right Sample Task Output Channel Input Channels θ 76 121, 122 0.946 81 122, 123, 125, 127 0.639 83 122, 123 0.944 85 123, 125, 126, 127 0.949 88 121, 123 0.951 Table 3.2 : Mann-Whitney Statistics for final MIMO model for LNM, RNM and RS tasks As a final hypothesis test that the final MIMO models capture an existing temporal dependence, we form the null hypothesis for these MIMO models. We generate 57 N randomly generated independent Poisson spike trains (whose rate parameter follows a normal distribution with a fixed mean value r), which have no causal relationship between each other, and randomly pick one of them as the output and the remaining spike trains (N-1) as the inputs of the system. We build the MISO model for each case and estimate the theta values. We repeat the procedure 50 times and get averaged statistics, that are reported in Table 3.3, for different values of N and r. This test reveals the baseline of the predictor performance of the proposed modeling approach and establishes the causal relationship captured by the aforementioned models. Mean MFRs of Data r 2 spikes/sec 10 spikes/sec 20 spikes/sec Number of Inputs N-1 mean theta (std theta) 2 0.609 (0.019) 0.558 (0.01) 0.543 (0.007) 3 0.611 (0.022) 0.583 (0.011) 0.563 (0.009) 4 0.620 (0.022) 0.609 (0.01) 0.583 (0.007) Table 3.3: Averaged Mann-Whitney Statistics for MISO models estimated from randomly generated independent Poisson spike trains for different MFRs and number of inputs 58 Figure 3.6. Multiple Input Multiple Output model, Left Non-Match Case. (A) ROC Curve of each output, (B) First order kernels, (C) Second Order self and cross kernels, (D) Conceptual representation of the input neurons from CA3, the output neurons from CA1 and the connections considered in this model, (E) The predicted and the actual spatiotemporal activity (A) (B) (C) (D) (E) 59 Figure 3.7. Multiple Input Multiple Output model, Left Sample Case. (A) ROC Curve of each output, (B) First order kernels, (C) Second Order self and cross kernels, (D) Conceptual representation of the input neurons from CA3, the output neurons from CA1 and the connections considered in this model, (E) The predicted and the actual spatiotemporal activity (A) (B) (C) (D) (E) 60 Figure 3.8. Multiple Input Multiple Output model, Right Sample Case. (A) ROC Curve of each output, (B) First order kernels, (C) Second Order self and cross kernels, (D) Conceptual representation of the input neurons from CA3, the output neurons from CA1 and the connections considered in this model, (E) The predicted and the actual spatiotemporal activity (D) (E) (A) (B) (C) 61 Functional Connectivity Results The projection patterns in the hippocampus CA3-CA1 synapse revealed from the input/output selection algorithm using activity from 8 different animals are shown in Figs. 3.9. to 3.18. Note that the weight of the arrows is proportional to the number of sessions and behavioral events in which the causal connection was present, as well as the theta value of the connection at each session/behavioral event. For the complete connectivity figures of single animal cases (figs. 3.11-3.18), along the connectivity patterns, also measured and reported is the contribution of second order terms in the models used to assign functional connectivity as a percentage of connections that emerge only when nonlinearities are taken into account, divided by the overall number of connections found. Also the same measure is reported for the contribution of pairwise comparisons performed at the second step of our input selection algorithm (described in detail in page 24). Furthermore, a convergence ratio (number of input cells that connect to an output, normalized by the total number of active inputs) and a divergence ratio (number of output cells that are connected to an input, normalized by the total number of active outputs) is computed for each input and output respectively and their distribution, along with their median are shown in the convergence/divergence ratio histograms. Finally a histogram of the weights for all connections in a specific case is shown. As shown in Fig. 3.9, the projection patterns remain roughly the same across all behavioral events - especially for those connections that remain consistent throughout days (sessions) of experiments. As shown in Fig. 3.10, the projection patterns also remain roughly the same across days of the experiments, although the degree of consistency is 62 smaller this time compared to the consistency throughout behavioral events. These two facts suggest the presence of animal-specific, rather than event-specific connection patterns. Evident from Figs. 3.11-3.18 for the eight different animals analyzed, is, as expected, the fact that a subset of synaptic connections revealed by our method are consistent across datasets. This is also depicted in the connection weight histograms, where the number of connections drops exponentially in most cases as the weight of the connections increases. Also, all figures reveal both strong convergence (outputs receive projections from multiple input cells) and strong divergence (inputs project to multiple outputs). Nevertheless, a close inspection on both the connectivity schematics as well as the convergence and divergence ratio histograms will reveal a higher degree of divergence, a consistent fact across all animals. Finally, consistent across rats are strong spatial connection trends that direct connections massively from top areas of CA3 to bottom areas of CA1. 63 Figure. 3.9. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for each of the behavioral tasks in a single animal case (1052). 64 Figure. 3.10. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for each of the recording days for all behavioral events in a single animal case (1052) 65 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 50 100 150 200 250 300 350 400 450 Ratio Convergence/Divergence Ratio Histogram Convergence Median =0.71429 Divergence Median =0.89474 Figure. 3.11. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1052), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. Mean 2 nd Order Contribution 12% Mean Pairwise Contribution 7.9% 0 2 4 6 8 10 12 14 0 20 40 60 80 100 120 140 Connection Weight Connection Weights Histogram 66 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 50 100 150 200 250 Ratio Convergence/Divergence Ratio Histogram Convergence Median =0.8 Divergence Median =0.88235 Figure. 3.12. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1029), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. Mean 2 nd Order Contribution 10% Mean Pairwise Contribution 5.5% 0 5 10 15 0 20 40 60 80 100 120 140 160 180 Connection Weight Connection Weights Histogram 67 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 50 100 150 200 250 300 350 400 450 Ratio Convergence/Divergence Ratio Histogram Convergence Median =0.68421 Divergence Median =0.82353 Figure. 3.13. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1036), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. Mean 2 nd Order Contribution 18.6% Mean Pairwise Contribution 4.4% 0 2 4 6 8 10 12 0 50 100 150 200 250 300 Connection Weight Connection Weights Histogram 68 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 50 100 150 200 250 300 350 Ratio Convergence/Divergence Ratio Histogram Convergence Median =0.75 Divergence Median =1 Figure. 3.14. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1040), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. Mean 2 nd Order Contribution 22.4% Mean Pairwise Contribution 8.3% 0 5 10 15 0 5 10 15 20 25 Connection Weight Connection Weights Histogram 69 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 100 200 300 400 500 Ratio Convergence/Divergence Ratio Histogram Convergence Median =0.6 Divergence Median =0.82353 Figure. 3.15. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1042), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. Mean 2 nd Order Contribution 15.6% Mean Pairwise Contribution 2.4% 0 2 4 6 8 10 12 0 20 40 60 80 100 120 140 160 180 Connection Weight Connection Weights Histogram 70 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 100 200 300 400 500 600 Ratio Convergence/Divergence Ratio Histogram Convergence Median =0.71429 Divergence Median =1 Figure. 3.16. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1051), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. Mean 2 nd Order Contribution 14.2% Mean Pairwise Contribution 2.5% 0 2 4 6 8 10 12 14 0 20 40 60 80 100 120 Connection Weight Connection Weights Histogram 71 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 100 200 300 400 500 Ratio Convergence/Divergence Ratio Histogram Convergence Median =0.66667 Divergence Median =0.83333 Figure. 3.17. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1053), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. Mean 2 nd Order Contribution 19.1% Mean Pairwise Contribution 1.7% 0 1 2 3 4 5 6 7 8 9 10 0 50 100 150 200 250 300 350 Connection Weight Connection Weights Histogram 72 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 50 100 150 200 Ratio Convergence/Divergence Ratio Histogram Convergence Median =0.58333 Divergence Median =0.71429 Figure. 3.18. Schematics of the topology of the input and output neurons in the CA3 and CA1 area respectively and the causal connections for all recording days and all behavioral events in a single animal case (1054), along with the Convergence/Divergence Ratio Histogram, the percent contribution of the second order to our model connections (nonlinearity), the contribution of the pairwise connections (step 2 of the algorithm) and the histogram of the connection weights. Mean 2 nd Order Contribution 18.7% Mean Pairwise Contribution 2.5% 0 1 2 3 4 5 6 7 8 9 0 10 20 30 40 50 60 70 80 90 100 Connection Weight Connection Weights Histogram 73 Principal Dynamic Modes Results The principal dynamic modes of all four MIMO models shown in Figs. 3.5 - 3.8 regarding four behavioral events were computed and are shown in Fig. 3.19 – 3.22. The principal dynamic modes for these MIMO systems were computed by constructing a matrix for each output and including all kernels computed for all inputs to that output, first and second order, as well as cross-kernels (if present). By using Singular Value Decomposition (SVD) on this matrix, we come up with the singular vectors, which are in effect the PDMs of this output, as well as the singular values for each PDM. The resulting models, as depicted in all four figures, are significantly more compact and tractable compared to the complete MIMO models shown in figures 3.5 – 3.8. The smaller numbers of final descriptors is self-evident, although the PDM models do not include the necessary non-linearity in these figures in order to be a complete predicting model. Another interesting fact depicted in the PDM models is that the vast majority of the PDM triplets are the same for all outputs, across all behavioral events. As shown in figures 3.5 – 3.8, the PDMs that appear almost for all outputs are a negative exponential PMD that reaches zero around 50 to 100msec (green), a biphasic PDM that starts off negative, peaks at positive values around 50msec and converges to zero (blue) and finally a slower PDM that starts off positive, becomes briefly negative and then positive again to slowly converge to zero values around 300msec. Although these PDMs seem fairly consistent for all connections across behavioral events, the task of linking these PDMs to basic physiological mechanisms is not straightforward. The PDMs do not 74 necessarily represent distinct the physiological mechanisms of the system modeled but they are the model primitives, the modular form of the Volterra model. To link these PDMs to such mechanisms, further experiments and studies are needed, akin to the Berger & Sclabbasi studies during the 1980s [Berger et al. 1988, Sclabbasi et al. 1988] and other similar studies in renal autoregulation dynamics [Chon et al. 1994]. In such studies, experiments that selectively eliminate basic mechanisms of the system are perfomed and the effect of this to the computed PDMs of the new system can elucidate the relationship between the physiological mechanisms and the PDMs. The presented singular values that correspond to each PDM differ in range between outputs and even between PDMs for one specific output. For the case of a specific output, the three singular values represent the effect of each PDM to the intermediate output, before it crosses the non-linearity. For example, in Fig. 3.19, in the first output, the green PDM seems to have the larger effect (0.377) while the red PDM has a lower effect (0.302) and the blue PDM has the lowest effect (0.201). The difference in range between outputs is not related to the significance of each output or a specific PDM to the overall model (for example a PDM with a singular value of 14 is not necessarily more significant than a PDM with a singular value of 0.3) but it depends on the difference of the firing rates of the inputs of a MISO module and the output of that module, a difference that affects the value range of the Volterra kernels and thus of the singular values of the PDMs. 75 0 50 100 150 200 250 300 350 400 450 500 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 (msec) 0 50 100 150 200 250 300 350 400 450 500 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 (msec) 0 50 100 150 200 250 300 350 400 450 500 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 (msec) 0 50 100 150 200 250 300 350 400 450 500 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 (msec) 0 50 100 150 200 250 300 350 400 450 500 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 (msec) Figure 3.19: The PDMs of the MIMO model for the Left Sample behavioral event. Each triplet of PDMs corresponds to an output in the model. Included in the plots are the singular values for each PDM. 0.068 (red) 0.009 (blue) 0.258 (green) 0.481 (red) 0.193 (blue) 0.594 (green) 0.036 (red) 0.012 (blue) 0.189 (green) 3.973 (red) 1.832 (blue) 2.52 (green) 0.302 (red) 0.201 (blue) 0.377 (green) 76 0 50 100 150 200 250 300 350 400 450 500 -1 -0.5 0 0.5 (msec) 0 50 100 150 200 250 300 350 400 450 500 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 (msec) 0 50 100 150 200 250 300 350 400 450 500 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 (msec) 0 50 100 150 200 250 300 350 400 450 500 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 (msec) Figure 3.20: The PDMs of the MIMO model for the Left Non Match behavioral event. Each triplet of PDMs corresponds to an output in the model. Included in the plots are the singular values for each PDM. 0.426 (red) 0.059 (blue) 0.006 (green) 0.117 (red) 0.666 (blue) 0.304 (green) 1.626 (red) 2.969 (blue) 9.755 (green) 4.504 (red) 3.948 (blue) 12.1 (green) 77 0 50 100 150 200 250 300 350 400 450 500 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 (msec) 0 50 100 150 200 250 300 350 400 450 500 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 (msec) 0 50 100 150 200 250 300 350 400 450 500 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 (msec) 0 50 100 150 200 250 300 350 400 450 500 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 (msec) 0 50 100 150 200 250 300 350 400 450 500 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 (msec) Figure 3.21: The PDMs of the MIMO model for the Right Sample behavioral event. Each triplet of PDMs corresponds to an output in the model. Included in the plots are the singular values for each PDM. 1.583 (red) 6.777 (blue) 7.675(green) 0.627 (red) 1.338 (blue) 14.46 (green) 3.299 (red) 1.232 (blue) 1.737 (green) 0.044 (red) 0.268 (blue) 0.095 (green) 0.795 (red) 2.28 (blue) 12.75 (green) 78 0 50 100 150 200 250 300 350 400 450 500 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 (msec) 0 50 100 150 200 250 300 350 400 450 500 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 (msec) 0 50 100 150 200 250 300 350 400 450 500 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 (msec) 0 50 100 150 200 250 300 350 400 450 500 -1 -0.5 0 0.5 (msec) 0 50 100 150 200 250 300 350 400 450 500 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 (msec) Figure 3.22: The PDMs of the MIMO model for the Right Non Match behavioral event. Each triplet of PDMs corresponds to an output in the model. Included in the plots are the singular values for each PDM. 0.963 (red) 1.514 (blue) 3.356 (green) 1.75 (red) 7.241 (blue) 5.438(green) 0.304 (red) 1.535 (blue) 11.5 (green) 0.958 (red) 2.813 (blue) 0.26 (green) 0.453 (red) 0.598 (blue) 1.58 (green) 79 CHAPTER FIVE: CONCLUSIONS The increasing availability of multi-unit recordings gives new urgency to the need for effective analysis of such data that are derived from the recorded activity of neuronal ensembles. A general methodological framework has been presented for analyzing the possible causal links among multi-unit recordings from neuronal ensembles by modeling them as systems with multiple inputs and outputs. This methodological framework is broadly applicable (because of the minimal assumptions regarding the underlying system structure) and robust in the presence of noise or errors in the data. It is computationally efficient and applicable in a practical context (e.g. it can be effective with relatively short data records). It is also scalable and flexible to an arbitrary number of inputs and outputs. The functional connectivity method proposed in this proposal seeks to identify the inputs that are causally-linked to each output and affect significantly the predictive capability of the model (in a statistical sense embodied in the Mann-Whitney statistic). The proposed approach explores all possible second-order interactions among the inputs and exceeds the capabilities of linear approaches, such as cross-correlation, coherence function, Granger causality and partial directed coherence [Sameshima & Baccala 1999] while keeping the problem computationally tractable size and requiring relatively small amounts of experimental data. As shown with the computer-simulated examples (where ground truth is available), as well as in the application to real hippocampal data, the method succeeds in capturing the causal relationships wherever they are present. The method was also shown with computer simulations to be remarkably robust in the presence of noise (spurious spikes in the inputs and outputs), jitter in the recorded spike 80 location and deleted or misassigned spikes that may occur in actual recordings and represent serious impediments in the use of other methods in actual applications. Let us recap the fundamental rationale of the proposed modeling approach. We posit the existence of an internal/intermediate continuous variable u(n) that is generated by the inputs being transformed through a multi-input Volterra model (the operator NVT in Fig. 1) operating in discrete time n = t/Dt, where Dt denotes the time-binwidth which is selected slightly larger than the effective refractory period of the specific neuron (10 ms in this case). We view this intermediate variable u(n) as representing the transmembrane potential at the axon hillock of the output neuron and we posit further that an output spike is generated if, and only if, this variable exceeds a threshold value T. This threshold value and the parameters of the aforementioned multi-input/single-output (MISO) discrete-time Volterra model of the operator NVT must be estimated from the input-output data in this deterministic formulation (using currently the least-squares method). It is evident that this modeling problem is comprised of two critical steps: (1) estimation of the kernels of the NVT model, and (2) the selection of the appropriate threshold for the operator TT. The former is currently accomplished via the least-squares method in an approximate manner (see further comments below) and the latter through the use of ROC curves – which also provide the quantitative means for evaluating the model performance. The use of the Volterra kernels as the basis of the proposed modeling methodology assures excellent predictive capabilities for the timing of output spikes in response to arbitrary inputs and provides reliable quantitative descriptors of the 81 underlying system dynamics – linear and nonlinear – in a manner that advances scientific knowledge and allows the rigorous testing of scientific hypotheses. Application of the proposed methodology to the real hippocampal data has revealed spatio-temporal connectivity patterns between neuronal populations in the CA3 (input) and CA1 (output) region of the rat hippocampus for each of four different behavioral tasks. The false- negatives (misses) that are seen in the output predictions by the model may be due either to the selection of the “optimum” threshold that balances the true-positives with the false- positives, thus resulting to some false-negatives (misses), or to the presence of spurious spikes in the actual output that are not related to the inputs and, therefore, are not predicted by the model. Among these “spurious” spikes there may be spikes due to “spontaneous activity” of the neurons that has been modeled in the past by incorporating in the threshold operator a monotonic function of time independent of the input (see, for instance, [Brillinger 1988]). This input-independent function leads to neuron firing even in the absence of input (spontaneous activity). Our model does not currently incorporate such a feature, but it can easily accommodate it as an additive term to the intermediate variable u(n) – provided that its form can be reasonably postulated. There are three known drawbacks of the proposed methodology. The first is model complexity. The obtained MIMO model generally has considerable complexity (many kernels for multiple, nonlinearly interacting inputs). However, the use of Laguerre expansions of the kernels and the input selection algorithm reduce the model complexity significantly and make the modeling problem tractable, requiring only modest computational effort for the ambitious task of MIMO modeling. Nonetheless, the large 82 number of kernels, as shown in Fig. 3.5, also makes the interpretation of the model difficult and comprehension of the model overwhelming for most investigators. This has provided the motivation for exploring the use of Principal Dynamic Modes to reduce this model complexity [Marmarelis 2004] and the results of this model reduction method are evident in Figs. 3.19-3.22. The second drawback is that our methodology is deterministic in nature and cannot incorporate the stochastic aspects of biological phenomena. One such limitation in the context of the present application is the inability to incorporate the likely stochastic variations in the spike-triggering threshold. The third drawback is the possible bias introduced in the kernel estimates of the NVT component of each MISO module by the employed least-squares method. Note that the optimality of the least- squares estimation method is guaranteed only when the model prediction errors follow a Gaussian distribution (when it becomes equivalent to the optimal maximum-likelihood estimation). In this application, the model prediction errors are not expected to follow a Gaussian distribution because of the binary nature of the output point-process. Therefore, the least-squares method of estimation is not optimal – i.e. it does not generally yield unbiased estimates with minimum variance. The extent of this estimation bias depends on the specific characteristics of each problem and relates to the issue of the truncation of the Volterra model and its implications for the obtained kernel estimates. The proposed methodology is based on the premise that the resulting estimation bias will generally be small and the obtained kernel estimates will be satisfactory approximations of the actual kernels. The validity of this premise can be examined after the model estimation is 83 completed by evaluating the model predictive performance. If the latter is satisfactory, then the validity of this premise can be accepted. This issue has been examined through computer simulations (where ground truth is available) and, as expected, the answer is that the severity of the estimation bias depends on the particular characteristics of each estimation problem (system structure and threshold statistics). Nonetheless, the results of the computer-simulated example presented in the manuscript (which seek to emulate the kernel forms found in this application), demonstrate that the obtained least-square estimates of the kernels of NVT are very close to the actual kernels used in the simulation. This provides some reassurance that the least-squares estimates of the kernels are likely to be reasonable approximations, but it does not guarantee the extent of possible estimation biases. Because of the potential importance of this issue, we will be exploring alternative estimation methods that will guarantee unbiased kernel estimates in the case of point- process outputs – including the maximum-likelihood approach using the binomial distribution of the point-process output -- akin to the Brillinger’s formulation of the problem in the linear case [Brillinger 1988] and iterative estimation methods that avoid estimation biases by employing a cost function consistent with the binary nature of the output. The proposed input selection method seeks to determine which of the pre- designated inputs are statistically significant in terms of their effects of the respective output in the MISO context. The proposed approach does not seek the “minimum set of inputs” in the sense of a maximally reduced model and, therefore, it does not exclude 84 “relay-type” of neurons (i.e. neurons that receive exclusive input from other pre- designated input neurons and simply relay it – perhaps transformed -- to the output neuron). If one is only interested in the overall input-output mapping, then these “relay- type” neurons can be excluded in the interest of model parsimony. However, the purpose of the proposed approach is simply to determine whether there is a significant causal link between each of the pre-designated inputs and the specific output of each MISO model. Nonetheless, the minimum set of inputs can be determined in a subsequent step of analysis whereby one input is excluded on a rotating basis and the effect on the output prediction is assessed statistically (using the MWS) so that this input can be excluded if the effect is insignificant. This problem was addressed in the past through the use of “partial coherence” measures in the linear case [Brillinger 1988]. Finally, it must be emphasized that the identified ‘causal’ relationships between input and output variables (point-processes) in a real neuronal system are only temporal statistical dependencies and cannot be fully substantiated as manifestations of causality without further evidence about the internal workings of the system. For instance, the common effect of a “third party” on both input and output may be viewed mistakenly as the result of an input-output causal relationship. An excellent example of this was presented in the aforementioned study of the three Aplysia neurons, where two of them were affected by the third in a manner that gave the false impression of a causal link between them, although none existed as ascertained by use of the “partial coherence” measures [Brillinger 1988]. It is important to keep this possibility in mind, although it should not discourage us from seeking the “statistical dependencies” in actual multi-unit 85 data, since the latter can provide important insight into the system function when cast in the proper modeling framework. 86 ALPHABETIZED BIBLIOGRAPHY Berger, T. W., J. L. Eriksson, et al. (1988). “Nonlinear-Systems Analysis of the Hippocampal Perforant Path-Dentate Projection .2. Effects of Random Impulse Train Stimulation.” Journal of Neurophysiology 60(3): 1077-1094. Berger, T. W., M. Baudry, et al. (2001). “Brain-implantable biomimetic electronics as the next era in neural prosthetics.” Proceedings of the IEEE 89(7): 993-1012. Berger, T.W. and Glanzman, D.L. (2005). “Toward Replacement Parts for the Brain: Implantable Biomimetic Electronics as the Next Era in Neural Prosthetics.” Cambridge, MA: MIT Press, 2005. Berger, T. W., A. Ahuja, et al. (2005). “Restoring lost cognitive function.” IEEE Eng Med Biol Mag 24(5): 30-44. Bernasconi, C. and P. Konig (1999). “On the directionality of cortical interactions studied by structural analysis of electrophysiological recordings.” Biol Cybern 81(3): 199-210. D. R. Brillinger, “The identification of point process systems,” Annals of Prob., vol. 3, pp. 909-924, 1975. D. R. Brillinger, H. L. Bryant and J. P. Segundo, “Identification of synaptic interactions,” Biological Cybernetics, vol. 22, pp. 213-228, 1976. D. R. Brillinger, “The maximum likelihood approach to the identification of neuronal firing systems,” Ann Biomed Eng, vol. 16, no. 1, pp.3-16, 1988. M. Carandini, D. J. Heeger and J. A. Movshon, “Linearity and normalization in simple cells of the macaque primary visual cortex,” J. Neurosci. vol. 17, pp. 8621–44, 1997. M. Carandini, J. B. Demb, V. Mante, D. J. Tolhurst, Y. Dan, B. A. Olshausen, J. L. Gallant and N. C. Rust, “Do we know what the early visual system does?,” J. Neurosci. vol. 25, pp. 10577–97, 2005 Chon, K.H., N.H. Holstein-Rathlou, D.J. Marsh and V.Z. Marmarelis: Parametric and nonparametric nonlinear modeling of renal autoregulation dynamics, In: Advanced Methods of Physiological System Modeling, Vol. III, edited by V.Z. Marmarelis, pp. 195-210, Plenum Press, New York, 1994. M. C. Citron, J. P. Kroeker and G. D. McCann, “Nonlinear interactions in ganglion cell receptive fields,” J. Neurophysiology, vol. 46, pp. 1161-1176, 1981. 87 M. C. Citron, R. C. Emerson and W. R. Levick, “Nonlinear measurement and classification of receptive fields in cat retinal ganglion cells,” Annals of Biomedical Eng., vol. 16, pp. 65-77, 1988. Courellis, S. H., G. Gholmieh, et al. (2004). “Modeling Nonlinear Neural Dynamics with Volterra-Poisson Kernels.” 2004 IEEE International Joint Conference on Neural Networks 4: 3219-3222. Courellis, S. H., T. P. Zanos, et al. (2006). “Modeling Hippocampal Nonlinear Dynamic Transformations with Principal Dynamic Modes.” Proc. of IEEE-EMBS 2006. D. R. Cox and P. A. W. Lewis, The Statistical Analysis of Series of Event, Ed. London: Nethuen, 1966. Dayan P. and L. F. Abbott (2001). “Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems”, Cambridge, MA: MIT Press Dimoka, A., S. H. Courellis, et al. (2003). “Identification of the Lateral and Medial Perforant Path of the Hippocampus Using Single and Dual Random Impulse Train Stimulation.” Proc. of IEEE-EMBS 2003. A. Dimoka, S. H. Courellis, G. Gholmieh, V. Z. Marmarelis and T. W. Berger, “Modeling the nonlinear properties of the in vitro hippocampal perforant path-dentate system using multi-electrode array technology,” IEEE Trans. on Biomedical Engineering 55(2): 693-702. R. C. Emerson, M. C. Citron, W. J. Vaughn and S. A. Klein, “Nonlinear directionally selective subunits in complex cells of the cat striate cortex,” J. Neurophysiology, vol. 58, pp. 33-65, 1987. Fellous, J. M., M. Rudolph, et al. (2003). “Synaptic background noise controls the input/output characteristics of single cells in an in vitro model of in vivo activity.” Neuroscience 122(3): 811-829. Friston, K.J., C.D. Frith, et al. (1993). “Functional connectivity: the principal component analysis of large (PET) data sets”. J. Cereb. Blood Flow Metab. 13: 5–14. Friston, K. J. (2001). “Brain function, nonlinear coupling, and neuronal transients.” Neuroscientist 7(5): 406-418. Gholmieh, G., S. Courellis, et al. (2002). “An efficient method for studying short-term plasticity with random impulse train stimuli.” Journal of Neuroscience Methods 121(2): 111-127. 88 Goldman, M. S., J. Golowasch, et al. (2001). “Global structure, robustness, and modulation of neuronal models.” Journal of Neuroscience 21(14): 5229-5238. Gruart, A., M. D. Munoz, et al. (2006). “Involvement of the CA3-CA1 synapse in the acquisition of associative learning in behaving mice.” Journal of Neuroscience 26(4): 1077-1087. Hampson, R. E., J. D. Simeral, et al. (1999). “Distribution of spatial and nonspatial information in dorsal hippocampus.” Nature 402(6762): 610-614. Hampson, R. E., J. D. Simeral, et al. (2002). “"Keeping on track": Firing of hippocampal neurons during delayed-nonmatch-to-sample performance.” Journal of Neuroscience 22(2) Hoeffding, W. (1948). “A class of statistics with asymptotically normal distribution.” Annals of Mathematical Statistics 19: 293-325. M. Ito, H. Tamura, I. Fujita and K. Tanaka, “Size and position invariance of neuronal responses in monkey inferotemporal cortex,” J. Neurophysiol., vol. 73, pp. 218–26, 1995. Kaminski, M., M. Ding, et al. (2001). “Evaluating causal relations in neural systems: granger causality, directed transfer function and statistical assessment of significance.” Biol Cybern 85(2): 145-57. Kamitani, Y. and F. Tong (2005). “Decoding the visual and subjective contents of the human brain.” Nature Neuroscience 8(5): 679-685. C. Koch and I. Segev (Eds.), Methods in Neuronal Modeling, Cambridge: MIT Press, 1989 Koch, M. A., D. G. Norris, et al. (2002). “An investigation of functional and anatomical connectivity using magnetic resonance imaging.” Neuroimage 16(1): 241-250. Korenberg, M.J. (1989). “Fast orthogonal algorithms for nonlinear system identification and time series analysis.” Advanced Methods of Physiological System Modeling, Volume II, V.Z. Marmarelis (Ed.), Plenum, New York, pp. 165-178 Lachaux, J. P., E. Rodriguez, et al. (1999). “Measuring phase synchrony in brain signals.” Hum Brain Mapp 8(4): 194-208. Lao, R., O. V. Favorov, et al. (2001). “Nonlinear dynamical properties of a somatosensory cortical model.” Information Sciences 132(1-4): 53-66. 89 Lebedev, M. A. and M. A. Nicolelis (2006). “Brain-machine interfaces: past, present and future.” Trends Neurosci 29(9): 536-46. Levy, W. B. (1996). “A sequence predicting CA3 is a flexible associator that learns and uses context to solve hippocampal-like tasks.” Hippocampus 6(6): 579-590. London, M. and M. Hausser (2005). “Dendritic computation.” Annual Review of Neuroscience 28: 503-532. W. Maas, T. Natschlager and H. Markram, “Real-time computing without stable states: a new framework for neural computation based on perturbations,” Neural Computation, vol. 14, no. 11, pp. 2531-60, 2002. Makarov, V. A., F. Panetsos, et al. (2005). “A method for determining neural connectivity and inferring the underlying network dynamics using extracellular spike recordings.” Journal of Neuroscience Methods 144(2): 265-279. P. Z. Marmarelis and K. I. Naka, “Nonlinear analysis and synthesis of receptive-field responses in the catfish retina. III. Two-input white-noise analysis,” J. Neurophysiology, vol. 36, pp. 634-648, 1973. P. Z. Marmarelis and K. I. Naka, “Identification of multi-input biological systems,” IEEE Trans. Biomedical Engineering, vol. 21, pp. 88-101, 1974. Marmarelis P.Z. and V.Z. Marmarelis. (1978). “Analysis of Physiological Systems: The White-Noise Approach.” New York, NY: Plenum Marmarelis V.Z. (1993). “Identification of Nonlinear Biological Systems Using Laguerre Expansions of Kernels”, Annals of Biomedical Engineering, 21: 573-589 Marmarelis V.Z. (2004). “Nonlinear Dynamic Modeling of Physiological Systems.” New York, NY: Wiley V. Z. Marmarelis and T. W. Berger, “General methodology for nonlinear modeling of neural systems with Poisson point-process inputs,” Mathematical Biosciences, vol. 196, pp. 1-13, 2005. D. McAlpine “Creating a sense of auditory space,” J. Physiol. vol. 566, pp. 21–28, 2005. Menschik, E. D., S. C. Yen, et al. (1999). “Model- and scale-independent performance of a hippocampal CA3 network architecture.” Neurocomputing 26-7: 443-453. 90 Mormann, F., K. Lehnertz, et al. (2000). “Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients.” Physica D- Nonlinear Phenomena 144(3-4): 358-369. Naselaris, T., H. Merchant, et al. (2006). “Large-Scale Organization of Preferred Directions in the Motor Cortex. I. Motor Cortical Hyperacuity for Forward Reaching.” J Neurophysiol 96(6): 3231-3236. Y. Ogata and H. Akaike, “On linear intensity models for mixed doubly stochastic Poisson and self-exciting point processes,” J. Roy. Statist. Soc. Ser. B, vol. 44, pp. 102-107, 1982. M. Okatan, M. A. Wilson and E. N. Brown, “Analyzing functional connectivity using a network likelihood model of ensemble neural spiking activity,” Neural Comput., vol. 17, no. 9, pp.1927-1961, 2005. Pack, C. C., B. R. Conway, et al. (2006). “Spatiotemporal structure of nonlinear subunits in macaque visual cortex.” Journal of Neuroscience 26(3): 893-907. L. Paninski, M. R. Fellows, N. G. Hatsopoulos and J. P. Donoghue, “Spatiotemporal tuning of motor neurons for hand position and velocity,” J. Neurophysiology, vol. 91, pp. 515–532, 2004. L. Paninski, J. W. Pillow and E. P. Simoncelli, “Maximum likelihood estimation of a stochastic integrate-and-fire neural encoding model,” Neural Comput., vol. 16, no. 12, pp. 2533-2561, 2004. Poirazi, P., T. Brannon, et al. (2003). “Pyramidal neuron as two-layer neural network.” Neuron 37(6): 989-999. R. R. Poznanski (Ed.), Modeling in the Neurosciences, Amsterdam: Harwood Academic Publishers, 1999 Prinz, A. A., C. P. Billimoria, et al. (2003). “"Alternative to hand-tuning conductance- based models: Construction and analysis of databases of model neurons.” Journal of Neurophysiology 90(6): 3998-4015. Quiroga, Q.R., A. Kraskov, et al. (2002). “Performance of different synchronization measures in real data: a case study on electroencephalographic signals.” Phys Rev E Stat Nonlin Soft Matter Phys 65(4 Pt 1): 041903. Rapela, J., J. M. Mendel, et al. (2006). “Estimating nonlinear receptive fields from natural images.” J Vis 6(4): 441-74. 91 K. Sameshima and L. A. Baccala, “Using partial directed coherence to describe neuronal ensemble interactions”. J. Neurosci. Meth., 15, 93–103, 1999.T. Schreiber, “Measuring information transfer,” Phys Rev Lett, vol. 85, no. 2, pp. 461-464, 2000. Sanchez, J. C., J. M. Carmena, et al. (2004). “Ascertaining the importance of neurons to develop better brain-machine interfaces.” IEEE Transactions on Biomedical Engineering 51(6): 943-953. Sargsyan, A., A. Melkonyan, et al. (2004). “A computer model of field potential responses for the study of short-term plasticity in hippocampus.” Journal of Neuroscience Methods 135(1-2): 175-191. Schreiber, T. (2000). “Measuring information transfer.” Phys Rev Lett 85(2): 461-4. Schubert, D., J. F. Staiger, et al. (2001). “Layer-specific intracolumnar and transcolumnar functional connectivity of layer V pyramidal cells in rat barrel cortex.” Journal of Neuroscience 21(10): 3580-3592. Schwartz, A. B., X. T. Cui, et al. (2006). “Brain-controlled interfaces: movement restoration with neural prosthetics.” Neuron 52(1): 205-20. Sclabassi, R. J., J. L. Eriksson, et al. (1988). “Nonlinear-Systems Analysis of the Hippocampal Perforant Path-Dentate Projection .1. Theoretical and Interpretational Considerations.” Journal of Neurophysiology 60(3): 1066-1077. E. P. Simoncelli, D. J. Heeger, “A model of neuronal responses in visual area MT,” Vis. Res. vol. 38, pp.743–761, 1998. Snider, R. K., J. F. Kabara, et al. (1998). “Burst firing and modulation of functional connectivity in cat striate cortex.” Journal of Neurophysiology 80(2): 730-744. D. L. Snyder , Random Point Processes, Ed. New York: John Wiley & Sons, 1975. D. Song, R. H. M. Chan, V. Z. Marmarelis, R. E. Hampson, S. A. Deadwyler and T. W. Berger, “Nonlinear dynamic modeling of spike train transformation for hippocampal- cortical prostheses,” IEEE Trans. on Biomedical Engineering,, vol. 54, pp. 1053-1065, 2007. L. Stark, Neurological Control Systems: Studies in Bioengineering, New York: Plenum Press, 1968. Touretzky, D. S. and A. D. Redish (1996). “Theory of rodent navigation based on interacting representations of space.” Hippocampus 6(3): 247-270. 92 Treves, A. and E. T. Rolls (1992). “Computational constraints suggest the need for two distinct input systems to the hippocampal CA3 network.” Hippocampus 2(2): 189-99. W. Truccolo, U. T. Eden, M. R. Fellows, J. P. Donoghue and E. N. Brown, (2005). “A point process framework for relating neural spiking activity to spiking history, neural ensemble, and extrinsic covariate effects,” J. Neurophysiology, vol. 93, no. 2, pp. 1074- 89, 2005. VanKleef, J., A. C. James, et al. (2005). “A spatiotemporal white noise analysis of photoreceptor responses to UV and green light in the dragonfly median ocellus.” Journal of General Physiology 126(5): 481-497. Weiland J.D. and M.S. Humayun (2006) “Intraocular Retinal Prosthesis”, IEEE Engineering in Medicine and Biology Magazine, 25(1): 60-66. Westwick D.T. and R.E. Kearney (1998) “Nonparametric identification of nonlinear biomedical systems, Part I: Theory”, Crit. Rev. biomed. Eng., 26: 153-226 Westwick, D. T., E. A. Pohlmeyer, et al. (2006). “Identification of multiple-input systems with highly coupled inputs: application to EMG prediction from multiple intracortical electrodes.” Neural Comput 18(2): 329-55. Wilson B. (2000). “Cochlear Implant Technology”, Philadelphia, PA: Lippincott Williams & Wilkins Wolfart J., D. Debay, et al. (2005). “Synaptic background activity controls spike transfer from thalamus to cortex” Nature Neuroscience 8(12): 1760-67 Wu, M. C., S. V. David, et al. (2006). “Complete functional characterization of sensory neurons by system identification.” Annu Rev Neurosci 29: 477-505. T. Yamazaki and S. Tanaka, “The cerebellum as a liquid state machine,” Neural Networks vol. 20, pp. 290-297, 2007. Young, E. D. and B. M. Calhoun (2005). “Nonlinear modeling of auditory-nerve rate responses to wideband stimuli.” J Neurophysiol 94(6): 4441-54. Zanos, T. P., S. H. Courellis, et al. (2006). “A multi-input modeling approach to quantify hippocampal nonlinear dynamic transformations.” Proc. of IEEE-EMBS 2006. Zanos, T. P., S. H. Courellis, et al. (2008a). “Nonlinear Modeling of Causal Interrelationships in Neuronal Ensembles”, IEEE Trans. Neural Systems Rehabilitation Engineering, 16(4): 336-352. 93 Zanos, T. P., R. E. Hamson, et al. (2008b). “Functional Connectivity through nonlinear modeling: An application to the rat hippocampus”, Proc. Of IEEE-EMBS 2008. 94 APPENDICES Appendix A: Kernel Computation using Laguerre Basis Functions Explicit estimation of Volterra kernels is computationally intensive whether it is performed indirectly through Wiener kernels (Marmarelis and Marmarelis 1978) or directly through least squares optimization (Korenberg 1989). Unless long datasets are employed in computing the kernels, explicit estimation tends to overfit the data and increase the estimation variance (Marmarelis and Marmarelis 1978). A frequent solution to address such estimation issues is the use of orthogonal basis function sets { } ( ), 0,1, 2,..., 1 l b n l L = − . The kernels to be estimated are expanded on a set of appropriately chosen orthogonal functions, reducing the kernel estimation to computing the coefficients of the basis functions. Consider, for example, the single-input/single-output case with a Volterra model representation as follows: 1 2 1 2 3 1 1 1 0 1 2 1 2 1 2 0 0 0 1 1 1 3 1 2 3 1 2 3 0 0 0 ( ) ( ) ( ) ( , ) ( ) ( ) ( , , ) ( ) ( ) ( ) M M M m m m M M M m m m y n k k m x n m k m m x n m x n m k m m m x n m x n m x n m higher order terms − − − = = = − − − = = = = + − + − − + + − − − + ∑ ∑ ∑ ∑ ∑ ∑ (A.1) where x(n) is the input dataset, y(n) is the output dataset, n is the time index, k 0 , k 1 , k 2 , and k 3 are respectively the zero-, first-, second-, and third-order kernels, and M is the kernel memory. Expanding the kernels on the set { } ( ), 0,1, 2,..., 1 l b n l L = − of orthogonal basis functions, we obtain: 95 1 2 1 2 1 2 1 2 3 1 2 3 1 2 3 1 1 1 (1) (2) 1 2 1 2 1 2 0 0 0 1 1 1 (3) 3 1 2 3 1 2 3 0 0 0 ( ) ( ) ( , ) ( ) ( ) ( , , ) ( ) ( ) ( ) L L L l l l l l l l l l L L L l l l l l l l l l k m c b m k m m c b m b m k m m m c b m b m b m − − − = = = − − − = = = = = = ∑ ∑∑ ∑∑∑ (A.2) where { } 1 1 2 1 2 3 (1) (2) (3) , , l l l l l l c c c is the set of the expansion coefficients. Substituting equations (A.2) into equation (A.1), we obtain: 1 2 1 2 1 2 1 2 1 2 3 1 2 3 2 3 1 2 3 1 1 1 1 1 1 (1) (2) 0 1 2 1 2 0 0 0 0 0 0 1 1 1 1 (3) 1 2 3 1 2 3 0 0 0 0 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( .3) ( ) ( ) ( ) ( ) ( ) ( ) M L M M L L l l l l l l m l m m l l M L L L l l l l l l m m l l l y n k c b m x n m c b m b m x n m x n m A c b m b m b m x n m x n m x n m − − − − − − = = = = = = − − − − = = = = = + − + − − + + − − − ∑∑ ∑ ∑∑∑ ∑∑∑∑ 1 1 1 0 0 M M m higher order terms − − = = + ∑ ∑ Rearranging the sequence of summations in (A.3), we get: 1 2 1 2 1 2 1 2 1 2 3 1 2 1 2 3 1 1 1 1 1 1 1 (1) (2) 0 1 1 2 2 0 0 0 0 0 0 1 1 1 1 (3) 1 1 2 0 0 0 0 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( .4) ( ) ( ) ( L M L L M M l l l l l l l m l l m m L L L M l l l l l l l l m y n k c b m x n m c b m x n m b m x n m A c b m x n m b m − − − − − − = = = = = = − − − − = = = =      = + − + − − +             + −     ∑ ∑ ∑∑ ∑ ∑ ∑∑∑ ∑ 3 2 3 1 1 2 3 3 0 0 ) ( ) ( ) ( ) M M l m m x n m b m x n m higher order terms − − = =     − − +        ∑ ∑ which reduces to: 1 2 1 2 1 2 3 1 2 3 1 2 1 2 3 1 1 1 1 1 1 (1) (2) (3) 0 0 0 0 0 0 0 ( ) ( ) ( ) ( ) ( ) ( ) ( ) L L L L L L l l l l l l l l l l l l l l l l l l y n k c v n c v n v n c v n v n v n higher order terms − − − − − − = = = = = = = + + + + ∑ ∑∑ ∑∑∑ (A.5) with ( ) i v n denoting the convolutions inside the square brackets between each basis function and the input. Equation (A.5) is equivalent to equation (A.1) and the kernel estimation problem has been reduced to the computation of the expansion coefficients. Employing least squares, we compute the expansion coefficients in (A.5), we substitute them in (A.2) and, thus, obtain the Volterra kernels. Similarly, higher order kernels can be expanded over the set of orthogonal functions, when higher order Volterra terms are included in equation (A.1). 96 We chose the associated Laguerre basis of functions to expand the Volterra kernels on (Marmarelis 1993) with the l th –order Laguerre function being mathematically expressed as: 1 1 ( ) 2 2 0 1 ( ) (1 ) ( 1) (1 ) m l k l k k l k n L m k k α α α α − − =    = − − −       ∑ (A.6) where α (alpha) is the decay factor also referred to as the “Laguerre pole” (0<α<1) and m is the m th -lag. Small values of α result on quickly decaying Laguerre functions while values of α approaching 1 increase the “temporal spread” of the Laguerre functions. The l th –order Laguerre function crosses the zero-axis l-times before it asymptotically returns to the zero line. Figure A1 illustrates the effect of different values of α on the “temporal spread” of the Laguerre functions and the dependency of the number of zero-crossings on the order of a Laguerre function. L=2 L=0 L=4 α = 0.94 α α α α = 0.94 α α α α = 0.95 α α α α = 0.90 L = 2 (A) (B) Figure A.1: Graphs of associated Laguerre functions: (A) zero (L=0), second (L=2), and fourth (L=4) order, for fixed α=0.94, and (B) second order (L=2) for three values of α. Note that the order of the Laguerre function is equal to the number of its zero crossings (A) and that as α approaches 1 the “temporal stretch” of a Laguerre function increases (B). 97 Expansion of the kernels on the set of associated Laguerre functions reduces significantly the number of parameters to be estimated. The number of Laguerre coefficients in equation (A.5) is markedly less than the number of values of the kernels in equation (A.1) if explicit kernel computation is to be pursued. However, computation of each Laguerre function at different lags according to equation (A.6) may reduce the performance improvement achieved. Fortunately, there is a recursive relationship among the associated Laguerre functions, mathematically expressed as follows: [ ] 1 1 ( ) ( 1) ( ) ( 1) l l l l L m L n L n L n α − − = − + − − (A.7) which simplifies the computation of Laguerre function values at different lags and for different orders. Numerical differences that may be observed between the values of Laguerre functions computed via equation (A.6) versus values computed via equations (A.7) are generally negligible for all practical purposes. As the set of orthogonal basis functions to expand the kernels on, the associated Laguerre functions provide a number of advantages. They offer a mechanism (based on α) to model kernels with different memory lengths. They can accommodate kernels with different degrees of variability due to the increased variability of a Laguerre function as its order increases. They improve kernel computation even more because of the existence of the recursive relationship of equation (A.7) that reduces considerably the computation overhead associated with each Laguerre function defined by equation (A.6). This recursive relationship can also be used to reduce computational and implementational complexity of the model when it is realized in software or in hardware. 98 Appendix B: Principal Dynamic Modes Modeling with Principal Dynamic Modes (PDMs) can be viewed as the expression of each kernel in a Volterra series in terms of a set of orthogonal functions that are kernel-specific instead of common for all the kernels, as is the case when the associated Laguerre functions to expand the kernels on (see Appendix A). The expression of a Volterra model in terms of its principal dynamic modes has it roots in the Wiener- Bose model and an approximation of the equivalent matrix-vector representation of equation (A.1) by approximating matrices by their most significant eigenvectors. Figure B.1 shows an example of a single input / single output N-th order Volterra model (Fig. B.1A) and its equivalent model based on principal dynamic modes (Fig. B.1B). (A) (B) Figure B.1: A single input / single output N-th order Volterra model (A) and its equivalent model based on principal dynamic modes (B). 99 The PDM model convolves the input x(n) with each PDM μ p (n) and feeds the result of each convolution u p (n) to a static nonlinearity F(u 0 , u 1 , …, u P-1 ), where P is the number of PDMs in the model. The output of this nonlinearity y(n) is the response of the PDM model. Mathematically, the PDM model is expressed as: ( ) ( ) ∑ M -1 p p m=0 0 1 P -1 u n = μ m x(n - m), p = 0,1, y(n) = f(u ,u ,...,u ..., P - 1 ) (B.1) One way to compute the principal dynamic modes is to compute the kernels first, then perform eigendecomposition on the matrices of the equivalent matrix-vector representation of equation (A.1), and consider the eigenvectors corresponding to the most significant eigenvalues as the principal dynamic modes (PDMs). F(u) then is computed from the u(n)=(u 1 (n), u 2 (n), …, u P (n)) T vectors and the corresponding output values r(n). As an example, consider the second order Volterra model shown below in equation (B.2): 1 2 1 1 1 0 1 2 1 2 1 2 0 0 0 ( ) ( ) ( ) ( , ) ( ) ( ) M M M m m m y n k k m x n m k m m x n m x n m − − − = = = = + − + − − ∑ ∑ ∑ (B.2) where, k 0 , k 1 , and k 2 are respectively the zero-, first-, and second-order kernels, and M is the kernel memory. In a matrix-vector form, equation (B.2) can be written as follows: 100 1 2 2 2 1 2 2 2 0 1 2 2 1 (0) ( ) ( ) (0,0) (0,1 ) (0 , 1 ) ( 1 ) ( 1 ) ( 1 ) ( 1 ,0) ( 1 ,1 ) ( 1 , 1 ) ( ) (2) ( 2) ( 2) ( 1 ,0) ( 1 , ( 1 ) ( 1 ) ( 1 ) T T k x n x n k k k M k x n x n k k k M y n k k x n x n k M k M k M x n M x n M       −       − −       −       = + + − −             − −       − − + − +       L L M M O M M M M 2 0 1 2 ( ) ( 1 ) ( 2) 1 ) (0 , 1 ) ( 1 ) ( ) ( ) ( ), T T x n x n x n k M x n M k K X n X n K X n       −        = −         −     − +   = + + M L (B.3) Performing eigendecomposition on K 2 as follows: 2 2 2 0 1 , , ( ,..., ) T T T K K K U U U U I and diag λ λ Μ− = = Λ = Λ = (B.4) and substituting in equation (B.3), we obtain: 0 1 0 1 0 1 0 1 1 ( ) ( ) ( ) ( ) ( ) [ ( )] [ ( )] ( ) ( ) ( ) ( ) [ ( ), ( ),..., ( )] ( ) ( ), : T T T T T T T T T T T M m m m y n k K X n X n U U X n k K X n U X n U X n k K X n V n V n with V n v n v n v n and v n u X n u m th eigenvector − = + + Λ = + + Λ = + + Λ = = − (B.5) Thus, 1 2 0 1 0 ( ) ( ) ( ) M T m m m y n k K X n v n λ − = = + + ∑ (B.6) In equation (B.6) we chose the P most significant (P<M) eigenvalues {λ 0 , λ 1 ,..., λ P-1 } and the corresponding eigenvectors {u p , p=0,1,…P-1} reflected in v p 2 (n). These P most significant eigenvectors comprise the part of the set of the principal dynamic modes. If the first order kernel can not be written as a linear combination of these PDMs then it becomes part of the PDM set. The PDM model corresponding to equation (B.6) is shown in Figure B.2. 101 Figure B.2: PDM model approximating a second order Volterra model. The PDM based models are significantly simpler than their Volterra counterparts. Each mode requires as much storage as a first order kernel and generally a small number of modes is required. The static nonlinearity can be implemented either explicitly through a lookup table, through an analytical approximation, or via an Artificial Neural Network. 102 Appendix C: ROC Curve / Mann – Whitney Statistic The Receiver Operating Characteristic curve (ROC) was developed in the 50’s as a by-product of research on radio signals contaminated by noise. The ROC curve is a graphical plot of the sensitivity versus the specificity of a binary classifier, as its trigger threshold is varied. The terms sensitivity and specificity are often referred to as the fraction of true positives (TP) versus the fraction of false positives (FP). The two quantities are computed as (number of correctly predicted positives) (total number of positives) (number of incorrectly predicted positives) (total number of non-positives) TP FP = = (C.1) The Mann-Whitney statistic estimates the probability θ that a randomly selected observation from a population will be less than, or equal to a randomly selected observation from another population. For our case, the two populations are the output of the kernel subsystem of our model (u p (n)) whenever we have a spike at the true output and whenever we have a zero at the true output. It is computed as the average of a kernel Ψ as 1 1 1 ˆ ( , ) n m i i j i X Y mn θ = = = Ψ ∑∑ (C.2) where 1 ( , ) 0 i i Y X X Y c Y X Y X <   Ψ = =   >  (C.3) 103 and c is a constant. An expression for the variance and the assumption of asymptotic normality can be derived from the fact that the Mann-Whitney statistic is a U-statistic and according to the theory developed by Hoeffding for generalized U-statistics, we can define 2 10 2 01 2 11 ( , ) ( , ) , ( , ) ( , ) , ( , ) ( , ) i j i k i j k j i j i j X Y X Y j k X Y X Y i k X Y X Y ξ θ ξ θ ξ θ   = Ε Ψ Ψ − ≠     = Ε Ψ Ψ − ≠     = Ε Ψ Ψ −   (C.4) and 10 01 11 ( 1) ( 1) ˆ var( ) n m mn ξ ξ ξ θ − + − + = (C.5)

Abstract (if available)

Abstract The increasing availability of multiunit recordings gives new urgency to the need for effective analysis of "multidimensional" time-series data that are derived from the recorded activity of neuronal ensembles in the form of multiple sequences of action potentials -- treated mathematically as point-processes and computationally as spike-trains. Whether in conditions of spontaneous activity or under conditions of external stimulation, the objective is the identification and quantification of possible causal links among the neurons generating the observed binary signals. Nonparametric, data driven models with predictive capabilities are excellent candidates for these purposes. When modeling input-output relations in multi-input neuronal systems, it is important to select the subset of inputs that are functionally and causally related to the output. Inputs that do not convey information about the actual transformation not only increase the computational burden but also affect the generalization of the model. Moreover, a reliable functional connectivity measure can provide patterns of information flow that can be linked to physiological and anatomical properties of the system. This Doctoral Thesis presents a multiple-input/multiple-output (MIMO) modeling methodology that can be used to quantify the neuronal dynamics of causal functional relationships in neuronal ensembles using spike-train data recorded from individual neurons. Part of this methodological framework is a novel functional connectivity algorithm. Results from simulated systems and for actual applications using multiunit data recordings from the hippocampus of behaving rats are presented.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Nonlinear dynamical modeling of single neurons and its application to analysis of long-term potentiation (LTP)

PDF

Functional consequences of network architecture in rat hippocampus: a computational study

PDF

Data-driven modeling of the hippocampus & the design of neurostimulation patterns to abate seizures

PDF

Decoding memory from spatio-temporal patterns of neuronal spikes

PDF

Nonlinear models for hippocampal synapses for multi-scale modeling

PDF

Parametric and non‐parametric modeling of autonomous physiologic systems: applications and multi‐scale modeling of sepsis

PDF

Experimental and computational models for seizure prediction

PDF

Computational investigation of glutamatergic synaptic dynamics: role of ionotropic receptor distribution and astrocytic modulation of neuronal spike timing

PDF

Simulating electrical stimulation and recording in a multi-scale model of the hippocampus

PDF

A million-plus neuron model of the hippocampal dentate gyrus: role of topography, inhibitory interneurons, and excitatory associational circuitry in determining spatio-temporal dynamics of granul...

PDF

Applications of graph theory to brain connectivity analysis

PDF

Investigating the role of muscle physiology and spinal circuitry in sensorimotor control

PDF

Dynamic neuronal encoding in neuromorphic circuits

PDF

Emulating variability in the behavior of artificial central neurons

PDF

A biomimetic approach to non-linear signal processing in ultra low power analog circuits

PDF

Causality and consistency in electrophysiological signals

PDF

Excitatory-inhibitory interactions in pyramidal neurons

PDF

Dealing with unknown unknowns

PDF

Intelligent knowledge acquisition systems: from descriptive to predictive models

PDF

Studies into computational intelligence approaches for the identification of complex nonlinear systems

Asset Metadata

Creator Zanos, Theodoros (author)

Core Title Nonlinear modeling of causal interrelationships in neuronal ensembles: an application to the rat hippocampus

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Biomedical Engineering

Publication Date 05/11/2009

Defense Date 10/31/2008

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag CA1,CA3,computational neuroscience,functional connectivity,hippocampus,Modeling,nonlinear,OAI-PMH Harvest,point process,system identification,Volterra

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Marmarelis, Vasilis Z. (committee chair), Berger, Theodore W. (committee member), Courellis, Spyros (committee member), D'Argenio, David (committee member), Schaal, Stefan (committee member)

Creator Email theozanos@gmail.com,zanos@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m2215

Unique identifier UC188765

Identifier etd-Zanos-2727 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-240402 (legacy record id),usctheses-m2215 (legacy record id)

Legacy Identifier etd-Zanos-2727.pdf

Dmrecord 240402

Document Type Dissertation

Rights Zanos, Theodoros

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu