PURPOSE: Anesthesiologists simultaneously manage several aspects of patient care during general anesthesia. Automating administration of hypnotic agents could enable more precise control of a patient's level of unconsciousness and enable anesthesiologists to focus on the most critical aspects of patient care. Reinforcement learning (RL) algorithms can be used to fit a mapping from patient state to a medication regimen. These algorithms can learn complex control policies that, when paired with modern techniques for promoting model interpretability, offer a promising approach for developing a clinically viable system for automated anesthestic drug delivery. METHODS: We expand on our prior work applying deep RL to automated anesthetic dosing by now using a continuous-action model based on the actor-critic RL paradigm. The proposed RL agent is composed of a policy network that maps observed anesthetic states to a continuous probability density over propofol-infusion rates and a value network that estimates the favorability of observed states. We train and test three versions of the RL agent using varied reward functions. The agent is trained using simulated pharmacokinetic/pharmacodynamic models with randomized parameters to ensure robustness to patient variability. The model is tested on simulations and retrospectively on nine general anesthesia cases collected in the operating room. We utilize Shapley additive explanations to gain an understanding of the factors with the greatest influence over the agent's decision-making. RESULTS: The deep RL agent significantly outperformed a proportional-integral-derivative controller (median episode median absolute performance error 1.9% ± 1.8 and 3.1% ± 1.1). The model that was rewarded for minimizing total doses performed the best across simulated patient demographics (median episode median performance error 1.1% ± 0.5). When run on real-world clinical datasets, the agent recommended doses that were consistent with those administered by the anesthesiologist. CONCLUSIONS: The proposed approach marks the first fully continuous deep RL algorithm for automating anesthestic drug dosing. The reward function used by the RL training algorithm can be flexibly designed for desirable practices (e.g. use less anesthetic) and bolstered performances. Through careful analysis of the learned policies, techniques for interpreting dosing decisions, and testing on clinical data, we confirm that the agent's anesthetic dosing is consistent with our understanding of best-practices in anesthesia care.

Short-time Fourier transform (STFT) is the most common window-based approach for analyzing the spectrotemporal dynamics of time series. To mitigate the effects of high variance on the spectral estimates due to finite-length, independent STFT windows, state-space multitaper (SSMT) method used a state-space framework to introduce dependency among the spectral estimates. However, the assumed time-invariance of the state-space parameters makes the spectral dynamics difficult to capture when the time series is highly nonstationary. We propose an adaptive SSMT (ASSMT) method as a time-varying extension of SSMT. ASSMT tracks highly nonstationary dynamics by adaptively updating the state parameters and Kalman gains using a heuristic, computationally efficient exponential smoothing technique. In analyses of simulated data and real human electroencephalogram (EEG) recordings, ASSMT showed improved denoising and smoothing properties relative to standard multitaper and SSMT approaches.

To capture the slowly time-varying spectral content of real-world time-series, a common paradigm is to partition the data into approximately stationary intervals and perform inference in the time-frequency domain. However, this approach lacks a corresponding nonstationary time-domain generative model for the entire data and thus, time-domain inference occurs in each interval separately. This results in distortion/discontinuity around interval boundaries and can consequently lead to erroneous inferences based on any quantities derived from the posterior, such as the phase. To address these shortcomings, we propose the Piecewise Locally Stationary Oscillation (PLSO) model for decomposing time-series data with slowly time-varying spectra into several oscillatory, piecewise-stationary processes. PLSO, as a nonstationary time-domain generative model, enables inference on the entire time-series without boundary effects and simultaneously provides a characterization of its time-varying spectral properties. We also propose a novel two-stage inference algorithm that combines Kalman theory and an accelerated proximal gradient algorithm. We demonstrate these points through experiments on simulated data and real neural data from the rat and the human brain.

Assessing directional influences between neurons is instrumental to understand how brain circuits process information. To this end, Granger causality, a technique originally developed for time-continuous signals, has been extended to discrete spike trains. A fundamental assumption of this technique is that the temporal evolution of neuronal responses must be due only to endogenous interactions between recorded units, including self-interactions. This assumption is however rarely met in neurophysiological studies, where the response of each neuron is modulated by other exogenous causes such as, for example, other unobserved units or slow adaptation processes. Here, we propose a novel point-process Granger causality technique that is robust with respect to the two most common exogenous modulations observed in real neuronal responses: within-trial temporal variations in spiking rate and between-trial variability in their magnitudes. This novel method works by explicitly including both types of modulations into the generalized linear model of the neuronal conditional intensity function (CIF). We then assess the causal influence of neuron i onto neuron j by measuring the relative reduction of neuron j's point process likelihood obtained considering or removing neuron i. CIF's hyper-parameters are set on a per-neuron basis by minimizing Akaike's information criterion. In synthetic data sets, generated by means of random processes or networks of integrate-and-fire units, the proposed method recovered with high accuracy, sensitivity and robustness the underlying ground-truth connectivity pattern. Application of presently available point-process Granger causality techniques produced instead a significant number of false positive connections. In real spiking responses recorded from neurons in the monkey pre-motor cortex (area F5), our method revealed many causal relationships between neurons as well as the temporal structure of their interactions. Given its robustness our method can be effectively applied to real neuronal data. Furthermore, its explicit estimate of the effects of unobserved causes on the recorded neuronal firing patterns can help decomposing their temporal variations into endogenous and exogenous components.

Ketamine is an NMDA receptor antagonist commonly used to maintain general anesthesia. At anesthetic doses, ketamine causes high power gamma (25-50 Hz) oscillations alternating with slow-delta (0.1-4 Hz) oscillations. These dynamics are readily observed in local field potentials (LFPs) of non-human primates (NHPs) and electroencephalogram (EEG) recordings from human subjects. However, a detailed statistical analysis of these dynamics has not been reported. We characterize ketamine's neural dynamics using a hidden Markov model (HMM). The HMM observations are sequences of spectral power in seven canonical frequency bands between 0 to 50 Hz, where power is averaged within each band and scaled between 0 and 1. We model the observations as realizations of multivariate beta probability distributions that depend on a discrete-valued latent state process whose state transitions obey Markov dynamics. Using an expectation-maximization algorithm, we fit this beta-HMM to LFP recordings from 2 NHPs, and separately, to EEG recordings from 9 human subjects who received anesthetic doses of ketamine. Our beta-HMM framework provides a useful tool for experimental data analysis. Together, the estimated beta-HMM parameters and optimal state trajectory revealed an alternating pattern of states characterized primarily by gamma and slow-delta activities. The mean duration of the gamma activity was 2.2s([1.7,2.8]s) and 1.2s([0.9,1.5]s) for the two NHPs, and 2.5s([1.7,3.6]s) for the human subjects. The mean duration of the slow-delta activity was 1.6s([1.2,2.0]s) and 1.0s([0.8,1.2]s) for the two NHPs, and 1.8s([1.3,2.4]s) for the human subjects. Our characterizations of the alternating gamma slow-delta activities revealed five sub-states that show regular sequential transitions. These quantitative insights can inform the development of rhythm-generating neuronal circuit models that give mechanistic insights into this phenomenon and how ketamine produces altered states of arousal.

Electrodermal activity (EDA) is a direct read-out of sweat-induced changes in the skin's electrical conductance. Sympathetically-mediated pulsatile changes in skin sweat measured as EDA resemble an integrate-and-fire process, which yields an inverse Gaussian model as the inter-pulse interval distribution. We have previously showed that the inter-pulse intervals in EDA follow an inverse Gaussian distribution. However, the statistical structure of EDA pulse amplitudes has not yet been characterized based on the physiology. Expanding upon the integrate-and-fire nature of sweat glands, we hypothesized that the amplitude of an EDA pulse is proportional to the excess volume of sweat produced compared to what is required to just reach the surface of the skin. We modeled this as the difference of two inverse Gaussian models for each pulse, one which represents the time required to produce just enough sweat to rise to the surface of the skin and one which represents the time requires to produce the actual volume of sweat. We proposed and tested a series of four simplifications of our hypothesis, ranging from a single difference of inverse Gaussians to a single simple inverse Gaussian. We also tested four additional models for comparison, including the lognormal and gamma distributions. All models were tested on EDA data from two subject cohorts, 11 healthy volunteers during 1 hour of quiet wakefulness and a different set of 11 healthy volunteers during approximately 3 hours of controlled propofol sedation. All four models which represent simplifications of our hypothesis outperformed other models across all 22 subjects, as measured by Akaike's Information Criterion (AIC), as well as mean and maximum distance from the diagonal on a quantile-quantile plot. Our broader model set of four simplifications offered a useful framework to enhance further statistical descriptions of EDA pulse amplitudes. Some of the simplifications prioritize fit near the mode of the distribution, while others prioritize fit near the tail. With this new insight, we can summarize the physiologically-relevant amplitude information in EDA with at most four parameters. Our findings establish that physiologically based probability models provide parsimonious and accurate description of temporal and amplitude characteristics in EDA.

Emery Brown, M.D., Ph.D. is an American statistician, neuroscientist, and anesthesiologist. He is the Warren M. Zapol Professor of Anesthesia at Harvard Medical School and at Massachusetts General Hospital(MGH), and a practicing anesthesiologist at MGH. At MIT he is the Edward Hood Taplin Professor of Medical Engineering and professor of computational neuroscience, the Associate Director of the Institute for Medical Engineering and Science, and the Director of the Harvard-MIT Program in Health Sciences and Technology. Brown is one of only 19 individuals who has been elected to all three branches of the National Academies of Sciences, Engineering, and Medicine, Brown is also the first African American and first anesthesiologist to be elected to all three National Academies.

High quality spectra are crucial in anesthesia related procedures where Electroencephalogram (EEG) frequency content can drastically help track different brain states. To this end, an adaptive autoregressive model framework to fit non-stationary EEG data using hybrid Kalman Filtering (HKF) is developed. In this setup, a state-space formulation is adopted. Hybridity arises from the fact that the state vector, which includes the autoregressive parameters, evolves in continuous time, while the observation equation is discrete, to account for the fact that the observations, i.e. EEG data, are discrete points in time. As a smoothing constraint, the parameters are modeled to follow a continuous multivariate random walk. As shown in this work, their adaptive estimation by means of HKF and expectation-maximization (EM) algorithm yields smoother estimation of frequency spectra, outperforming other current purely discrete parametric methods as well as various non-parametric approaches. Examples of dynamic EEG data were taken from patients under gradually varying doses of propofol. The suitability of our method for online use in combination with its ability to smoothly track frequency changes in human EEG signals suggests that it can be used for real time brain state tracking under general anesthesia, facilitating the design of closed loop systems for automatic and precise control of brain states.

Spectral properties of the electroencephalogram (EEG) are commonly analyzed to characterize the brain's oscillatory properties in basic science and clinical neuroscience studies. The spectrum is a function that describes power as a function of frequency. To date inference procedures for spectra have focused on constructing confidence intervals at single frequencies using large sample-based analytic procedures or jackknife techniques. These procedures perform well when the frequencies of interest are chosen before the analysis. When these frequencies are chosen after some of the data have been analyzed, the validity of these conditional inferences is not addressed. If power at more than one frequency is investigated, corrections for multiple comparisons must also be incorporated. To develop a statistical inference approach that considers the spectrum as a function defined across frequencies, we combine multitaper spectral methods with a frequency-domain bootstrap (FDB) procedure. The multitaper method is optimal for minimizing the bias-variance tradeoff in spectral estimation. The FDB makes it possible to conduct Monte Carlo based inferences for any part of the spectrum by drawing random samples that respect the dependence structure in the EEG time series. We show that our multitaper FDB procedure performs well in simulation studies and in analyses comparing EEG recordings of children from two different age groups receiving general anesthesia.

Clinical researchers, historians, educators and field researchers alike still regularly capture data on paper spreadsheets. In the case of health care and education, data will often contain sensitive personal information, further complicating the process of transcribing paper-based archives into digital form. In this work, we describe a tool that utilizes machine learning and crowd intelligence to automatically transcribe images of paper-based spreadsheets into electronic form while protecting sensitive personal information. Our solution consists of four high-level stages: (1) the extraction of cell-level images from the spreadsheet grid, (2) machine recognition of digits within the cells, (3) human transcription of cell contents that the machine was uncertain of and (4) feedback of human transcription results to the machine to improve future classification performance. We test the algorithm on a novel data-set of 135 heterogeneous clinical flow-sheet images collected from the Massachusetts General Hospital (MGH), 2 hand-drawn spreadsheets, one chalk-board drawing, and one printed table. we demonstrate that our algorithm provides a generalized solution for spreadsheet transcription that maintains privacy, is up to 10 times faster and twice as cost effective than existing alternatives. Our work is valuable both as a tool and as a starting point for the development of better algorithms.

We consider the problem of estimating the power spectral density of the neural covariates underlying the spiking of a neuronal population. We assume the spiking of the neuronal ensemble to be described by Bernoulli statistics. Furthermore, we consider the conditional intensity function to be the logistic map of a second-order stationary process with sparse frequency content. Using the binary spiking data recorded from the population, we calculate the maximum a posteriori estimate of the power spectral density of the process while enforcing sparsity-promoting priors on the estimate. Using both simulated and clinically recorded data, we show that our method outperforms the existing methods for extracting a frequency domain representation from the spiking data of a neuronal population.

The anesthetic propofol elicits many different spectral properties on the EEG, including alpha oscillations (8-12 Hz), Slow Wave Oscillations (SWO, 0.1-1.5 Hz), and dose-dependent phase-amplitude coupling (PAC) between alpha and SWO. Propofol is known to increase GABAA inhibition and decrease H-current strength, but how it generates these rhythms and their interactions is still unknown. To investigate both generation of the alpha rhythm and its PAC to SWO, we simulate a Hodgkin-Huxley network model of a hyperpolarized thalamus and corticothalamic inputs. We find, for the first time, that the model thalamic network is capable of independently generating the sustained alpha seen in propofol, which may then be relayed to cortex and expressed on the EEG. This dose-dependent sustained alpha critically relies on propofol GABAA potentiation to alter the intrinsic spindling mechanisms of the thalamus. Furthermore, the H-current conductance and background excitation of these thalamic cells must be within specific ranges to exhibit any intrinsic oscillations, including sustained alpha. We also find that, under corticothalamic SWO UP and DOWN states, thalamocortical output can exhibit maximum alpha power at either the peak or trough of this SWO; this implies the thalamus may be the source of propofol-induced PAC. Hyperpolarization level is the main determinant of whether the thalamus exhibits trough-max PAC, which is associated with lower propofol dose, or peak-max PAC, associated with higher dose. These findings suggest: the thalamus generates a novel rhythm under GABAA potentiation such as under propofol, its hyperpolarization may determine whether a patient experiences trough-max or peak-max PAC, and the thalamus is a critical component of propofol-induced cortical spectral phenomena. Changes to the thalamus may be a critical part of how propofol accomplishes its effects, including unconsciousness.

Common spatial patterns (CSP) is a well-known spatial filtering algorithm for multichannel electroencephalogram (EEG) analysis. In this paper, we cast the CSP algorithm in a probabilistic modeling setting. Specifically, probabilistic CSP (P-CSP) is proposed as a generic EEG spatio-temporal modeling framework that subsumes the CSP and regularized CSP algorithms. The proposed framework enables us to resolve the overfitting issue of CSP in a principled manner. We derive statistical inference algorithms that can alleviate the issue of local optima. In particular, an efficient algorithm based on eigendecomposition is developed for maximum a posteriori (MAP) estimation in the case of isotropic noise. For more general cases, a variational algorithm is developed for group-wise sparse Bayesian learning for the P-CSP model and for automatically determining the model size. The two proposed algorithms are validated on a simulated data set. Their practical efficacy is also demonstrated by successful applications to single-trial classifications of three motor imagery EEG data sets and by the spatio-temporal pattern analysis of one EEG data set recorded in a Stroop color naming task.

Likelihood-based encoding models founded on point processes have received significant attention in the literature because of their ability to reveal the information encoded by spiking neural populations. We propose an approximation to the likelihood of a point-process model of neurons that holds under assumptions about the continuous time process that are physiologically reasonable for neural spike trains: the presence of a refractory period, the predictability of the conditional intensity function, and its integrability. These are properties that apply to a large class of point processes arising in applications other than neuroscience. The proposed approach has several advantages over conventional ones. In particular, one can use standard fitting procedures for generalized linear models based on iteratively reweighted least squares while improving the accuracy of the approximation to the likelihood and reducing bias in the estimation of the parameters of the underlying continuous-time model. As a result, the proposed approach can use a larger bin size to achieve the same accuracy as conventional approaches would with a smaller bin size. This is particularly important when analyzing neural data with high mean and instantaneous firing rates. We demonstrate these claims on simulated and real neural spiking activity. By allowing a substantive increase in the required bin size, our algorithm has the potential to lower the barrier to the use of point-process methods in an increasing number of applications.

Although the existence of correlated spiking between neurons in a population is well known, the role such correlations play in encoding stimuli is not. We address this question by constructing pattern-based encoding models that describe how time-varying stimulus drive modulates the expression probabilities of population-wide spike patterns. The challenge is that large populations may express an astronomical number of unique patterns, and so fitting a unique encoding model for each individual pattern is not feasible. We avoid this combinatorial problem using a dimensionality-reduction approach based on regression trees. Using the insight that some patterns may, from the perspective of encoding, be statistically indistinguishable, the tree divisively clusters the observed patterns into groups whose member patterns possess similar encoding properties. These groups, corresponding to the leaves of the tree, are much smaller in number than the original patterns, and the tree itself constitutes a tractable encoding model for each pattern. Our formalism can detect an extremely weak stimulus-driven pattern structure and is based on maximizing the data likelihood, not making a priori assumptions as to how patterns should be grouped. Most important, by comparing pattern encodings with independent neuron encodings, one can determine if neurons in the population are driven independently or collectively. We demonstrate this method using multiple unit recordings from area 17 of anesthetized cat in response to a sinusoidal grating and show that pattern-based encodings are superior to those of independent neuron models. The agnostic nature of our clustering approach allows us to investigate encoding by the collective statistics that are actually present rather than those (such as pairwise) that might be presumed.

Medically-induced coma is a drug-induced state of profound brain inactivation and unconsciousness used to treat refractory intracranial hypertension and to manage treatment-resistant epilepsy. The state of coma is achieved by continually monitoring the patient's brain activity with an electroencephalogram (EEG) and manually titrating the anesthetic infusion rate to maintain a specified level of burst suppression, an EEG marker of profound brain inactivation in which bursts of electrical activity alternate with periods of quiescence or suppression. The medical coma is often required for several days. A more rational approach would be to implement a brain-machine interface (BMI) that monitors the EEG and adjusts the anesthetic infusion rate in real time to maintain the specified target level of burst suppression. We used a stochastic control framework to develop a BMI to control medically-induced coma in a rodent model. The BMI controlled an EEG-guided closed-loop infusion of the anesthetic propofol to maintain precisely specified dynamic target levels of burst suppression. We used as the control signal the burst suppression probability (BSP), the brain's instantaneous probability of being in the suppressed state. We characterized the EEG response to propofol using a two-dimensional linear compartment model and estimated the model parameters specific to each animal prior to initiating control. We derived a recursive Bayesian binary filter algorithm to compute the BSP from the EEG and controllers using a linear-quadratic-regulator and a model-predictive control strategy. Both controllers used the estimated BSP as feedback. The BMI accurately controlled burst suppression in individual rodents across dynamic target trajectories, and enabled prompt transitions between target levels while avoiding both undershoot and overshoot. The median performance error for the BMI was 3.6%, the median bias was -1.4% and the overall posterior probability of reliable control was 1 (95% Bayesian credibility interval of [0.87, 1.0]). A BMI can maintain reliable and accurate real-time control of medically-induced coma in a rodent model suggesting this strategy could be applied in patient care.

We consider the problem of recovering a sequence of vectors, $(x_k)_{k=0}^K$, for which the increments $x_k-x_{k-1}$ are $S_k$-sparse (with $S_k$ typically smaller than $S_1$), based on linear measurements $(y_k = A_k x_k + e_k)_{k=1}^K$, where $A_k$ and $e_k$ denote the measurement matrix and noise, respectively. Assuming each $A_k$ obeys the restricted isometry property (RIP) of a certain order---depending only on $S_k$---we show that in the absence of noise a convex program, which minimizes the weighted sum of the $\ell_1$-norm of successive differences subject to the linear measurement constraints, recovers the sequence $(x_k)_{k=1}^K$ \emph{exactly}. This is an interesting result because this convex program is equivalent to a standard compressive sensing problem with a highly-structured aggregate measurement matrix which does not satisfy the RIP requirements in the standard sense, and yet we can achieve exact recovery. In the presence of bounded noise, we propose a quadratically-constrained convex program for recovery and derive bounds on the reconstruction error of the sequence. We supplement our theoretical analysis with simulations and an application to real video data. These further support the validity of the proposed approach for acquisition and recovery of signals with time-varying sparsity.

The moving bar experiment is a classic paradigm for characterizing the receptive field (RF) properties of neurons in primary visual cortex (V1). Current approaches for analyzing neural spiking activity recorded from these experiments do not take into account the point-process nature of these data and the circular geometry of the stimulus presentation. We present a novel analysis approach to mapping V1 receptive fields that combines point-process generalized linear models (PPGLM) with tomographic reconstruction computed by filtered-back projection. We use the method to map the RF sizes and orientations of 251 V1 neurons recorded from two macaque monkeys during a moving bar experiment. Our cross-validated goodness-of-fit analyses show that the PPGLM provides a more accurate characterization of spike train data than analyses based on rate functions computed by the methods of spike-triggered averages or first-order Wiener-Volterra kernel. Our analysis leads to a new definition of RF size as the spatial area over which the spiking activity is significantly greater than baseline activity. Our approach yields larger RF sizes and sharper orientation tuning estimates. The tomographic reconstruction paradigm further suggests an efficient approach to choosing the number of directions and the number of trials per direction in designing moving bar experiments. Our results demonstrate that standard tomographic principles for image reconstruction can be adapted to characterize V1 RFs and that two fundamental properties, size and orientation, may be substantially different from what is currently reported.

Precise spike coordination between the spiking activities of multiple neurons is suggested as an indication of coordinated network activity in active cell assemblies. Spike correlation analysis aims to identify such cooperative network activity by detecting excess spike synchrony in simultaneously recorded multiple neural spike sequences. Cooperative activity is expected to organize dynamically during behavior and cognition; therefore currently available analysis techniques must be extended to enable the estimation of multiple time-varying spike interactions between neurons simultaneously. In particular, new methods must take advantage of the simultaneous observations of multiple neurons by addressing their higher-order dependencies, which cannot be revealed by pairwise analyses alone. In this paper, we develop a method for estimating time-varying spike interactions by means of a state-space analysis. Discretized parallel spike sequences are modeled as multi-variate binary processes using a log-linear model that provides a well-defined measure of higher-order spike correlation in an information geometry framework. We construct a recursive Bayesian filter/smoother for the extraction of spike interaction parameters. This method can simultaneously estimate the dynamic pairwise spike interactions of multiple single neurons, thereby extending the Ising/spin-glass model analysis of multiple neural spike train data to a nonstationary analysis. Furthermore, the method can estimate dynamic higher-order spike interactions. To validate the inclusion of the higher-order terms in the model, we construct an approximation method to assess the goodness-of-fit to spike data. In addition, we formulate a test method for the presence of higher-order spike correlation even in nonstationary spike data, e.g., data from awake behaving animals. The utility of the proposed methods is tested using simulated spike data with known underlying correlation dynamics. Finally, we apply the methods to neural spike data simultaneously recorded from the motor cortex of an awake monkey and demonstrate that the higher-order spike correlation organizes dynamically in relation to a behavioral demand.

Characterizing neural spiking activity as a function of intrinsic and extrinsic factors is important in neuroscience. Point process models are valuable for capturing such information; however, the process of fully applying these models is not always obvious. A complete model application has four broad steps: specification of the model, estimation of model parameters given observed data, verification of the model using goodness of fit, and characterization of the model using confidence bounds. Of these steps, only the first three have been applied widely in the literature, suggesting the need to dedicate a discussion to how the time-rescaling theorem, in combination with parametric bootstrap sampling, can be generally used to compute confidence bounds of point process models. In our first example, we use a generalized linear model of spiking propensity to demonstrate that confidence bounds derived from bootstrap simulations are consistent with those computed from closed-form analytic solutions. In our second example, we consider an adaptive point process model of hippocampal place field plasticity for which no analytical confidence bounds can be derived. We demonstrate how to simulate bootstrap samples from adaptive point process models, how to use these samples to generate confidence bounds, and how to statistically test the hypothesis that neural representations at two time points are significantly different. These examples have been designed as useful guides for performing scientific inference based on point process models.

The ability to identify directional interactions that occur among multiple neurons in the brain is crucial to an understanding of how groups of neurons cooperate in order to generate specific brain functions. However, an optimal method of assessing these interactions has not been established. Granger causality has proven to be an effective method for the analysis of the directional interactions between multiple sets of continuous-valued data, but cannot be applied to neural spike train recordings due to their discrete nature. This paper proposes a point process framework that enables Granger causality to be applied to point process data such as neural spike trains. The proposed framework uses the point process likelihood function to relate a neuron's spiking probability to possible covariates, such as its own spiking history and the concurrent activity of simultaneously recorded neurons. Granger causality is assessed based on the relative reduction of the point process likelihood of one neuron obtained excluding one of its covariates compared to the likelihood obtained using all of its covariates. The method was tested on simulated data, and then applied to neural activity recorded from the primary motor cortex (MI) of a Felis catus subject. The interactions present in the simulated data were predicted with a high degree of accuracy, and when applied to the real neural data, the proposed method identified causal relationships between many of the recorded neurons. This paper proposes a novel method that successfully applies Granger causality to point process data, and has the potential to provide unique physiological insights when applied to neural spike trains.

A common feature in biological neuromuscular systems is the redundancy in joint actuation. Understanding how these redundancies are resolved in typical joint movements has been a long-standing problem in biomechanics, neuroscience and prosthetics. Many empirical studies have uncovered neural, mechanical and energetic aspects of how humans resolve these degrees of freedom to actuate leg joints for common tasks like walking. However, a unifying theoretical framework that explains the many independent empirical observations and predicts individual muscle and tendon contributions to joint actuation is yet to be established. Here we develop a computational framework to address how the ankle joint actuation problem is resolved by the neuromuscular system in walking. Our framework is founded upon the proposal that a consideration of both neural control and leg muscle-tendon morphology is critical to obtain predictive, mechanistic insight into individual muscle and tendon contributions to joint actuation. We examine kinetic, kinematic and electromyographic data from healthy walking subjects to find that human leg muscle-tendon morphology and neural activations enable a metabolically optimal realization of biological ankle mechanics in walking. This optimal realization (a) corresponds to independent empirical observations of operation and performance of the soleus and gastrocnemius muscles, (b) gives rise to an efficient load-sharing amongst ankle muscle-tendon units and (c) causes soleus and gastrocnemius muscle fibers to take on distinct mechanical roles of force generation and power production at the end of stance phase in walking. The framework outlined here suggests that the dynamical interplay between leg structure and neural control may be key to the high walking economy of humans, and has implications as a means to obtain insight into empirically inaccessible features of individual muscle and tendons in biomechanical tasks.

Point process generalized linear models (GLMs) have been widely used for neural spike trains analysis. Statistical inference for GLMs include maximum likelihood and Bayesian estimation. Variational Bayesian (VB) methods provide a computationally appealing means to infer the posterior density of unknown parameters, in which conjugate priors are designed for the regression coefficients in logistic and Poisson regression. In this paper, we develop and apply VB inference for point process GLMs in neural spike train analysis. The hierarchical Bayesian framework allows us to tackle the variable selection problem. We assess and validate our methods with ensemble neuronal recordings from rat's hippocampal place cells and entorhinal cortical cells during foraging in an open field environment.

Assessing the causal relationship among multivariate time series is a crucial problem in many fields. Granger causality has been widely used to identify the causal interactions between continuous-valued time series based on multivariate autoregressive models in the Gaussian case. In order to extend the application of the Granger causality concept to non-Gaussian time series, we propose a general statistical framework for assessing the causal interactions. In this study, the Granger causality from a time series x2to a time series x1is assessed based on the relative reduction of the likelihood of x1by the exclusion of x2compared to the likelihood obtained using all the time series. Simulation results indicated that the proposed algorithm accurately predicted nature of interactions between discrete-valued time series as well as between continuous-valued time series.

We develop a state space approach to multimodal integration of simultaneously recorded EEG and fMRI. The EEG is represented with a distributed current source model using realistic MRI-based forward models, whose temporal evolution is governed by a linear state space model. The fMRI signal is similarly modeled by a linear state space model describing the hemodynamic response to underlying EEG current activity. We explore the feasibility of high dimensional dynamic estimation of simultaneous EEG/fMRI using simulation studies of the alpha wave.

Brain machine interfaces work by mapping the relevant neural activity to the intended movement known as `decoding'. Here, we develop a recursive Bayesian decoder for goal-directed movements from neural observations, which exploits the optimal feedback control model of the sensorimotor system to build better prior state-space models. These controlled state models depend on the movement duration that is not known a priori. We thus consider a discretization of the task duration and develop a decoder consisting of a bank of parallel point-process filters, each combining the neural observation with the controlled state model of a discretization point. The final reconstruction is made by optimally combining these filter estimates. Using very coarse discretization and hence only a few parallel branches, our decoder reduces the root mean square (RMS) error in trajectory reconstruction in reaches made by a rhesus monkey by approximately 40%.

In numerous neuroscience studies, multichannel EEG data are often recorded over multiple trial periods under the same experimental condition. To date, little effort is aimed to learn spatial patterns from EEG data to account for trial-to-trial variability. In this paper, a hierarchical Bayesian framework is introduced to model inter-trial source variability while extracting common spatial patterns under multiple experimental conditions in a supervised manner. We also present a variational Bayesian algorithm for model inference, by which the number of sources can be determined effectively via automatic relevance determination (ARD). The efficacy of the proposed learning algorithm is validated with both synthetic and real EEG data. Using two brain-computer interface (BCI) motor imagery data sets we show the proposed algorithm consistently outperforms the common spatial patterns (CSP) algorithm while attaining comparable performance with a recently proposed discriminative approach.

One approach for understanding the encoding of information by spike trains is to fit statistical models and then test their goodness of fit. The time-rescaling theorem provides a goodness-of-fit test consistent with the point process nature of spike trains. The interspike intervals (ISIs) are rescaled (as a function of the model's spike probability) to be independent and exponentially distributed if the model is accurate. A Kolmogorov-Smirnov (KS) test between the rescaled ISIs and the exponential distribution is then used to check goodness of fit. This rescaling relies on assumptions of continuously defined time and instantaneous events. However, spikes have finite width, and statistical models of spike trains almost always discretize time into bins. Here we demonstrate that finite temporal resolution of discrete time models prevents their rescaled ISIs from being exponentially distributed. Poor goodness of fit may be erroneously indicated even if the model is exactly correct. We present two adaptations of the time-rescaling theorem to discrete time models. In the first we propose that instead of assuming the rescaled times to be exponential, the reference distribution be estimated through direct simulation by the fitted model. In the second, we prove a discrete time version of the time-rescaling theorem that analytically corrects for the effects of finite resolution. This allows us to define a rescaled time that is exponentially distributed, even at arbitrary temporal discretizations. We demonstrate the efficacy of both techniques by fitting generalized linear models to both simulated spike trains and spike trains recorded experimentally in monkey V1 cortex. Both techniques give nearly identical results, reducing the false-positive rate of the KS test and greatly increasing the reliability of model evaluation based on the time-rescaling theorem.

Evaluation of baroreflex control of heart rate (HR) has important implications in clinical practice of anesthesia and postoperative care. In this paper, we present a point process method to assess the dynamic baroreflex gain using a closed-loop model of the cardiovascular system. Specifically, the inverse Gaussian probability distribution is used to model the heartbeat interval, whereas the instantaneous mean is identified by a linear or bilinear bivariate regression on the previous R-R intervals and blood pressure (BP) measures. The instantaneous baroreflex gain is estimated in the feedback loop with a point process filter, while the RR→BP feedforward frequency response is estimated by a Kalman filter. In addition, the instantaneous cross-spectrum and cross-bispectrum (as well as their ratio) can also be estimated. All statistical indices provide a valuable quantitative assessment of the interaction between heartbeat dynamics and hemodynamics during general anesthesia.

A state-space method for simultaneously estimating time-dependent rate and higher-order correlation underlying parallel spike sequences is proposed. Discretized parallel spike sequences are modeled by a conditionally independent multivariate Bernoulli process using a log-linear link function, which contains a state of higher-order interaction factors. A nonlinear recursive filtering formula is derived from a log-quadratic approximation to the posterior distribution of the state. Together with a fixed-interval smoothing algorithm, time-dependent log-linear parameters are estimated. The smoothed estimates are optimized via EM-algorithm such that their prior covariance matrix maximizes the expected complete data log-likelihood. In addition, we perform model selection on the hierarchical log-linear state-space models to avoid over-fitting. Application of the method to simultaneously recorded neuronal spike sequences is expected to contribute to uncover dynamic cooperative activities of neurons in relation to behavior.

UP and DOWN states, the periodic fluctuations between increased and decreased spiking activity of a neuronal population, are a fundamental feature of cortical circuits. Understanding UP-DOWN state dynamics is important for understanding how these circuits represent and transmit information in the brain. To date, limited work has been done on characterizing the stochastic properties of UP-DOWN state dynamics. We present a set of Markov and semi-Markov discrete- and continuous-time probability models for estimating UP and DOWN states from multiunit neural spiking activity. We model multiunit neural spiking activity as a stochastic point process, modulated by the hidden (UP and DOWN) states and the ensemble spiking history. We estimate jointly the hidden states and the model parameters by maximum likelihood using an expectation-maximization (EM) algorithm and a Monte Carlo EM algorithm that uses reversible-jump Markov chain Monte Carlo sampling in the E-step. We apply our models and algorithms in the analysis of both simulated multiunit spiking activity and actual multi- unit spiking activity recorded from primary somatosensory cortex in a behaving rat during slow-wave sleep. Our approach provides a statistical characterization of UP-DOWN state dynamics that can serve as a basis for verifying and refining mechanistic descriptions of this process.

In this paper, we compare and validate different probabilistic models of human heart beat intervals for assessment of the electrocardiogram data recorded with varying conditions in posture and pharmacological autonomic blockade. The models are validated using the adaptive point process filtering paradigm and Kolmogorov-Smirnov test. The inverse Gaussian model was found to achieve the overall best performance in the analysis of autonomic control. We further improve the model by incorporating the respiratory covariate measurements and present dynamic respiratory sinus arrhythmia (RSA) analysis. Our results suggest the instantaneous RSA gain computed from our proposed model as a potential index of vagal control dynamics.

Electrophysiological recordings of brain activity include point process spike trains as well as continuous valued signals such as electroencephalograms (EEG), electrocorticograms (ECoG), and local field potentials (LFP). The brain represents information about the outside world in neural spiking activity, which is reflected in each of these signal modalities. An important problem in neuroscience data analysis involves estimating dynamic biological and behavioral signals from neural recordings. Here, we develop an adaptive filtering paradigm for estimating dynamic state processes from mixed observation processes that contain both point process and continuous valued observations. In our analysis of these filtering algorithms, we draw analogies to well-studied linear estimation algorithms such as the Kalman and Extended Kalman filters. We demonstrate the application of this mixed filtering paradigm to the problem of estimating a reaching movement trajectory from simulated simultaneously recorded motor cortical spiking and LFP activity. We demonstrate that the mixed filter is better able to capture information about the movement trajectory than are filters based on the spiking activity or LFPs alone.

The execution of reaching movements involves the coordinated activity of multiple brain regions that relate variously to the desired target and a path of arm states to achieve that target. These arm states may represent positions, velocities, torques, or other quantities. Estimation has been previously applied to neural activity in reconstructing the target separately from the path. However, the target and path are not independent. Because arm movements are limited by finite muscle contractility, knowledge of the target constrains the path of states that leads to the target. In this letter, we derive and illustrate a state equation to capture this basic dependency between target and path. The solution is described for discrete-time linear systems and gaussian increments with known target arrival time. The resulting analysis enables the use of estimation to study how brain regions that relate variously to target and path together specify a trajectory. The corresponding reconstruction procedure may also be useful in brain-driven prosthetic devices to generate control signals for goal-directed movements.

Analyzing the dependencies between spike trains is an important step in understanding how neurons work in concert to represent biological signals. Usually this is done for pairs of neurons at a time using correlation-based techniques. Chornoboy, Schramm, and Karr (1988) proposed maximum likelihood methods for the simultaneous analysis of multiple pair-wise interactions among an ensemble of neurons. One of these methods is an iterative, continuous-time estimation algorithm for a network likelihood model formulated in terms of multiplicative conditional intensity functions. We devised a discrete-time version of this algorithm that includes a new, efficient computational strategy, a principled method to compute starting values, and a principled stopping criterion. In an analysis of simulated neural spike trains from ensembles of interacting neurons, the algorithm recovered the correct connectivity matrices and interaction parameters. In the analysis of spike trains from an ensemble of rat hippocampal place cells, the algorithm identified a connectivity matrix and interaction parameters consistent with the pattern of conjoined firing predicted by the overlap of the neurons' spatial receptive fields. These results suggest that the network likelihood model can be an efficient tool for the analysis of ensemble spiking activity.

In order to construct spatial activation plots from functional magnetic resonance imaging (fMRI) data, a complex spatio-temporal modeling problem must be solved. A crucial part of this process is the estimation of the hemodynamic response (HR) function, an impulse response relating the stimulus signal to the measured noisy response. The estimation of the HR is complicated by the presence of low frequency colored noise. The standard approach to modeling the HR is to use simple parametric models, although FIR models have been used. We offer two contributions. First, we pursue a nonparametric approach using orthonormal causal Laguerre polynomials which have become popular in the system identification literature. It also happens that the shape of the basis elements is similar to that of a typical HR. We thus expect to achieve a compact, and so bias reduced, and low noise representation of the HR. Additionally, we develop a procedure for providing confidence intervals for the whole HR function. This feature is completely lacking in all previous work.

In order to construct spatial activation plots from functional magnetic resonance imaging (fMRI) data, a complex spatio-temporal modeling problem must be solved. A crucial part of this process is the estimation of the hemodynamic response (HR) function, an impulse response relating the stimulus signal to the measured noisy response. The estimation of the HR is complicated by the presence of low frequency colored noise. The standard approach to modeling the HR is to use simple parametric models, although FIR models have been used. We pursue a nonparametric approach using orthonormal causal Laguerre polynomials which have become popular in the system identification literature. It also happens that the shape of the basis elements is similar to that of a typical HR. We thus expect to achieve a compact and so bias reduced and low noise representation of the HR. This is not the case in FIR modeling, because a low FIR order is unable to cover the whole length of the HR over its region of support while a high FIR order results in overestimation of signal and underestimation of noise leading to misleading interpretations.

Neural spike train decoding algorithms and techniques to compute Shannon mutual information are important methods for analyzing how neural systems represent biological signals. Decoding algorithms are also one of several strategies being used to design controls for brain-machine interfaces. Developing optimal strategies to design decoding algorithms and compute mutual information are therefore important problems in computational neuroscience. We present a general recursive filter decoding algorithm based on a point process model of individual neuron spiking activity and a linear stochastic state-space model of the biological signal. We derive from the algorithm new instantaneous estimates of the entropy, entropy rate, and the mutual information between the signal and the ensemble spiking activity. We assess the accuracy of the algorithm by computing, along with the decoding error, the true coverage probability of the approximate 0.95 confidence regions for the individual signal estimates. We illustrate the new algorithm by reanalyzing the position and ensemble neural spiking activity of CA1 hippocampal neurons from two rats foraging in an open circular environment. We compare the performance of this algorithm with a linear filter constructed by the widely used reverse correlation method. The median decoding error for Animal 1 (2) during 10 minutes of open foraging was 5.9 (5.5) cm, the median entropy was 6.9 (7.0) bits, the median information was 9.4 (9.4) bits, and the true coverage probability for 0.95 confidence regions was 0.67 (0.75) using 34 (32) neurons. These findings improve significantly on our previous results and suggest an integrated approach to dynamically reading neural codes, measuring their properties, and quantifying the accuracy with which encoded information is extracted.

Neural receptive fields are dynamic in that with experience, neurons change their spiking responses to relevant stimuli. To understand how neural systems adapt their representations of biological information, analyses of receptive field plasticity from experimental measurements are crucial. Adaptive signal processing, the well-established engineering discipline for characterizing the temporal evolution of system parameters, suggests a framework for studying the plasticity of receptive fields. We use the Bayes' rule Chapman-Kolmogorov paradigm with a linear state equation and point process observation models to derive adaptive filters appropriate for estimation from neural spike trains. We derive point process filter analogues of the Kalman filter, recursive least squares, and steepest-descent algorithms and describe the properties of these new filters. We illustrate our algorithms in two simulated data examples. The first is a study of slow and rapid evolution of spatial receptive fields in hippocampal neurons. The second is an adaptive decoding study in which a signal is decoded from ensemble neural spiking activity as the receptive fields of the neurons in the ensemble evolve. Our results provide a paradigm for adaptive estimation for point process observations and suggest a practical approach for constructing filtering algorithms to track neural receptive field dynamics on a millisecond timescale.

Within the context of a properly formulated spatio-temporal model for functional magnetic resonance imaging (fMRI), we develop a procedure for wavelet-based activation estimation that explicitly accounts for the temporally correlated 'brain' noise. This leads to a kind of resolution level - dependent wavelet thresholding, a novelty in the spatial setting. The method is illustrated on a visual A-B fMRI experiment and is seen to give significantly better results compared with those obtainable via a more conventional smoothing step involving preprocessing with a spatially homogeneous filter.

A widely used signal processing paradigm is the state-space model. The state-space model is defined by two equations: an observation equation that describes how the hidden state or latent process is observed and a state equation that defines the evolution of the process through time. Inspired by neurophysiology experiments in which neural spiking activity is induced by an implicit (latent) stimulus, we develop an algorithm to estimate a state-space model observed through point process measurements. We represent the latent process modulating the neural spiking activity as a gaussian autoregressive model driven by an external stimulus. Given the latent process, neural spiking activity is characterized as a general point process defined by its conditional intensity function. We develop an approximate expectation-maximization (EM) algorithm to estimate the unobservable state-space process, its parameters, and the parameters of the point process. The EM algorithm combines a point process recursive nonlinear filter algorithm, the fixed interval smoothing algorithm, and the state-space covariance algorithm to compute the complete data log likelihood efficiently. We use a Kolmogorov-Smirnov test based on the time-rescaling theorem to evaluate agreement between the model and point process data. We illustrate the model with two simulated data examples: an ensemble of Poisson neurons driven by a common stimulus and a single neuron whose conditional intensity function is approximated as a local Bernoulli process.

Measuring agreement between a statistical model and a spike train data series, that is, evaluating goodness of fit, is crucial for establishing the model's validity prior to using it to make inferences about a particular neural system. Assessing goodness-of-fit is a challenging problem for point process neural spike train models, especially for histogram-based models such as perstimulus time histograms (PSTH) and rate functions estimated by spike train smoothing. The time-rescaling theorem is a well-known result in probability theory, which states that any point process with an integrable conditional intensity function may be transformed into a Poisson process with unit rate. We describe how the theorem may be used to develop goodness-of-fit tests for both parametric and histogram-based point process models of neural spike trains. We apply these tests in two examples: a comparison of PSTH, inhomogeneous Poisson, and inhomogeneous Markov interval models of neural spike trains from the supplementary eye field of a macque monkey and a comparison of temporal and spatial smoothers, inhomogeneous Poisson, inhomogeneous gamma, and inhomogeneous inverse gaussian models of rat hippocampal place cell spiking activity. To help make the logic behind the time-rescaling theorem more accessible to researchers in neuroscience, we present a proof using only elementary probability theory arguments. We also show how the theorem may be used to simulate a general point process model of a spike train. Our paradigm makes it possible to compare parametric and histogram-based neural spike train models directly. These results suggest that the time-rescaling theorem can be a valuable tool for neural spike train data analysis.

We consider signal estimation for functional MRI studies on multiple subjects. There are two major issues; alignment or registration of images across subjects, and using the multisubject information to capture covariance information; we discuss only the latter. Capturing this covariance information properly can lead to great improvements in statistical efficiency beyond what simple averaging can offer as well as compact description of group features.

In the last half decade, fast methods of magnetic resonance imaging have led to the possibility, for the first time, of non-invasive dynamic brain imaging. This has led to an explosion of work in the Neurosciences. From a signal processing viewpoint the problems are those of nonlinear spatio-temporal system identification. In this paper, we develop new methods of identification using novel spatial regularization. We also develop a new model comparison technique and use that to compare our method with existing techniques on some experimental data.

We consider a model comparison based on a new model selection criterion. We treat fMRI as a spatio-temporal system identification problem and compare our model fitting method based on spatial regularization methods with the so-called statistical parametric map technique currently popular in the fMRI literature. We illustrate our results with data from the brain obtained during a combined visual and motor experiment.

The authors consider spatio-temporal modelling of functional MRI data from an inverse problems point of view. Most of the modelling to date has been on a pixel by pixel basis with no acknowledgment given to spatial smoothness when it exists. The authors discuss regularization methods that address that issue and illustrate their results with experimental data.

Discusses modelling of functional magnetic resonance imaging (fMRI) data as an ill-conditioned inverse problem. The model involved is non-linear and requires both temporal and spatial regularization.