PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Conf Proc IEEE Eng Med Biol Soc. Author manuscript; available in PMC 2012 August 1.
Published in final edited form as:
PMCID: PMC3293406
NIHMSID: NIHMS353910

Robust Time-Varying Multivariate Coherence Estimation: Application to Electroencephalogram Recordings During General Anesthesia

Abstract

Coherence analysis characterizes frequency-dependent covariance between signals, and is useful for multivariate oscillatory data often encountered in neuroscience. The global coherence provides a summary of coherent behavior in high-dimensional multivariate data by quantifying the concentration of variance in the first mode of an eigenvalue decomposition of the cross-spectral matrix. Practical application of this useful method is sensitive to noise, and can confound coherent activity in disparate neural populations or spatial locations that have a similar frequency structure. In this paper we describe two methodological enhancements to the global coherence procedure that increase robustness of the technique to noise, and that allow characterization of how power within specific coherent modes change through time.

I. INTRODUCTION

Coherence analysis characterizes frequency-dependent covariance between signals, and is useful for multivariate oscillatory data often encountered in neuroscience. In high-dimensional data sets, such as high-density electroencephalogram (EEG), the large number of channels, ranging from 64 to 256, makes pairwise coherence analysis cumbersome in practice. Mitra and Bokil have proposed a ”global coherence” method for summarizing coherence relationships within large multivariate data sets [1]. This procedure consists of an eigenvalue decomposition of the cross-spectral matrix, producing a summary statistic that quantifies the extent to which variance at a particular frequency can be concentrated within the first mode of the decomposition. Recently, this method has been applied to analyze EEG data recorded during general anesthesia [2]. A number of challenges exist in the practical application of this useful procedure. Noise from subject movement or muscle activity, which can introduce large-amplitude transient signals, can pose problems by introducing broad-band artifacts covering a range of frequencies of interest. This can occur frequently in neuroscience experiments, and in clinical settings such as sleep laboratories or the operating room. If temporal averaging is required to estimate the cross-spectral matrix, relatively brief instances of noise can contaminate global coherence estimates from a much longer segment of otherwise artifact-free data. These noisy segments can be removed manually, but doing so can be tedious and subjective, and is unsuitable for real-time clinical monitoring applications. For time-varying signals, such as those recorded during sleep or induction of general anesthesia, coherent oscillations within a frequency band of interest may adopt different spatial distributions, which may be difficult to appreciate from examining the global coherence summary statistic alone.

In this paper we describe two methodological enhancements to the global coherence procedure that increase robustness of the technique to noise, and that allow characterization of how activity within specific coherent spatial modes change through time. We demonstrate the feasibility of this approach by analyzing high-density EEG recorded during induction of general anesthesia.

II. METHODS

A. Global Coherence

The global coherence quantifies the concentration of variance in the first mode of an eigenvalue decomposition of the cross spectral matrix, computed using a singular value decomposition [1]:

equation M1
(1)

where S(k) (f) is the k-th largest entry of the singular value matrix S (f) = diag (S(1) (f), S(2) (f), …, S(N) (f)), which is unitarily similar to the cross spectral power density matrix P (f) and satisfies

equation M2
(2)

for some unitary matrix U (f) at frequency f.

In [2], the cross spectral matrix is computed by averaging cross-spectral estimates over an ensemble of neighboring windows:

equation M3
(3)

where P (f, t) is the outer product of the tapered Fourier transform of a data window with itself, and J represents the number of windows in the interval [t − Δt, t]. A single discrete prolate spheroidal (Slepian) sequence is used for tapering. Estimates are allowed to vary in time by sliding the ensemble of windows in a non-overlapping fashion.

B. Robust Estimation of the Cross Spectral Matrix

Estimates of the cross-spectral matrix as described above are sensitive to noise. To make cross-spectral estimation more robust, we use the median in place of the mean. The median is a robust estimator of centrality, and is much less sensitive to outliers than the mean.

equation M4
(4)

equation M5
(5)

is the set of Slepian tapered cross spectral density matrices. Here the median is computed elementwise and also separately over the real and imaginary components. Estimates are allowed to vary in time by sliding in single-window increments. This overlapping window arrangement allows the cross-spectral and global coherence estimates to track dynamics more rapidly than non-overlapping windows.

C. Projection Onto Global Coherent Spatial Modes

Coherent oscillations within a frequency band of interest may adopt different spatial distributions, which may be difficult to appreciate from the global coherence summary statistic alone. To analyze how power within a specific coherent spatial mode varies as a function of time, we can compute the variance of the projection onto that mode. In particular, we can take a coherent spatial mode estimated from data at a particular moment in time using the eigenvalue decomposition in (2), and apply the projection over the full data record to see how the power within that mode changes in time throughout the course of the experiment.

equation M6
(6)

Normalizing equation M7 by the total power at frequency f and time t, we obtain a spatial mode projection power ratio:

equation M8
(7)

D. EEG Recordings During General Anesthesia

We induced general anesthesia in normal volunteers with the approval of the Human Research Committee at the Massachusetts General Hospital, and followed all hospital safety regulations for administration of general anesthesia. We show data from one subject to illustrate application of the robust estimation methods described above. For the anesthetic induction, we increased the effect-site concentration in step-wise levels of 0, 1, 2, 3, 4, and 5 μg/ml every 14 minutes using a computer controlled infusion [3], [4]. The subject listened to a series of pre-recorded auditory stimuli consisting of their name, affectively-neutral words, or a train of 40 and 84 Hz clicks, presented once every 4 seconds. Subjects were instructed to press a button to differentiate between these sounds. We used the loss of button responses as an indicator of loss of consciousness (LOC). We recorded 64-channel EEG continuously during this time (BrainAmp MRPlus, BrainProducts, GMBH). Prior to EEG recording, we acquired structural MRI for each subject and digitized scalp electrode positions (Polhemus Fastrack 3D). Prior to data analysis, we co-registered the electrode positions with MRI-based scalp surface reconstructions obtained with Freesurfer [5], and re-montaged EEG signals to a Laplacian reference, using distances along the scalp surface to weigh neighboring electrode contributions.

For the average-based global coherence procedure described in II-A, we used ten 2-second data windows to estimate the cross-spectral matrix, tapering with a single Slepian taper with time-frequency bandwidth product of NW = 2. For the median-based global coherence estimate described in II-B, we again used 10 2-second data windows, but used the first 3 Slepian tapers to construct the sample from which to calculate the median.

III. RESULTS AND DISCUSSIONS

In Figure 1 we show the global coherence calculated in three different ways: using the original method described in [2] (top), using that same method but with sliding overlapping windows in time (middle), and using the robust median estimates of the cross-spectral matrix (bottom). Loss of consciousness (loss of response to auditory stimuli) occurs at approximately 65 minutes. We observe that the original global coherence calculation is highly sensitive to noise, taking the form of broad-band artifacts (red vertical lines, Figure 1, top). Application of a sliding window appears to improve temporal resolution (Figure 1, middle), but noise contamination remains. Application of the robust cross-spectral median estimate removes the majority of the noise artifacts (Figure 1, bottom).

Fig. 1
Global coherence estimates using time-averaged non-overlapping windows (top), time-averaged overlapping windows (middle), and robust median estimates of the cross-spectral matrix (bottom) during induction of general anesthesia. The y-axis frequency range ...

From this analysis we observe a high degree of global coherence within the alpha band (8–12 Hz). When consciousness is lost, the global coherence appears to shift upwards in frequency by a small amount (~1 Hz). To analyze the spatial distributions and dynamics of these coherent modes, we apply the spatial projection described in section II-C above, first using the coherent spatial mode estimated from the baseline period (Figure 2, middle), and then using one estimated from the unconscious period (Figure 2, bottom). The coherent spatial mode estimated from the baseline period shows a strong posterior dominance in the alpha band, consistent with the normal alpha rhythm. The spatial mode coming from the unconscious period, however, shows a strong frontal distribution in the alpha band. If we examine the variance within this mode over the full data record, we observe that power in the posterior mode decreases as the experiment proceeds and the level of propofol is increased (Figure 2, top left), whereas the frontal mode has little power during the conscious period, but turns on after loss of consciousness (Figure 2, top right).

Fig. 2
Spatial distribution of the alpha-band (8–12 Hz) coherent modes and power ratio of spatial mode projection. The spatial distribution, an interpolation of the modulus of the 64 complex values in the projection vectors U(1) (equation M9, equation M10) which corresponds ...

In this paper we focused on methods for robust estimation, but did not discuss statistical inference, which will be addressed in future work. Although the cross spectral matrices follow a multivariate Chi-squared distribution (Wishart), the global coherence does not have a closed-form distribution, because the singular value decomposition lacks a fixed functional form. Markov Chain Monte Carlo (MCMC) methods may provide a means of calculating an empirical distribution for the global coherence which can be used to calculate confidence intervals or for hypothesis testing.

In additional to coherence analysis, results with other indicators of the state of consciousness for this experiment can be found in [6] and [7].

IV. CONCLUSIONS

The global coherence is a useful method for summarizing patterns of frequency-dependent co-variation in multivariate time series. In this paper we have developed a robust estimation method that improves the performance of this procedure in noisy data that are frequently encountered in neuroscience experiments, and in clinical settings such as sleep laboratories or the operating room. Furthermore, we have introduced a simple projection procedure that can be used to assess the dynamics of specific oscillatory modes discovered using global coherence analysis. We have illustrated the application of this method on EEG data recorded during induction of general anesthesia, but the method is applicable to a wide variety of neuroscience applications featuring multivariate neurophysiological recordings.

Acknowledgments

This work was supported by NIH Grants DP2-OD006454 (Purdon), K25-NS057580 (Purdon), T32-NS048005 (Harrell), DP1-OD003646 (Brown), R01-EB006385 (Brown, Purdon), and R01-MH071847 (Brown).

References

[1] Mitra P, Bokil H. Observed Brain Dynamics. Oxford University Press; New York: 2007.
[2] Cimenser A, Pierce ET, Salazar-Gomez A, Walsh JL, Harrell PG, Tavares CL, Habeeb K, Purdon PL, Brown EN. Developing new neurophysiological signatures of general anesthesia induced loss of consciousness. BMC Neuroscience. 2009;10:79. [PubMed]
[3] Shafer SL, Gregg KM. Algorithms to rapidly achieve and maintain stable drug concentrations at the sit of drug effect with a computer-controlled infusion pump. Journal of Pharmacokinetics and Pharmacodynamics. 1992;20:147–169. [PubMed]
[4] Schnider TW, Minto CF, Shafer SL, Gambus PL, Goodale DB, Youngs EJ. The influence of age on propofol pharmacodynamics. Anesthesiology. 1999;90:1502–1516. [PubMed]
[5] Dale AM, Fischl B, Sereno MI. Cortical surface-based analysis. I. segmentation and surface reconstruction. Neuroimage. 1999;9:179–194. [PubMed]
[6] Mukamel EA, Wong KFK, Prerau MJ, Brown EN, Purdon PL. Phase-based measures of cross-frequency coupling in brain electrical dynamics under general anesthesia. Conference Proceedings IEEE Engineering in Medicine and Biology Society.2011. [PMC free article] [PubMed]
[7] Wong KFK, Smith AC, Pierce ET, Harrell PG, Walsh JL, Salazar AF, Tavares CL, Cimenser A, Prerau MJ, Mukamel EA, Sampson A, Purdon PL, Brown EN. Bayesian analysis of trinomial data in behavioral experiments and its application to human studies of general anesthesia. Conference Proceedings IEEE Engineering in Medicine and Biology Society.2011. [PMC free article] [PubMed]