Measuring functional connectivity from fMRI recordings is important in understanding processing in cortical networks. However, because the brain's connection pattern is complex, currently used methods are prone to producing false functional connections. We introduce differential covariance analysis, a new method that uses derivatives of the signal for estimating functional connectivity. We generated neural activities from dynamical causal modeling and a neural network of Hodgkin-Huxley neurons and then converted them to hemodynamic signals using the forward balloon model. The simulated fMRI signals, together with the ground-truth connectivity pattern, were used to benchmark our method with other commonly used methods. Differential covariance achieved better results in complex network simulations. This new method opens an alternative way to estimate functional connectivity.

Epilepsy is a neurological disorder characterized by the sudden occurrence of unprovoked seizures. There is extensive evidence of significantly altered brain connectivity during seizure periods in the human brain. Research on analyzing human brain functional connectivity during epileptic seizures has been limited predominantly to the use of the correlation method. However, spurious connectivity can be measured between two brain regions without having direct connection or interaction between them. Correlations can be due to the apparent interactions of the two brain regions resulting from common input from a third region, which may or may not be observed. Hence, researchers have recently proposed a sparse-plus-latent-regularized precision matrix (SLRPM) when there are unobserved or latent regions interacting with the observed regions. The SLRPM method yields partial correlations of the conditional statistics of the observed regions given the latent regions, thus identifying observed regions that are conditionally independent of both the observed and latent regions. We evaluate the performance of the methods using a spring-mass artificial network and assuming that some nodes cannot be observed, thus constituting the latent variables in the example. Several cases have been considered, including both sparse and dense connections, short-range and long-range connections, and a varying number of latent variables. The SLRPM method is then applied to estimate brain connectivity during epileptic seizures from human ECoG recordings. Seventy-four clinical seizures from five patients, all having complex partial epilepsy, were analyzed using SLRPM, and brain connectivity was quantified using modularity index, clustering coefficient, and eigenvector centrality. Furthermore, using a measure of latent inputs estimated by the SLRPM method, it was possible to automatically detect 72 of the 74 seizures with four false positives and find six seizures that were not marked manually.

High-density electrocorticogram (ECoG) electrodes are capable of recording neurophysiological data with high temporal resolution with wide spatial coverage. These recordings are a window to understanding how the human brain processes information and subsequently behaves in healthy and pathologic states. Here, we describe and implement delay differential analysis (DDA) for the characterization of ECoG data obtained from human patients with intractable epilepsy. DDA is a time-domain analysis framework based on embedding theory in nonlinear dynamics that reveals the nonlinear invariant properties of an unknown dynamical system. The DDA embedding serves as a low-dimensional nonlinear dynamical basis onto which the data are mapped. This greatly reduces the risk of overfitting and improves the method's ability to fit classes of data. Since the basis is built on the dynamical structure of the data, preprocessing of the data (e.g., filtering) is not necessary. We performed a large-scale search for a DDA model that best fit ECoG recordings using a genetic algorithm to qualitatively discriminate between different cortical states and epileptic events for a set of 13 patients. A single DDA model with only three polynomial terms was identified. Singular value decomposition across the feature space of the model revealed both global and local dynamics that could differentiate electrographic and electroclinical seizures and provided insights into highly localized seizure onsets and diffuse seizure terminations. Other common ECoG features such as interictal periods, artifacts, and exogenous stimuli were also analyzed with DDA. This novel framework for signal processing of seizure information demonstrates an ability to reveal unique characteristics of the underlying dynamics of the seizure and may be useful in better understanding, detecting, and maybe even predicting seizures.

With our ability to record more neurons simultaneously, making sense of these data is a challenge. Functional connectivity is one popular way to study the relationship of multiple neural signals. Correlation-based methods are a set of currently well-used techniques for functional connectivity estimation. However, due to explaining away and unobserved common inputs (Stevenson, Rebesco, Miller, & Körding, 2008 ), they produce spurious connections. The general linear model (GLM), which models spike trains as Poisson processes (Okatan, Wilson, & Brown, 2005 ; Truccolo, Eden, Fellows, Donoghue, & Brown, 2005 ; Pillow et al., 2008 ), avoids these confounds. We develop here a new class of methods by using differential signals based on simulated intracellular voltage recordings. It is equivalent to a regularized AR(2) model. We also expand the method to simulated local field potential recordings and calcium imaging. In all of our simulated data, the differential covariance-based methods achieved performance better than or similar to the GLM method and required fewer data samples. This new class of methods provides alternative ways to analyze neural signals.

In estimating the frequency spectrum of real-world time series data, we must violate the assumption of infinite-length, orthogonal components in the Fourier basis. While it is widely known that care must be taken with discretely sampled data to avoid aliasing of high frequencies, less attention is given to the influence of low frequencies with period below the sampling time window. Here, we derive an analytic expression for the side-lobe attenuation of signal components in the frequency domain representation. This expression allows us to detail the influence of individual frequency components throughout the spectrum. The first consequence is that the presence of low-frequency components introduces a 1/f component across the power spectrum, with a scaling exponent of . This scaling artifact could be composed of diffuse low-frequency components, which can render it difficult to detect a priori. Further, treatment of the signal with standard digital signal processing techniques cannot easily remove this scaling component. While several theoretical models have been introduced to explain the ubiquitous 1/f scaling component in neuroscientific data, we conjecture here that some experimental observations could be the result of such data analysis procedures. Abstract In estimating the frequency spectrum of real-world time series data, we must violate the assumption of infinite-length, orthogonal components in the Fourier basis. While it is widely known that care must be taken with discretely sampled data to avoid aliasing of high frequencies, less attention is given to the influence of low frequencies with period below the sampling time window. Here, we derive an analytic expression for the side-lobe attenuation of signal components in the frequency domain representation. This expression allows us to detail the influence of individual frequency components throughout the spectrum. The first consequence is that the presence of low-frequency components introduces a 1/f component across the power spectrum, with a scaling exponent of . This scaling artifact could be composed of diffuse low-frequency components, which can render it difficult to detect a priori. Further, treatment of the signal with standard digital signal processing techniques cannot easily remove this scaling component. While several theoretical models have been introduced to explain the ubiquitous 1/f scaling component in neuroscientific data, we conjecture here that some experimental observations could be the result of such data analysis procedures. Abstract In estimating the frequency spectrum of real-world time series data, we must violate the assumption of infinite-length, orthogonal components in the Fourier basis. While it is widely known that care must be taken with discretely sampled data to avoid aliasing of high frequencies, less attention is given to the influence of low frequencies with period below the sampling time window. Here, we derive an analytic expression for the side-lobe attenuation of signal components in the frequency domain representation. This expression allows us to detail the influence of individual frequency components throughout the spectrum. The first consequence is that the presence of low-frequency components introduces a 1/f component across the power spectrum, with a scaling exponent of . This scaling artifact could be composed of diffuse low-frequency components, which can render it difficult to detect a priori. Further, treatment of the signal with standard digital signal processing techniques cannot easily remove this scaling component. While several theoretical models have been introduced to explain the ubiquitous 1/f scaling component in neuroscientific data, we conjecture here that some experimental observations could be the result of such data analysis procedures.

The correlation method from brain imaging has been used to estimate functional connectivity in the human brain. However, brain regions might show very high correlation even when the two regions are not directly connected due to the strong interaction of the two regions with common input from a third region. One previously proposed solution to this problem is to use a sparse regularized inverse covariance matrix or precision matrix (SRPM) assuming that the connectivity structure is sparse. This method yields partial correlations to measure strong direct interactions between pairs of regions while simultaneously removing the influence of the rest of the regions, thus identifying regions that are conditionally independent. To test our methods, we first demonstrated conditions under which the SRPM method could indeed find the true physical connection between a pair of nodes for a spring-mass example and an RC circuit example. The recovery of the connectivity structure using the SRPM method can be explained by energy models using the Boltzmann distribution. We then demonstrated the application of the SRPM method for estimating brain connectivity during stage 2 sleep spindles from human electrocorticography (ECoG) recordings using an electrode array. The ECoG recordings that we analyzed were from a 32-year-old male patient with long-standing pharmaco-resistant left temporal lobe complex partial epilepsy. Sleep spindles were automatically detected using delay differential analysis and then analyzed with SRPM and the Louvain method for community detection. We found spatially localized brain networks within and between neighboring cortical areas during spindles, in contrast to the case when sleep spindles were not present.

Our understanding of neural population coding has been limited by a lack of analysis methods to characterize spiking data from large populations. The biggest challenge comes from the fact that the number of possible network activity patterns scales exponentially with the number of neurons recorded ( ). Here we introduce a new statistical method for characterizing neural population activity that requires semi-independent fitting of only as many parameters as the square of the number of neurons, requiring drastically smaller data sets and minimal computation time. The model works by matching the population rate (the number of neurons synchronously active) and the probability that each individual neuron fires given the population rate. We found that this model can accurately fit synthetic data from up to 1000 neurons. We also found that the model could rapidly decode visual stimuli from neural population data from macaque primary visual cortex about 65 ms after stimulus onset. Finally, we used the model to estimate the entropy of neural population activity in developing mouse somatosensory cortex and, surprisingly, found that it first increases, and then decreases during development. This statistical model opens new options for interrogating neural population data and can bolster the use of modern large-scale in vivo Ca and voltage imaging tools. Abstract Our understanding of neural population coding has been limited by a lack of analysis methods to characterize spiking data from large populations. The biggest challenge comes from the fact that the number of possible network activity patterns scales exponentially with the number of neurons recorded ( ). Here we introduce a new statistical method for characterizing neural population activity that requires semi-independent fitting of only as many parameters as the square of the number of neurons, requiring drastically smaller data sets and minimal computation time. The model works by matching the population rate (the number of neurons synchronously active) and the probability that each individual neuron fires given the population rate. We found that this model can accurately fit synthetic data from up to 1000 neurons. We also found that the model could rapidly decode visual stimuli from neural population data from macaque primary visual cortex about 65 ms after stimulus onset. Finally, we used the model to estimate the entropy of neural population activity in developing mouse somatosensory cortex and, surprisingly, found that it first increases, and then decreases during development. This statistical model opens new options for interrogating neural population data and can bolster the use of modern large-scale in vivo Ca and voltage imaging tools.

Inhibition-stabilized networks (ISNs) are neural architectures with strong positive feedback among pyramidal neurons balanced by strong negative feedback from inhibitory interneurons, a circuit element found in the hippocampus and the primary visual cortex. In their working regime, ISNs produce damped oscillations in the -range in response to inputs to the inhibitory population. In order to understand the properties of interconnected ISNs, we investigated periodic forcing of ISNs. We show that ISNs can be excited over a range of frequencies and derive properties of the resonance peaks. In particular, we studied the phase-locked solutions, the torus solutions, and the resonance peaks. Periodically forced ISNs respond with (possibly multistable) phase-locked activity, whereas networks with sustained intrinsic oscillations respond more dynamically to periodic inputs with tori. Hence, the dynamics are surprisingly rich, and phase effects alone do not adequately describe the network response. This strengthens the importance of phase-amplitude coupling as opposed to phase-phase coupling in providing multiple frequencies for multiplexing and routing information.

We propose a time-domain approach to detect frequencies, frequency couplings, and phases using nonlinear correlation functions. For frequency analysis, this approach is a multivariate extension of discrete Fourier transform, and for higher-order spectra, it is a linear and multivariate alternative to multidimensional fast Fourier transform of multidimensional correlations. This method can be applied to short and sparse time series and can be extended to cross-trial and cross-channel spectra (CTS) for electroencephalography data where multiple short data segments from multiple trials of the same experiment are available. There are two versions of CTS. The first one assumes some phase coherency across the trials, while the second one is independent of phase coherency. We demonstrate that the phase-dependent version is more consistent with event-related spectral perturbation analysis and traditional Morlet wavelet analysis. We show that CTS can be applied to short data windows and yields higher temporal resolution than traditional Morlet wavelet analysis. Furthermore, the CTS can be used to reconstruct the event-related potential using all linear components of the CTS.

Nonlinear dynamical system analysis based on embedding theory has been used for modeling and prediction, but it also has applications to signal detection and classification of time series. An embedding creates a multidimensional geometrical object from a single time series. Traditionally either delay or derivative embeddings have been used. The delay embedding is composed of delayed versions of the signal, and the derivative embedding is composed of successive derivatives of the signal. The delay embedding has been extended to nonuniform embeddings to take multiple timescales into account. Both embeddings provide information on the underlying dynamical system without having direct access to all the system variables. Delay differential analysis is based on functional embeddings, a combination of the derivative embedding with nonuniform delay embeddings. Small delay differential equation (DDE) models that best represent relevant dynamic features of time series data are selected from a pool of candidate models for detection or classification. We show that the properties of DDEs support spectral analysis in the time domain where nonlinear correlation functions are used to detect frequencies, frequency and phase couplings, and bispectra. These can be efficiently computed with short time windows and are robust to noise. For frequency analysis, this framework is a multivariate extension of discrete Fourier transform (DFT), and for higher-order spectra, it is a linear and multivariate alternative to multidimensional fast Fourier transform of multidimensional correlations. This method can be applied to short or sparse time series and can be extended to cross-trial and cross-channel spectra if multiple short data segments of the same experiment are available. Together, this time-domain toolbox provides higher temporal resolution, increased frequency and phase coupling information, and it allows an easy and straightforward implementation of higher-order spectra across time compared with frequency-based methods such as the DFT and cross-spectral analysis.

Data sets with high dimensionality such as natural images, speech, and text have been analyzed with methods from condensed matter physics. Here we compare recent approaches taken to relate the scale invariance of natural images to critical phenomena. We also examine the method of studying high-dimensional data through specific heat curves by applying the analysis to noncritical systems: 1D samples taken from natural images and 2D binary pink noise. Through these examples, we concluded that due to small sample sizes, specific heat is not a reliable measure for gauging whether high-dimensional data are critical. We argue that identifying order parameters and universality classes is a more reliable way to identify criticality in high-dimensional data.

The analysis of natural images with independent component analysis (ICA) yields localized bandpass Gabor-type filters similar to receptive fields of simple cells in visual cortex. We applied ICA on a subset of patches called position-centered patches, selected for forming a translation-invariant representation of small patches. The resulting filters were qualitatively different in two respects. One novel feature was the emergence of filters we call double-Gabor filters. In contrast to Gabor functions that are modulated in one direction, double-Gabor filters are sinusoidally modulated in two orthogonal directions. In addition the filters were more extended in space and frequency compared to standard ICA filters and better matched the distribution in experimental recordings from neurons in primary visual cortex. We further found a dual role for double-Gabor filters as edge and texture detectors, which could have engineering applications.

When subjects adapt their reaching movements in the setting of a systematic force or visual perturbation, generalization of adaptation can be assessed psychophysically in two ways: by testing untrained locations in the work space at the end of adaptation (slow postadaptation generalization) or by determining the influence of an error on the next trial during adaptation (fast trial-by-trial generalization). These two measures of generalization have been widely used in psychophysical studies, but the reason that they might differ has not been addressed explicitly. Our goal was to develop a computational framework for determining when a two-state model is justified by the data and to explore the implications of these two types of generalization for neural representations of movements. We first investigated, for single-target learning, how well standard statistical model selection procedures can discriminate two-process models from single-process models when learning and retention coefficients were systematically varied. We then built a two-state model for multitarget learning and showed that if an adaptation process is indeed two-rate, then the postadaptation generalization approach primarily probes the slow process, whereas the trial-by-trial generalization approach is most informative about the fast process. The fast process, due to its strong sensitivity to trial error, contributes predominantly to trial-by-trial generalization, whereas the strong retention of the slow system contributes predominantly to postadaptation generalization. Thus, when adaptation can be shown to be two-rate, the two measures of generalization may probe different brain representations of movement direction.

Neurons in sensory systems convey information about physical stimuli in their spike trains. In vitro, single neurons respond precisely and reliably to the repeated injection of the same fluctuating current, producing regions of elevated firing rate, termed events. Analysis of these spike trains reveals that multiple distinct spike patterns can be identified as trial-to-trial correlations between spike times (Fellous, Tiesinga, Thomas, & Sejnowski, 2004 ). Finding events in data with realistic spiking statistics is challenging because events belonging to different spike patterns may overlap. We propose a method for finding spiking events that uses contextual information to disambiguate which pattern a trial belongs to. The procedure can be applied to spike trains of the same neuron across multiple trials to detect and separate responses obtained during different brain states. The procedure can also be applied to spike trains from multiple simultaneously recorded neurons in order to identify volleys of near-synchronous activity or to distinguish between excitatory and inhibitory neurons. The procedure was tested using artificial data as well as recordings in vitro in response to fluctuating current waveforms.

The pattern of spikes recorded from place cells in the rodent hippocampus is strongly modulated by both the spatial location in the environment and the theta rhythm. The phases of the spikes in the theta cycle advance during movement through the place field. Recently intracellular recordings from hippocampal neurons (Harvey, Collman, Dombeck, & Tank, 2009 ) showed an increase in the amplitude of membrane potential oscillations inside the place field, which was interpreted as evidence that an intracellular mechanism caused phase precession. Here we show that an existing network model of the hippocampus (Tsodyks, Skaggs, Sejnowski, & McNaughton, 1996 ) can equally reproduce this and other aspects of the intracellular recordings, which suggests that new experiments are needed to distinguish the contributions of intracellular and network mechanisms to phase precession.

Convolutive mixtures of signals, which are common in acoustic environments, can be difficult to separate into their component sources. Here we present a uniform probabilistic framework to separate convolutive mixtures of acoustic signals using independent vector analysis (IVA), which is based on a joint distribution for the frequency components originating from the same source and is capable of preventing permutation disorder. Different gaussian mixture models (GMM) served as source priors, in contrast to the original IVA model, where all sources were modeled by identical multivariate Laplacian distributions. This flexible source prior enabled the IVA model to separate different type of signals. Three classes of models were derived and tested: noiseless IVA, online IVA, and noisy IVA. In the IVA model without sensor noise, the unmixing matrices were efficiently estimated by the expectation maximization (EM) algorithm. An online EM algorithm was derived for the online IVA algorithm to track the movement of the sources and separate them under nonstationary conditions. The noisy IVA model included the sensor noise and combined denoising with separation. An EM algorithm was developed that found the model parameters and separated the sources simultaneously. These algorithms were applied to separate mixtures of speech and music. Performance as measured by the signal-to-interference ratio (SIR) was substantial for all three models.

Perceiving and identifying an object is improved by prior exposure to the object. This perceptual priming phenomenon is accompanied by reduced neural activity. But whether suppression of neuronal activity with priming is responsible for the improvement in perception is unclear. To address this problem, we developed a rate-based network model of visual processing. In the model, decreased neural activity following priming was due to stimulus-specific sharpening of representations taking place in the early visual areas. Representation sharpening led to decreased interference of representations in higher visual areas that facilitated selection of one of the competing representations, thereby improving recognition. The model explained a wide range of psychophysical and physiological data observed in priming experiments, including antipriming phenomena, and predicted two functionally distinct stages of visual processing.

Pharmacologically isolated GABAergic irregular spiking and stuttering interneurons in the mouse visual cortex display highly irregular spike times, with high coefficients of variation ≈0.9–3, in response to a depolarizing, constant current input. This is in marked contrast to cortical pyramidal cells, which spike quite regularly in response to the same current injection. We applied time-series analysis methods to show that the irregular behavior of the interneurons was not a consequence of low-dimensional, deterministic processes. These methods were also applied to the Hindmarsh and Rose neuronal model to confirm that the methods are adequate for the types of data under investigation. This result has important consequences for the origin of fluctuations observed in the cortex in vivo.

Gaussian scale mixture models offer a top-down description of signal generation that captures key bottom-up statistical characteristics of filter responses to images. However, the pattern of dependence among the filters for this class of models is prespecified. We propose a novel extension to the gaussian scale mixturemodel that learns the pattern of dependence from observed inputs and thereby induces a hierarchical representation of these inputs. Specifically, we propose that inputs are generated by gaussian variables (modeling local filter structure), multiplied by a mixer variable that is assigned probabilistically to each input from a set of possible mixers. We demonstrate inference of both components of the generative model, for synthesized data and for different classes of natural images, such as a generic ensemble and faces. For natural images, the mixer variable assignments show invariances resembling those of complex cells in visual cortex; the statistics of the gaussian components of the model are in accord with the outputs of divisive normalization models. We also show how our model helps interrelate a wide range of models of image statistics and cortical processing.

Missing data are common in real-world data sets and are a problem for many estimation techniques. We have developed a variational Bayesian method to perform independent component analysis (ICA) on high-dimensional data containing missing entries. Missing data are handled naturally in the Bayesian framework by integrating the generative density model. Modeling the distributions of the independent sources with mixture of gaussians allows sources to be estimated with different kurtosis and skewness. Unlike the maximum likelihood approach, the variational Bayesian method automatically determines the dimensionality of the data and yields an accurate density model for the observed data without overfitting problems. The technique is also extended to the clusters of ICA and supervised classification framework.