Oscillations are observed at various frequency bands in continuous-valued neural recordings like the electroencephalogram (EEG) and local field potential (LFP) in bulk brain matter, and analysis of spike-field coherence reveals that spiking of single neurons often occurs at certain phases of the global oscillation. Oscillatory modulation has been examined in relation to continuous-valued oscillatory signals, and independently from the spike train alone, but behavior or stimulus triggered firing-rate modulation, spiking sparseness, presence of slow modulation not locked to stimuli and irregular oscillations with large variability in oscillatory periods, present challenges to searching for temporal structures present in the spike train. In order to study oscillatory modulation in real data collected under a variety of experimental conditions, we describe a flexible point-process framework we call the Latent Oscillatory Spike Train (LOST) model to decompose the instantaneous firing rate in biologically and behaviorally relevant factors: spiking refractoriness, event-locked firing rate non-stationarity, and trial-to-trial variability accounted for by baseline offset and a stochastic oscillatory modulation. We also extend the LOST model to accommodate changes in the modulatory structure over the duration of the experiment, and thereby discover trial-to-trial variability in the spike-field coherence of a rat primary motor cortical neuron to the LFP theta rhythm. Because LOST incorporates a latent stochastic auto-regressive term, LOST is able to detect oscillations when the firing rate is low, the modulation is weak, and when the modulating oscillation has a broad spectral peak.
Neural oscillations have generated considerable interest for their roles in cognition and as indicators for disease [1–14]. Electroencephalograms (EEGs) and local field potential (LFPs) recordings have revealed transient oscillations in many cortical and subcortical structures, to which neurons both near and far from the LFP recording site show phase preferences in spiking. Groups of neurons that are locked to a common oscillation, and are therefore active in tightly confined temporal windows, may define a cell assembly whose synchronous spiking could select and activate afferent structures [15–17]. Such transient cell assemblies may form and disband as objects are attended to in the visual scene [3,6–8,11], as preparation for movements are made [2,13], or during choice points for reward [9]. Prominent slow oscillations are observed as patients undergo anesthesia [10], and exhibit changes in dynamics during transitions in unconscious brain state [14]. Furthermore, strong oscillatory signals are prominent in neurological disorders such as Parkinson’s disease, and can be used to characterize pathophysiology and its relation to behavior [12,18,19]. The oscillatory structure of individual neurons can thus provide insight into the dynamics of cognitive processing and motor planning, and their implementation in health and disease. However, identifying oscillations in neural spike trains, particularly within individual experimental trials, can be challenging because they occur in conjunction with both non-oscillatory fluctuations in neural firing rate and extraneous noise. To tease apart multiple factors that affect spiking activity, statisticians have suggested the use of point process models together with modern regression methods associated with generalized linear model (GLM) technology [20,21]. In this paper we develop a method that can find oscillations in spike trains even when the signal is comparatively weak, and can track changes in oscillatory behavior across trials. The starting point for our approach is the observation by Smith and Brown [22] that, for many purposes, evolving neural firing rates can be described using state-space models imbedded in point processes. We modify the models used by Smith and Brown, replacing their first-order autoregressive processes with higher-order processes [23] that can fit oscillations found in spiking neurons. This requires careful attention to the form of the higher-order autoregressive process. We take advantage of a Bayesian time-series decomposition introduced by Huerta and West [24], using prior probability distributions to constrain the fit so that it has appropriate power spectral content. We also use a Gibbs sampling method developed recently in a different context [25] to compute posterior distributions efficiently. Our approach is superficially related to that of Allcroft et al. [26], who used autoregressive moving average models to estimate spectral content from censored data, but their method does not accommodate easily the special situation we face with spiking neurons, including history effects that can account for hard and soft refractory periods. Point process models employing the GLM framework [21,27] have been employed to explain the spiking behavior by taking into account spiking refractoriness, behavioral and stimulus-induced non-stationarities, and spiking of neighboring neurons [20]. A series of investigations [12,18,19] have added long-term history effects to capture oscillation in the history dependence and have also used a state-space smoothing algorithm to track changes in the oscillatory dynamics over time. Eden et al used a long history to capture both, but as we shall see, LOST can better capture the irregularities present in realistic biological oscillations. Our approach is based on what we call the Latent Oscillatory Spike Train (LOST) model, which not only includes terms to describe the spiking refractoriness, event (behavior or stimulus) -locked effects, and trial-to-trial variability in firing rate, but also explicitly models the dynamics of the modulating oscillation itself. However, the latent state cannot capture the discontinuous change in the spiking probability that occurs after every spike, which the aforementioned GLM models do well in characterizing, so we utilize a short-term spiking history to capture the refractory period, while also introducing the latent state to model the oscillatory dynamics. Oscillatory structure is often considerably degraded when it is observed as a spike train. On the one hand, if the firing rate is low, or the oscillatory modulation is weak, noise associated with erratic spiking will dominate. On the other hand, even if the modulation is strong and the firing rate is high, the oscillatory signal may itself be unsteady, with substantial variation in period from cycle to cycle. LOST accommodates both spiking noise and oscillatory variation by imbedding into the point process intensity function a latent auto-regressive process, and then imposing a suitable soft constraint (in the form of a prior probability distribution) on the oscillatory dynamics. We provide details of the model, along with the posterior sampling scheme and a discussion on how to interpret the posterior samples, and study its effectiveness with leaky integrate-and-fire neural simulations. We then demonstrate the ability of the LOST framework to find interesting trial-to-trial variation in motor cortical neuron during a lever-pulling task.
Gibbs sampling [28] and data augmentation [29,30] are used to infer the model parameters Θ = [F1: p, σ2, α, μ, v] and the oscillation x. The spike history and TAE are functions of time relative to the last spike of an event or behavior, respectively, and they are estimated by parameterizing a subset of continuous functions with splines. The estimation requires us to choose a set of basis splines, which we choose heuristically, using some prior belief about their general shape to choose the knot locations. In the current work, the number and locations of the knots are not found with the Gibbs sampling procedure, so we briefly describe how they are determined before we describe the Gibbs sampling procedure. We jointly sample model parameters and latent states from the joint posterior distribution, Eq (4), using Gibbs sampling with data augmentation. The conditional posterior distribution can be read off from the full joint posterior distribution by considering all parameters fixed, except the parameter whose conditional posterior we are interested in. Our problem deviates from the standard Gibbs sampling by the missing data x in Eq (4), and we utilize the data augmentation strategy [29,30], a scheme of simplifying analysis by augmenting the observed data with missing values, to generate samples of X from the predictive distribution p (X|y, Θ) in between sampling parameters from the standard conditional posteriors. Even with samples of X in hand, the posterior distribution does not yield conditional posteriors that can be sampled easily. The recently developed Pólya-Gamma data augmentation scheme [25] allows sampling from simple Gaussian conditional posteriors at the cost of introducing the M × N matrix of Pólya-Gamma variables ω through the following identity for the spike probability: p (y m n | Θ, H n, l, x m n) = (e x m n + μ m + λ l (m, n) R + f n) y m n 1 + e x m n + μ m + λ l (m, n) R + f n = 1 2 e κ m n (x m n + f n + μ m + λ l (m, n) R) × ∫ 0 ∞ e - ω m n (x m n + f n + μ m + λ l (m, n) R) 2 2 P G (ω m n | b = 1, z = 0) d ω m n ∝ ∫ 0 ∞ N (x m n + f n + μ m + λ l (m, n) R | κ m n ω m n, 1 ω m n) P G (ω m n | b = 1, z = 0) d ω m n (8) where PG (ωmn|b = 1, z = 0) is the Pólya-Gamma distribution with parameters b = 1 and z = 0, and κ ≡ y - 1 2. With the introduction of the Pólya-Gamma variables, the joint posterior Eq (4) becomes p (Θ | y) ∝ p (Θ) ∫ X [ ∏ m = 1 M ∏ n ′ = 1 N p (X m n ′ | X m, n ′ − 1, Θ) p (X m 0 | Θ) × ∏ n = 0 N ∫ 0 ∞ N (x m n + f n + μ m + λ l (m, n) R | κ m n ω m n, 1 ω m n) P G (ω m n | b = 1, z = 0) d ω m n ] d X, (9) from which we derive the conditional posteriors for the parameters and the predictive distributions for the augmented variables. We refer to the set of augmented variables as V+ = {X, ω}, and denote by Θ\x to be the set of all parameters except parameter x. Gibbs samplding was implemented in the Python programming language using the Numpy, SciPy, Matplotlib, StatsModels, and Patsy toolboxes. Numerical routines were written in C/C++ and Cython. Software for the Gibbs sampler and a Python wrapper for the Pólya-Gamma routine based on the original code by Jesse Windle, available at https: //github. com/jwindle, is available at https: //github. com/AraiKensuke/PP-AR and https: //github. com/AraiKensuke/pyPG, respectively.
The stochastic oscillation for the mth trial wm is generated by first generating a stochastic phase tm and stochastic amplitude Am as t m, n + 1 = t m n + Δ t (1 + ξ m n C ξ) w m n = 1 + A m n C A sin (2 π ν t m n). (44) where ξm and Am are generated by AR (1) processes, and tm0 ∈ [0,1] a uniform random number. These are then used to generate the irregular oscillation Cξ and CA control the size of deviation from uniformity, while the timescale of the AR (1) processes controls the timescale of variation of amplitude and phase in the oscillation. AR (1) s with timescales on the order of the oscillation period T = 1/ν produce fluctuations that causes considerable variability of period and amplitude in w, which we quantify with the oscillator irregularity, the oscillatory coefficient of variation (OCV) of the instantaneous periods T of oscillation, defined as the ratio of the standard deviation to the mean of time intervals between phase 0 crossings of an oscillatory signal. Fig 1 shows 2 example oscillations generated by this procedure, their OCVs and their power spectral densities. Irregular oscillations have a higher OCV and a less peaked spectral density. To generate spikes, we employ an inhomogeneous renewal process or the LIF model, driven by a stochastic model of oscillation to modulate the firing of spikes. For the inhomogeneous renewal process, for each trial m and time n, if for a uniform random number rmn ∈ [0,1], r m n < Δ t (e μ m + w m n + λ l (m, n) G), we generated a spike, with λ l (m, n) G the user-defined ground truth refractory history function. For the LIF model, spikes were generated using Δ V m n = (- V m n τ + f n + μ m + w m n + χ m n) Δ t, (45) where the Gaussian random variable χ m n ∼ N (0, σ b 2) represents irregular arrival times of afferent spikes, τ = 0. 2 the membrane time constant, μm a baseline DC current and σ b 2 = 210 the amplitude of background fluctuation. Spikes are elicited when Vmn ≥ 1 passes the threshold value of 1, which then causes a reset to Vm, n+1 = 0. LIF have a refractory period that depends on the model parameters, and therefore the CIF naturally is dependent on the last spike time. We analyze neurons from M1 of rat performing a self-paced lever push-hold-pull task [13,41] that have been found to be significantly modulated to the theta rhythm in the LFP. These 2 neurons were recorded on separate electrodes, with “neuron 1” in Figs 10 and 11 recorded and spike sorted from a tetrode on a siliconprobe, and “neuron 2” recorded using juxtacellular (cell-attached) recording, where the spike trains did not need to be sorted. Both neurons were located in layer 1/2, and “neuron 2” identified morphologically as an interneuron, while “neuron 1” is likely to be an interneuron based on spike shape. We identified trials by selecting lever hold periods lasting more than 1 second followed by a large-amplitude pull, and analyzed the 1. 2 second period encompassing the hold-pull period.
We have defined the LOST model, together with accompanying posterior simulation technology, in order to detect the presence of oscillatory firing-rate modulation in a spike train, and infer its phase of oscillation in the presence of a variety of non-stationarities in firing rate that are present in experimental data. Previous methods have assessed the oscillatory content in spike trains by comparing the spiking to a known oscillatory signal like a band-passed LFP [42,43], by detecting oscillation directly from the spike train [12,18,44–46], or by point process regression using the oscillatory signal as a covariate [47]. The LOST model not only detects oscillatory modulation from the spike train itself, but does so in the presence of both spiking noise and oscillatory irregularity, and also allows extensions to inferring additional latent structure, such as non-stationarities in the modulational strength across trials. In addition, LOST is able to separately account for modulational signals of different frequency band, which has proven to be vital in the analysis of real spike trains. Structural priors on the latent state dynamics together with explicit consideration of spiking noise and oscillatory irregularity allows LOST to uncover oscillatory structure even when the firing rate is low, the modulation weak or when the oscillation is irregular, compared to methods that directly regressed the spiking probability on the raw spiking history itself without an intermediary latent state [12,18,19]. LOST uses the intermediary latent state to addresses the irregularity characteristic of neural oscillations. Another approach to modeling the spectral features of time series using Gaussian processes has recently been developed by Wilson et al [48], where different classes of covariance kernels, the SE and SM, respectively, model non-oscillatory and oscillatory structures, analogous to the real and imaginary roots of the characteristic polynomials. It would be interesting to compare the two approaches in future work. The increased sensitivity and flexibility in specifying latent structures, should allow LOST to be used in investigating the role of oscillations in cognitive functions in the cortex. Investigators are increasingly interested in characterizing the change in neural responses to time-varying stimulus or behavior [49,50]. Theta and gamma oscillations in hippocampus change their coupling structure and prevalence during learning and memory acquisition [51,52]. The inclusion of trial-specific structure independently of the LFP in the LOST model may also allow detection of increased recruitment of a given neuron into cell assemblies during periods when oscillations in the LFP are changing. Further, the variability of the timing and presence of oscillations in the LFP seen in many areas of the cortex and hippocampus [2,13,53], suggests oscillatory modulation in single neurons may likewise exhibit finer structure on a per-trial basis. Use of the LOST model in the analysis of such systems may reveal richer dynamics of recruitment into cell assemblies, and a better understanding of the role of single neurons in cognition.
Oscillatory modulation of neural activity in the brain is widely observed under conditions associated with a variety of cognitive tasks and mental states. Within individual neurons, oscillations may be uncovered in the moment-to-moment variation in neural firing rate. This, however, is often challenging because many factors may affect fluctuations in neural firing rate and, in addition, neurons fire irregular sets of action potentials, or spike trains, due to an unknown combination of meaningful signals and extraneous noise. We have devised a statistical Latent Oscillatory Spike Train (LOST) model with accompanying model-fitting technology, that is able to detect subtle oscillations in spike trains by taking into account both spiking noise and temporal variation in the oscillation itself. The method couples two techniques developed for other purposes in the literature on Bayesian analysis. Using data simulated from theoretical neurons and real data recorded from cortical motor neurons, we demonstrate the method’s ability to track changes in the modulatory structure of the oscillation across experimental trials.
action potentials
medicine and health sciences
engineering and technology
signal processing
applied mathematics
membrane potential
random variables
electrophysiology
neuroscience
covariance
motor neurons
simulation and modeling
algorithms
mathematics
computational neuroscience
waves
research and analysis methods
animal cells
kalman filter
probability theory
physics
cellular neuroscience
square waves
cell biology
physiology
neurons
single neuron function
biology and life sciences
cellular types
physical sciences
computational biology
neurophysiology
modulation
Inferring oscillatory modulation in neural spike trains