Multivariate crossfrequency coupling via generalized eigendecomposition
 Cited
 1
 Views
 997
 Comments
 0
Abstract
This paper presents a new framework for analyzing crossfrequency coupling in multichannel electrophysiological recordings. The generalized eigendecompositionbased crossfrequency coupling framework (gedCFC) is inspired by sourceseparation algorithms combined with dynamics of mesoscopic neurophysiological processes. It is unaffected by factors that confound traditional CFC methods—such as nonstationarities, nonsinusoidality, and nonuniform phase angle distributions—attractive properties considering that brain activity is neither stationary nor perfectly sinusoidal. The gedCFC framework opens new opportunities for conceptualizing CFC as network interactions with diverse spatial/topographical distributions. Five specific methods within the gedCFC framework are detailed, these are validated in simulated data and applied in several empirical datasets. gedCFC accurately recovers physiologically plausible CFC patterns embedded in noise that causes traditional CFC methods to perform poorly. The paper also demonstrates that spikefield coherence in multichannel local field potential data can be analyzed using the gedCFC framework, which provides significant advantages over traditional spikefield coherence analyses. Nullhypothesis testing is also discussed.
https://doi.org/10.7554/eLife.21792.001Introduction
Crossfrequency coupling (CFC) refers to the phenomenon in which dynamics across two or more frequency bands are related. CFC has been reported in many cortical and subcortical brain regions in multiple species, and has inspired computational models and theories (Canolty and Knight, 2010; Fell and Axmacher, 2011; Lisman and Jensen, 2013). A core idea that permeates CFC theories is that the phase of a slower brain rhythm coordinates a temporal sequence of faster processes that represent specific items in memory or in sensory space. There are several manifestations of CFC (Canolty and Knight, 2010; Hyafil et al., 2015; Jensen et al., 2007; Jirsa and Müller, 2013), and the most widely used CFC analysis methods involve the application of Euler’s formula, phase synchronization, or distribution analyses, with the goal of determining whether highfrequency power values are nonuniformly distributed over lowfrequency phase (Tort et al., 2010).
Nevertheless, the standard corpus of CFC measures is increasingly criticized, with valid concerns about spurious or biased CFC estimates resulting from nonstationarities, sharp transients, nonuniform phase distributions, and other issues (Aru et al., 2015; Kramer et al., 2008; LozanoSoldevilla et al., 2016; van Driel et al., 2015). These potential biases are not fatal flaws, and in some cases can be alleviated by using alternative methods, shufflingbased permutation testing, or debiasing terms (Cohen, 2014; PittmanPolletta et al., 2014; van Driel et al., 2015; Voytek et al., 2013). However, these concerns limit the ubiquitous applicability of standard CFC measures.
The source of these potential biases comes from violating the assumption that brain oscillations are sinusoidal with a stable frequency. This assumption stems from the use of sinusoidal narrowband filters, such as Morlet wavelets or narrow FIR filters. On the other hand, it is increasingly becoming clear that neural oscillations are not always sinusoidal (Jones, 2016; Mazaheri and Jensen, 2008). The nonsinusoidal waveform shape of neural oscillations may provide important insights into underlying biophysical processes, and yet, for CFC analyses, the nonsinusoidal waveform shapes are potential confounds that must be avoided or corrected. Therefore, analysis methods that do not assume sinusoidality are important 1) for methodological reasons to examine CFC in a wider range of datasets and 2) for theoretical reasons to facilitate discovery of CFC regardless of waveform shape (or explicitly incorporating waveform shape), which is an important step towards uncovering the neurobiological origins and functional significances of neural oscillations.
Multichannel recordings provide new opportunities for the discovery of physiologically interpretable patterns of CFC. For example, theories of CFC suggest that the lowfrequency rhythm regulates the timing of different networks (Lisman and Jensen, 2013), implying that the neural circuits that produce the low and highfrequency rhythms should have different electrical projections onto multielectrode recordings. Furthermore, because volume conduction carries activity from a single neuroelectric source to multiple electrodes, incorporating weighted averages from multiple electrodes can boost signaltonoise ratio and help in isolating patterns of activity (Parra et al., 2005), including in CFC analyses (Canolty et al., 2012; Soto et al., 2016). This is particularly relevant for noninvasive recordings, in which highfrequency activity can have low power and can be contaminated by muscle artifacts. Finally, multichannel recordings allow the application of advanced matrix analysis methods, including dimensionalityreduction and source separation algorithms, which can identify spatiotemporal patterns of activity that might be difficult to uncover when considering activity from a single electrode at a time (Grootswagers et al., 2016; Orekhova et al., 2011).
The purpose of this paper is to present a new hypothesisdriven framework for conceptualizing and quantifying CFC in multichannel datasets. This framework is based on using generalized eigendecomposition (GED) of multichannel covariance matrices, and is therefore termed generalized eigendecompositionbased crossfrequency coupling (gedCFC). The backbone of gedCFC for source separation and dimensionality reduction is grounded in decades of mathematics and problemsolving in engineering, with applications in neuroscience (Parra et al., 2005; Särelä and Valpola, 2005). This paper will show that gedCFC has several advantages over the commonly used CFC methods, including increasing signaltonoise characteristics and avoiding confounds of nonsinusoidal oscillations. Perhaps the most important advantages, however, are its potential for use in hypothesisdriven network discovery and its increased flexibility and robustness for uncovering CFC patterns in nonstationary multichannel data.
Five specific methods are derived from this framework (a brief overview is provided in Table 1 for reference), with each method designed for different assumptions about underlying neural processes and different analysis goals. Each method is validated in simulated data, and proofofprinciple applications are shown in empirical data recorded in humans (EEG, MEG, ECoG) and in rodents (LFP and spikeLFP). Furthermore, it is also demonstrated that spikefield coherence in multichannel datasets can be conceptualized as a special case of CFC, permitting the application of gedCFC in this context. Advantages, assumptions, and limitations underlying the different methods are discussed, along with suggestions for practical aspects of data analysis and inferential statistics.
Generalized eigendecomposition
This paper describes the gedCFC framework in a manner that is approachable to scientists with diverse backgrounds; deeper mathematical discussions that justify using GED to optimize linear spatiotemporal filters according to userspecified objectives are presented in many other publications (Blankertz et al., 2008; de Cheveigné and Parra, 2014; Nikulin et al., 2011; Parra et al., 2005; Tomé, 2006).
The goal of gedCFC is to create a component, formed from a weighted sum of all electrodes, that optimizes the ratio between userspecified minimization and maximization criteria. For example, a component might maximize the difference between highfrequency activity that appears during peaks of a lowfrequency oscillation, versus highfrequency activity that is unrelated to the lowfrequency oscillation. This component is created through a matrix decomposition procedure called eigendecomposition.
Eigendecomposition (also called eigenvalue or eigenvector decomposition) involves finding certain vectors that are associated with square matrices. The basic eigenvalue equation is Sw = wλ, where S is a square matrix, w is a vector, and λ is a single number called an eigenvalue. This equation means that multiplying the eigenvector w by matrix S has the same effect as multiplying w by a single number λ. In other words, matrix S does not change the direction of w, S merely stretches or shrinks it. If S is a channelbychannel covariance matrix formed by multiplying the MbyN (channels by time points) data matrix by its transpose, then the eigenvector w points in the direction of maximal covariance. The full set of eigenvectors produces matrix W (with corresponding eigenvalues in the diagonal matrix Λ) that spans the space of S using basis vectors that each account for maximal variance in S while selected to be pairwise orthogonal. This is known as a principal components analysis.
The eigenvalue equation is generalized to two square matrices R and S, as SW = RWΛ. This equation can be more intuitively understood by moving R to the left side via multiplication of its inverse: (R^{−1}S)W = WΛ. R^{−1}S is the matrix analog of S divided by R (think of ⅔ as being equivalent to 3^{−1}× 2). If S is a covariance matrix of a ‘signal’ dataset and R is a covariance matrix of a ‘reference’ dataset, then GED can be understood to produce eigenvectors that identify directions of maximal power ratio (highest gain) in the matrix product R^{−1}S, in other words, directions that best differentiate matrices S from R. The columns of W are called spatial filters (spatiotemporal filters for delayembedded matrices) or unmixing coefficients, and the column with the largest associated eigenvalue is the one that maximally differentiates the two matrices. The filters can be applied to the data to compute components, and they can be visualized by inspecting the columns of W^{−T}, which is also called the forward model of the filter or the activation pattern (Haufe et al., 2014).
gedCFC identifies multichannel CFCrelated networks by contrasting covariance matrices computed from tobemaximized data features (matrix S) against tobeminimized data features (matrix R). The two covariance matrices should be similar enough to suppress CFCunrelated activity, while being different enough to isolate the neural networks that exhibit CFC. For example, one could identify a thetaband (~6 Hz) rhythm and then compute matrix S from peripeak data and matrix R from peritrough data; the GED of S and R would reveal a component that maximally differentiates peakrelated from troughrelated activity, which is a manifestation of networklevel phaseamplitude coupling.
The spatial filter comprises values for each electrode that are used to compute a weighted sum of activity from all electrodes. The resulting time series is the gedCFC component. Note that although defining the spatial filter involves two covariance matrices that might be computed from discrete time windows, the gedCFC component is a "continuous" signal to which timedomain, frequencydomain, or timefrequency domain analyses can be applied.
In addition to being able to specify both maximization and minimization criteria, another important distinction between GED of two covariance matrices and eigendecomposition of one matrix is that eigenvectors are selected to form an orthogonal set only when the matrix S is symmetric. Although both S and R are symmetric covariance matrices, R^{−1}S is generally not symmetric. In other words, the eigenvectors are independent but not orthogonal. This is an important advantage over principal components analysis, which often performs poorly when used as a source separation method (Delorme et al., 2012) (see also Figure 1—figure supplement 1).
Using GED to identify the lowfrequency rhythm
The GED framework is not limited to CFC; it is also recommended that GED is used to define a spatial filter that maximizes power at the lower frequency (Blankertz et al., 2008; Nikulin et al., 2011; Parra and Sajda, 2003; Särelä and Valpola, 2005). There are several ways to define the S and R covariance matrices (see Cohen, 2016, for comparisons); the approach taken here is to compute S from narrowband filtered data and R from the broadband data (de Cheveigné and Arzounian, 2015). The eigenvector with the largest eigenvalue is taken as a spatial filter optimized for activity in that frequency band. The advantages of GED for source separation of narrowband activity include increased signaltonoise ratio, more accurate recovery of source time course dynamics, and eliminated necessity of electrode selection.
The MATLAB code used to simulate data and to implement the five methods presented here are available online (mikexcohen.com/data). Readers are encouraged to reproduce the findings, modify the code to evaluate success and failure scenarios, and adapt and extend the code to specific data analysis applications.
Results
Method 1: gedCFC on troughlocked covariance vs. entire time series covariance
Method 1 is designed for a situation in which the activity of one network fluctuates as a function of the phase of a lowerfrequency network. To implement Method 1, a component (or single electrode) from which to compute the lowerfrequency rhythm is first defined and then the troughs are identified (or peaks, or rising slopes, or any theoretically relevant phase position). The S covariance matrix is computed from data surrounding each trough, and the R covariance matrix is computed from the entire time series. Figure 1—figure supplement 2 provides a graphical overview of Method 1.
Simulated EEG data to illustrate this method were constructed by generating activity in three dipoles in occipital cortex: one to provide a theta rhythm, and two to provide gamma oscillations. The topographical projections of these dipoles were overlapping but not identical (Figure 1a). The theta rhythm (peak frequency of 6 Hz) contained amplitude and frequency nonstationarities, one gamma dipole had a peak frequency at 40 Hz and an amplitude modulated by theta phase, and the other gamma dipole had a peak frequency at 50 Hz and an amplitude that was modulated independently of theta phase. This second dipole served as a ‘distractor’ to test whether it would produce an artifact. These dipole time series, along with random correlated 1/f noise at 2,001 additional brain dipoles, were projected to 64 scalp EEG channels (see Materials and methods for additional details about simulations).
GED was used to extract a component that maximized power in the theta band. The forward model of the filter closely matched the dipole projection (Figure 1a), indicating accurate reconstruction of the source activity. The S covariance matrix was formed using broadband data from ¼ of a thetafrequency cycle surrounding each theta trough (⅛ of a cycle before and ⅛ after the trough), and the R covariance matrix was formed using the entire broadband time series. The MATLAB code used to obtain the weights is [V,D]=eig(covTrough, covTot). The eigenvector in matrix V with the largest eigenvalue (diagonals of matrix D) is the spatial filter that can be used to obtain the ‘trough component'. This component is a weighted sum of all electrodes that maximally differentiates activity during theta troughs from activity during all theta phases.
Method 1 accurately reconstructed both the topographical map and the relationship between gamma power and theta phase (Figure 1). It is noteworthy that the channel data were not bandpass filtered; the 40 Hz gamma component was identified because of the separability of the covariance matrices. It is also noteworthy that the 50 Hz gamma component had twice the power of the 40 Hz gamma component (Figure 1b) and yet was fully suppressed by the spatial filter. This occurred because its dynamics were uncorrelated with theta phase, meaning that activity from the 50 Hz dipole contributed equally to both S and R covariance matrices.
CFC strength can be quantified in a number of ways. One approach is to apply narrowband temporal filters to the thetatrough component and to compute the average peak minus trough amplitude differences. This procedure revealed a maximum at 40 Hz and smaller peaks at 40 ± 6 Hz, corresponding to the nonstationarities induced from the theta rhythmicity of the 40 Hz amplitude modulations (Figure 1b). Other quantitative methods, including computing timefrequency power spectra surrounding troughs, R^{2} fits between the lowfrequency time series and the CFC component power time series, and phase synchronization between the two components, will be illustrated in later sections. Time courses of the reconstructed theta component and power time series of the thetatrough component revealed strong corhythmicity, as was constructed in the simulation (Figure 1b).
A traditional phaseamplitude coupling analysis based on Euler’s formula from POz (the electrode with maximum power) failed to capture the simulated pattern of CFC due to low signaltonoise characteristics of the simulation. Furthermore, the identical temporal filters applied to electrode POz showed little suggestion of CFC (Figure 1c). Together, these results show that even in simple simulations, gedCFC can accurately identify weak patterns of CFC in noisy data while traditional CFC methods can produce uninspiring results.
To demonstrate that gedCFC is robust to topographical distribution, this simulation was repeated by having several different electrodes (not dipoles projected to electrodes) that formed a network of theta phasecoupled gamma oscillators. The neurophysiological interpretation of this simulation is that a single theta generator regulates the timing of a spatially distributed synchronous gammaoscillating network. The results remained robust, as seen in Figure 1—figure supplement 3. Although this is a useful demonstration, it is not a surprising finding, considering that gedCFC makes no assumptions about topographical characteristics or spatial smoothing.
Nonsinusoidal oscillations produce spurious phaseamplitude coupling, either due to nonuniform distribution of phase angles or to fast derivatives causing highfrequency spikes at certain phases (Kramer et al., 2008; van Driel et al., 2015). Therefore, this simulation was repeated, replacing the theta sine wave with repeated Gaussians or van der Pol oscillators at theta frequency (these two waveforms are known to produce spurious phaseamplitude coupling). gedCFC was completely robust to these nonsinusoidal time series, in that (1) the simulated pattern of CFC was accurately recovered regardless of nonstationarities or a nonsinusoidal shape of the lower frequency rhythm, and (2) when gamma oscillations were not amplitudemodulated by the lower frequency rhythm, no spurious CFC was identified. The latter is illustrated in Figure 1b, and additional simulations are shown in the online MATLAB code.
As a proofofprinciple application, Method 1 was applied to restingstate MEG data taken from the Human Connectome Project (Van Essen et al., 2013). The first step was to identify an alpha component by comparing covariance matrices from data filtered in the alpha range (peak 11 Hz; FWHM 5 Hz; these parameters were selected on the basis of visual inspection of power spectra from posterior sensors). There were several physiologically plausible alpha components, as expected on the basis of previous work (Haegens et al., 2014; van der Meij et al., 2016; Walsh, 1958). The alpha component with the largest eigenvalue was selected for subsequent analyses; it is possible that additional crossfrequency coupling dynamics could emerge when examining other alpha components. A covariance matrix of broadband data was constructed around alpha peaks (only peaks corresponding to >1 standard deviation above the mean alpha peak amplitudes were included; this limits the analysis to periods of high alpha power) and compared against the broadband covariance from the entire resting period. The three gedCFC components with the largest eigenvalues were selected for analysis and visualization. Results show that, in this dataset, bursts of broadband and highfrequency activity were timelocked to alpha phase, ranging up to 150 Hz for components 2 and 3 (Figure 2), which appears consistent with previous reports (Osipova et al., 2008). Note that the different topographies of the alpha and alphapeak components rule out the possibility that the crossfrequency coupling resulted from an artifact of the nonsinusoidal shape of alpha.
The primary assumptions of Method 1 are that there is only one highfrequency network of interest, and that its spatiotemporal characteristics (and therefore its covariance matrix) are consistent across each trough. The lowerfrequency rhythm is assumed to be sufficiently rhythmic to be able to identify peaks and troughs, although frequency stationarity is not required; it would be valid to use, for example, empirical mode decomposition (a timefrequency decomposition method that can estimate nonsinusoidal oscillations) to obtain the lowerfrequency rhythm. The primary limitation of Method 1 is that it is valid only when the singlenetwork assumption is valid.
Method 2: gedCFC on covariances of peripeak vs. peritrough
Method 2 is an extension of Method 1, and is designed for a situation in which a lowfrequency rhythm regulates the timing of two different networks that are activated at distinct lowfrequency phases. This might be the case, for example, if theta phase controls the timing of different populations of neurons that represent different items held in memory (Heusser et al., 2016; Zheng et al., 2016).
The implementation is similar to that of Method 1, except that the two covariance matrices are formed from data around the trough and around the peak (or any other two phase regions). The MATLAB code is [V,D]=eig(covTrough, covPeak). The eigenvector in matrix V with the largest eigenvalue (diagonals of matrix D) is the spatial filter used to obtain the ‘trough component,’ while the eigenvector with the smallest eigenvalue is the spatial filter used to obtain the ‘peak component.’ Figure 3—figure supplement 1 provides a graphical overview of Method 2.
EEG data were simulated as for Method 1, except that activity from one dipole had a peak frequency at 40 Hz and an amplitude modulated by theta phase, while activity from the second dipole had a peak frequency at 45 Hz and an amplitude modulated by the inverse of theta phase. However, different frequencies are not necessary, and the two dipoles could oscillate at the same frequency. The key to their separability is topographical projections that produce differentiable covariance matrices, and different amplitude time courses with respect to the lowerfrequency phase.
gedCFC recovered the crossfrequency dynamics of the network, accurately identifying the topographical distributions, frequency ranges, and time courses (Figure 3). Timefrequency power spectra timelocked to the theta troughs and peaks accurately captured the simulated dynamics. Note that this result emerged despite the absence of visually prominent gamma peaks in the channel power spectra. In other words, consistent spatiotemporal patterns that allow component extraction and CFC identification do not necessarily require visually compelling spectral peaks, as has been suggested by Aru et al. (2015); (see He et al. (2010)for another counterexample). This is because the power spectrum computed over an extended period of time can fail to reveal oscillatory activity that is temporally brief or that contains frequency nonstationarities.
The standard phaseamplitude coupling measure relying on Euler’s formula failed to identify the pattern of crossfrequency coupling. The difficulty in this case was the sourcelevel mixing in combination with weak gamma power relative to the noise spectrum. This simulation therefore also demonstrates how gedCFC can increase signaltonoise characteristics, which is particularly important for noninvasilve measurements such as EEG.
As a proofofprinciple illustration, Method 2 was applied to human EEG data taken from a previously published study (Cohen and van Gaal, 2013), in which we reported that thetaalpha phaseamplitude coupling was modulated by error adaptation. In the present reanalysis, GED was applied to define a spatial filter that maximized thetaband power (peak 4 Hz, FWHM 4 Hz, based on the empirical peak theta frequency from channel FCz), which was then used to identify theta peaks and troughs. Next, the 64channel data were filtered in the alpha band (peak 10.5 Hz, FWHM 4 Hz, based on the empirical alpha peak frequency from channel POz), and covariance matrices were computed using data surrounding ¼ of a theta cycle centered on peaks and troughs of the theta component (data were bandpass filtered here because of a priori hypotheses about thetaalpha coupling). The two components with the largest and smallest eigenvalues were selected as spatial filters. The power of those components was computed as a function of 30 theta phase bins, and the distributions are shown in Figure 4c. The relationship between alpha power and theta phase bin was shuffled 1,000 times to produce a nullhypothesis distribution. The 95% range of this distribution is shown as a gray patch in Figure 4c. Data values outside this range can be considered unlikely to occur by chance. Additional possibilities for statistical evaluation include fitting a sine wave to the alpha power distribution, or performing a test against a null hypothesis of a uniform distribution (e.g., a KolmogorovSmirnov test).
The main assumption underlying Method 2 is that the neural populations activated during peaks and troughs of the lowfrequency rhythm produce field potential fluctuations that have separable projections onto multichannel recordings. Other assumptions are similar to those of Method 1: consistent spatiotemporal characteristics of each network, and a lower frequency dynamic that is rhythmic enough to be able to define peaks and troughs.
The main limitation of Method 2 is that the assumption of topographically separable networks according to different lowfrequency phases is critical. If Method 2 were applied to a single network that varies only in amplitude according to lowerfrequency phase, it would produce two identical (plus noise) covariance matrices. gedCFC would then return uninterpretable results (Method 1 would be appropriate in this case).
Method 3: lowfrequency waveform shape as a bias filter on sphered data
Methods 1 and 2 are robust to lowfrequency nonsinusoidal oscillations precisely because the waveform shape is ignored, except to identify specific time points for timelocking the covariance matrices. On the other hand, very different waveform shapes could produce peaks at the same times. Given that waveform shape is the result of biophysical processes, one may wish to incorporate the lowerfrequency waveform shape into the analysis. In this case, Methods 1 and 2 are inappropriate. Therefore, the goal of Method 3 is to use the lowerfrequency time series as a continuous regressor without necessitating a priori specification of the relevant phase regions.
In Method 3, a lowfrequency time series is used as a ‘bias filter’ that is applied to the multichannel data in order to determine the linear combination of electrodes that best matches the bias filter. The multichannel data are first sphered (i.e. whitened; after sphering, the covariance matrix is diagonal, meaning that interchannel covariances are zero), which allows the filter to highlight only the spatiotemporal patterns that covary with the bias filter without influence from endogenous oscillatory activity.
The mechanics of Method 3 closely follow the ‘joint decorrelation’ method presented in de Cheveigné and Parra (2014). A brief overview is provided here and depicted in Figure 5—figure supplement 1; interested readers are directed to the 2014 paper for more details. The first step is to extract the time series of a lowfrequency component, as was done in Methods 1 and 2, to be used as a bias filter. The next step is to narrowband filter the multichannel data and extract the amplitude envelope (e.g., from filterHilbert or from Morlet wavelet convolution). This is an important step: the goal is to identify power fluctuations that vary according to the lowerfrequency phase; therefore, it is the highfrequency power time series, not the bandpass filtered signal, that is of primary interest. The next step is to sphere the multichannel power fluctuations, which is implemented by scaling their eigenvalues. If the channelsbytime power time series data X has eigenvectors and eigenvalues V and D, then the sphered data are defined as Y=X^{T}VD^{½}, where ^{T} indicates the matrix transpose and ^{½} indicates the square root of the matrix inverse. Next, the bias filter is expanded into a Toeplitz matrix (B) that leftmultiplies the data as BY. The eigendecomposition of the covariance matrix of BY provides a new set of eigenvectors W, which are used to rotate the original eigenspace of the data as VD^{½}W. The spatial filter with the largest eigenvalue is applied to the nonsphered narrowband multichannel power time series data to obtain the component. This procedure is repeated over a range of higher frequencies, using the same B matrix for all frequencies.
Data for Method 3 were simulated like those for Method 1. Figure 5 shows the fit of the filtered component to the simulated dipole gamma power time series (R^{2}). Fits were equivalently good (up to 5) when using correlations and phase synchronization measures, and were poor (<0.1) when using data from the gammamaximum electrode instead of the component. The theta peaklocked power spectra from the component showed a precise spectraltemporal localization consistent with how the data were simulated.
Empirical data for a proofofprinciple demonstration were taken from recordings in a human epilepsy patient. Eleven channels were simultaneously recorded, ten along the medial temporal lobe axis (including the amygdala and hippocampus, locallyaverage reference) and one from surface channel Cz (on the vertex on the scalp, referenced to mastoids). Lowfrequency midline cortical activity (peak frequency 4 Hz, FWHM 3 Hz) was extracted from Cz and was used as the bias filter in the Toeplitz matrix B. To reduce computational load and increase signaltonoise characteristics, the continuous data were cut into 50 equally sized epochs of 4 s each, and the covariance matrix of BY was computed separately per epoch and then averaged together.
The gedCFC was applied iteratively using power time series from the medial temporal lobe in 70 frequencies, ranging from 10 to 120 Hz. CFC was quantified using R^{2} fit between the component power time series and the Cz lowfrequency rhythm, and as phase synchronization (which ignores the amplitude fluctuations). Figure 6a reveals a peak in this coupling at around 65 Hz. Inspection of timefrequency power spectra timelocked to theta peaks (Figure 6b) indicates that highfrequency power built up prior to the Cz theta peak. This component was more strongly driven by anterior channels (Figure 6c). Statistical evaluation of the couplingbyfrequency analysis was implemented by shifting the lowfrequency time series by a random amount relative to the power time series, and by recomputing the correlation and phase synchronization coupling measures. This procedure was repeated 1,000 times to generate a null hypothesis distribution. The 99% values of these distributions are shown as dotted lines in Figure 6a.
Interactions between the medial temporal lobe (and hippocampus in particular) and prefrontal cortex are widely implicated in memory formation (Preston and Eichenbaum, 2013). A multivariate technique like gedCFC might be useful in providing insights into how different hippocampalprefrontal networks are involved in different aspects of memory (Shin and Jadhav, 2016). Incorporating the waveform shape into the filter might also prove insightful considering that rat hippocampal theta is nonsinusoidal (Belluscio et al., 2012).
The key assumption of Method 3 is that the lowerfrequency waveform shape is important (as opposed to being used simply to identify peaks and troughs). The main limitation is that this method relies crucially on proper specification of the bias filter. The bias filter could be misspecified if its spatial and spectral characteristics are not known or are suboptimally estimated. Thus, Method 3 should be used only when there is a strong a priori motivation to justify using an appropriate bias filter.
Method 4: gedCFC using timedelay embedded matrices
In Methods 1–3, the higherfrequency components are computed using covariance matrices, but their temporal and spectral characteristics are computed using sinusoidbased filters, which impose sinusoidality and temporal symmetry on the results, particularly when using narrowband filters. Therefore, the purpose of Method 4 is to empirically compute the higherfrequency waveform shape, as well as its spatial distribution, directly from the data, without relying on narrowband filters. This allows the identification of CFCrelated activity that is not necessarily sinusoidal or even rhythmic. And if the activity is sinusoidal, the spatiotemporal filter will empirically identify the sinusoidality without the confound of using a sinusoidal filter.
Method 4 is accomplished through timedelay embedding, which means that additional rows are added to the data matrix, which correspond to timeshifted versions of the channel time series. For example, a firstorder timedelay embedded matrix of 64 electrodes would have 128 rows; the first 64 rows are taken from the 64 electrodes and time points 1 to end1, while the next 64 rows are taken from the same 64 electrodes and time points two to end. Importantly, because GED doesn’t ‘know’ that the data come from timedelayed versions of the same electrodes, the weights assigned to timedelayed rows act as temporal weights from which a temporal filter is empirically computed (de Cheveigné, 2010). For example, a filter that finds the temporal derivative would have weights of 1 and −1 for the original and timedelayed rows. In practice, more than one delay embedding is useful. Timedelay embedding was developed in dynamical systems analysis, and has been successfully applied to psychometric and neural data (Brunton et al., 2016; de Cheveigné and Simon, 2007; Lainscsek and Sejnowski, 2015; Tome et al., 2004; von Oertzen and Boker, 2010). After creating this timedelay embedded data matrix, Method 4 proceeds similarly as Method 1: two covariance matrices are computed, one from data surrounding troughs and one from the entire time series, and GED is applied to those two matrices. The eigenvector with the largest eigenvalue is a spatiotemporal filter of size 1 ×MN, where M is the number of channels and N is the number of embeddings. To interpret this filter, its forward model can be reshaped to M×N and visualized as a timebychannels matrix.
The number and spacing of delay embeds are parameters of the timedelayembedded matrix. If one wishes to use the Fourier transform to determine the filter’s frequency characteristics, the delays should be linear (that is, embedding 1, 2, 3, 4, instead of 1, 2, 4, 8), and the number of delay embeds determines the frequency resolution. Here, a 60^{th} order matrix was used, meaning the data matrix contains 3,840 rows (60 embeddings times 64 channels). At 1,024 Hz, this provided a frequency resolution of 17 Hz, which was then zeropadded to a frequency resolution of 2 Hz. There are other methods to determine the embedding dimension (Cao, 1997; Maus and Sprott, 2011); it is beyond the scope of this paper to include an exhaustive discussion of these methods.
Data were simulated by adding a burst of 75 Hz activity to one dipole at each trough of theta from another dipole (Figure 7a), then adding 1/f noise and projecting to 64 EEG electrodes. The gamma power bursts were timelocked to the troughs, but the phase of the gamma burst varied randomly across bursts, thus producing a nonphaselocked (i.e., ‘induced’) response. A highpass temporal filter (lower edge 20 Hz) was applied to the data to prevent the component from simply reflecting the shape of the lowerfrequency time series around the peak.
Figure 7 shows that the theta and gamma components were accurately recovered. This result can be compared with the GED results to the troughtriggered averages (this is an eventrelated potential, where the ‘event’ is the trough), which failed to detect the simulated crossfrequency coupling. This occurred even without adding noise to the data, and is attributable to the nonphaselocked nature of the gamma bursts. Note that, as with previous simulations, the gamma bursts were generally too small and transient to be visually detected in the channel power spectrum, even at the electrode with the maximal thetatroughcomponent projection (Figure 7D).
Method 4 was applied to empirical data taken from the rat hippocampus (data downloaded from crcns.org; Diba and Buzsáki, 2008; Mizuseki et al., 2013; dataset ec013.156). A theta component was obtained by comparing covariance matrices between 8 Hz and the average of 4 Hz and 12 Hz (a broadband reference was not used because it was dominated by theta). As with the simulation, the data were highpass filtered at 20 Hz to prevent the largest components from reflecting theta phase. Peripeak covariance matrices were computed, and compared against the total covariance. The eigenvector with the largest eigenvalue identified a component with a spectral peak in the gamma band that had a similar spatial projection onto all channels (Figure 8). By contrast, the peaklocked average had a complex waveform shape without a clear spectral concentration.
As implemented here, Method 4 assumes a single network that varies according to lowfrequency phase. The method could be adapted to identify two networks, as with Method 2. Parameterizing the matrix is also important, as insufficient delay embeds will reduce the temporal precision and frequency resolution.
The primary limitation is that working with delayembedded matrices requires sufficient computing resources. For example, two minutes of 256 channels at 1 kHz and a 40fold embedding would produce a matrix of size 10,240by120,000. Matrix operations can be slow and can produce numerical inaccuracies. Two potential solutions are to temporally downsample the data or to use a datareduction technique, such as principal components analysis, to reduce the data, to 40 dimensions for example., makng it possible to perform analyses in the 40dimensional subspace.
Method 5: gedCFC for spikefield coherence
Spikefield coherence refers to local field potential (LFP) correlates of singlecell spiking. It is often quantified either through spiketriggered averages, in which the LFP traces are averaged around each spike (similar to an eventrelated potential, where the event is the spike), or through phaseclusteringbased methods, in which the phase angles are extracted from frequencyspecific analytic time series (e.g., from complex wavelet convolution).
Although spikefield coherence is not typically conceptualized as a manifestation of CFC, in the gedCFC framework, spikefield coherence is simply a special case in which multichannel data are timelocked to the action potential instead of to a lowfrequency trough. In this sense, any of the methods presented above (or combinations thereof) can be adapted to work for multichannel spikefield coherence. In this section, an adaptation of Method 4 is demonstrated. Method 4 is highlighted here because it allows unconstrained empirical discovery of asymmetric spikerelated LFP activity; it can be used, for example, if the LFP waveform shape differs before and after the spike.
Simulated data were modeled after a 16channel linear probe that is often used to study laminar dynamics in the cortex. Spikes were generated at random times, and bursts of 70 Hz gamma were centered on each spike, with a ‘rugged’ laminar profile (similar to that in Figure 1—figure supplement 3, which illustrates that gedCFC makes no assumptions about spatial smoothness) (Figure 9a). In addition, a 75 Hz gamma rhythm that was uncorrelated with spike timing was added as a ‘distractor.’ Random 1/f noise was also added.
As expected, Method 5 recovered the simulated spatiotemporal LFP pattern (Figure 9c–d). With few spikes (N = 40), only a small amount of noise will prevent the spiketriggered average from revealing the true effect. Even when no noise was added, the spiketriggered average failed to capture the true effect if the gamma burst was timelocked but not phaselocked to the spike.
Method 5 was applied to an empirical dataset from recordings of rat medial prefrontal cortex and hippocampus (data downloaded from CRCNS.org; Fujisawa et al., 2008, 2015; dataset EE.042). Data were recorded from 64 channels in the medial prefrontal cortex and 32 channels in the hippocampus. Spiking data were taken from a single neuron in the prefrontal cortex. The first several components appeared to capture a spike artifact in the LFP (Figure 10a). Two later components were selected on the basis of visual inspection. Neither component appeared to suffer from spike artifacts, and both revealed rhythmic LFP dynamics surrounding the spike with phase offsets between the prefrontal and hippocampal electrodes. Temporal and spectral plots of the filter forward model, as well as the spiketriggered average component time series (Figure 10b), revealed that these two spikefield networks had different characteristic frequencies and temporal dynamics. Note the temporal asymmetries (before vs. after the spike) and nonstationarities in the component time series; these are readily visible because no filters were applied that would artificially impose sinusoidality or acausality on the data; instead, the filter kernel was empirically estimated on the basis of broadband data.
Components that capture spike artifacts can be identified by visual inspection or by examining the power spectrum of the filter forward model: artifacts have higher energy in a broad highfrequency range relative to lower frequencies, whereas true effects should have higher energy below ~200 Hz compared to above. Nonartifact components can be individually selected on the basis of prespecified criteria such as frequency band, or can be summarized using principal components analysis.
Figure 10 shows one example of how to apply the gedCFC framework to investigate multichannel spikefield coherence. There are several other possibilities. One could identify the spike artifact components and reconstruct the data with the artifact removed (de Cheveigné, 2010) in order to apply traditional spikefield coherence analyses. Another possibility is to adapt Method 2 to identify LFP components corresponding to different temporal patterns of spiking, such as singlets vs. trains of action potentials. An example of this approach is presented in Figure 10—figure supplement 1.
Overfitting and nullhypothesis testing
Components analyses that are guided by minimization or maximization criteria (such as GED) entail a risk of overfitting. Essentially, one is searching through a highdimensional space for a particular feature; even in pure noise, the components will identify some pattern that best fits the criteria. There are at least three appropriate statistical methods for evaluating the results of gedCFC.
One approach is to compare the resulting components against a nullhypothesis distribution. This nullhypothesis distribution can be derived by applying the method but using randomly selecting time points to be ‘troughs’ (the total number of randomized ‘troughs’ should be the same as the number of real troughs). When repeated, for example 1,000 times with different sets of random ‘troughs,’ the spatiotemporal characteristics, including time courses and power spectra, of the true components can be compared against the empirical null hypothesis distribution. This procedure follows the nonparametric permutation framework that is widely used in neuroimaging and electrophysiology (Cohen, 2014; Maris and Oostenveld, 2007). A difficulty with this approach is that the resulting components may capture the lowfrequency component. One could apply a highpass filter or ignore any coupling close to the frequency used to extract the lowfrequency component.
A second approach is to apply gedCFC to data pooled across all experiment conditions, and then to apply statistical comparisons of component time series or power spectra across different conditions. In this case, although the componentfitting procedure itself may include overfitting, differences across conditions are not biased because the filter was defined orthogonal to condition differences. A related approach is to apply the method to a range of frequencies and then to test whether the coupling strength at different frequencies is statistically significant, as illustrated in Figure 6.
A third approach is to compute the gedCFC filter from independent data and then apply the filter to the experimental data. Independent data could be drawn from different trials from the main experiment (as in crossvalidation) or from an independent dataset such as a localizer. In the case of crossvalidation, confidence intervals can be computed.
Discussion
There is spatiotemporal structure in the mesoscale electrical activity of the brain, which reflects the dynamic interactions within and across populations of neurons. Because this spatiotemporally structured activity is larger than that which can be recorded by a single electrode, singleelectrode analyses will miss or underestimate these patterns. Therefore, as the number of simultaneous electrodes used in neural recordings increases (Stevenson and Kording, 2011), so does the need to consider methods for extracting the neural patterns that underlie those recordings. Componentextraction and dimensionalityreduction methods have proven useful in neuroscience, both in terms of managing largescale data sets and in terms of providing insights into neural mechanics that are difficult to obtain only from visualizing topographical maps (Cunningham and Yu, 2014; Onton et al., 2006).
Crossfrequency coupling in neural time series data
CFC is a fascinating observation that has garnered interest amongst empiricists, modelers, theoreticians, and mathematicians. It has also been the subject of considerable debate, as there seems to be a growing number of publications criticizing CFC methods (Aru et al., 2015; Hyafil, 2015; LozanoSoldevilla et al., 2016). Many of these methodological concerns affect only one framework for analyzing CFC using individual electrodes via narrow bandpass filtering and applying Euler’s formula, phase synchronization, or other distribution analyses. These methods provide accurate results under the assumption of uniform phase angle distribution; CFC biases result from violations of this assumption (van Driel et al., 2015). The framework introduced here evades several of these problems. This new framework does not invalidate or antiquate existing CFC measures; traditional CFC methods are wellsuited for situations in which the lower and higher frequency dynamics are sinusoidal and are produced by neural circuits that have the same spatial projections to multichannel recordings.
The advantages of gedCFC go beyond evading potential biases in traditional CFC measures. This multivariate framework allows detection of diverse spatiotemporal manifestations of CFC, some of which may be difficult to detect with other methods. The inspiration for this framework came from considering the physics of M/EEG and LFP field propagation, by theoretical predictions that the lowfrequency rhythm should regulate anatomically diverse neural populations, and by developments in spatial filtering methods used in source separation and dimensionality reduction analyses.
Advantages of gedCFC
Several specific advantages were highlighted in the Results subsections; a few general remarks are made here. One advantage of the gedCFC framework is that there are no assumptions about underlying generative models that may have produced the data (other than the plausible assumption that electromagnetic fields propagate simultaneously and linearly to multiple channels). It is not necessary to assume Gaussian or Poisson processes for example, nor is it necessary to employ complicated statistical or biophysical models. Instead, the components are extracted directly from the empirical data.
A second advantage is flexibility. From the single equation SW=RWΛ, five applications were derived. By carefully selecting parameters, frequency ranges, experiment conditions, time windows, and so on, this framework is easily extended and tailored to specific hypotheses. One should be cognizant that componentsbased and dimensionalityreductionbased analyses will perform well when there are clean data and when the researcher has clear and physiologically motivated objectives (van Ede and Maris, 2016). Excessively noisy data or poorly specified analysis goals (i.e., from poorly specific hypotheses) may produce misleading or uninterpretable results.
A third advantage is increased signaltonoise characteristics (and therefore, increased sensitivity to detecting CFC), which come from analyzing a weighted combination of data from all electrodes instead of a single electrode. The key is to create those weights appropriately and in an hypothesisdriven manner. The presence of a visually robust peak in the FFTderived power spectrum is not a hard constraint on the success of gedCFC (as illustrated in Figures 3 and 5). An oscillation may exist in a data set and yet have a small peak in a static power spectrum due to nonstationarities in power and phase characteristics.
A fourth advantage is computation time. Many traditional CFC analyses can be prohibitively timeconsuming. gedCFC is fast because software programs like MATLAB use efficient libraries for estimating eigenvalues. On most computers and for most covariance matrices, the MATLAB eig function takes a few milliseconds. The exception is Method 4, which can be processorintensive for very large matrices (the Results section provided a few suggestions for improving computation time, including temporal downsampling and dimensionality reduction). Bandpass filtering, when used, is also fast, and the narrowband filter used here involves only one FFT and one IFFT.
Limitations of gedCFC
Specific limitations were discussed in each Results subsection; a few additional remarks are made here. First, signflipping of components produces uncertainty of peak vs. trough in the lowfrequency component. The signing strategy used here is discussed in the Materials and methods section. Relatedly, bursts of highfrequency activity might be locked to the rising slope or falling slope of the lowfrequency rhythm instead of the peak or trough (Fujisawa and Buzsáki, 2011). Careful data examination and a priori theoretical guidance should be used to determine the appropriate phase values for timelocking.
A second limitation is that the component with the largest eigenvalue is not necessarily the only theoretically relevant component (as illustrated in Figure 10, the largest eigenvalued component may contain artifacts). It is possible that the first few components define a subspace of physiological CFC dynamics. In other words, a single component does not necessarily correspond to a single brain dynamic (de Cheveigné and Parra, 2014). It is advisable to inspect several components for physiological interpretability and modulation by experimental conditions, as was done in Figures 2 and 10.
Third, the gedCFC framework is designed for two covariance matrices. To identify more than two networks, one could iteratively apply, for example, Method 1 to six different phase bins. The six spatial filters could be applied the data to reconstruct six different component time courses. This approach is valid if there really are six distinct networks; otherwise, the same spatial filter might be recreated multiple times, with apparent differences attributable to overfitting noise.
Finally, gedCFC should not be treated as a ‘blackbox’ analysis procedure in which results are interpreted without careful inspection of the data, analysis procedures, and parameters. Misleading or uninterpretable results can occur if excessively noisy data or inappropriate parameters are applied. Instead, gedCFC is best conceptualized as an analysis strategy that mixes physiologically inspired hypothesis testing (e.g., when defining the lowerfrequency dynamic) with blind network discovery. The discovery aspect stems from not needing to specify which topographical regions or frequency bands will manifest CFC.
Nonstationarities and condition differences
The intrinsic nonstationarities of the brain can be problematic for analyses that rely on stationarity over long periods of time, and are also problematic for traditional CFC analyses. Sinusoidal stationarity is not assumed in the gedCFC framework. gedCFC is also not affected by phase lags across different electrodes, as long as the lags are consistent. The core assumption is that the neural patterns maintain a consistent spatiotemporal signature, for example, that each troughlocked covariance matrix is a representative sample. If narrowband filtering is applied to the higherfrequency activity, then an additional assumption is that the activity can be reasonably approximated by sinusoidal basis functions (this is generally a reasonable assumption at the scale of hundreds of ms).
Taskrelated experiments often have multiple conditions. This leads to the decision to compute covariance matrices separately for each condition, or from conditionpooled data to which the spatial filter is then applied to each condition separately. If the covariance structure is expected to be qualitatively different across conditions, separate covariance matrices should be preferred. On the other hand, if the network structure is expected to be the same and only the strength of the modulations are expected to vary, it is preferable to pool data across all conditions when computing the GED.
The quality of gedCFC components is related to the quality of the covariance matrices, and a sufficient amount of clean data ensures highquality covariance matrices. Determining the quality of a covariance matrix can be difficult. One metric is the condition number, which is the ratio between the largest and smallest singular values of the matrix. However, there is no hard threshold for considering a matrix to be illconditioned, and GED can be successfully applied to singular matrices (which have a condition number of infinity). When working with a small amount of data, using alternative methods of estimating covariance, such as shrinkage estimators, may be helpful (Daniels and Kass, 2001).
gedCFC for spikefield coherence provides new insights and avoids potential artifacts
One of the main confounds in spikefield coherence analyses is that the brief largeamplitude spike can cause an artifact in the lowpass filtered LFP (Ray, 2015; Zanos et al., 2011). This is a wellknown phenomenon in signal processing that is often called an edge artifact.
Spikefield coherence here was investigated by adapting Method 4 (Figures 9 and 10) to give rise to Method 5 (Figure 10—figure supplement 1), but other situations might call for different approaches. For example, Method 2 could be further adapted to identify multivariate components that differentiate spikefield coherence patterns between two cognitive or behavioral states (e.g., attend vs. ignore, or rest vs. run). A primary advantage of using gedCFC for spikefield coherence is that traditional spikefield coherence analyses can detect only phaselocked activity; gedCFC will additionally detect nonphaselocked responses, such as a spikelocked burst of nonphaselocked gamma.
Further advances
GED is a powerful technique, but it is not the only sourceseparation method. Many other componentsbased and dimensionalityreductionbased methods exist, including iterative and nonlinear algorithms (Jutten and Karhunen, 2004). For example, phaseamplitude coupling can be identified using parallel factor analysis or tensor decomposition (van der Meij et al., 2012). A primary difference between gedCFC and such techniques is that gedCFC is a hypothesisdriven approach, whereas other techniques are blind decomposition methods that are wellsuited for exploratory analyses. Another important advantage of gedCFC is that it is straightforward both conceptually and in implementation, and therefore it will facilitate analyses for researchers with diverse backgrounds and expertise in mathematics and programming. Furthermore, linear methods tend to produce robust and easily interpretable results (Parra et al., 2005). Nonetheless, the goal of this paper was not to argue that gedCFC is the only sourceseparation algorithm applicable to CFC. Instead, the goal was to highlight the significant insights into the neural mechanisms and the implications of CFC that can be gained by expanding the repertoire of data analyses, as well as theoretical conceptualizations, from a singleelectrodesinusoidbased framework to a multivariatecomponentsbased framework.
Materials and methods
Simulated EEG data
Simulated EEG data were created by generating time series in 2,004 dipoles in the brain according to different assumptions of CFC dynamics as described in the Results, and then projecting those activities to virtual scalp electrodes. Random data were generated by taking the inverse Fourier transform of random complex numbers sampled from a uniform distribution. A 1/f shape was imposed by tapering the spectrum by a sigmoidal curve and then concatenating a mirrored version of the tapered spectrum to produce the negative frequencies. This procedure was done separately for each voxel, thereby producing 2,004 uncorrelated dipoles.
Next, crossvoxel correlations were imposed across all dipoles by creating a random dipoletodipole correlation matrix with a maximum correlation of .8, and computing the new data as Y = X^{T}VD^{½}, where V and D are matrices of eigenvectors and eigenvalues of the correlation matrix, ^{½} indicates the square root, and X is the data matrix. Finally, time series data from selected dipoles were replaced with sinusoidlike time series in the theta or gamma bands, as described in the Results section.
Dipole locations were based on a standard MRI brain. The forward model to project three cardinal directions at each dipole location to the scalp EEG channels was created using algorithms developed by openmeeg (Gramfort et al., 2010) and implemented in the Brainstorm toolbox in MATLAB (Tadel et al., 2011). MATLAB code to generate the simulations, and to apply all five methods, can be found online (mikexcohen.com/data).
Empirical datasets
The following text is copied verbatim from the Human Connectome Project, as requested: MEG data in Figure 2 ‘were provided by the Human Connectome Project, WUMinn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.’ The dataset used here was restingstate data from subject 104012_MEG, session 3.
EEG data in Figure 4 were from a single subject randomly selected from Cohen and van Gaal (Cohen and van Gaal, 2013). The human intracranial EEG data in Figure 6 were taken from 200 s of restingstate data recorded in a patient with epilepsy, who had electrodes implanted as part of presurgical mapping. Data were acquired from the Department of Epileptology at the University Hospital of Bonn, and informed consent was acquired for the recording.
Rodent LFP and singleunit recording data were downloaded from crcns.org. The appropriate references and dataset identifiers are cited in the Results section.
Visualizing gedCFC components
The topographical projections of the components were obtained from columns of the inverse of the transpose of the eigenvectors matrix (Haufe et al., 2014). For the MEG data, projections to the brain were computed using a procedure described in Hild and Nagarajan (2007), and adapted by Cohen and Gulbinaite, 2016. The forward model for this subject was computed using the Brainstorm toolbox. Brain voxels with values lower than the median value across all 15,002 voxels were not colored.
Signflipping of components
One issue that arises in eigenvectorbased components analysis is that the sign of a vector is often not meaningful — the eigenvector points along a dimension, and stretching, compressing, or flipping the sign of the vector does not change the dimension. When extracting power, the sign of the timedomain signal is irrelevant. For topographical projections, the sign also doesn’t matter, although for visual clarity, the sign was adjusted so that the electrode with the largest magnitude was forced to be positive (this is a common procedure in principal components analysis).
However, identifying ‘peaks’ vs. ‘troughs’ clearly requires the correct sign. The solution used here was to correlate the lowfrequency GEDbased component with the bandpass filtered data at the electrode with the maximal topographical projection. If the correlation coefficient had a negative sign, the sign of the lowfrequency component was flipped. This solution is not guaranteed to produce the correct answer, particularly with noisy data. If the interpretation of ‘peak’ vs. ‘trough’ is crucial, it is advisable to run gedCFC separately for peaks and for troughs.
Narrowband temporal filtering
Temporal filtering was done by circular convolution. The Fourier transform of the data was computed using the MATLAB function fft, and was pointwise multiplied by a frequencydomain Gaussian, defined as exp(−0.5[x/s]^{2}), where exp indicates the natural exponential, x is a vector of frequencies shifted such that the desired peak frequency of the filter corresponds to 0, and s is σ(2π−1)/(4π), where σ is the full width at half maximum of the filter, specified in Hz. This filter is advantageous over equivalently narrow FIR or IIR filters because it contains fewer parameters and has no sharp edges in the frequency domain, and therefore does not require careful inspection of the filter kernel’s temporal and spectral response profiles. The inverse Fourier transform was then applied to the tapered frequency domain signal to get back to the time domain. When power time series were required, the Hilbert transform was applied and its magnitude was taken. MATLAB code to implement this filter is provided with the online code.
References

1
Untangling crossfrequency coupling in neuroscienceCurrent Opinion in Neurobiology 31:51–61.https://doi.org/10.1016/j.conb.2014.08.002

2
Crossfrequency phasephase coupling between θ and γ oscillations in the hippocampusJournal of Neuroscience 32:423–435.https://doi.org/10.1523/JNEUROSCI.412211.2012

3
Optimizing spatial filters for Robust EEG SingleTrial analysisIEEE Signal Processing Magazine 25:41–56.https://doi.org/10.1109/MSP.2008.4408441

4
Extracting spatialtemporal coherent patterns in largescale neural recordings using dynamic mode decompositionJournal of Neuroscience Methods 258:1–15.https://doi.org/10.1016/j.jneumeth.2015.10.010

5
Multivariate phaseamplitude crossfrequency coupling in neurophysiological signalsIEEE Transactions on Biomedical Engineering 59:8–11.https://doi.org/10.1109/TBME.2011.2172439

6
The functional role of crossfrequency couplingTrends in Cognitive Sciences 14:506–515.https://doi.org/10.1016/j.tics.2010.09.001

7
Practical method for determining the minimum embedding dimension of a scalar time seriesPhysica D: Nonlinear Phenomena 110:43–50.https://doi.org/10.1016/S01672789(97)001188
 8
 9

10
Analyzing Neural Time Series Data: Theory and PracticeCambridge, Massachusetts: MIT Press.

11
Comparison of linear spatial filters for identifying oscillatory activity in multichannel dataJournal of Neuroscience Methods 278:1–12.https://doi.org/10.1016/j.jneumeth.2016.12.016

12
Dimensionality reduction for largescale neural recordingsNature Neuroscience 17:1500–1509.https://doi.org/10.1038/nn.3776

13
Shrinkage estimators for covariance matricesBiometrics 57:1173–1184.https://doi.org/10.1111/j.0006341X.2001.01173.x

14
Scanning for oscillationsJournal of Neural Engineering 12:066020.https://doi.org/10.1088/17412560/12/6/066020
 15

16
Denoising based on timeshift PCAJournal of Neuroscience Methods 165:297–305.https://doi.org/10.1016/j.jneumeth.2007.06.003

17
Timeshift denoising source separationJournal of Neuroscience Methods 189:113–120.https://doi.org/10.1016/j.jneumeth.2010.03.002
 18

19
Hippocampal network dynamics constrain the time lag between pyramidal cells across modified environmentsJournal of Neuroscience 28:13448–13456.https://doi.org/10.1523/JNEUROSCI.382408.2008

20
The role of phase synchronization in memory processesNature Reviews Neuroscience 12:105–118.https://doi.org/10.1038/nrn2979

21
Behaviordependent shortterm assembly dynamics in the medial prefrontal cortexNature Neuroscience 11:823–833.https://doi.org/10.1038/nn.2134

22
Simultaneous Electrophysiological Recordings of Ensembles of Isolated Neurons in Rat Medial Prefrontal Cortex and Intermediate CA1 Area of the Hippocampus During a Working Memory TaskCollaborative Research in Computational Neuroscience.
 23

24
OpenMEEG: opensource software for quasistatic bioelectromagneticsBioMedical Engineering OnLine 9:45.https://doi.org/10.1186/1475925X945
 25
 26
 27
 28

29
Episodic sequence memory is supported by a thetagamma phase codeNature Neuroscience 19:1374–1380.https://doi.org/10.1038/nn.4374

30
Source localization of EEG/MEG Data by Correlating Columns of ICA Solution with Lead Field Matrix2007 IEEE International Conference On Acoustics Speech and Signal Processing  ICASSP’ 07.https://doi.org/10.1109/icassp.2007.367285

31
Neural CrossFrequency coupling: Connecting architectures, Mechanisms, and FunctionsTrends in Neurosciences 38:725–740.https://doi.org/10.1016/j.tins.2015.09.001

32
Misidentifications of specific forms of crossfrequency coupling: three warningsFrontiers in Neuroscience 9:370.https://doi.org/10.3389/fnins.2015.00370

33
Crossfrequency coupling between neuronal oscillationsTrends in Cognitive Sciences 11:267–269.https://doi.org/10.1016/j.tics.2007.05.003

34
Crossfrequency coupling in real and virtual brain networksFrontiers in Computational Neuroscience 7:78.https://doi.org/10.3389/fncom.2013.00078

35
When brain rhythms aren't 'rhythmic': implication for their mechanisms and meaningCurrent Opinion in Neurobiology 40:72–80.https://doi.org/10.1016/j.conb.2016.06.010

36
Advances in blind source separation (BSS) and independent component analysis (ICA) for nonlinear mixturesInternational Journal of Neural Systems 14:267–292.https://doi.org/10.1142/S012906570400208X

37
Sharp edge artifacts and spurious coupling in EEG frequency comodulation measuresJournal of Neuroscience Methods 170:352–357.https://doi.org/10.1016/j.jneumeth.2008.01.020

38
Delay differential analysis of time seriesNeural Computation 27:594–614.https://doi.org/10.1162/NECO_a_00706
 39

40
Neuronal oscillations with nonsinusoidal morphology produce spurious phasetoamplitude coupling and directionalityFrontiers in Computational Neuroscience 10:87.https://doi.org/10.3389/fncom.2016.00087

41
Nonparametric statistical testing of EEG and MEGdataJournal of Neuroscience Methods 164:177–190.https://doi.org/10.1016/j.jneumeth.2007.03.024

42
Neural network method for determining embedding dimension of a time seriesCommunications in Nonlinear Science and Numerical Simulation 16:3294–3302.https://doi.org/10.1016/j.cnsns.2010.10.030

43
Asymmetric amplitude modulations of brain oscillations generate slow evoked responsesJournal of Neuroscience 28:7781–7787.https://doi.org/10.1523/JNEUROSCI.163108.2008

44
Multiple single unit recordings from different rat hippocampal and entorhinal regions while the animals were performing multiple behavioral tasksCollaborative Research in Computational Neuroscience.
 45

46
Imaging human EEG dynamics using independent component analysisNeuroscience & Biobehavioral Reviews 30:808–822.https://doi.org/10.1016/j.neubiorev.2006.06.007

47
Unraveling superimposed EEG rhythms with multidimensional decompositionJournal of Neuroscience Methods 195:47–60.https://doi.org/10.1016/j.jneumeth.2010.11.010
 48

49
Blind source separation via generalized eigenvalue decompositionJournal of Machine Learning Research 4:1261–1269.
 50

51
Detecting phaseamplitude coupling with high frequency resolution using adaptive decompositionsJournal of Neuroscience Methods 226:15–32.https://doi.org/10.1016/j.jneumeth.2014.01.006

52
Interplay of hippocampus and prefrontal cortex in memoryCurrent Biology 23:R764–R773.https://doi.org/10.1016/j.cub.2013.05.041

53
Challenges in the quantification and interpretation of spikeLFP relationshipsCurrent Opinion in Neurobiology 31:111–118.https://doi.org/10.1016/j.conb.2014.09.004

54
Multiple modes of hippocampalprefrontal interactions in memoryguided behaviorCurrent Opinion in Neurobiology 40:161–169.https://doi.org/10.1016/j.conb.2016.07.015
 55

56
How advances in neural recording affect data analysisNature Neuroscience 14:139–142.https://doi.org/10.1038/nn.2731
 57

58
Brainstorm: a userfriendly application for MEG/EEG analysisComputational Intelligence and Neuroscience 2011:1–13.https://doi.org/10.1155/2011/879716

59
Blind source separation using timedelayed signalsIEEE International Joint Conference on Neural Networks (IJCNN). pp. 2187–2191.

60
The generalized eigendecomposition approach to the blind source separation problemDigital Signal Processing 16:288–302.https://doi.org/10.1016/j.dsp.2005.06.002

61
Measuring phaseamplitude coupling between neuronal oscillations of different frequenciesJournal of Neurophysiology 104:1195–1210.https://doi.org/10.1152/jn.00106.2010
 62
 63

64
Phaseclustering bias in phaseamplitude crossfrequency coupling and its removalJournal of Neuroscience Methods 254:60–72.https://doi.org/10.1016/j.jneumeth.2015.07.014

65
Physiological plausibility can increase reproducibility in cognitive neuroscienceTrends in Cognitive Sciences 20:567–569.https://doi.org/10.1016/j.tics.2016.05.006
 66
 67
 68

69
Autonomy of alpha rhythm generators studied by multiple channel crosscorrelationElectroencephalography and Clinical Neurophysiology 10:121–130.https://doi.org/10.1016/00134694(58)901093

70
Removal of spurious correlations between spikes and local field potentialsJournal of Neurophysiology 105:474–486.https://doi.org/10.1152/jn.00642.2010
 71
Decision letter

Lila DavachiReviewing Editor; New York University, United States
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Multivariate crossfrequency coupling via generalized eigendecomposition" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and David Van Essen as the Senior Editor. The following individual involved in review of your submission has agreed to reveal his identity: Christopher J Honey (Reviewer #2).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
The paper nicely addresses an alternative method to measuring crossfrequency coupling, namely using a multivariatecomponentsbased framework. The paper notes a number of challenges for existing analytic approaches and proposes and examines a new analytic framework named generalized eigendecomposition CFC (gedCFC).
The reviewers all agreed that the basics of the paper are novel, important and sound. However, they also all agreed that the paper was in need of an extensive rewriting coupled with some additional analyses to bring this to a level appropriate for eLife readership. Although we typically at eLife will provide a condensed list of essential revisions, in this case, I am appending the actual reviews below because they all nicely describe the needed changes.
Reviewer #1:
This paper reports results of various methods of applying generalized eigendecomposition to perform crossfrequency correlation analysis. The discussion in the Introduction seems to be too cursory to be of much benefit to those not already familiar with the idea, so depending on the purpose and target audience, some elaboration may be warranted.
There seems to be an error in the statement of the equation: SW = LRW (where L stands for Λ) – I think SW = RWL is meant.
I would humbly suggest giving an intuitive illustrative example in the Introduction to demonstrate the "maximization and minimization" effect referred to, perhaps using the Rayleigh coefficient framework.
The Discussion seems to be somewhat lacking in structure. It may be helpful to move some of the discussion in the corresponding Methods sections to the general Discussion section. The identification of different methods is interesting, but it would be nice to have a prior list of "general CFC contexts" in the Introduction, in which one would like to apply the GED framework, perhaps with simple examples of each case, then the methods description referring to the corresponding contexts, focusing on detailed description of the simulation/data used and method employed, followed by a more general discussion in the Discussion section.
I believe the paper could benefit from a quantitative report of the resulting covariances obtained with the covariance matrices, e.g. the theta and gamma template vectors with each of the two covariance matrices used in the method, as well as the generalized eigenvalue. And it would be of use to know the relative power in the high and low frequency bands used in the simulations, with some corresponding discussion of minimal relative power necessary e.g. gamma, to obtain interpretable results, and the applicability of the required power levels assumptions to real data, e.g. EEG.
The paper is likely to be useful for specialists working in the area, but if a wider audience is intended, a more generous introductory section and structured Discussion section would be helpful.
Reviewer #2:
The goal of this paper is "to highlight that significant insights into the neural mechanisms and implications of CFC can be gained by expanding the repertoire of data analyses, as well as theoretical conceptualizations, from a singleelectrodesinusoidbased framework to a multivariatecomponentsbased framework." The author notes a number of challenges for existing CFC approaches include nonstationarities and sharp transients in signal, nonuniform phase distributions, and the fact that CFC sources and targets are spatially distributed and potentially overlapping. To address these challenges they propose a new analytic framework named generalized eigendecomposition CFC (gedCFC). Although the method proposed is very promising, the manuscript reads more like a proof of principle than a convincing demonstration. This is not because the framework is flawed, but because the range of applications discussed (5 different settings) leads to a lack of focus, and so the strengths and weaknesses of the framework are not clearly enough tested and reported within each setting.
1) The paper describes five methods for characterizing changes in spatiotemporal properties of brain data, using a common linear algebraic formalism. On the one hand, the range of methods is valuable because CFC measures are applied in a range of contexts to address a range of hypotheses. However, this also means that manuscript lacks conceptual unity and focus. For example, it is not entirely clear which of the 5 methods are solutions to which of the problems (e.g. nonstationarities, sharp transients, nonuniform phase distributions, etc.) with existing CFC methods. I would recommend:i) try to group the 5 methods into clusters with a common purpose (for example Methods 1 and 2 could be grouped), rather than simply using numeric labels;ii) assign names to the different approaches;iii) create a table which summarizes, for each method, the kind of hypothesis it is best suited for testing, which problems it solves, and which problems it does not solve;iv) for each of the claims above concerning the problems solved by gedCFC, please design a test that specifically demonstrates that claim. The manuscript already includes valuable demonstration examples that do perform this function, but it is not always explicitly specified which facet of gedCFC they are meant to test.
More generally, I found that the readability of the manuscript decreased from Method 3 onward, as more complex methods were introduced. In method 4 delaycoordinate embedding is introduced, just in passing, but this is a complicated and rich area of research. The investigation of how to choose the correct parametrization (e.g. number of time delays) in Method 4 could be a paper of its own. So, again, it is nice to see how the gedCFC framework can be applied, but the range of uses also limits the manuscript's depth and clarity.
2) Why compare gedCFC against methods, such as PACz, which are (according to the citations in the introduction) known to be flawed? The manuscript would be more compelling if gedCFC was compared against stateoftheart methods (some of which are mentioned in the Introduction) and/or against the most successful of the methods reviewed in Tort et al. (2010, J Neurophys). In addition, when existing methods fail, it would help to specify precisely why they are failing. For example, it could be that gedCFC has, effectively, greater statistical power when applied to multichannel data due to pooling of information across multiple electrodes, but it was unclear to me whether this was the core reason for its advantage over PACz.
3) "This alleviates the feature of principle components analysis that makes it suboptimal for brain data: components are forced to be orthogonal but neural dynamics (certainly at the meso to macroscopic scale of LFP and EEG) are not orthogonal."
It seems unlikely to me that this undesirable feature of PCA – the fact that, e.g., the second eigenvector's direction is constrained by the direction of the first eigenvector – is truly absent from the gedCFC framework. Generally, there is a distinction drawn between (i) dimensionality reduction methods, such as PCA, which rotate data so that it can be compactly reexpressed; and (ii) source separation methods, such as ICA, which aim to identify statistically distinct generators within a mixed signal. Although in many cases PCA does approximate the function of source separation, this is not its goal, and the orthogonality assumption is a testament to the fact that compact expression is prioritized over source separation.
My (admittedly intuitive) impression is that the gedCFC framework is closer to a dimensionality reduction method than to a source separation method. Specifically, even though the components in gedCFC are not, in general, constrained to be orthogonal, this does not imply that (e.g.) the direction of the eigenvector of secondlargest eigenvalue is unaffected by the direction of the eigenvector whose eigenvalue is largest. Consider the special case in which R^(1)S is indeed symmetric: in this case, the eigenvectors are again constrained to be orthogonal and the wellknown PCA problems return. So, should we expect PCAlike difficulties whenever R^(1)S is nearsymmetric? What is our guarantee that these issues do not arise even when R^(1)S is far from symmetric?
The manuscript already presents some examples of cases where gedCFC successfully extracts signals from a mixture, but in these cases the mixed signals have different frequencies, so the filtering operations may be doing much of the work. To demonstrate that gedCFC really does not run into the interpretational difficulties seen with PCA, I recommend setting up a situation in which multiple signals (say 3) with similar spatial profiles, similar frequencies, and overlapping functional response profiles, are mixed – a case in which PCA should fail – and then show that gedCFC does not suffer from similar problems.
4) Because the model works by comparing a reference covariance to a target covariance, it assumes a binary transition between two modes. But how well would the model perform if the underlying generative process is actually changing continuously? For example if the generative model of gamma power is something like probability(gamma burst) ~ phase(theta), then a method that assumes a continuous variation as a function of phase will have far greater statistical power than a method which simply contrasts peak and offpeak data? Now, Method 3 may present an answer this "binary assumption" issue, because a continuous timevarying signal is convolved and included the model. But in this case, does not one not reencounter the problem of nonstationarity and sharp transients, because now the estimation of the lowfrequency envelope matrix, B, could be affected by these factors? Similar to the point I made in (1) above, I would like to see a more careful discussion of the strengths and weaknesses of each method, so that the manuscript explicitly describes for each method, which of the common challenges for CFC are overcome, and which are not overcome.
5) The paper would be more practically useful if it provided a slightly more thorough treatment of at least one of the statistical methods for testing components, and noted briefly any problems that may arise in this setting.
Reviewer #3:
This paper describes a novel technique to calculate crossfrequency coupling using generalized eigenvalue decomposition. The technique is novel and has several significant advantages that have currently not been clearly addressed, such as its independence of selecting electrodes and not assuming a sinusoidal shape. Moreover, the technique can be used in several ways (5 methods are described in the paper). The paper is well written; as I was reading it and wrote down questions, a subsequent section would often already provide an answer. Overall, I think the paper is of high quality and I only have a few comments that can be easily addressed.
1) It is not quite clear to me whether the method assumes that multiple electrodes have the exact same phase for an underlying rhythm. For example, if the method were to be used to analyze LFP recordings across motor cortex, with traveling β oscillations, where there is a shift in phase across electrodes, how would that reflect upon component waves?
2) It is said that the quality of gedCFC components is related to the quality of the covariance matrices, and a sufficient amount of clean data ensures highquality covariance matrices. Can this be quantified: is there a way to compute from one of the data sets used in the paper how the amount of data used relates to the quality of the gedCFC components?
3) In method 3 example recordings from a human epilepsy patient are used with recordings of the medial temporal lobe and electrode Cz. Why would there be any coupling in these data? How does referencing affect the method? More anterior temporal lobe channels are also closer the face with potential muscle artifacts, which may not even be removed with bipolar referencing.
https://doi.org/10.7554/eLife.21792.025Author response
[…] Reviewer #1:
This paper reports results of various methods of applying generalized eigendecomposition to perform crossfrequency correlation analysis. The discussion in the Introduction seems to be too cursory to be of much benefit to those not already familiar with the idea, so depending on the purpose and target audience, some elaboration may be warranted.
There seems to be an error in the statement of the equation: SW = LRW (where L stands for Λ) – I think SW = RWL is meant.
Indeed, that was a typo (copied from the initial introduction of a single eigenvalue). Thanks for catching that.
I would humbly suggest giving an intuitive illustrative example in the Introduction to demonstrate the "maximization and minimization" effect referred to, perhaps using the Rayleigh coefficient framework.
I agree that an intuitive example would be useful. However, I would guess that anyone who is familiar with the Rayleigh quotient would also understand GED. I’ve opted for a visual example presented in Figure 1—figure supplement 1, (inspired by Figure 6 in Blankertz et al., 2008) that I hope is more intuitive to readers without a strong mathematical background.
Furthermore, I added several new supplemental figures that provide graphical illustrations of the methods.
The Discussion seems to be somewhat lacking in structure. It may be helpful to move some of the discussion in the corresponding Methods sections to the general Discussion section. The identification of different methods is interesting, but it would be nice to have a prior list of "general CFC contexts" in the Introduction, in which one would like to apply the GED framework, perhaps with simple examples of each case, then the methods description referring to the corresponding contexts, focusing on detailed description of the simulation/data used and method employed, followed by a more general discussion in the Discussion section.
I have tried to improve the readability and organization of the manuscript, including more discussions of the assumptions and goals of each specific method. The Discussion section is now more focused, and several methodspecific points have been incorporated into the relevant Results section. Furthermore, the revision includes a table that summarizes the key goals and assumptions of each method.
I believe the paper could benefit from a quantitative report of the resulting covariances obtained with the covariance matrices, e.g. the theta and gamma template vectors with each of the two covariance matrices used in the method, as well as the generalized eigenvalue.
I’m afraid I don’t completely understand this request. The covariance matrices themselves are not terribly insightful and they are dependent on the location of the simulated dipole, which is arbitrary and not relevant for the method. If the reviewer was referring to the use of GED to obtain withinfrequency components (e.g., the theta component to extract theta phase), then this is a validated method in neuroscience, and I have added additional references in the subsection “Using GED to identify the lowfrequency rhythm” for interested readers. If the reviewer is referring to assessing the quality of a covariance matrix, this is now discussed in the last paragraph of the subsection “Nonstationarities and condition differences”. If I have misunderstood this comment, then I apologize and ask the reviewer to clarify.
And it would be of use to know the relative power in the high and low frequency bands used in the simulations, with some corresponding discussion of minimal relative power necessary e.g. gamma, to obtain interpretable results, and the applicability of the required power levels assumptions to real data, e.g. EEG.
The relative power between high and low frequency bands is not important; indeed, they are not directly compared. Instead, it is only the covariance matrices within a frequency band from different time periods that are compared (or, in other cases, relative to the broadband power). For example, in the simulations shown in Figure 1, the gamma peak is present but weaker than the “distractor” (noncoupled) gamma; in the simulations shown in Figures 3 and 5, the gamma peaks are hardly visible in the power spectrum. This is now discussed in the Results section among other places.
The paper is likely to be useful for specialists working in the area, but if a wider audience is intended, a more generous introductory section and structured Discussion section would be helpful.
I hope the reviewer agrees that the revised version is more approachable. I welcome any additional specific suggestions for improvement.
Reviewer #2:
The goal of this paper is "to highlight that significant insights into the neural mechanisms and implications of CFC can be gained by expanding the repertoire of data analyses, as well as theoretical conceptualizations, from a singleelectrodesinusoidbased framework to a multivariatecomponentsbased framework." The author notes a number of challenges for existing CFC approaches include nonstationarities and sharp transients in signal, nonuniform phase distributions, and the fact that CFC sources and targets are spatially distributed and potentially overlapping. To address these challenges they propose a new analytic framework named generalized eigendecomposition CFC (gedCFC). Although the method proposed is very promising, the manuscript reads more like a proof of principle than a convincing demonstration. This is not because the framework is flawed, but because the range of applications discussed (5 different settings) leads to a lack of focus, and so the strengths and weaknesses of the framework are not clearly enough tested and reported within each setting.
Thank you for your evaluation and for your many comments. I agree with this assessment, but I would not consider it to be a weakness. The intention was not to propose a single method with a single set of parameters; instead, the goal was to provide the CFC field with an alternative perspective for thinking about and analyzing CFC in multichannel data. The specific implementation of the gedCFC framework should be left to the researcher, depending on the data and goals of the analysis. I agree that a paper presenting only one method would be more focused as a piece of written literature, but it would also be less useful as a scientific contribution.
Furthermore, because the different methods are generally similar to each other, I do not think it is necessary to provide a full complement of testing for each method; this was the reason for setting up the different simulations in different ways. I have tried to make the consistency across methods more apparent in the revised version.
As detailed below, many aspects of the revision have been clarified and tightened. I hope the reviewer agrees.
1) The paper describes five methods for characterizing changes in spatiotemporal properties of brain data, using a common linear algebraic formalism. On the one hand, the range of methods is valuable because CFC measures are applied in a range of contexts to address a range of hypotheses. However, this also means that manuscript lacks conceptual unity and focus. For example, it is not entirely clear which of the 5 methods are solutions to which of the problems (e.g. nonstationarities, sharp transients, nonuniform phase distributions, etc.) with existing CFC methods. I would recommend:i) try to group the 5 methods into clusters with a common purpose (for example Methods 1 and 2 could be grouped), rather than simply using numeric labels;ii) assign names to the different approaches;iii) create a table which summarizes, for each method, the kind of hypothesis it is best suited for testing, which problems it solves, and which problems it does not solve;iv) for each of the claims above concerning the problems solved by gedCFC, please design a test that specifically demonstrates that claim. The manuscript already includes valuable demonstration examples that do perform this function, but it is not always explicitly specified which facet of gedCFC they are meant to test.
Thank you for these suggestions. The summary table is a great idea, and is included in the revised version. I decided against explicit grouping of methods, as they do not clearly group together in a parsimonious way. The text now discusses better how the different Methods are related to each other. I also decided against giving each method a name, as the names would either be (1) as long as their descriptions or (2) difficult to remember.
I appreciate that the reviewer feels the paper “lacks conceptual unity”; however, the goal here is not to present a simple literary work, but instead to present a group of related data analysis techniques, linked by a single mathematical and conceptual framework, that will benefit the electrophysiology community. Readers will be able to determine the most appropriate method for their data and hypotheses, rather than relying on (and, typically, being disappointed by) one particular blackbox algorithm that may or may not be appropriate for their data.
The revision is more explicit about the underlying assumptions and limitations of each method, including extra text at the beginning and end of each of the subsections that introduce each method.
It is not really the case that there is a onetoone mapping between a CFC problem and a gedCFC method. All of the methods address the limitations faced by standard CFC analyses. But the purpose of gedCFC is not simply that it addresses limitations of existing methods; it is a totally new technique for measuring CFC. Therefore, it is not so simple as writing “method A has problem X, and gedCFCmethodY resolves it.” That would be analogous to saying that the primary purpose of EEG is that fMRI has poor temporal resolution.
More generally, I found that the readability of the manuscript decreased from Method 3 onward, as more complex methods were introduced. In method 4 delaycoordinate embedding is introduced, just in passing, but this is a complicated and rich area of research. The investigation of how to choose the correct parametrization (e.g. number of time delays) in Method 4 could be a paper of its own. So, again, it is nice to see how the gedCFC framework can be applied, but the range of uses also limits the manuscript's depth and clarity.
Indeed, parameterization of delayembedded matrices is a nontrivial problem. But without making this paper a fulllength book, it is not possible to go into great detail about each method. I appreciate that this necessarily requires sacrificing some focus, but as written above, the goal here was to present a framework, not a blackbox recipe for a specific analysis.
I have tried to make the paper more readable. However, I believe it simply will necessarily be the case that people will need some background in signal processing and matrix analysis to be able to apply the analyses described here.
More concretely, the revision includes more approachable introductions to each method, and pictorial overviews of the mechanics of each method.
2) Why compare gedCFC against methods, such as PACz, which are (according to the citations in the introduction) known to be flawed? The manuscript would be more compelling if gedCFC was compared against stateoftheart methods (some of which are mentioned in the Introduction) and/or against the most successful of the methods reviewed in Tort et al. (2010, J Neurophys). In addition, when existing methods fail, it would help to specify precisely why they are failing. For example, it could be that gedCFC has, effectively, greater statistical power when applied to multichannel data due to pooling of information across multiple electrodes, but it was unclear to me whether this was the core reason for its advantage over PACz.
PACz is a standard and widely used measure of CFC, probably the most commonly used method. It has its limitations, but I would not consider it to be “flawed.” Instead, it can overestimate or underestimate true PAC in the case of nonsinusoidal oscillations, depending on the nature of the nonstationarities (this is detailed in van Driel et al., 2014, J Neuro Methods). PACz will produce the same results as the distribution analysis described in Tort et al. 2010 when power has a monopolar distribution over lowfrequency phase (this is what is observed empirically) and when the lower frequency phase angles are uniformly distributed.
The reviewer is correct that the increased SNR of gedCFC is partly due to the weighted average of electrodes. However, the primary benefit is the robustness to nonstationarities in the lowfrequency and the highfrequency dynamics, which are the primary sources of bias in PACz, as described in the manuscript. These points are now detailed in the paper.
3) "This alleviates the feature of principle components analysis that makes it suboptimal for brain data: components are forced to be orthogonal but neural dynamics (certainly at the meso to macroscopic scale of LFP and EEG) are not orthogonal."
It seems unlikely to me that this undesirable feature of PCA – the fact that, e.g., the second eigenvector's direction is constrained by the direction of the first eigenvector – is truly absent from the gedCFC framework. Generally, there is a distinction drawn between (i) dimensionality reduction methods, such as PCA, which rotate data so that it can be compactly reexpressed; and (ii) source separation methods, such as ICA, which aim to identify statistically distinct generators within a mixed signal. Although in many cases PCA does approximate the function of source separation, this is not its goal, and the orthogonality assumption is a testament to the fact that compact expression is prioritized over source separation.
My (admittedly intuitive) impression is that the gedCFC framework is closer to a dimensionality reduction method than to a source separation method. Specifically, even though the components in gedCFC are not, in general, constrained to be orthogonal, this does not imply that (e.g.) the direction of the eigenvector of secondlargest eigenvalue is unaffected by the direction of the eigenvector whose eigenvalue is largest. Consider the special case in which R^(1)S is indeed symmetric: in this case, the eigenvectors are again constrained to be orthogonal and the wellknown PCA problems return. So, should we expect PCAlike difficulties whenever R^(1)S is nearsymmetric? What is our guarantee that these issues do not arise even when R^(1)S is far from symmetric?
The manuscript already presents some examples of cases where gedCFC successfully extracts signals from a mixture, but in these cases the mixed signals have different frequencies, so the filtering operations may be doing much of the work. To demonstrate that gedCFC really does not run into the interpretational difficulties seen with PCA, I recommend setting up a situation in which multiple signals (say 3) with similar spatial profiles, similar frequencies, and overlapping functional response profiles, are mixed – a case in which PCA should fail – and then show that gedCFC does not suffer from similar problems.
GED is often referred to as a source separation method, both in the EEG (e.g., Parra 2003) and in the mathematics and engineering literatures (e.g., Tomé 2006). It is also, as the reviewer suggests, referred to as a dimensionality reduction method, although source separation is also dimensionality reduction (the reverse, of course, is not true).
Although the paper provides only a theoretical remark about the superiority of GED over PCA, this has been empirically demonstrated in several papers in neuroscience, a few of which are now cited in the text. Given the many instances in the literature showing the superiority of GED over PCA for source separation, and given that PCA does not fit into the gedCFC framework, and that the purpose here is to provide a method for CFC and PCA is not a method for CFC, it seems redundant and beyond the scope of this paper to provide another formal demonstration of GED vs. PCA for source separation. PCA was mentioned only to introduce readers to GED, because in my experience, most electrophysiologists are familiar with PCA and this provides an easy scaffold for understanding GED. I believe the new Figure 1—figure supplement 1 also helps clarify the improvement of GED over PCA for source separation.
Furthermore, I’m not sure I follow the reviewer’s suggestion: If three signals have the same spatial and spectral profiles and overlapping functional response profiles, then it would seem to me that these are actually reflecting a single functional source. gedCFC is not a dipolelocalization technique; in fact, it works entirely indifferently to spatial configuration (as illustrated in Figure 1—figure supplement 2). The term “source” has varied interpretations in neuroscience; components analyses refer to statistical sources, as opposed to anatomical sources (the two sometimes overlap but often do not).
As for temporal filtering, the reviewer’s point about temporal filtering taking care of the separation is not correct. For one thing, many of the analyses reported here did not involve any temporal filtering to obtain the CFC components. Secondly, no spurious CFC was identified for strong gamma components that were in the data but not locked to the lowfrequency rhythm.
Regarding the symmetry: It is difficult to imagine a situation in which R^{1}S will just happen to be symmetric. Symmetric matrices are very rare, and basically only occur by specific construction (multiplying a matrix by its transpose). Even multiplying a symmetric matrix with the identity matrix that has a single offdiagonal element replaced by .1 will produce a nonsymmetric matrix. Furthermore, the law (it is mathematically required; it is not an assumption) that a symmetric matrix has orthogonal eigenvectors simply states that if the eigenvalues are different, the eigenvectors are orthogonal (if there are repeated eigenvalues, the vectors need not be orthogonal, but most algorithms select them to be orthogonal); there is nothing in the math about vectors approaching orthogonality as a matrix approaches symmetry.
Finally, regarding the point about the direction of the second vector being constrained by the direction of the first vector, this is true, but it is also true for ICA and any other source separation method, in the sense that two eigenvectors cannot point in exactly the same direction (analogously, two ICs cannot be exactly the same). For eigendecomposition, any two distinct eigenvalues will produce different eigenvectors. Therefore, any one eigenvector will constrain the direction of any other eigenvector in the senses that (1) their angles are nonzero for distinct eigenvalues and (2) they are orthogonal for symmetric matrices. Even for reducedrank matrices that lack a full spectrum of eigenvalues, the uncertainty about inplane eigenvectors affects only the smaller components that are not considered in gedCFC.
I’m not sure what to do about this comment. It is beyond the scope of this paper to provide a lengthy discussion about linear algebra concepts (readers can find such information elsewhere). I hope the supplemental figure illustrating GED in a 2D example is informative.
4) Because the model works by comparing a reference covariance to a target covariance, it assumes a binary transition between two modes. But how well would the model perform if the underlying generative process is actually changing continuously? For example if the generative model of gamma power is something like probability(gamma burst) ~ phase(theta), then a method that assumes a continuous variation as a function of phase will have far greater statistical power than a method which simply contrasts peak and offpeak data? Now, Method 3 may present an answer this "binary assumption" issue, because a continuous timevarying signal is convolved and included the model. But in this case, does not one not reencounter the problem of nonstationarity and sharp transients, because now the estimation of the lowfrequency envelope matrix, B, could be affected by these factors? Similar to the point I made in (1) above, I would like to see a more careful discussion of the strengths and weaknesses of each method, so that the manuscript explicitly describes for each method, which of the common challenges for CFC are overcome, and which are not overcome.
We should distinguish between two different stages of analysis: (1) designing the spatial filter; (2) using the spatial filter to extract a component that is then used for analysis. Only the first stage involves a binary discretization; the second stage provides a continuous measure of activity.
In fact, the situation the reviewer describes (gamma power varying continuously as a function of theta phase) is exactly what was implemented in the simulations for Methods 1 and 2. As shown in Figures 1B and 3E, the methods work quite well at recovering continuous gamma power as a function of continuous theta phase, even though the spatial filter was defined only using two covariance matrices taken from discrete time windows. As long as the spatial structure of the network remains the same, the method will work. In this sense, the reviewer is not correct that the discretization is a weakness. This distinction is now clarified on “Generalized eigendecomposition” and in Figure 3 legend.
With regards to the final point, the revised version now more clearly states the assumptions, advantages, goals, and limitations of each method. This information is listed at the start and end of each subsection of the Results.
5) The paper would be more practically useful if it provided a slightly more thorough treatment of at least one of the statistical methods for testing components, and noted briefly any problems that may arise in this setting.
Statistical evaluations are now presented in two figures, and the section about statistical evaluation is expanded.
Reviewer #3:
This paper describes a novel technique to calculate crossfrequency coupling using generalized eigenvalue decomposition. The technique is novel and has several significant advantages that have currently not been clearly addressed, such as its independence of selecting electrodes and not assuming a sinusoidal shape. Moreover, the technique can be used in several ways (5 methods are described in the paper). The paper is well written; as I was reading it and wrote down questions, a subsequent section would often already provide an answer. Overall, I think the paper is of high quality and I only have a few comments that can be easily addressed.
1) It is not quite clear to me whether the method assumes that multiple electrodes have the exact same phase for an underlying rhythm. For example, if the method were to be used to analyze LFP recordings across motor cortex, with traveling β oscillations, where there is a shift in phase across electrodes, how would that reflect upon component waves?
As long as the phase lags are consistent over repeated epochs, they are not detrimental to the analysis. This is now in the first paragraph of the subsection “Nonstationarities and condition differences”. Method 4 would reveal phase lags as delays in the filter time series for different channels, e.g., a slanted bar in Figure 7c. I have seen this in other datasets, but not in the data shown here.
2) It is said that the quality of gedCFC components is related to the quality of the covariance matrices, and a sufficient amount of clean data ensures highquality covariance matrices. Can this be quantified: is there a way to compute from one of the data sets used in the paper how the amount of data used relates to the quality of the gedCFC components?
This is an important point, and it was too ambiguous in the original submission. Unfortunately, there are no formal criteria for what makes a covariance matrix good enough for eigendecomposition. This is an important issue and it is now discussed in more detail in the paper (subsection “Nonstationarities and condition differences”, third paragraph).
3) In method 3 example recordings from a human epilepsy patient are used with recordings of the medial temporal lobe and electrode Cz. Why would there be any coupling in these data? How does referencing affect the method? More anterior temporal lobe channels are also closer the face with potential muscle artifacts, which may not even be removed with bipolar referencing.
There is considerable evidence for the role of hippocampalprefrontal interactions in memory consolidation and retrieval. This has been studied in rodents using electrophysiology and lesions, and in humans using invasive electrophysiology, MEG, and fMRI. A popular review paper is now cited in the last paragraph of the subsection “Method 3: Lowfrequency waveform shape as a bias filter on sphered data”.
The hippocampal electrode was locally referenced to its average while Cz was referenced to linked mastoids. Intracranial sources, particularly with local referencing, are not affected by facial muscles. Nonetheless, this N=1 dataset was included as proofofprinciple illustration in empirical data; strong claims about brain function are not made based on these data.
https://doi.org/10.7554/eLife.21792.026Article and author information
Author details
Funding
European Research Council (ERCStG 638589)
 Michael X Cohen
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Reviewing Editor
 Lila Davachi, Reviewing Editor, New York University, United States
Publication history
 Received: September 23, 2016
 Accepted: January 8, 2017
 Version of Record published: January 24, 2017 (version 1)
Copyright
© 2017, Cohen et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 997
 Page views

 252
 Downloads

 1
 Citations
Article citation count generated by polling the highest count across the following sources: PubMed Central, Scopus, Crossref.
Comments
Download links
Downloads (link to download the article as PDF)
Download citations (links to download the citations from this article in formats compatible with various reference manager tools)
Open citations (links to open the citations from this article in various online reference manager services)
Further reading

 Genes and Chromosomes