A procedure for testing across-condition rhythmic spike-field association change

Abstract

Many experiments in neuroscience have compared the strength of association between neural spike trains and rhythms present in local field potential (LFP) recordings. The measure employed in these comparisons, "spike-field coherence", is a frequency dependent measure of linear association, and is shown to depend on overall neural activity (Lepage et al., 2011). Dependence upon overall neural activity, that is, dependence upon the total number of spikes, renders comparison of spike-field coherence across experimental context difficult. In this paper, an inferential procedure based upon a generalized linear model is shown to be capable of separating the effects of overall neural activity from spike train-LFP oscillatory coupling. This separation provides a means to compare the strength of oscillatory association between spike train-LFP pairs independent of differences in spike counts. Following a review of the generalized linear modelling framework of point process neural activity a specific class of generalized linear models are introduced. This model class, using either a piece-wise constant link function, or an exponential function to relate an LFP rhythm to neural response, is used to develop hypothesis tests capable of detecting changes in spike train-LFP oscillatory coupling. The performance of these tests is validated, both in simulation and on real data. The proposed method of inference provides a principled statistical procedure by which across-context change in spike train-LFP rhythmic association can be directly inferred that explicitly handles between-condition differences in total spike count.

Fig. 1

Testing procedure and associated inference. The flow chart describes the key steps necessary to perform the proposed change-in-modulation testing procedure for a single frequency interval.

Fig. 2

Description of synthetic data used in Fig. (3). Upper left half: Example spike trains and LFP timeseries associated with the left hand side of Fig. (3). The link between the LFP and the spiking rate is constant, while the background rate, α_k, increases from left to right and from top to bottom. Rates are prevented from falling below zero by employing the piecewise linear link function. The preferred phase of spiking in all plots, ϕ_p, is zero. The dark thick lines indicate the spike times, and increase in occurrence frequency with increasing background rate. Note that the fraction of spike times that occur at the LFP peak decreases with increasing background rate. For k > 1, that is for all but the upper left plot, the thick non-positive gray line indicates the occurrence of spikes that have not been removed in the thinning procedure employed in Mitchell et al. (2009) and Gregoriou et al. (2009). Here spike trains are thinned pairwise between the spike trains corresponding to k = 2, 3, 4, and to the spike trains associated with k = 1. As expected, the times at which spikes occur at LFP peaks is greatly reduced in the thinned spike train depicted in the lower right hand plot for the k = 4 case. Upper right half: Example spike trains and LFP timeseries associated with the right hand side of Fig. (3). The link between the LFP and spiking rate decreases from left to right and top to bottom, while the background rate remains constant. Bottom left: One of twenty realizations of the spiking rate associated with the upper left plot in Fig. (2). Bottom right: The theoretical LFP spectrum associated with the synthetic rate plotted on left. The LFP is oscillatory with a spectral peak near 50 Hz. Realizations of this random process are used to generate all of the synthetic LFPs used in simulation.

Fig. 3

Proposed test breaks average rate/association confound. The squared magnitude of the spike field coherence computed between spike trains and LFP time series for synthetic data generated using the model specified by Eq. (29) and illustrated by Fig. 2. When the spike probability modulation due to the LFP is held constant the squared magnitude of the spike field coherence decreases with increasing background rate (left four plots). When the background rate, α_k is held constant while the modulation, in this case equal to β_k, decreases, the squared magnitude of the spike field coherence decreases. Increasing background rate and decreasing association between neural spike times and LFP rhythm are indistinguishable in spike-field coherence. Bar plots: Negative of the base 10 logarithm of the p-values of the proposed test using the piecewise-linear link (PL), see Section 3, for differences in LFP rhythm/spike train association. The test is computed between the data used to compute the upper left plot and the data used to create the other three plots for each of the two cases (constant background rate (left half) and non-constant background rate (right half)). The test disambiguates changes in background rate from changes in the association between LFP rhythm and neural spike times. Ninety-nine percent confidence intervals are depicted in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

Fig. 4

Impact of PL and log link functions on tests of association. Both simulations: The LFP is as used in the first simulation and is depicted in Fig. 2. Left-half (A–D): The modulation, ρ, is changed across conditions while the background α is held constant. The spike-train data are generated using the piecewise linear link (PL), example in plot C. The modulation changes across conditions and both tests are significant (D). Right-half (E–H): The spike-train data are generated using the log link, the preferred phase of spiking changes across conditions, the spike-triggered LFP phase densities are otherwise equivalent (F), and the modulation is much stronger as demonstrated by the response of the rate as a function of LFP phase (E). Only the test computed using the piecewise linear link function is significant (H). This behaviour is consistent with the description in Section 3.1 of the difference between the test computed using the piecewise linear link function and the log link function. The ninety-five percent phase density confidence intervals are computed by boostrapping the trials. They have not been corrected for multiple comparisons. The black horizontal line is the uniform probability density function over phase.

Fig. 5

Modulation, ρ, estimated from the data depicted in Fig. 4. All plots: LFP as described in Fig. 2. Left-half (A,B): Per-frequency rate modulation, ρ, as well as background rate, α, estimated from the spiking data depicted in Fig. 4, plot C. This spiking data is generated with use of the PL link function and is used to produce the plots, (A–D) of Fig. 4. Right-half (C,D): Per-frequency rate modulation, ρ, as well as background rate, α, estimated from the spiking data depicted in Fig. 4, plot G. This spiking data is generated with use of the log link function and is used to produce the plots, (E–H) of Fig. 4. Note that the interpretation of the ρ and α parameters depends upon the type of link function used. When the log link is employed, the modulation is multiplicative; e^ρ, multiplies a background rate of Δ⁻¹e^α Hz to yield the modulation due to LFP rhythm in units of Hz. The data generated with the log link is significantly modulated at a number of frequencies about 50 Hz (C and D) for both conditions, but only modulation ρ associated with the PL statistical model relating LFP rhythm to spiking is significantly different across conditions. The confidence intervals are approximate, 95 percent confidence intervals computed with the delta method (see Appendix D). These intervals are Bonferroni corrected for multiple comparisons. The results are consistent with the results of the proposed hypothesis test plotted in Fig. 4.

Fig. 6

Spike field coherence computed from the data generated in making (Fig. 4). A: Spike field coherence computed from the non-thinned spike-trains associated with the left-half of Fig. 4. B: Spike field coherence computed from the thinned spike-trains associated with the left-half of Fig. 4. C: Spike field coherence computed from the non-thinned spike-trains associated with the right-half of Fig. 4. D: Spike field coherence computed from the thinned spike-trains associated with the right-half of Fig. 4. The thinning procedure slightly reduces spike magnitude-squared coherence between the spikes and the LFP. Magnitude-squared coherence is more affected by the thinning procedure for the data generated with the log link function (C,D). This is due to the larger across-condition difference in the total number of spikes leading to larger probabilities of spike removal in the spike thinning procedure. The frequency resolution of these estimates, equal to twice the time-bandwidth parameter (NW = 1) divided by the observation duration of 200 ms, is 10 Hz.

Fig. 7

Example local field potential (black curves) recorded from Macaque visual cortex during the attend-in (left) and attend-out (right) experimental conditions. Black vertical bars indicate the times at which spikes occur in multi-cell recordings from the FEF. The across-condition test for associativity between LFP rhythm and neural spike train is applied to time-series such as these.

Fig. 8

Application of the proposed test (low-frequency) to real data. Top: The proposed test for spike-field coupling change applied, as a function of frequency, to LFP rhythm and neural spike trains recorded from one Macaca mulatta monkey trained in a covert attention task; as described in Gregoriou et al. (2009). Each frequency has 76 tests, two for each spike-train/LFP pair. Each test was computed using both the PL link function (blue) and the log link function (red). Each test is sensitive to a change across attentional condition but in different ways. The tests have been Bonferroni corrected for multiple comparisons. Bottom: Difference in the modulation, ρ, across attention conditions. At 4 Hz, spike coupling to LFP rhythm is reduced during attention and is accompanied by a spike-phase distribution that is less tuned to a preferred phase of spiking. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

Fig. 9

Application of the proposed test (high frequency) to real data. Top: The proposed test for spike-field coupling change applied, as a function of frequency, to LFP rhythm and neural spike trains recorded from one Macaca mulatta monkey trained in a covert attention task; as described in Gregoriou et al. (2009). Each frequency has 76 tests, two for each spike-train/LFP pair. Each test was computed using both the PL link function (blue) and the log link function (red). Each test is sensitive to a change across attentional condition but in different ways. The tests have been Bonferroni corrected for multiple comparisons. Bottom: Number of pairs with significant differences across conditions vs. frequency. Most changes in association occur for a frequency interval centered upon 50 Hz. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

Fig. 10

Modulation vs. experimental context for the Macaca monkey data as a function of frequency. The modulation is computed using the PL model (left, top), and using the log link model (left-bottom). Associated modulation uncertainty is plotted below, center. Magnitude-squared coherence (MSC) (right-top) and MSC computed from thinned spike-trains (right-bottom), behave in a similar fashion to modulation. In all cases, modulation (PL link or log link) and MSC behave similarly; during the attend-in condition it is larger than during the attend-out condition, and this is true for almost all frequencies. In addition, there is a linearly increasing trend from 0 Hz to 50 Hz followed by a sharp decrease. In all plots, each bar, for a given condition, represents the modulation, MSC, or modulation uncertainty for a single LFP/electrode pair for a specific experimental condition. The ordering of the bars is consistent from one condition to the other so that direct comparisons across condition can be made visually. Due to the large, spiking rates, and the small difference in spiking rates across condition, the MSC computed from the thinned and the non-thinned spike trains is nearly indistinguishable.

Fig. 11

Background rates (α): The height of each bar of the same color indicates the estimated background rate for one LFP/spike-train pair. There are 38 such pairs and the color indicates the experimental condition. The background rate estimates, computed with either the PL link function or the log link function are, as plotted, indistinguishable. Uncertainty is larger for larger rates, a property expected of count-type data. Background rates tend to be larger during the attend-in experimental condition. The uncertainty associated with the log link is e^σ_α ; specifying the multiplicative-modulation of the background rate resulting from one positive standard deviation in alpha when using the log link function.

Fig. 12

Difference of the across-condition average spiking rates plotted against the average of the across-condition average spiking rates. The average spiking rates tend to be large while the across-condition rate changes are approximately an order of magnitude less. The thinning operation has a no effect due to the typical pair-wise (i.e. across-condition) similarity of the spiking rates (Fig. 12); as well as the relatively large firing rates.

Fig. 13

Test comparisons. Summary: all three tests are correlated. The thinning procedure tends to reduce MSC difference across-condition. There exist a number of significant tests associated with the PL link function, that do not have significant counter-parts when using the log link function. This latter fact indicates that for some LFP/spike-train pairs a distributional change in the LFP phase of spiking is not accompanying a change in spike-field association. Left: The difference between the magnitude-squared coherence evaluated at 50 Hz for the attend-in condition is reduced by the magnitude-squared coherence at 50 Hz for the attend-out condition and is plotted against the – log ₁₀(p – values) computed using the proposed testing procedure with the model employing the piecewise linear link function. Circles indicate the use of non-thinned spike-trains and squares indicate that computations have occurred using spike-trains thinned to the minimum average firing rate of the two spike-trains. Each shape corresponds to a single spike-train/LFP pair and the size of each shape is proportional to the difference in the across-condition average spiking rates. Right: Comparison of p-values computed using the proposed testing procedure with differing functions linking the expected intensity to the covariates. A linear trend relating the transformed test p-values exists; with the test computed using the log link tending to have smaller p-values.

Fig. 14

A raster of significant detections. MSC refers to tests performed by comparing 99% magnitude-squared coherence bootstrap confidence intervals across conditions (a detection corresponds to non-overlapping confidence intervals), MSCt refers to the same test but computed using spike-trains randomly thinned to the minimum pairwise average spiking rate; while Log and Piecewise Linear refer to the proposed test computed using the different link functions. The MSC and MSCt detections are detected by both of the proposed tests and the proposed tests are similar with the exception of a few spike-train/LFP pairs. Here a significant p-value is taken to be 0.01 after Bonferroni correction. There are no tests computed using the log link function that are significant where the associated test computed using the PL link function is not significant. Tests computed using the PL link function that are significant where the associated test computed using the log link function is not significant indicate changes in the modulation due to LFP at 50 Hz that do not change the width of the distribution of the LFP phase associated with spikes.

Fig 15

Black curve: empirical probability density functions for ${\widehat{\beta }}_0$ (left column), ${\widehat{\beta }}_c$ (middle column) and ${\widehat{\beta }}_s$ (right column), for background rates increasing from top row to bottom row. Empirical probability density functions are computed from 150 estimates. Each estimate is computed from a single trial consisting of a one second recording generated according to Eqs. (8) and (9), with a background rate specified by the left-hand-side label for each row of the figure, and with a constant LFP-rate coupling constant equal to 100 Hz. The Gaussian approximation is good for all rates and is excellent for rates greater than or equal to 50 Hz. The non-linear link function is in a linear regime for background rates exceeding 100 Hz. As expected, in this regime estimator variance increases with background rate. When the background rate is 0 Hz, the rate is a rectified cosine, and though bias is evident; the estimates are sufficiently accurate to perform inference.